// exp01.c
//   > make 
//   printd exp01.dat

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nrlmsise-00.h"

// #define  TEST           1

#define  A0         120.0      // km     bottom altitude simulated
#define  A1         600.0      // km     top altitude simulated
#define  ASP        100.0      // km     space altitude

#define  RGEO  42164170.0      // m      geosynchronous radius
#define  RE     6378137.0      // m      equatorial radius
#define  MU    3.98600448e14   // m3/s2  gravitational parameter
#define  AM           5.0      // m2/kg  thinsat sail ratio
#define  AE           2.0      // m2/kg  experiment sail ratio, edge on?

#define  KM        1000.0      // m/km
#define  HOUR      3600.0      // seconds per hour
#define  DAY      86400.0      // seconds per day

int main() {
   struct nrlmsise_output output;
   struct nrlmsise_input input;
   struct nrlmsise_flags flags;
   struct ap_array aph;
   int i;

   double pi = 4.0*atan(1.0);  // compute pi

   // KHL loop
   // time of day 
   double hour ;
   double h0    =  0.0   ;
   double h1    = 24.0   ;
   double hs    = 0.25   ;
   double havg  = hs/(h1-h0) ;

   // altitude above equator
   double alt            ;
   double a0    =     A0 ;
   double a1    =     A1 ;
   double as    =    1.0 ;
   double dalt  =   0.01 ;

   // f07 
   double f107[3] ;
   f107[0] = 70.0 ;  f107[1] = 120.0 ;  f107[2] = 200.0 ;
   int    fi   ;

   double denm ;  //  density dalt km down
   double dens ;
   double denp ;  //  density dalt km up
// double temp ;
   double laps ;
   double tdec[3] ;

   /* input values */
   for (i=0;i<7;i++)
      aph.a[i]=100;

   for (i=0;i<24;i++)
      flags.switches[i]=1;

   input.doy  =   79  ;   // March 21 vernal equinox 
   input.year   = 0   ;   // without effect
   input.g_lat  = 0.0 ;
   input.g_long = 0.0 ;
   input.ap     = 4.0 ;   // magnetic index 

   // header

   printf( "#f107" );
   for( i=0 ; i < 3 ; i++) printf( "%6.0f%9.0f%9.0f ", f107[i],f107[i],f107[i] ) ; 
   printf( "\n#alt");
   for( i=0 ; i < 3 ; i++) printf( "  dens.   exper.    decay" ) ;
   printf(  "\n"   );

   // altitude loop

   tdec[0] = 0.0 ;      //  decay time to atmosphere 
   tdec[1] = 0.0 ;      //  decay time to atmosphere (100   km )
   tdec[2] = 0.0 ;      //  decay time to atmosphere (100   km )

   for( alt =a0 ; alt  <= a1 ; alt += as ) {
      printf( "%4.0f", alt );

      // decay from geo -------------
      // starting conditions experimental orbit

      double rp   = RE + alt * KM   ; // m perigee radius
      double ra   = RGEO            ; // m apogee  radius
      double sem  = 0.5*(ra+rp)     ; // m semimajor axis
      double ecc  = 0.5*(ra-rp)/sem ; //   eccentricity
      double v0   = sqrt( MU*((0.5/rp)+(0.5/ra))); // m/s 
      double vp   = (1.0+ecc)*v0    ; // m/s perigee velocity
      double va   = (1.0-ecc)*v0    ; // m/s apogee  velocity

#ifdef  TEST
      printf( "\n");
      printf( "rp  =%14.4e\n", rp  );
      printf( "ra  =%14.4e\n", ra  );
      printf( "sem =%14.4e\n", sem );
      printf( "ecc =%14.4e\n", ecc );
      printf( "v0  =%14.4e\n", v0  );
      printf( "vp  =%14.4e\n", vp  );
      printf( "va  =%14.4e\n", va  );
#endif

      for( fi  = 0 ; fi   <=  2 ; fi++ ) {
         input.f107A=f107[fi];
         input.f107 =f107[fi];

         // prepare to take average
         denm  = 0.0 ;
         denp  = 0.0 ;
         dens  = 0.0 ;

         for( hour=h0 ; hour < h1  ; hour += hs ) {
            input.sec = hour * HOUR ;
            input.lst = hour        ;

            input.alt = alt-dalt;
            gtd7(&input, &flags, &output);
            denm += output.d[5] ;

            input.alt = alt+dalt;
            gtd7(&input, &flags, &output);
            denp += output.d[5] ;

            input.alt = alt;
            gtd7(&input, &flags, &output);
            dens += output.d[5] ;
         }

         denm *= havg ;   // normalize average
         denp *= havg ;
         dens *= havg ;

         laps  = 2.0*dalt /(log(denm/denp)); // compute lapse rate, km

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

         // Compute decay times
         
         double texp = 0.0 ;      //  experiment time from geo (42164 km )
         double H    = laps*KM ;  // lapse rate in meters
         double vpt  = sqrt( MU*((2.0/rp)-(1.0/(rp+H)))); // m/s vp thresh

#ifdef  TEST
         printf( "\n\n");
         printf( "dens=%14.4e\n", dens);
	 printf( "H   =%14.4e\n", H   );
	 printf( "vpt =%14.4e\n", vpt );
#endif
         double   vpx = vp  ;                       ; // experimental vp

         // experimental decay loop 
         while( vpx > vpt ) {
            double rvv = rp*vpx*vpx ;
            double sq2 = rp / ( 2.0*MU - rvv ) ;
            double dte = 2.0*pi*MU*sq2*sqrt(sq2)    ; // orbit time
            texp      += dte                        ; // add orbit time
            double dvp = dens*AE*vpx*sqrt((pi*rp*H)/(1.0-MU/rvv));
            vpx       -= dvp                        ; // slow down
#ifdef  TEST
            printf( "dt=%14.4e t=%14.4e dv=%14.4e vpx=%14.4e\n",
                     dte, texp, dvp, vpx );
#endif
         } // integration loop

         // circular decay, tumbling --------------
	 // dr/dt = dens * AM/2 * sqrt( mu rp ) 
	 // 2.0 <- half the area: disk divided by two faces of sphere

         tdec[fi] += 2.0 * KM * as / ( dens * AM * sqrt( MU * rp ) ) ;
         printf( "%7.0e%9.2e%9.2e", dens, texp/DAY, tdec[fi]/DAY );

      } // F107 loop
      printf( "\n" );
   }
   exit(0);
}