// gcc -o ir1 ir1.c -lm ; ./ir1 > ir1.dat ; gp ir1.gp #define EFILT 0.3542 // eV high pass cutoff energy NOT USED #define SATTEMP 391.92810 // Kelvin Thinsat temperature STARTING #define SUNTEMP 5777.992911 // Kelvin Sun temperature #define BG 1.35 // ev InP bandgap // MKS UNITS #define SB 5.670373e-8 // W/m²K Stefan Boltzmann constant #define KB 1.3806488e-23 // J/K Boltzmann constant #define H 6.62606957e-34 // J-s #define C 299792458.0 // m/s speed of light #define Q 1.602176565e-19 // Joules per electron volt #define RS 695500.0 // km Radius of Sun #define DE 149600000.0 // km Average distance to sun #define E0 0.0 // eV plot starting energy #define E1 4.0 // eV plot ending energy #define NMAX 5000 // eV maximum sun energy integration #define NPLOT 1000 // number of plot points #define NCALC 256 // number of calc steps per point #define EFILT0 0.1 // eV minimum high pass cutoff energy #define EFILT1 0.6 // eV maximum high pass cutoff energy #include #include #include double sunpower[NPLOT] ; // W/m²-eV power absorbed by thinsat double suntotal[NPLOT] ; // W/m² power absorbed by thinsat double heatmin[ NPLOT] ; // W/m²-eV satellite heating at high cutoff double heatmax[ NPLOT] ; // W/m²-eV satellite heating at low cutoff double thrust[ NPLOT] ; // W/m² satellite thrust power double sattemp[ NPLOT] ; // K satellite temperature // black body radiation: // dP/dEv (W/m²-eV) = 2(qEv/h)³ / ( c² exp(qEv/kT) - 1 ) double pi , pi2, pih, pis ; int main() { // pi constants pih = 2.0*atan( 1.0 ); pi = 2.0*pih ; pi2 = 4.0*pih ; pis = pi*pi ; int ie0, ie1 ; double dcalc = (double) NCALC ; // double dplot = (double) NPLOT ; // double dev = (E1-E0)/(dcalc*dplot) ; // plot step double sunatten = ( RS/DE )*( RS/DE ) ; // distance attenuation of sun double ktsun = KB*SUNTEMP ; int ia = (int)( dplot *(EFILT0-E0)/(E1-E0) ) ; // start filter sweep int ib = (int)( dplot *(EFILT1-E0)/(E1-E0) )+1 ; // end filter sweep // FIXME TEST THESE SCALINGS double escale = dev * dcalc ; double pscale1 = pi * 0.001 / dcalc ; // to kW/m²eV graph double pscale2 = pi * dev ; // to W/m² totals double isat = 2.0 * Q / ( H * H * H * C * C ); double isun = isat * sunatten ; double powtotal = 0.0 ; // compute sunpower graph for( ie0 = NMAX ; ie0 >= 0 ; ie0-- ) { double de0 = (double)( ie0 * NCALC ); double psunavg = 0.0 ; for( ie1 = 0 ; ie1 < NCALC ; ie1++ ) { double de1 = 0.5 + (double) ie1 ; double ev = dev * ( de0 + de1 ); double eng = Q * ev ; double eng3 = eng * eng * eng ; psunavg += isun * eng3 / ( exp( eng / ktsun ) - 1.0 ) ; } if( ie0 < NPLOT ) { sunpower[ ie0 ] = pscale1*psunavg ; suntotal[ ie0 ] = powtotal ; } powtotal += pscale2 * psunavg ; } // powtotal contains total sun power integrated from high energy to zero // suntotal contains the total power down to the associated energy // sunpower contains the incremental power per electron volt double ea = dev * (double) ( ia * NCALC ); double eb = dev * (double) ( ib * NCALC ); double powa = suntotal[ ia ] ; double powb = suntotal[ ib ] ; printf( "// powtotal=%12.6f\n", powtotal ); printf( "// powa =%12.6f at %10.4feV\n", powa, ea ); printf( "// powb =%12.6f at %10.4feV\n", powb, eb ); double ktsat = KB*SATTEMP; // double powA = 0.0 ; // W/m² power absorbed ( 1.0 thrust ) double powR = 0.0 ; // W/m² power reflected ( 1.0 thrust ) double powF = 0.0 ; // W/m² power emitted frontside ( 2/3 thrust) double powB = 0.0 ; // W/m² power emitted backside (-2/3 thrust) double powT, powS, powE ; double powthr = 1e-3 ; int ie ; // energy loop fprintf( stderr, "sun finished\n" ); fflush( stderr ); for( ie = ia ; ie <= ib ; ie++ ) { // add binary search loop to seek ktsat ; double eavg = escale * ((double)ie) ; double sattemp = 400.0 ; double dsattemp = 200.0 ; double psun = suntotal[ie] ; int convg = 0 ; double pfront ; double pback ; while( convg < 2 ) { double ktsat = KB * sattemp ; pfront = 0.0 ; pback = 0.0 ; for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) { double de0 = (double)( ie0 * NCALC ); double psat = 0.0 ; for( ie1 = 0 ; ie1 < NCALC ; ie1++ ) { double eng = Q * dev * ( 0.5 + de0 + (double) ie1 ); double eng3 = eng * eng * eng ; psat += eng3 / ( exp( eng / ktsat ) - 1.0 ) ; } double pheat = psat * isat * pscale2 ; psat *= isat * pscale1 ; // kW/m²-eV if( convg == 1 ) { if( ie == ia ) { heatmin[ie0] = psat ; } if( ie == ib ) { heatmax[ie0] = psat ; } } pback += pheat ; if( ie0 > ie ) { pfront += pheat ; } } double dpow = pfront + pback - psun ; if( dpow < -powthr ) { sattemp += dsattemp ; } else if( dpow > powthr ) { sattemp -= dsattemp ; } else { convg++ ; } dsattemp *= 0.5 ; fprintf( stderr, "%6.3f%10.4f%12.3f%12.3f%12.3f\r", eavg, sattemp, pfront, pback, psun ); fflush( stderr ); } thrust[ ie ] = (2.0/3.0)*(pfront-pback)+2.0*powtotal-psun ; } fprintf( stderr, "filter loop done\n" ); fflush( stderr ); for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) { double eavg = escale * (0.5+(double)ie0) ; if( ( ie0 < ia ) || ( ie0 > ib ) ) { printf("%7.4f%12.4e%12.4e%12.4e\n", eavg, heatmin[ie0], heatmax[ie0], sunpower[ie0] ); } else { printf("%7.4f%12.4e%12.4e%12.4e%12.2f\n", eavg, heatmin[ie0], heatmax[ie0], sunpower[ie0],thrust[ie0] ); } } fprintf( stderr, "printout done\n" ); }