// gcc -o ir0 ir0.c -lm ; ./ir0 > ir0.dat ; gp ir0.gp

// #define EFILT    0.00                // eV      high pass cutoff energy
// #define SATTEMP  331.28555           // Kelvin  Thinsat temperature
//  1366.003A  1366.003T   2.38084  inf

// #define EFILT    0.30                // eV      high pass cutoff energy
// #define SATTEMP  391.05045           // Kelvin  Thinsat temperature
//  1353.854A   512.724T   0.89364  4.133μ

// #define EFILT    0.34                // eV      high pass cutoff energy
// #define SATTEMP  391.84015           // Kelvin  Thinsat temperature
//  1339.437A   503.855T   0.87156  3.647μ

// #define EFILT    0.35                // eV      high pass cutoff energy
// #define SATTEMP  391.91000           // Kelvin  Thinsat temperature
//  1339.437A   503.855T   0.87017  3.542μ

#define EFILT    0.3542              // eV      high pass cutoff energy
#define SATTEMP  391.92810           // Kelvin  Thinsat temperature
//  1349.857A   499.131T   0.86995  3.500μ

// #define EFILT    0.36                // eV      high pass cutoff energy
// #define SATTEMP  391.94315           // Kelvin  Thinsat temperature
//  1345.998A   499.144T   0.86997  3.444μ

// #define EFILT    0.40                // eV      high pass cutoff energy
// #define SATTEMP  391.80445           // Kelvin  Thinsat temperature
//  1339.437A   503.855T   0.87819  3.100μ

// #define EFILT    0.50                // eV      high pass cutoff energy
// #define SATTEMP  390.45285           // Kelvin  Thinsat temperature
//  1318.119A   535.392T   0.93315  2.480μ

#define SUNTEMP  5778.0              // Kelvin  Sun temperature
#define BG       1.35                // InP bandgap eV

//  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 RE       149600000.0         // km      Average distance to sun

#define E0       0.0                 // eV      plot starting energy
#define E1       10.0                // eV      plot ending energy
#define NPLOT    1000                //         number of plot points
#define NCALC    256                 //         number of calc steps per point

#include <stdio.h>
#include <math.h>

// black body radiation: 
//  dP/dEv (W/m²-eV) =  2(qEv/h)³ / ( c² exp(qEv/kT) - 1 )  

double  pi , pi2, pih, pis ;

main() {
   // pi constants
   pih = 2.0*atan( 1.0 );
   pi  = 2.0*pih ;
   pi2 = 4.0*pih ;
   pis = pi*pi   ;

   double dcalc    = (double) NCALC         ; //
   double dplot    = (double) NPLOT         ; //
   double dev      = (E1-E0)/(dcalc*dplot)  ; // plot step
   double sunatten = ( RS/RE )*( RS/RE )    ; // distance attenuation of sun
   double ktsun    = KB*SUNTEMP ;
   double ktsat    = KB*SATTEMP ;

   double lambda   = 1E6*C*H / (Q*EFILT)    ; // filter wavelength in micrometers

   // 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
   
   int ie0, ie1 ;
    
   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 powP   = 0.0 ;      // W/m²  solar cell power

   double isat   = 2.0 * Q / ( H * H * H * C * C );
   double isun   = isat * sunatten ;


   // powA  power absorbed  = powF + powB  --- adjust TSat until true
   // powT  thrust power      powA + 2.0*powR + 2/3 ( powF - powB )

   for( ie0 = 0 ; ie0 < NPLOT ; ie0++ ) {
      double de0     = (double)( ie0 * NCALC );
      double psatavg = 0.0 ;
      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 ;
         double psat = isat * eng3 / ( exp( eng / ktsat ) - 1.0 ) ;
         double psun = isun * eng3 / ( exp( eng / ktsun ) - 1.0 ) ;
         psatavg += psat ;   // black body curve for incoming sun
         psunavg += psun ;   // black body curve or outgoing IR
         powB    += psat ;   // IR out back

         if( ev < EFILT ) {
            powR += psun ;   // sun power reflected
         } else {
            powA += psun ;   // sun power absorbed
	    powF += psat ;   // IR out front also
         }

         if( ev > BG ) {
            powP += (BG/ev)*psun ;
         }
      }
      psatavg *= pscale1 ;   // kW/m²-eV
      psunavg *= pscale1 ;   // kW/m²-eV
      double  eavg = escale * (0.5+(double)ie0) ;
      printf( "%12.4f%16.6e%16.6e\n", eavg, psatavg, psunavg );
   }
   
   powP *= pscale2 ;
   powA *= pscale2 ;
   powR *= pscale2 ;
   powF *= pscale2 ;
   powB *= pscale2 ;
   powS  = powA + powR ;
   powE  = powF + powB ;
   powT  = powA + 2.0*powR + (2.0/3.0) * (powF - powB );

   // print totals
   printf( "//%10.3f W/m² - sun power\n", powS );
   printf( "//%10.3f W/m² - emitted power\n", powE );
   printf( "//%10.3f W/m² - absorbed power\n", powA  );
   printf( "//%10.3f W/m² - reflected power\n", powR );
   printf( "//%10.3f W/m² - IR power emitted frontside\n", powF );
   printf( "//%10.3f W/m² - IR power emitted backside\n",  powB );
   printf( "//%10.3f W/m² - forward thrust power\n", powT );
   printf( "//%10.3f W/m² - photovoltaic illumination power\n", powP );
   printf( "//%10.5f      - thrust/photovoltaic ratio\n", powT/powP );
   printf( "//%10.3f μm   - filter wavelength\n", lambda );
}