// hc04.c
// compute heating and cooling 
//
// gcc -o hc04 hc04.c -lm

#define  NAME   "hc04"
#define  CMD    "hc"

#define  RATE       12   //          Flash frame rate ( 60x )
#define  TIME     3600   //          Total time
#define  NPLOT     720   //          Number of plots
#define  NCOMP      70   //          Computation passes per plot
#define  PPLOT    2400   //          Plot width (may be ugly!)
#define  OPLOT    -200   //

#define  TEMP0   315.0   // K        Starting temperature
#define  TEMPA   250.0   // K        Ambient temperature
#define  PV    10800.0   // m/s      Payload maximum velocity
#define  PA       30.0   // m/s2     Payload acceleration
#define  PT      360.0   // s        Payload acceleration time
#define  PM     5000.0   // kg       Payload mass
#define  PL     1944.0   // km       Payload acceleration distance

#define  PS       24.0   // s        Payload spacing in time
#define  PN         15   //          Number of computed payloads

#define  RL       5180   // km       Rotor Length
#define  RM        3.0   // kg/m     Rotor mass/length
#define  RR      0.025   // m        Rotor radius in meters
#define  RV    14000.0   // m/s      Rotor velocity

#define  EM        0.8   //          Emissivity
#define  SB    5.67E-8   // W/m2K4   Blackbody radiation
#define  HC        600   // J/kg-K   Rotor heat capacity
#define  PI 3.14159265 

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

double  rtemp[RL]                  ; // rotor temperature K

// modulus function
int rmax = RL * ( 500000000 / RL ) ; // large offset multiple of RL
int rmod( int arg ) { return ( ( rmax + arg ) % RL ) ; }

//-----------------------------------------------------------------
int main() {

   FILE    *datfile                     ; //
   char    rmmkdir[200]                 ; //
   char    png2swf[200]                 ; //
   char    gnuplot[200]                 ; //
   char    filename[80]                 ; //

   double  ptime                        ; // acceleration time of payload
   double  vp                           ; // velocity of payload
   double  xp                           ; // position of payload

   int     ipoint                       ; // point index for plotting, cooling
   int     index = 0                    ; // movement of rotor since t0
   int     iplot                        ; // plot number
   int     icomp                        ; // computation index
   int     ipay                         ; // payload index
   double  rtp                          ; // temporary rtemp
   double  time                         ; // simulating time
   double  tdelta = 1000.0/RV           ; // timestep in seconds
   double  xpoint                       ; //

   double  dtp = tdelta                   // payload heating factor 
               * PA * PM                  // payload force
               * 0.001                    // kilometer of mass
               / ( RM * HC )            ; // rotor heat capacity  J/K

   double  heat                         ; // added heat
   int     ip                           ; // integer position of payload
   int     ip0                          ; // payload position 0,  rotor-coords
   int     ip1                          ; // payload position 1,  rotor-coords
   double  xpi                          ; // integer part of payload position
   double  xpf                          ; // fractional part of position

   double  dt0 = tdelta                   // rotor cooling factor
               * EM * SB                  // blackbody radiation
               * 2.0 * PI * RR            // surface area
               / ( RM * HC )            ; // rotor heat capacity  J/K

   double  dt1 = dt0                      // rotor ambient 
               * TEMPA*TEMPA              //
               * TEMPA*TEMPA            ; // 

   // --------------------------------------------------------------------

   printf ( "tdelta = %14.4e\n", tdelta );
   printf ( "dt0    = %14.4e\n", dt0    );
   printf ( "dt1    = %14.4e\n", dt1    );
   printf ( "dtp    = %14.4e\n", dtp    );

   sprintf( rmmkdir, "rm -rf %sdir ; mkdir %sdir", NAME, NAME );
   system( rmmkdir )           ; // make directory

   // initialize
   for( ipoint=0 ; ipoint < RL ; ipoint++ ) { rtemp[ipoint] = TEMP0 ; }

   // plot loop
   for( iplot=0 ; iplot < NPLOT ; iplot++ ) { 

      printf( "plot %4d\r", iplot ); fflush( stdout );

      // compute loop per plot
      for( icomp=0 ; icomp < NCOMP ; icomp++ ) {

         index = iplot*NCOMP + icomp ;
         time  = tdelta * ((double) index );

         //  cooling
         for( ipoint=0 ; ipoint < RL ; ipoint++ ) {
            rtp  = rtemp[ ipoint ] ;
            rtp *= rtp ;
            rtemp[ ipoint ] -= dt0*rtp*rtp - dt1 ;
         }

         // heating by payloads - start a little later
         for( ipay = 1 ; ipay <= PN ; ipay++ ) {
            ptime = time - PS * ((double)ipay) ;
            if( ( ptime >= 0.0 ) && ( ptime <= PT ) ) {

               vp   = PA * ptime           ;
               heat = dtp * ( RV - vp )    ;

               xp   = 0.0005 * vp * ptime  ;
               xpi  = floor( xp )          ;
               ip   = (int) xpi            ;
               xpf  = xp-xpi               ;
               ip0  = rmod( ip   - index ) ;
               ip1  = rmod( ip+1 - index ) ;

               rtemp[ ip0 ] += heat*(1.0-xpf) ;
               rtemp[ ip1 ] += heat*(    xpf) ;
            }
         }
      } // end of icomp compute loop 

      index   = (iplot+1)*NCOMP ;
      datfile = fopen( "tmp.dat", "wb" );

      for( ipoint = 0 ; ipoint < PPLOT ; ipoint++ ) {
         xpoint = (double)( OPLOT+ipoint) ;  // starts before 0
         ip     = rmod( OPLOT+ipoint-index );
         fprintf( datfile, "%16.4f%16.4f\n", xpoint, rtemp[ ip ] ) ;
      }
      fclose( datfile );

      sprintf( gnuplot, "/usr/local/bin/gp %s.cmd", CMD );
      system( gnuplot );

      sprintf( filename, "%sdir/b%04d.png", NAME, iplot );
      rename( "tmp.png", filename );
   }

   // make flash video
   sprintf( png2swf, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME );
   system( png2swf );

   return(0);
}