attachment:irin01.c of IRtelescope

   1 // irin01.c
   2 // gcc -o irin01 irin01.c -lm ; ./irin01 
   3 // percentile wavelengths
   4 
   5 #define	RS	149.6E11   	// m	 100AU Radius of shell
   6 #define LY	9.46E+15	// m/LY	 Light year
   7 #define R100    (100.0*LY)      // m     100 light years
   8 
   9 #define Q	1.60217E-19	// Coul	 Electronic Charge
  10 #define K 	1.3806488e-23	// J/K   Boltzmann constant
  11 #define	H	6.62607E-34	// J-s   Planck constant
  12 #define C	299792458.0	// m/s	 Speed of light
  13 #define WEIN    2897.7721       // K-μm  Wein displacement law constant
  14 #define N       1000000         //       number of integration steps
  15 #define PSCALE  6.598811205712441e+15 // KLUDGE power scaling
  16 
  17 #include <stdlib.h>
  18 #include <stdio.h>
  19 #include <math.h>
  20 
  21 int main(int argc, char *argv[] ) {
  22    double T                                   ; // K   temperature
  23    double wl                                  ; // μm  mean wavelength 
  24    double wr                                  ; // λ/Δλ  wavelength ratio
  25 
  26    if( argc !=4 ) {
  27      fprintf( stderr, "usage: irin0 Temp(K) λ(μm) λ/Δλ\n");
  28      return 1 ;
  29    }
  30 
  31    sscanf( argv[1], "%lf", &T  );             ; // temperature
  32    sscanf( argv[2], "%lf", &wl );             ; // mean wavelength
  33    sscanf( argv[3], "%lf", &wr );             ; // mean wavelength
  34 
  35    double  hckt = 1e6*H*C/(K*T)               ; // μm exponent
  36    double  rad  = sqrt( 4.0*wr*wr + 1.0 )     ; // eqn square root
  37    double  mlt  = 0.5 * wl / wr               ; // eqn scale
  38    double  wl1  = mlt*(rad-1.0)               ; // μm  min wavelength
  39    double  wl2  = mlt*(rad+1.0)               ; // μm  max wavelength
  40    double  ns   = (double)N                   ; //     number of steps
  41    double  ws   = (wl2-wl1)/ns                ; // μm  integration step
  42    double  w0   = wl1 + 0.5*ws                ; // μm  first integration wl
  43 
  44    double  scal = ws*PSCALE/(T*T*T*T)         ;
  45 
  46    double p                                   ; // W/m² integrated BB power
  47    double w                                   ;
  48    int    n                                   ;
  49 
  50    for(	 n = 0 ; n < N ; n++ ) 
  51    {
  52       w  = w0 + ((double)n)*ws                ;
  53       p += 1.0/(w*w*w*w*w*(exp(hckt/w)-1.0 )) ;
  54    }
  55 
  56    p *= scal                                  ;
  57    double wein  = WEIN/T                      ;
  58 
  59    printf("T    = %12.5e K \n", T    )        ;
  60    printf("peak = %12.5e μm\n", wein )        ;
  61    printf("hckt = %12.5e   \n", hckt )        ;
  62    printf("wl   = %12.5e μm\n", wl   )        ;
  63    printf("wr   = %12.5e μm\n", wr   )        ;
  64    printf("ws   = %12.5e μm\n", ws   )        ;
  65    printf("wl1  = %12.5e μm\n", wl1  )        ;
  66    printf("wl2  = %12.5e μm\n", wl2  )        ;
  67    printf("frac = %12.5e   \n", p    )        ;
  68    printf("frac = %12.2e   \n", p    )        ;
  69 
  70 /*
  71    printf("scal = %12.3e\n", scal   )         ;
  72    printf("1/p  = %20.12e\n", 1.0/p )         ;
  73    printf("pm1  = %20.12e\n", p-1.0 )         ;
  74    printf("power fraction =%20.12e\n", p )    ;
  75    printf("PSC  = %24.16e\n", PSCALE  )       ;
  76    printf("PS   = %24.16e\n", PSCALE/p)       ;
  77 */
  78 
  79    return 0 ;
  80 }