// sge1.c space grain elevator v1
// cc -o sge1 sge1.c -lm
// takes about 5 minutes to run on a 2.2 GHz Intel Core Duo, 32 bit mode

#define  IREDUCE   0.953639070     // reduced angular velocity
#define  RSTEP     200.0           // m  integration step

#define  AU        149597870700.0  // m   astronomical unit
#define  C         299792458.0     // m/s speed of light
#define  MU        3.986400414e14  // m3s2 earth grav. param., unused
#define  SDAY      86164.09        // seconds
#define  NGRAIN    3.0e7           // grains per meter
#define  DGRAIN    700.0           // kg/m3 density of grain, initial guess
#define  GEO       4.2164e7        // m  GEO radius
#define  RPRINT    1.0E8           // print step

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

//-------------------------------------------------------------------
// pressure in Pascal, density in kg/m^3

double dens( double pres ) {
   double pr = pres/1e12 ;
   if( pr < 1.0 ) {  return 700.0 + 1700.0*pr ; }
   else { return  2363.0 + 37.0*pow( pr, 0.6 ); }
}

//-------------------------------------------------------------------
int main() {
   double  pi2    = 8.0*atan(1.0)     ; 
   double  area   = 0.004             ; // m2 cross section area
   double  omega0 = pi2 / SDAY        ; // rad/s  initial angular velocity
   double  iearth = 8.034e37          ; // kg-m2 earth moment of inertia
   double  ngrain = pow( 2.0, 63.0 )  ; // number of grains of wheat
   double  vgrain = ngrain / NGRAIN   ; // m3 volume of grain, zero pressure
   double  mgrain = vgrain * DGRAIN   ; // kg mass of grain
   double  max    = mgrain / area     ; // terminate loop

   double  rprint = 0.0               ; // print step

   // initial guess
   double  omega  = IREDUCE*omega0    ; // rad/s angular velocity, tweaked
   double  om2    = omega*omega       ; // squared angular velocity
   double  iela   = 0.0               ; // kg elevator moment of inertia over area
   double  r      = GEO               ; // meters radius
   double  p      = 0.0               ; // Pa pressure
   double  ma     = 0.0               ; // kg/m2 column density
   double  dma                        ; // kg/m2 mass step

   // core calculation loop
   while( ma < max ) {

      if( r >= rprint ) {
         printf( "%14.6e"  , r       );
         printf( "%10.3e"  , p       );
         printf( "%10.3e"  , dens(p) );
	 printf( "%10.6f\r", ma/max  );
	 rprint = r + RPRINT ;
      }

      dma    = RSTEP*dens(p)          ; // kg/m2 additional mass
      ma    += dma                    ; // kg/m2 add to total mass
      iela  += dma * r*r              ; // kg    add to moment of inertia
      p     += dma * r * om2          ; // Pa    add to pressure
      r     += RSTEP                  ; // m     add to distance
   }

   // results 
   double v    = r * omega            ; // m/s  rotational velocity at end
   double m    = ma * area            ; // kg   total mass
   double d    = dens(p)              ; // kg/m3 density at end
   double iel  = iela * area          ; // kg-m2 elevator moment of inertia
   double ifac = iel/iearth           ; // fraction of earth moment of inertia
   double inew = 1.0/(1.0+ifac)       ; // reduced angular velocity
   
   printf( "%16.9f  Reduced angular velocity factor\n", inew );
   printf( "%16.9f  Elevator fraction of earth moment of inertia\n", ifac );
   printf( "%16.6e  meters radius of elevator\n", r );
   printf( "%16.6f    AU radius of elevator\n", r/AU );
   printf( "%16.6e  kg/m3 density of grain at end\n", d );
   printf( "%16.6e  Pa pressure on grain at end\n", p );
   printf( "%16.6e  m/s velocity of end\n", v );
   printf( "%16.6f    fraction of speed of light\n", v/C );
   printf( "%16.6e  kg mass of elevator\n", m );
   printf( "%16.6e    kg target mass of elevator\n", mgrain );
   printf( "%16.2f    m integration step\n", RSTEP );

   return 0 ;
}