Attachment 'ir1.c'
Download 1 // gcc -o ir1 ir1.c -lm ; ./ir1 > ir1.dat ; gp ir1.gp
2
3 #define EFILT 0.3542 // eV high pass cutoff energy NOT USED
4 #define SATTEMP 391.92810 // Kelvin Thinsat temperature STARTING
5 #define SUNTEMP 5777.992911 // Kelvin Sun temperature
6 #define BG 1.35 // ev InP bandgap
7
8 // MKS UNITS
9 #define SB 5.670373e-8 // W/m²K Stefan Boltzmann constant
10 #define KB 1.3806488e-23 // J/K Boltzmann constant
11 #define H 6.62606957e-34 // J-s
12 #define C 299792458.0 // m/s speed of light
13 #define Q 1.602176565e-19 // Joules per electron volt
14 #define RS 695500.0 // km Radius of Sun
15 #define DE 149600000.0 // km Average distance to sun
16
17 #define E0 0.0 // eV plot starting energy
18 #define E1 4.0 // eV plot ending energy
19 #define NMAX 5000 // eV maximum sun energy integration
20 #define NPLOT 1000 // number of plot points
21 #define NCALC 256 // number of calc steps per point
22 #define EFILT0 0.1 // eV minimum high pass cutoff energy
23 #define EFILT1 0.6 // eV maximum high pass cutoff energy
24
25 #include <stdlib.h>
26 #include <stdio.h>
27 #include <math.h>
28
29 double sunpower[NPLOT] ; // W/m²-eV power absorbed by thinsat
30 double suntotal[NPLOT] ; // W/m² power absorbed by thinsat
31 double heatmin[ NPLOT] ; // W/m²-eV satellite heating at high cutoff
32 double heatmax[ NPLOT] ; // W/m²-eV satellite heating at low cutoff
33 double thrust[ NPLOT] ; // W/m² satellite thrust power
34 double sattemp[ NPLOT] ; // K satellite temperature
35
36 // black body radiation:
37 // dP/dEv (W/m²-eV) = 2(qEv/h)³ / ( c² exp(qEv/kT) - 1 )
38
39 double pi , pi2, pih, pis ;
40
41 int main() {
42 // pi constants
43 pih = 2.0*atan( 1.0 );
44 pi = 2.0*pih ;
45 pi2 = 4.0*pih ;
46 pis = pi*pi ;
47
48 int ie0, ie1 ;
49
50 double dcalc = (double) NCALC ; //
51 double dplot = (double) NPLOT ; //
52 double dev = (E1-E0)/(dcalc*dplot) ; // plot step
53 double sunatten = ( RS/DE )*( RS/DE ) ; // distance attenuation of sun
54 double ktsun = KB*SUNTEMP ;
55
56 int ia = (int)( dplot *(EFILT0-E0)/(E1-E0) ) ; // start filter sweep
57 int ib = (int)( dplot *(EFILT1-E0)/(E1-E0) )+1 ; // end filter sweep
58
59 // FIXME TEST THESE SCALINGS
60 double escale = dev * dcalc ;
61 double pscale1 = pi * 0.001 / dcalc ; // to kW/m²eV graph
62 double pscale2 = pi * dev ; // to W/m² totals
63 double isat = 2.0 * Q / ( H * H * H * C * C );
64 double isun = isat * sunatten ;
65 double powtotal = 0.0 ;
66
67 // compute sunpower graph
68 for( ie0 = NMAX ; ie0 >= 0 ; ie0-- ) {
69 double de0 = (double)( ie0 * NCALC );
70 double psunavg = 0.0 ;
71 for( ie1 = 0 ; ie1 < NCALC ; ie1++ ) {
72 double de1 = 0.5 + (double) ie1 ;
73 double ev = dev * ( de0 + de1 );
74 double eng = Q * ev ;
75 double eng3 = eng * eng * eng ;
76 psunavg += isun * eng3 / ( exp( eng / ktsun ) - 1.0 ) ;
77 }
78 if( ie0 < NPLOT ) {
79 sunpower[ ie0 ] = pscale1*psunavg ;
80 suntotal[ ie0 ] = powtotal ;
81 }
82 powtotal += pscale2 * psunavg ;
83 }
84
85 // powtotal contains total sun power integrated from high energy to zero
86 // suntotal contains the total power down to the associated energy
87 // sunpower contains the incremental power per electron volt
88
89 double ea = dev * (double) ( ia * NCALC );
90 double eb = dev * (double) ( ib * NCALC );
91 double powa = suntotal[ ia ] ;
92 double powb = suntotal[ ib ] ;
93
94 printf( "// powtotal=%12.6f\n", powtotal );
95 printf( "// powa =%12.6f at %10.4feV\n", powa, ea );
96 printf( "// powb =%12.6f at %10.4feV\n", powb, eb );
97
98 double ktsat = KB*SATTEMP; //
99 double powA = 0.0 ; // W/m² power absorbed ( 1.0 thrust )
100 double powR = 0.0 ; // W/m² power reflected ( 1.0 thrust )
101 double powF = 0.0 ; // W/m² power emitted frontside ( 2/3 thrust)
102 double powB = 0.0 ; // W/m² power emitted backside (-2/3 thrust)
103 double powT, powS, powE ;
104
105 double powthr = 1e-3 ;
106
107 int ie ; // energy loop
108
109 fprintf( stderr, "sun finished\n" );
110 fflush( stderr );
111
112 for( ie = ia ; ie <= ib ; ie++ ) {
113 // add binary search loop to seek ktsat ;
114 double eavg = escale * ((double)ie) ;
115 double sattemp = 400.0 ;
116 double dsattemp = 200.0 ;
117 double psun = suntotal[ie] ;
118 int convg = 0 ;
119 double pfront ;
120 double pback ;
121 while( convg < 2 ) {
122 double ktsat = KB * sattemp ;
123 pfront = 0.0 ;
124 pback = 0.0 ;
125 for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) {
126 double de0 = (double)( ie0 * NCALC );
127 double psat = 0.0 ;
128 for( ie1 = 0 ; ie1 < NCALC ; ie1++ ) {
129 double eng = Q * dev * ( 0.5 + de0 + (double) ie1 );
130 double eng3 = eng * eng * eng ;
131 psat += eng3 / ( exp( eng / ktsat ) - 1.0 ) ;
132 }
133 double pheat = psat * isat * pscale2 ;
134 psat *= isat * pscale1 ; // kW/m²-eV
135
136 if( convg == 1 ) {
137 if( ie == ia ) { heatmin[ie0] = psat ; }
138 if( ie == ib ) { heatmax[ie0] = psat ; }
139 }
140 pback += pheat ;
141 if( ie0 > ie ) { pfront += pheat ; }
142 }
143 double dpow = pfront + pback - psun ;
144 if( dpow < -powthr ) { sattemp += dsattemp ; }
145 else if( dpow > powthr ) { sattemp -= dsattemp ; }
146 else { convg++ ; }
147
148 dsattemp *= 0.5 ;
149 fprintf( stderr, "%6.3f%10.4f%12.3f%12.3f%12.3f\r",
150 eavg, sattemp, pfront, pback, psun );
151 fflush( stderr );
152 }
153 thrust[ ie ] = (2.0/3.0)*(pfront-pback)+2.0*powtotal-psun ;
154 }
155 fprintf( stderr, "filter loop done\n" );
156 fflush( stderr );
157
158 for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) {
159 double eavg = escale * (0.5+(double)ie0) ;
160 if( ( ie0 < ia ) || ( ie0 > ib ) ) {
161 printf("%7.4f%12.4e%12.4e%12.4e\n",
162 eavg, heatmin[ie0], heatmax[ie0], sunpower[ie0] );
163 } else {
164 printf("%7.4f%12.4e%12.4e%12.4e%12.2f\n",
165 eavg, heatmin[ie0], heatmax[ie0], sunpower[ie0],thrust[ie0] );
166 }
167 }
168 fprintf( stderr, "printout done\n" );
169 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.