// nl02.c night illumination // random, full, partial, and non-illuminating orientations // circular orbits - we will do elliptical later, as well as solar inclination // // cc -o nl02 nl02.c -lm // // This produces two graphs: // 1) Illumination versus hour of the night, in W/m^4^ // multiply by m^2^ of thinsat // 2) Power versus orbital angle // // It produces three numbers per orientation: // 1) Peak equatorial illumination at 7pm // 2) Average power around orbit // 3) // // sunwards is 0=noon, orbit 0 degrees // // note: I tried to use a min() function, that demands a -std=c99 flag // (which is very silly) // // NOTE - - INCLUDE FAR SIDE ILLUMINATION // #define DEBUG 1 #define NAME "nl02" #define HOURSTART 6.0 // Start of night #define HOUREND 12.0 // Middle of night (use symmetry) #define DHOUR 0.0625 // Time steps through the night #define DORB 0.01 // Orbital angle degree steps // #define DHOUR 1.0 // Time steps through the night // #define DORB 5.0 // Orbital angle degree steps #define ALBEDO 0.5 // Worst case thinsat albedo for night light #define EFF 0.123 // efficiency of solar cell #define FILL 0.7 // solar cell fill #define SUNPOWER 1367.0 // W/m^2 Sun illumination #define TW 1E12 // Terawatt #define AXIAL 22.3 // Axial tilt degrees #define FULLMOON 27E-3 // W/m^2 Full moon illumination #define C 2.997E8 // Speed of light KM per sec #define KM 1000.0 // 1000m/km #define RE (6378.0*KM) // Earth equatorial radius kilometers #define M288 (12789.0*KM) // Server Sky orbit radius kilometers #define RORB (M288) #define PI 3.1415926535897932 #define PI2 (PI/2.0) #define RAD (PI/180.0) //---------------------------------------------------------------- // full moon equivalents per terawatt #define SATAREA (1E12/(FILL*EFF*SUNPOWER)) #define REFLECT (SATAREA*ALBEDO*SUNPOWER) #define FET ((DORB/360.0)*(REFLECT/FULLMOON)) //----------------------------------------------------------------- #include #include #include #include //----------------------------------------------------------------- int main() { FILE *datfile ; // char gnuplot[200] ; // char filename[80] ; // double hour ; // double orbit ; // double asoe = asin( RE / RORB ) ; // shadow angle of earth on orbit double acoe = acos( RE / RORB ) ; // shadow angle of earth on orbit double adorb = RAD*DORB ; // orbital increment angle // first, compute power level ---------------------------------- double pmax = 90.0/DORB ; double pmin = 90.0/DORB ; double pzero = 90.0/DORB ; sprintf( filename, "%spow.dat", NAME ); // power datfile = fopen( filename, "wb" ); for( orbit=90.0+DORB/2.0 ; orbit < 180.0 ; orbit += DORB ) { double aorbit = RAD*orbit ; double xorbit = RORB*sin(aorbit); double yorbit = -RORB*cos(aorbit); double imax = 0.0 ; double imin = 0.0 ; double izero = 0.0 ; if( xorbit > RE ) { imax = 1.0 ; imin = sin( aorbit ); izero = cos( atan2( yorbit, xorbit-RE )); } pmax += imax ; pmin += imin ; pzero += izero ; fprintf( datfile, "%8.3f%9.5f%9.5f%9.5f\n", orbit, imax, imin, izero ); } fclose( datfile ); // normalize power pmax *= DORB/180.0 ; pmin *= DORB/180.0 ; pzero *= DORB/180.0 ; printf( "power max=%6.4f min=%6.4f zero=%6.4f\n", pmax, pmin, pzero ); // Scaling factor for full moon equivalents. // this provides backend scaling for the night // time scaling integtrated segment of orbit. // // FET = SATAREA // for 1 terawatt of thinsats, // = 1E12/( fillfactor // * efficiency // * sunpower ) // * ALBEDO * sunpower // * DORB/360 // scaling for integration // / FULLMOON // full moon power, small // // see the #define statements above double femax = FET/pmax ; double femin = FET/pmin ; double fezero = FET/pzero ; double ferand = FET/pzero ; // pzero smallest printf( "%12.3e=SATAREA%12.3e=ALBEDO%12.3e=DORB360%12.3e=FULLMOON\n", SATAREA, ALBEDO, DORB/360.0, FULLMOON ); printf( "%12.3e=femax%12.3e=femin%12.3e=fezero%12.3e=ferand\n", femax, femin, fezero, ferand ); // second, compute night time illumination ---------------------- sprintf( filename, "%snlp.dat", NAME ); // power datfile = fopen( filename, "wb" ); for( hour=HOURSTART ; hour <= HOUREND ; hour += DHOUR ) { double ahour = RAD*15.0*hour ; double xhour = RE*sin(ahour); // 0 radians towards sun double yhour = -RE*cos(ahour); double nlrand = 0.0; double nlmax = 0.0; double nlmin = 0.0; double nlzero = 0.0; double aostart = 0.5*adorb + ahour - acoe ; double aostop = ahour + acoe ; // integrate over visible thinsats double aorbit ; for( aorbit=aostart ; aorbit < aostop ; aorbit += adorb ) { double xorbit = RORB*sin(aorbit); // 0 radians towards sun double yorbit = -RORB*cos(aorbit); double distx = xorbit-xhour ; double disty = yorbit-yhour ; double dist2 = distx*distx + disty*disty ; double dist = sqrt( dist2 ); double adist = atan2( distx, -disty ); // 0 radians towards sun // assume in shadow, make perpendicular (fake) for 0 light double amax = RAD*90.0 ; double amin = RAD*90.0 ; double azero = RAD*90.0 ; double irand = 0.00 ; if( yorbit < 0.0 ) { // in front of earth, daylight sky amax = 0.0 ; amin = 0.0 ; azero = 0.0 ; irand = 0.25 ; } else if( xorbit > RE ) { // westward of earth amax = 0.0 ; amin = aorbit-RAD*90.0 ; azero = atan2( yorbit, xorbit-RE ); irand = 0.25 ; } else if( xorbit < -RE ) { // eastward of earth amax = 0.0 ; amin = aorbit-RAD*270.0 ; azero = atan2( -yorbit, -xorbit-RE ); irand = 0.25 ; } // thinsat illumination double imax = cos( amax ); double imin = cos( amin ); double izero = cos( azero ); // distance attenuation, receive angle // cos( ahour-adist -pi/2 ) = sin( ahour-adist ) // Lambert, transmit angle double cadmax = -cos(adist - amax ) ; double cadmin = -cos(adist - amin ) ; double cadzero = -cos(adist - azero ) ; // receive angle double cp = cos(ahour-adist)/dist2 ; if( cp > 0.0 ) { // add increment, unnormalized if( cadmax > 0.0 ) { nlmax += cp * imax * cadmax ; } if( cadmin > 0.0 ) { nlmin += cp * imin * cadmin ; } if( cadzero > 0.0 ) { nlzero += cp * izero * cadzero ; } nlrand += cp * irand ; } #ifdef DEBUG printf( "X%9.4f=hour %9.4f=hstart %9.4f=hstop %9.4f=horbit\n", hour, aostart/(RAD*15.0), aostop/(RAD*15.0), aorbit/(RAD*15.0) ); printf( "X%9.4f=aorbit%9.4f=amax%11.4f=amin", aorbit/RAD, amax/RAD, amin/RAD ); printf( "%11.3f=azero%9.3f=ahour\n", azero/RAD, ahour/RAD ); printf( "X%10.1f=xorbit%10.1f=yorbit", xorbit/KM, yorbit/KM ); printf( "%9.1f=distx%10.1f=disty%9.1f=dist\n", distx/KM, disty/KM, dist/KM ); printf( "X%9.4f=adist %9.4f=cadmax%9.4f=cadmin%9.4f=cadzero\n", adist/RAD, cadmax, cadmin, cadzero ); printf( "X%12.3e= irand %12.3e= imax%13.3e= imin%13.3e= izero\n", irand, imax, imin, izero ); printf( "X%11.5f=cos(recv)%13.3e=dist2%13.3e=cp\n", cp*dist2, dist2, cp ); printf( "X%12.3e=nlrand%13.3e=nlmax%13.3e=nlmin%13.3e=nlzero\n\n", nlrand, nlmax, nlmin, nlzero ); #endif // DEBUG } // end of orbit loop nlmax *= femax ; nlmin *= femin ; nlzero *= fezero ; nlrand *= ferand ; double nlmoon = -cos( ahour ); fprintf( datfile, "%8.4f %11.4e %11.4e %11.4e %11.4e %11.4e\n", hour, nlrand, nlmax, nlmin, nlzero, nlmoon ); } fclose( datfile ); }