attachment:irn1.c of IRfilter

   1 // gcc -o irn1 irn1.c -lm ; ./irn1 > irn1.dat ; gp irn1.gp
   2 // 2013/08/09
   3 // infrared filtered thinsat,  night orientation and temperature
   4 //   circular orbit approximation, no light pressure eccentricity 
   5 
   6 #define VMULT    1.16625404993432468 //          Multiply starting velocity
   7 
   8 #define EFILT    0.3542              // eV       high pass cutoff energy
   9 #define ETEMP    260.0               // Kelvin   Earth temperature
  10 
  11 //  MKS UNITS
  12 #define SB       5.670373e-8         // W/m²K    Stefan Boltzmann constant
  13 #define KB       1.3806488e-23       // J/K      Boltzmann constant
  14 #define H        6.62606957e-34      // J-s      Planck's constant
  15 #define C        299792458.0         // m/s      speed of light 
  16 #define Q        1.602176565e-19     // J/eV     Joules per electron volt
  17 #define RE       6378.0              // km       Average radius of earth
  18 #define RSAT     12789.0             // km       Average distance to earth
  19 #define ANG0    -60.0                // degrees  Starting earth relative angle
  20 #define ANG1     60.0                // degrees  Ending earth relative angle
  21 #define OANG0    150.0               // degrees  Starting orbit angle from noon
  22 #define OANG1    210.0               // degrees  Ending orbit angle
  23 #define NPLOT    401                 //          number of plot points
  24 #define NCALC    10000               //          number of calculation points
  25 #define OTIME    14400.0             // seconds  Orbit period
  26 
  27 #include <stdio.h>
  28 #include <math.h>
  29 
  30 double  pi , pi2, pih, pis, d2r, r2d ;
  31 
  32 main() {
  33    // pi constants
  34    pih = 2.0*atan( 1.0 );
  35    pi  = 2.0*pih ;
  36    pi2 = 4.0*pih ;
  37    pis = pi*pi   ;
  38    d2r = pi/180.0 ;
  39    r2d = 180.0/pi ;
  40 
  41    double dplot1   = (double) (NPLOT-1)     ; // 
  42    double dcalc    = (double) NCALC         ; // 
  43    double dintrvl  = dplot1*dcalc           ; // number of calculation intervals
  44    double angchng  = ANG1-ANG0              ; // angular change in degrees
  45    double rangchng = d2r*angchng            ; // angular change in radians
  46    double dang     = angchng/dintrvl        ; // plot calculation step
  47    double t0       = OTIME*OANG0/360.0      ; // start time in seconds
  48    double t1       = OTIME*OANG1/360.0      ; // end time in seconds
  49    double tchng    = t1-t0                  ; // simulation time, seconds
  50    double dt       = tchng/dintrvl          ; // time step, seconds
  51    double dtplot   = tchng/dplot1           ; // plot time step, seconds
  52    double ktearth  = KB*ETEMP               ; // earth IR temperature
  53    double lambda   = 1E6*C*H / (Q*EFILT)    ; // filter wavelength in micrometers
  54    double omega    = pi2/OTIME              ; // orbit angular frequency
  55    double rav0     = VMULT*rangchng/tchng   ; // starting angular velocity
  56    double ascale   = 1.5*omega*omega        ; // scale for angular acceleration
  57    double rang0    = d2r * ANG0             ; //
  58 
  59    double rhoe     = asin( RE / RSAT )      ; // earth apparent angle
  60    double sang0    = 1.0-cos(2.0*rhoe)      ; // earth solid angle
  61    double temp0    = ETEMP*sqrt(sqrt(0.5*sang0)) ; // temperature scaling
  62 
  63    int i0, i1 ;
  64 
  65    // orbit loop 
  66    
  67    double t        = t0                     ; // time in seconds
  68    double rang     = rang0                  ; // thinsat angle in radians
  69    double rav      = rav0                   ; // angular velocity rad/sec
  70    double rangcel                           ; // angular acceleration
  71 
  72    for( i0 = 0 ; i0 < NPLOT ; i0++ ) {
  73       for( i1 = 0 ; i1 < NCALC ; i1++ ) {
  74          rangcel = ascale * sin( 2.0*rang + rav*dt ) ; // avg ang accel
  75 	 rang   += (rav + 0.5*rangcel*dt)*dt ;
  76 	 rav    += rangcel * dt ;
  77       }
  78       
  79       double eangle  = 360.0*t/OTIME;
  80       double sattemp = temp0*sqrt(sqrt(cos(rang))) ;
  81 
  82       //       Time,  Eangle,  Angle, Velocity, Temperature
  83       printf( "%10.4f%10.4f%10.4f%12.3e%10.4f\n",
  84                    t,  eangle, rang*r2d, rav*r2d, sattemp );
  85       t += dtplot ;
  86    }
  87 
  88    double   ang1 = r2d*rang ;
  89 
  90    printf( "# %12.9f starting angular velocity to omega ratio\n", 1.0+rav0/omega );
  91    
  92    printf( "# angular velocity%14.6e =start %14.6e =end%14.6e =diff\n",
  93                             rav0,       rav,      rav-rav0        );
  94  
  95    printf( "# end angle       %14.6e =expect%14.6e =end%14.6e =diff\n",
  96                             ANG1,       ang1,     ang1-ANG1       );
  97 }