Attachment 'ir0.c'
Download 1 // gcc -o ir0 ir0.c -lm ; ./ir0 > ir0.dat ; gp ir0.gp
2
3 // #define EFILT 0.00 // eV high pass cutoff energy
4 // #define SATTEMP 331.28555 // Kelvin Thinsat temperature
5 // 1366.003A 1366.003T 2.38084 inf
6
7 // #define EFILT 0.30 // eV high pass cutoff energy
8 // #define SATTEMP 391.05045 // Kelvin Thinsat temperature
9 // 1353.854A 512.724T 0.89364 4.133μ
10
11 // #define EFILT 0.34 // eV high pass cutoff energy
12 // #define SATTEMP 391.84015 // Kelvin Thinsat temperature
13 // 1339.437A 503.855T 0.87156 3.647μ
14
15 // #define EFILT 0.35 // eV high pass cutoff energy
16 // #define SATTEMP 391.91000 // Kelvin Thinsat temperature
17 // 1339.437A 503.855T 0.87017 3.542μ
18
19 #define EFILT 0.3542 // eV high pass cutoff energy
20 #define SATTEMP 391.92810 // Kelvin Thinsat temperature
21 // 1349.857A 499.131T 0.86995 3.500μ
22
23 // #define EFILT 0.36 // eV high pass cutoff energy
24 // #define SATTEMP 391.94315 // Kelvin Thinsat temperature
25 // 1345.998A 499.144T 0.86997 3.444μ
26
27 // #define EFILT 0.40 // eV high pass cutoff energy
28 // #define SATTEMP 391.80445 // Kelvin Thinsat temperature
29 // 1339.437A 503.855T 0.87819 3.100μ
30
31 // #define EFILT 0.50 // eV high pass cutoff energy
32 // #define SATTEMP 390.45285 // Kelvin Thinsat temperature
33 // 1318.119A 535.392T 0.93315 2.480μ
34
35 #define SUNTEMP 5778.0 // Kelvin Sun temperature
36 #define BG 1.35 // InP bandgap eV
37
38 // MKS UNITS
39 #define SB 5.670373e-8 // W/m²K Stefan Boltzmann constant
40 #define KB 1.3806488e-23 // J/K Boltzmann constant
41 #define H 6.62606957e-34 // J-s
42 #define C 299792458.0 // m/s speed of light
43 #define Q 1.602176565e-19 // Joules per electron volt
44 #define RS 695500.0 // km Radius of Sun
45 #define RE 149600000.0 // km Average distance to sun
46
47 #define E0 0.0 // eV plot starting energy
48 #define E1 10.0 // eV plot ending energy
49 #define NPLOT 1000 // number of plot points
50 #define NCALC 256 // number of calc steps per point
51
52 #include <stdio.h>
53 #include <math.h>
54
55 // black body radiation:
56 // dP/dEv (W/m²-eV) = 2(qEv/h)³ / ( c² exp(qEv/kT) - 1 )
57
58 double pi , pi2, pih, pis ;
59
60 main() {
61 // pi constants
62 pih = 2.0*atan( 1.0 );
63 pi = 2.0*pih ;
64 pi2 = 4.0*pih ;
65 pis = pi*pi ;
66
67 double dcalc = (double) NCALC ; //
68 double dplot = (double) NPLOT ; //
69 double dev = (E1-E0)/(dcalc*dplot) ; // plot step
70 double sunatten = ( RS/RE )*( RS/RE ) ; // distance attenuation of sun
71 double ktsun = KB*SUNTEMP ;
72 double ktsat = KB*SATTEMP ;
73
74 double lambda = 1E6*C*H / (Q*EFILT) ; // filter wavelength in micrometers
75
76 // FIXME TEST THESE SCALINGS
77 double escale = dev * dcalc ;
78 double pscale1 = pi * 0.001 / dcalc ; // to kW/m²eV graph
79 double pscale2 = pi * dev ; // to W/m² totals
80
81 int ie0, ie1 ;
82
83 double powA = 0.0 ; // W/m² power absorbed ( 1.0 thrust )
84 double powR = 0.0 ; // W/m² power reflected ( 1.0 thrust )
85 double powF = 0.0 ; // W/m² power emitted frontside ( 2/3 thrust)
86 double powB = 0.0 ; // W/m² power emitted backside (-2/3 thrust)
87 double powT, powS, powE ;
88 double powP = 0.0 ; // W/m² solar cell power
89
90 double isat = 2.0 * Q / ( H * H * H * C * C );
91 double isun = isat * sunatten ;
92
93
94 // powA power absorbed = powF + powB --- adjust TSat until true
95 // powT thrust power powA + 2.0*powR + 2/3 ( powF - powB )
96
97 for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) {
98 double de0 = (double)( ie0 * NCALC );
99 double psatavg = 0.0 ;
100 double psunavg = 0.0 ;
101 for( ie1 = 0 ; ie1 < NCALC ; ie1++ ) {
102 double de1 = 0.5 + (double) ie1 ;
103 double ev = dev * ( de0 + de1 );
104 double eng = Q * ev ;
105 double eng3 = eng * eng * eng ;
106 double psat = isat * eng3 / ( exp( eng / ktsat ) - 1.0 ) ;
107 double psun = isun * eng3 / ( exp( eng / ktsun ) - 1.0 ) ;
108 psatavg += psat ; // black body curve for incoming sun
109 psunavg += psun ; // black body curve or outgoing IR
110 powB += psat ; // IR out back
111
112 if( ev < EFILT ) {
113 powR += psun ; // sun power reflected
114 } else {
115 powA += psun ; // sun power absorbed
116 powF += psat ; // IR out front also
117 }
118
119 if( ev > BG ) {
120 powP += (BG/ev)*psun ;
121 }
122 }
123 psatavg *= pscale1 ; // kW/m²-eV
124 psunavg *= pscale1 ; // kW/m²-eV
125 double eavg = escale * (0.5+(double)ie0) ;
126 printf( "%12.4f%16.6e%16.6e\n", eavg, psatavg, psunavg );
127 }
128
129 powP *= pscale2 ;
130 powA *= pscale2 ;
131 powR *= pscale2 ;
132 powF *= pscale2 ;
133 powB *= pscale2 ;
134 powS = powA + powR ;
135 powE = powF + powB ;
136 powT = powA + 2.0*powR + (2.0/3.0) * (powF - powB );
137
138 // print totals
139 printf( "//%10.3f W/m² - sun power\n", powS );
140 printf( "//%10.3f W/m² - emitted power\n", powE );
141 printf( "//%10.3f W/m² - absorbed power\n", powA );
142 printf( "//%10.3f W/m² - reflected power\n", powR );
143 printf( "//%10.3f W/m² - IR power emitted frontside\n", powF );
144 printf( "//%10.3f W/m² - IR power emitted backside\n", powB );
145 printf( "//%10.3f W/m² - forward thrust power\n", powT );
146 printf( "//%10.3f W/m² - photovoltaic illumination power\n", powP );
147 printf( "//%10.5f - thrust/photovoltaic ratio\n", powT/powP );
148 printf( "//%10.3f μm - filter wavelength\n", lambda );
149 }
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.