attachment:sge1.c of RwsSpaceGrainElevator

   1 // sge1.c space grain elevator v1
   2 // cc -o sge1 sge1.c -lm
   3 // takes about 5 minutes to run on a 2.2 GHz Intel Core Duo, 32 bit mode
   4 
   5 #define  IREDUCE   0.953639070     // reduced angular velocity
   6 #define  RSTEP     200.0           // m  integration step
   7 
   8 #define  AU        149597870700.0  // m   astronomical unit
   9 #define  C         299792458.0     // m/s speed of light
  10 #define  MU        3.986400414e14  // m3s2 earth grav. param., unused
  11 #define  SDAY      86164.09        // seconds
  12 #define  NGRAIN    3.0e7           // grains per meter
  13 #define  DGRAIN    700.0           // kg/m3 density of grain, initial guess
  14 #define  GEO       4.2164e7        // m  GEO radius
  15 #define  RPRINT    1.0E8           // print step
  16 
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 
  21 //-------------------------------------------------------------------
  22 // pressure in Pascal, density in kg/m^3
  23 
  24 double dens( double pres ) {
  25    double pr = pres/1e12 ;
  26    if( pr < 1.0 ) {  return 700.0 + 1700.0*pr ; }
  27    else { return  2363.0 + 37.0*pow( pr, 0.6 ); }
  28 }
  29 
  30 //-------------------------------------------------------------------
  31 int main() {
  32    double  pi2    = 8.0*atan(1.0)     ; 
  33    double  area   = 0.004             ; // m2 cross section area
  34    double  omega0 = pi2 / SDAY        ; // rad/s  initial angular velocity
  35    double  iearth = 8.034e37          ; // kg-m2 earth moment of inertia
  36    double  ngrain = pow( 2.0, 63.0 )  ; // number of grains of wheat
  37    double  vgrain = ngrain / NGRAIN   ; // m3 volume of grain, zero pressure
  38    double  mgrain = vgrain * DGRAIN   ; // kg mass of grain
  39    double  max    = mgrain / area     ; // terminate loop
  40 
  41    double  rprint = 0.0               ; // print step
  42 
  43    // initial guess
  44    double  omega  = IREDUCE*omega0    ; // rad/s angular velocity, tweaked
  45    double  om2    = omega*omega       ; // squared angular velocity
  46    double  iela   = 0.0               ; // kg elevator moment of inertia over area
  47    double  r      = GEO               ; // meters radius
  48    double  p      = 0.0               ; // Pa pressure
  49    double  ma     = 0.0               ; // kg/m2 column density
  50    double  dma                        ; // kg/m2 mass step
  51 
  52    // core calculation loop
  53    while( ma < max ) {
  54 
  55       if( r >= rprint ) {
  56          printf( "%14.6e"  , r       );
  57          printf( "%10.3e"  , p       );
  58          printf( "%10.3e"  , dens(p) );
  59 	 printf( "%10.6f\r", ma/max  );
  60 	 rprint = r + RPRINT ;
  61       }
  62 
  63       dma    = RSTEP*dens(p)          ; // kg/m2 additional mass
  64       ma    += dma                    ; // kg/m2 add to total mass
  65       iela  += dma * r*r              ; // kg    add to moment of inertia
  66       p     += dma * r * om2          ; // Pa    add to pressure
  67       r     += RSTEP                  ; // m     add to distance
  68    }
  69 
  70    // results 
  71    double v    = r * omega            ; // m/s  rotational velocity at end
  72    double m    = ma * area            ; // kg   total mass
  73    double d    = dens(p)              ; // kg/m3 density at end
  74    double iel  = iela * area          ; // kg-m2 elevator moment of inertia
  75    double ifac = iel/iearth           ; // fraction of earth moment of inertia
  76    double inew = 1.0/(1.0+ifac)       ; // reduced angular velocity
  77    
  78    printf( "%16.9f  Reduced angular velocity factor\n", inew );
  79    printf( "%16.9f  Elevator fraction of earth moment of inertia\n", ifac );
  80    printf( "%16.6e  meters radius of elevator\n", r );
  81    printf( "%16.6f    AU radius of elevator\n", r/AU );
  82    printf( "%16.6e  kg/m3 density of grain at end\n", d );
  83    printf( "%16.6e  Pa pressure on grain at end\n", p );
  84    printf( "%16.6e  m/s velocity of end\n", v );
  85    printf( "%16.6f    fraction of speed of light\n", v/C );
  86    printf( "%16.6e  kg mass of elevator\n", m );
  87    printf( "%16.6e    kg target mass of elevator\n", mgrain );
  88    printf( "%16.2f    m integration step\n", RSTEP );
  89 
  90    return 0 ;
  91 }