// gcc -o irn1 irn1.c -lm ; ./irn1 > irn1.dat ; gp irn1.gp
// 2013/08/09
// infrared filtered thinsat,  night orientation and temperature
//   circular orbit approximation, no light pressure eccentricity 

#define VMULT    1.16625404993432468 //          Multiply starting velocity

#define EFILT    0.3542              // eV       high pass cutoff energy
#define ETEMP    260.0               // Kelvin   Earth temperature

//  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      Planck's constant
#define C        299792458.0         // m/s      speed of light 
#define Q        1.602176565e-19     // J/eV     Joules per electron volt
#define RE       6378.0              // km       Average radius of earth
#define RSAT     12789.0             // km       Average distance to earth
#define ANG0    -60.0                // degrees  Starting earth relative angle
#define ANG1     60.0                // degrees  Ending earth relative angle
#define OANG0    150.0               // degrees  Starting orbit angle from noon
#define OANG1    210.0               // degrees  Ending orbit angle
#define NPLOT    401                 //          number of plot points
#define NCALC    10000               //          number of calculation points
#define OTIME    14400.0             // seconds  Orbit period

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

double  pi , pi2, pih, pis, d2r, r2d ;

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

   double dplot1   = (double) (NPLOT-1)     ; // 
   double dcalc    = (double) NCALC         ; // 
   double dintrvl  = dplot1*dcalc           ; // number of calculation intervals
   double angchng  = ANG1-ANG0              ; // angular change in degrees
   double rangchng = d2r*angchng            ; // angular change in radians
   double dang     = angchng/dintrvl        ; // plot calculation step
   double t0       = OTIME*OANG0/360.0      ; // start time in seconds
   double t1       = OTIME*OANG1/360.0      ; // end time in seconds
   double tchng    = t1-t0                  ; // simulation time, seconds
   double dt       = tchng/dintrvl          ; // time step, seconds
   double dtplot   = tchng/dplot1           ; // plot time step, seconds
   double ktearth  = KB*ETEMP               ; // earth IR temperature
   double lambda   = 1E6*C*H / (Q*EFILT)    ; // filter wavelength in micrometers
   double omega    = pi2/OTIME              ; // orbit angular frequency
   double rav0     = VMULT*rangchng/tchng   ; // starting angular velocity
   double ascale   = 1.5*omega*omega        ; // scale for angular acceleration
   double rang0    = d2r * ANG0             ; //

   double rhoe     = asin( RE / RSAT )      ; // earth apparent angle
   double sang0    = 1.0-cos(2.0*rhoe)      ; // earth solid angle
   double temp0    = ETEMP*sqrt(sqrt(0.5*sang0)) ; // temperature scaling

   int i0, i1 ;

   // orbit loop 
   
   double t        = t0                     ; // time in seconds
   double rang     = rang0                  ; // thinsat angle in radians
   double rav      = rav0                   ; // angular velocity rad/sec
   double rangcel                           ; // angular acceleration

   for( i0 = 0 ; i0 < NPLOT ; i0++ ) {
      for( i1 = 0 ; i1 < NCALC ; i1++ ) {
         rangcel = ascale * sin( 2.0*rang + rav*dt ) ; // avg ang accel
	 rang   += (rav + 0.5*rangcel*dt)*dt ;
	 rav    += rangcel * dt ;
      }
      
      double eangle  = 360.0*t/OTIME;
      double sattemp = temp0*sqrt(sqrt(cos(rang))) ;

      //       Time,  Eangle,  Angle, Velocity, Temperature
      printf( "%10.4f%10.4f%10.4f%12.3e%10.4f\n",
                   t,  eangle, rang*r2d, rav*r2d, sattemp );
      t += dtplot ;
   }

   double   ang1 = r2d*rang ;

   printf( "# %12.9f starting angular velocity to omega ratio\n", 1.0+rav0/omega );
   
   printf( "# angular velocity%14.6e =start %14.6e =end%14.6e =diff\n",
                            rav0,       rav,      rav-rav0        );
 
   printf( "# end angle       %14.6e =expect%14.6e =end%14.6e =diff\n",
                            ANG1,       ang1,     ang1-ANG1       );
}