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.
  • [get | view] (2013-08-14 23:22:23, 9.0 KB) [[attachment:g200.c]]
  • [get | view] (2013-08-14 21:02:48, 814.8 KB) [[attachment:g200_1024.swf]]
  • [get | view] (2013-08-14 23:22:33, 814.9 KB) [[attachment:g200_384.swf]]
  • [get | view] (2013-08-14 21:03:12, 814.9 KB) [[attachment:g200_512.swf]]
  • [get | view] (2013-08-29 18:50:31, 0.2 KB) [[attachment:gridhex1.ini]]
  • [get | view] (2013-08-29 08:11:48, 117.3 KB) [[attachment:gridhex1.png]]
  • [get | view] (2013-08-29 18:50:17, 1.3 KB) [[attachment:gridhex1.pov]]
  • [get | view] (2013-08-08 05:34:27, 118.4 KB) [[attachment:ir-reflect0.png]]
  • [get | view] (2013-08-08 05:23:38, 5.8 KB) [[attachment:ir0.c]]
  • [get | view] (2013-08-08 05:23:49, 1.9 KB) [[attachment:ir0.gp]]
  • [get | view] (2013-08-08 05:24:06, 16.9 KB) [[attachment:ir0.png]]
  • [get | view] (2013-08-22 15:39:57, 6.6 KB) [[attachment:ir1.c]]
  • [get | view] (2013-08-22 15:40:06, 2.3 KB) [[attachment:ir1.gp]]
  • [get | view] (2013-08-22 15:40:15, 16.5 KB) [[attachment:ir1.png]]
  • [get | view] (2013-08-14 22:55:55, 4.2 KB) [[attachment:irn1.c]]
  • [get | view] (2013-08-14 22:56:15, 1.8 KB) [[attachment:irn1.gp]]
  • [get | view] (2013-08-14 22:56:26, 17.9 KB) [[attachment:irn1.png]]
 All files | Selected Files: delete move to page

You are not allowed to attach a file to this page.