// in01.c  incline
//
// gcc -o in01 in01.c -lm

#define  NAME   "in01"
#define  CMD    "in"          //          tmp.dat -> tmp.png

#define  SANGLE     20.0      // degrees  Ground start angle
#define  ZTOP    80000.0      // m        Top altitude
#define  DX          0.1      // m        Increment
#define  MS         10.0      // kg       Sheath mass
#define  MRTOP       3.0      // kg       Rotor mass
#define  VRTOP   14000.0      // m/s      Rotor speed at top  (reference)
#define  VRBOT   14055.2537   // km/s     Rotor speed at surface

// cables
#define  Z0      20000.0      // m        Cable "start" height
#define  HC      80000.0      // m        Support height
#define  ANGLE      45.0      //          Side Angle

#define  RE    6378000.0      // m        Earth Radius
#define  AG0         9.81     // m/s2     Gravity at surface

#define  PI 3.14159265 

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

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

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

   int     points= 0                    ;
   double  dzdx0 = tan( PI*SANGLE/180.0);
   double  cang  = cos( PI/4.0 )        ; // 45 degrees cable angle
   double  ag    = AG0                  ; // gravity
   double  ms    = MS                   ;
   double  dx    = DX                   ;
   double  dz                           ;
   double  x     = 0.0                  ;
   double  z     = 0.0                  ;
   double  kx    = 0.0                  ;
   double  kz    = 0.0                  ;
   double  vr    = VRBOT                ;
   double  dzdx  = dzdx0                ;
   double  dzdx1                        ;
   double  dzdx2                        ;
   double  mrbot = MRTOP*VRTOP/vr       ;
   double  mr    = mrbot                ;
   double  tmp                          ;
   double  gz                           ;
   double  angle                        ; // force angle radians
   double  dangle                       ; // force angle degrees
   double  force                        ; // force kiloNewtons
   double  forcex                       ; // horizontal force kiloNewtons
   double  forcez                       ; // vertical   force kiloNewtons
   double  forcec                       ; // cable      force kiloNewtons
   double  forces                       ; // station    force kiloNewtons
   double  smass                        ; // station mass metric tons
   double  cmass                        ; // main cable mass metric tons
   double  cinc = 0.0                   ; // incline cable mass

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

   datfile = fopen( "tmp.dat", "wb" );

   while( z < ZTOP ) {
      tmp   = 1.0/( 1.0 + z/RE )       ;
      ag    = AG0*tmp*tmp              ;
      mr    = MRTOP*VRTOP/vr           ;

      gz    = sqrt(exp(2.0*(z+Z0)/HC)-1.0)/cang ;

      dzdx1 = 1.0 + dzdx*dzdx          ;
      dzdx2 = dzdx1/(RE+z)
            - ag*(ms+mr)*dzdx1*dzdx1
	    / (mr*vr*vr*(1.0-gz*dzdx)) ;
      dzdx += dzdx2 *dx                ;
      dz    = dzdx * dx                ;
      vr   -= ag*dz/vr                 ; 
      z    += dz                       ;
      x    += dx                       ;
      cinc += (dz*ms/cang)
            * (  sqrt(exp(2.0*(z+Z0)/HC)-1.0)
               - sqrt(exp(2.0*Z0/HC)-1.0) ) ; // cable mass kg

      kx    = 0.001*x                  ;
      kz    = 0.001*z                  ;
     

      if( ((points++)%5000) == 0 ) {
         fprintf( datfile, "%14.4f%14.4f\n", kx, kz ) ;
      }
   }
   fprintf( datfile, "%14.4f%14.4f\n", kx, kz ) ;
   fclose( datfile );

   printf( "--------------------------------------------------\n"     );

   angle  = atan( dzdx0 )                                 ; // radians
   dangle = (180.0/PI)*angle                              ; // degrees
   force  = 2*mrbot*VRBOT*VRBOT*sin(0.5*angle)/1000.0     ; // kiloNewtons
   forcez = force * cos( 0.5*angle )                      ; //
   forcex = force * sin( 0.5*angle )                      ; //
    
   printf( "             bottom slope =%11.4f vert/horiz\n" ,  dzdx0 );
   printf( "             bottom angle =%11.4f degrees\n"    , dangle );
   printf( "             bottom force =%11.0f kiloNewtons\n",  force );
   printf( "    bottom force vertical =%11.0f kiloNewtons\n", forcez );
   printf( "  bottom force horizontal =%11.0f kiloNewtons\n", forcex );
   printf( "      bottom mass density =%11.4f kg/meter\n"   ,  mrbot );
   printf( "          bottom velocity =%11.4f meter/sec\n"  ,  VRBOT );
   printf( "           bottom gravity =%11.4f meter^2/sec\n",    AG0 );
   printf( "--------------------------------------------------\n"    );

   angle  = atan( dzdx )                                  ; // radians
   dangle = (180.0/PI)*angle                              ; // degrees
   force  = 2*mr*vr*vr*sin(0.5*angle)/1000.0              ; // kiloNewtons
   forcez = force * cos( 0.5*angle )                      ; //
   forcex = force * sin( 0.5*angle )                      ; //
   forcec = forcex * sqrt( 1.0 + gz*gz )                  ; // cable force kN
   forces = forcez-gz*forcex                              ; // station support
   smass  = forces/ag                                     ; // station mass tons
   cmass  = (forcex/(ag*cang))
          * (  sqrt(exp(2.0*(ZTOP+Z0)/HC)-1.0)
             - sqrt(exp(2.0*Z0/HC)-1.0) )                 ; // cable mass tons
   cinc  /= 1000.0                                        ; // incline cable tons

   printf( "      horizontal distance =%11.4f kilometers\n" ,     kx );
   printf( "        vertical distance =%11.4f kilometers\n" ,     kz );
   printf( "                top slope =%11.4f vert/horiz\n" ,   dzdx );
   printf( "                top angle =%11.4f degrees\n"    , dangle );
   printf( "                top force =%11.0f kiloNewtons\n",  force );
   printf( "       top force vertical =%11.0f kiloNewtons\n", forcez );
   printf( "     top force horizontal =%11.0f kiloNewtons\n", forcex );
   printf( "         top mass density =%11.4f kg/meter\n"   ,     mr );
   printf( "             top velocity =%11.4f meter/sec\n"  ,     vr );
   printf( "              top gravity =%11.4f meter^2/sec\n",     ag );
   printf( "--------------------------------------------------\n"    );
   printf( "     cable force ratio gz =%11.4f vert/horiz\n" ,     gz );
   printf( "    top cable axial force =%11.4f kiloNewtons\n", forcec );
   printf( "    station support force =%11.4f kiloNewtons\n", forces );
   printf( "  station + elevator mass =%11.4f metric tons\n", smass  );
   printf( " main diagonal cable mass =%11.4f metric tons\n", cmass  );
   printf( "       incline cable mass =%11.4f metric tons\n", cinc   );
   printf( "--------------------------------------------------\n"    );

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

   sprintf( filename, "%s.png", NAME );
   rename( "tmp.png", filename );

   return(0);
}