// irin01.c
// gcc -o irin01 irin01.c -lm ; ./irin01 
// percentile wavelengths

#define	RS	149.6E11   	// m	 100AU Radius of shell
#define LY	9.46E+15	// m/LY	 Light year
#define R100    (100.0*LY)      // m     100 light years

#define Q	1.60217E-19	// Coul	 Electronic Charge
#define K 	1.3806488e-23	// J/K   Boltzmann constant
#define	H	6.62607E-34	// J-s   Planck constant
#define C	299792458.0	// m/s	 Speed of light
#define WEIN    2897.7721       // K-μm  Wein displacement law constant
#define N       1000000         //       number of integration steps
#define PSCALE  6.598811205712441e+15 // KLUDGE power scaling

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

int main(int argc, char *argv[] ) {
   double T                                   ; // K   temperature
   double wl                                  ; // μm  mean wavelength 
   double wr                                  ; // λ/Δλ  wavelength ratio

   if( argc !=4 ) {
     fprintf( stderr, "usage: irin0 Temp(K) λ(μm) λ/Δλ\n");
     return 1 ;
   }

   sscanf( argv[1], "%lf", &T  );             ; // temperature
   sscanf( argv[2], "%lf", &wl );             ; // mean wavelength
   sscanf( argv[3], "%lf", &wr );             ; // mean wavelength

   double  hckt = 1e6*H*C/(K*T)               ; // μm exponent
   double  rad  = sqrt( 4.0*wr*wr + 1.0 )     ; // eqn square root
   double  mlt  = 0.5 * wl / wr               ; // eqn scale
   double  wl1  = mlt*(rad-1.0)               ; // μm  min wavelength
   double  wl2  = mlt*(rad+1.0)               ; // μm  max wavelength
   double  ns   = (double)N                   ; //     number of steps
   double  ws   = (wl2-wl1)/ns                ; // μm  integration step
   double  w0   = wl1 + 0.5*ws                ; // μm  first integration wl

   double  scal = ws*PSCALE/(T*T*T*T)         ;

   double p                                   ; // W/m² integrated BB power
   double w                                   ;
   int    n                                   ;

   for(	 n = 0 ; n < N ; n++ ) 
   {
      w  = w0 + ((double)n)*ws                ;
      p += 1.0/(w*w*w*w*w*(exp(hckt/w)-1.0 )) ;
   }

   p *= scal                                  ;
   double wein  = WEIN/T                      ;

   printf("T    = %12.5e K \n", T    )        ;
   printf("peak = %12.5e μm\n", wein )        ;
   printf("hckt = %12.5e   \n", hckt )        ;
   printf("wl   = %12.5e μm\n", wl   )        ;
   printf("wr   = %12.5e μm\n", wr   )        ;
   printf("ws   = %12.5e μm\n", ws   )        ;
   printf("wl1  = %12.5e μm\n", wl1  )        ;
   printf("wl2  = %12.5e μm\n", wl2  )        ;
   printf("frac = %12.5e   \n", p    )        ;
   printf("frac = %12.2e   \n", p    )        ;

/*
   printf("scal = %12.3e\n", scal   )         ;
   printf("1/p  = %20.12e\n", 1.0/p )         ;
   printf("pm1  = %20.12e\n", p-1.0 )         ;
   printf("power fraction =%20.12e\n", p )    ;
   printf("PSC  = %24.16e\n", PSCALE  )       ;
   printf("PS   = %24.16e\n", PSCALE/p)       ;
*/

   return 0 ;
}