/*
 * RASP89 - Rocket Altitude Simulation Program
 * Based on the simple BASIC program written by
 * Harry Stine in 1979, modified and converted
 * into C by Mark A. Storin, Sept/Oct 1989.
 * Rev 2 by Kent Hoult, June 1990.
 *
 * 10-20-89  modified to allow multiple stages up to 3.
 * 07-01-90  modified to use linear segment engine thrust tables from data file
 *           smaller time increments, and crude supersonic transition drag.
 * 6/7/93    Modified to support staging delays.
 *           Made changes to the supersonic transition drag equations.
 *           Migrated lots of global variables to local (with meaninful names)
 *           Rev #3  Stu Barrett
 * 12/29/93  Added support for VAX/VMS - Rusty Whitman
 * 02/15/94  Decrease air density with altitude/temp - added air_density()
 *           Make output gnuplot compatible 
 *           Rev 3.2  - Mark C. Spiegl (spiegl@rtsg.mot.com) 
 * 02/18/94  Added MAC support
 *           Revised GNUplot support
 *           Added input error checking and input defaults 
 *           Rev 3.3 - Warren Massey (massey@travis.llnl.gov)
 */

/* #define MSDOS  1 */
/* #define UNIX   1 */
/* #define VMS    1 */
#define MAC    1 

#include <stdio.h>
#include <ctype.h>

#define VERSION   "3.3" 
#define TRUE      1
#define FALSE     0
#define PI        3.1415927
#define G         9.80621
#define IN2M      0.0254
#define GM2LB     0.00220462
#define OZ2KG     0.028349523
#define M2FT      3.280840
#define MACH0_8   265.168     /* Mach 0.8 in Meters/sec  870 ft/sec */
#define MACH1     331.460     /* Mach 1.0 in Meters/sec  1086 ft/sec*/
#define MACH1_2   397.752     /* Mach 1.2 in Meters/sec  1305 ft/sec */
#define TTYmax    20          /* Max length of an input string from stdin */
#define MCmax     12          /* Max length of engine designation string, e.g. G80 */
#define MAXT      64          /* Max array size used to hold time/thrust 
                                 datapairs read from RASP.ENG file. A datapoint
                                 consists of time & thrust datums so this array
                                 can hold MAXT/2 datapoints. */

#define ROD        1.2192     /* Length of launch rod in meters (48") */
#define MINENG     1          /* Minimum number of engines per stage */
#define MAXENG    12          /* Maximum number of engines per stage */
#define MINSTAGE   1          /* Minimum number of stages allowed */
#define MAXSTAGE   5          /* Maximum number of stages allowed */
#define MINWGHT    0.001      /* Minimum weight per stage in ounces */
#define MAXWGHT 1600.0        /* Maximum weight per stage in ounces */
#define MINDIA     0.1        /* Minimum diameter of stage in inches */
#define MAXDIA    24.0        /* Maximum diameter of stage in inches */
#define MAXDELAY  60.0        /* Maximum interstage delay in seconds */
#define MINDELAY   0.0        /* Minimum interstage delay in seconds */
#define MAXDRAGC 100.0        /* Maximum drag coefficient */
#define MINDRAGC   0.0        /* Minimum drag coefficient */
#define MAXBASET 120.0        /* Maximum base temperature (Deg F) */
#define MINBASET   0.0        /* Minimum base temperature (Deg F) */

#define GNUcmd  1              /* Enable production of GNUplot command file */

/*
 * Global Variables
 */
int verbose  = TRUE;
int GNUplot  = FALSE;
int Error    = FALSE;

double rocketwt;
double drag_coff;
double base_temp = 72.0;	/* Initial default */
double base_temp_K;

int
  stagenum = 1,     /* number of stages */
  nexteng,          /* next free engine index */
  destnum = 1,      /* index into array of devices for output - dest */
  ianswer,
  minval,
  maxval;
  
double danswer;

char noGNUprefix = ' ';
char   GNUprefix = '#';
char prefix_char;
char Bell = '\7';
char Mcode[MCmax];   /* array to hold engine code */
char fname[64];      /* output file name for simulation */
#if GNUcmd
char fnam2[64];      /* output file name for GNUplot commands */
#endif
char fnam_prefix[64];
char sanswer[64];
char tty_input[TTYmax];

char GNUplot_Dat_Ext[] = {".dat"};    /* file name extension for GNUplot data file */
char GNUplot_Cmd_Ext[] = {".cmd"};    /* file name extension for GNUplot command file */

FILE *stream,  *efile, *fopen();
#if GNUcmd
FILE *stream2;
#endif

#ifdef MAC
#include <console.h>
#include <fonts.h>
FILE *fopenc();
int MsgSent = FALSE;
#endif

char efile_name[] = {"rasp.eng"}; /* name of engine database file */

struct engine {
    double t2;           /* thrust duration    (seconds)                 */
    double m2;           /* propellant wt.     (kilograms)               */
    double wt;           /* initial engine wt. (kilograms)               */
    double thrust[MAXT]; /* thrust curve     (Time-sec,Thrust-Newtons)   */
    double navg;         /* average thrust     (newtons)                 */
    double ntot;         /* total impulse      (newton-seconds)          */
    double npeak;        /* peak thrust        (newtons)                 */
    int    dia;          /* engine diameter    (millimeters)             */
    int    len;          /* engine length      (millimeters)             */
    int    delay[8];     /* ejection delays available (sec)              */
    char   mfg[32];      /* manufacturer info                            */
    char   Mcode[MCmax]; /* engine name                                  */
   };

struct stage_info {
    int engcnum;         /* index into e_info */
    int numeng;          /* # of engines in stage */
    double start_burn;   /* Start of engine (includes previous stage) */
    double end_burn;     /* End of engine burn  (includes previous stage) */
    double end_stage;    /* End of stage (includes previous stage) */
    double weight_Oz;    /* stage wt w/o engine (oz)*/
    double weight_Kg;    /* stage wt w/o engine (Kg)*/
    double totalw;       /* wt of stage w/ engine(s) */
    double maxdia;       /* max diameter of stage  */
    double delay;        /* interstage delay */
    char Mcode[MCmax];   /* array to hold engine code */
   };


    
struct stage_info stages[MAXSTAGE];
struct engine     e_info[MAXENG];


int *get_input( input, input_max )
    int input_max;
    char *input;
{
    int i, c;
    c = getc( stdin );
    if( (c != '\n') && (c != EOF) ) {
        ungetc( c, stdin );
        fgets( input, input_max, stdin );
        sscanf( input, "%s%*c", input );
    }
    else {
        strcpy( input, "") ;
    }
}        


/*   "dump_buffer()" must be called after every "scanf()" to avoid screwing up
 *   "get_input()"
 */   
int dump_buffer() {
  for (;;)
     switch( getchar() ) {
     case EOF: return EOF;
     case '\n': return 0;
     default: ;
  }
}


print_engine_info(e)
     int e;
{
  int d;
  fprintf(stream,"%c %-8s  %-9s  %8.2lf   %8.2lf  %8.2lf  %6.2lf  ", prefix_char,
         e_info[e].Mcode, e_info[e].mfg, e_info[e].navg, e_info[e].ntot,
         e_info[e].npeak, e_info[e].t2);
  for(d = 0; e_info[e].delay[d] >= 0; d++)
    fprintf(stream,"%d ", e_info[e].delay[d]);
  fprintf(stream,"\n");
}

print_engine_header() {
  fprintf(stream,"%c\n", prefix_char);
  fprintf(stream,"%c                       Average    Total      Peak      Burn\n",prefix_char);
  fprintf(stream,"%c                       Thrust    Impulse    Thrust     Time    Delays\n",prefix_char);
  fprintf(stream,"%c   Type     Mfg/Opt     (Nt)     (Nt-Sec)    (Nt)     (Sec)    (Sec)\n",prefix_char);
  fprintf(stream,"%c --------  ---------  --------   --------  --------  ------  ----------\n",prefix_char);

}

load_engine(e) int e; {
  int k, dia, len, more, x;
  char s[256], *c, sdelay[32];
  double sum,peak, t1, t2, f1, f2;
  double *p;

  for(more = TRUE; more;) {
    if((c = fgets(s, 256, efile)) != NULL) {
      if(*c != ';') {
        x = sscanf(s, "%s %d %d %s %lf %lf %s", e_info[e].Mcode, &e_info[e].dia,
                   &e_info[e].len, sdelay, &e_info[e].m2, &e_info[e].wt,
                   e_info[e].mfg);
        if( x == EOF) 
        {
          return(-1);
        }
        else if( x != 7) {
          printf("\n*** Error - Bad line in rasp.eng\n%c", Bell);
          Error = TRUE;
          return(-1);
        }
        more = FALSE;
      }
    }
    else
      return(-1);
  }
  k = 0;
  for(more = TRUE; more;) {
    if((c = fgets(s, 256, efile)) != NULL) {
      if(*c != ';') {
        if(k >= MAXT) {
          printf("\n*** Error - Thrust table overflow in rasp.eng\n%c", Bell);
          Error = TRUE;
          return(-1);
        }

        x = sscanf(s, "%lf %lf", &e_info[e].thrust[k], &e_info[e].thrust[k+1]);
        if( x != 2 ) {
          printf("\n*** Error - Bad thrust table line in rasp.eng\n%c", Bell);
          Error = TRUE;
          return(-1);
        }
        if(e_info[e].thrust[k+1] == 0) {
          e_info[e].t2 = e_info[e].thrust[k];
          more = FALSE;
        }
        else
          k += 2;
      }
    }
    else {
      printf("\n*** Error - unexpected EOF while loading thrust curve.\n%c", Bell);
      Error = TRUE;
      return(-1);
    }
  }
  t1 = f1 = sum = peak = 0.0;
  p = &e_info[e].thrust[0];
  do {
    t2 = *p++;
    f2 = *p++;
    if(f2 > peak) peak = f2;
    sum += (t2 - t1) * (f1 + f2) * 0.5;
    t1 = t2;
    f1 = f2;
  } while (f1 != 0.0);
  e_info[e].ntot = sum;
  e_info[e].npeak = peak;
  if( t1 > 0.0 )
     e_info[e].navg = sum / t1;
  else
     e_info[e].navg = 0.0;
  c = sdelay;
  for(k = 0; k < 7; k++) {
    e_info[e].delay[k] = (int) strtol(c, &c, 10);
    if(*c != '-') {
      e_info[e].delay[k+1] = -1;
      break;
    }
    else
      c++;

  }
  return(e);
}

int getmotor(st) int st; 
  {
    int i, match, e, ReTry;
    char elist[80], stemp[14];
    
    match = FALSE;
    e = nexteng++;

    do{
       if (stagenum > 1)
          printf("\n\n\n ********** Stage %d ********** \n", st+1 );

       printf("\nMotor Code? [ %s ]:  ",stages[st].Mcode);          
       get_input( tty_input, TTYmax );
       if( sscanf( tty_input, "%s", Mcode) == EOF ) 
              strcpy( Mcode, stages[st].Mcode );
       else               
          for (i = 0; i < strlen(Mcode); i++)
              if(islower(Mcode[i]))
                 Mcode[i] = toupper(Mcode[i]);
   
       if((efile = fopen(efile_name,"r")) == NULL) {
         printf("\n*** Error - unable to open engine database %s ***\n%c",
                efile_name, Bell);
         Error = TRUE;
         return(-1);
       }

       while(load_engine(e) >= 0)
       {
         if (strcmp( Mcode, e_info[e].Mcode ) == 0) 
         {
           strcpy( stages[st].Mcode, Mcode );
           match = TRUE;
           break;
         }
       }
       if( (st == 0) && ( strcmp( Mcode, "NONE") == 0) )
       {
          printf("\n*** Error - \"NONE\" is not a valid choice for a first stage engine\n\n%c",
                 Bell );
          strcpy( stages[st].Mcode, "");
          match = FALSE;
       }
       else if( !match )
       {
          printf("\n*** Error - Engine \"%s\" not found in file RASP.ENG\n%c", Mcode, Bell);
          strcpy( stages[st].Mcode, "");
          printf("Valid choices are: \n\n");
          rewind(efile);
          strcpy( elist, "" );
          i = 0;
          while(load_engine(e) >= 0)
          {
             if( (st >0) || ( strcmp( e_info[e].Mcode, "NONE") != 0) )
             {
                sprintf( stemp, "%13.12s", e_info[e].Mcode);          
                if( i <= 5 )                /* Columns 0 - thru - 5 on each line */
                {
                   i++;
                   strcat( elist, stemp );
                 }
                else
                {
                   printf("%s\n", elist);
                   strcpy( elist, stemp );
                   i = 1;
                }
             }
          }
          printf("%s\n\n", elist);
        }
        fclose(efile);
    } while( !match );
    fclose(efile);
    if( Error ) return(-1);
    return(e);
  }


choices() {
    int i, jj;
    double wt;

    for(;;)
    {
       printf("\nNumber of stages? [ %d ]:  ", stagenum );
       get_input( tty_input, TTYmax );
       if( sscanf(tty_input, "%d", &ianswer) != EOF )
          stagenum = ianswer;
       if( stagenum < MINSTAGE )
       { 
          printf("\n*** Error - \"%d\" is too small! Minimum allowed is \"%d\" ***\n%c",
                 stagenum, MINSTAGE, Bell );
          stagenum = MINSTAGE;
       }
       else if( stagenum > MAXSTAGE )
       { 
          printf("\n*** Error - \"%d\" is too large! Maximum allowed is \"%d\" ***\n%c",
                 stagenum, MAXSTAGE, Bell );
          stagenum = MAXSTAGE;
       }
       else break;
    }


    for (i = 0; i < stagenum; i++) {
        stages[i].engcnum = getmotor(i);
        if( Error )
            return(-1);

      for(;;)
      {
        if (stagenum > 1)
            printf("\nWeight of stage %d w/o engine in ounces? [ %.3f ]:  ",
                    i+1, stages[i].weight_Oz);
        else
            printf("\nWeight of rocket w/o engine in ounces? [ %.3f ]:  ",
                   stages[i].weight_Oz);
        get_input( tty_input, TTYmax );
        if( sscanf(tty_input, "%lE", &danswer) != EOF )
           stages[i].weight_Oz = danswer;
        stages[i].weight_Kg = stages[i].weight_Oz * OZ2KG;
        if( stages[i].weight_Oz < MINWGHT )
        { 
           printf("\n*** Error - \"%.3f\" is too small! Minimum allowed is \"%.3f\" ***\n%c",
                  stages[i].weight_Oz, MINWGHT, Bell );
           stages[i].weight_Oz = MINWGHT;
        }
        else if( stages[i].weight_Oz > MAXWGHT )
        { 
           printf("\n*** Error - \"%.3f\" is too large! Maximum allowed is \"%.3f\" ***\n%c",
                  stages[i].weight_Oz, MAXWGHT, Bell );
           stages[i].weight_Oz = MAXWGHT;
        }
        else break;
      }
        

      for(;;)
      {
        if (stagenum > 1)
            printf("\nNumber of engines in stage %d? [ %d ]:  ", i+1, stages[i].numeng);
        else
            printf("\nNumber of engines? [ %d ]:  ", stages[i].numeng);
        get_input( tty_input, TTYmax );
        if( sscanf(tty_input, "%d", &ianswer) != EOF )
           stages[i].numeng = ianswer;
        if( stages[i].numeng < MINENG )
        { 
           printf("\n*** Error - \"%d\" is too small! Minimum allowed is \"%d\" ***\n%c",
                  stages[i].numeng, MINENG, Bell );
           stages[i].numeng = MINENG;
        }
        else if( stages[i].numeng > MAXENG )
        { 
           printf("\n*** Error - \"%d\" is too large! Maximum allowed is \"%d\" ***\n%c",
                  stages[i].numeng, MAXENG, Bell );
           stages[i].numeng = MAXENG;
        }
        else break;
      }


      for(;;)
      {
        if (stagenum > 1)
            printf("\nMaximum diameter of stage %d in inches? [ %.3f ]:  ",
                   i+1, stages[i].maxdia);
        else
            printf("\nMaximum body diameter in inches? [ %.3f ]:  ", stages[i].maxdia);
        get_input( tty_input, TTYmax );
        if( sscanf(tty_input, "%lE", &danswer) != EOF )
           stages[i].maxdia = danswer;
        if( stages[i].maxdia < MINDIA )
        { 
           printf("\n*** Error - \"%.3f\" is too small! Minimum allowed is \"%.3f\" ***\n%c",
                  stages[i].maxdia, MINDIA, Bell );
           stages[i].maxdia = MINDIA;
        }
        else if( stages[i].maxdia > MAXDIA )
        { 
           printf("\n*** Error - \"%.3f\" is too large! Maximum allowed is \"%.3f\" ***\n%c",
                  stages[i].maxdia, MAXDIA, Bell );
           stages[i].maxdia = MAXDIA;
        }
        else break;
      }
        

        if (i + 1 < stagenum) 
        {
          for(;;)
          {
            printf("\nStaging delay for stage %d in seconds? [ %.3f ]:  ",
                    i+2, stages[i].delay);
            get_input( tty_input, TTYmax );
            if( sscanf(tty_input, "%lE", &danswer) != EOF )
               stages[i].delay = danswer;
            if( stages[i].delay < MINDELAY )
            { 
               printf("\n*** Error - \"%.3f\" is too small! Minimum allowed is \"%.3f\" ***\n%c",
                      stages[i].delay, MINDELAY, Bell );
               stages[i].delay = MINDELAY;
            }
            else if( stages[i].delay > MAXDELAY )
            { 
               printf("\n*** Error - \"%.3f\" is too large! Maximum allowed is \"%.3f\" ***\n%c",
                      stages[i].delay, MAXDELAY, Bell );
               stages[i].delay = MAXDELAY;
            }
            else break;
          }
        }

               /* figure start and stop times for motor burn  and stage */
        if (i == 0)
            stages[i].start_burn = 0;
         else
            stages[i].start_burn = stages[i-1].end_stage;

            stages[i].end_burn = stages[i].start_burn + e_info[stages[i].engcnum].t2;
            stages[i].end_stage = stages[i].end_burn + stages[i].delay;

                  /* figure weight for stage and total rocket */


        stages[i].totalw = stages[i].weight_Kg +
                           (e_info[stages[i].engcnum].wt * stages[i].numeng);
        rocketwt += stages[i].totalw;

    }


    for(;;)
    {
       printf("\n\n\nDrag coefficient? [ %.3f ]:  ", drag_coff);
       get_input( tty_input, TTYmax );
       if( sscanf(tty_input, "%lE", &danswer) != EOF )
       drag_coff = danswer;
       if( drag_coff < MINDRAGC )
       { 
          printf("\n*** Error - \"%.3f\" is too small! Minimum allowed is \"%.3f\" ***\n%c",
                 drag_coff, MINDRAGC, Bell );
          drag_coff = MINDRAGC;
       }
       else if( drag_coff > MAXDRAGC )
       { 
          printf("\n*** Error - \"%.3f\" is too large! Maximum allowed is \"%.3f\" ***\n%c",
                 drag_coff, MAXDRAGC, Bell );
          drag_coff = MAXDRAGC;
       }
       else break;
    }
            
    for(;;)
    {
       printf("\nBase temperature (Deg F)? [ %.3f ]:  ", base_temp);
       get_input( tty_input, TTYmax );
       if( sscanf(tty_input, "%lE", &danswer) != EOF )
       base_temp = danswer;
       if( base_temp < MINBASET )
       { 
          printf("\n*** Error - \"%.3f\" is too small! Minimum allowed is \"%.3f\" ***\n%c",
                 base_temp, MINBASET, Bell );
          base_temp = MINBASET;
       }
       else if( base_temp > MAXBASET )
       { 
          printf("\n*** Error - \"%.3f\" is too large! Maximum allowed is \"%.3f\" ***\n%c",
                 base_temp, MAXBASET, Bell );
          base_temp = MAXBASET;
       }
       else
       {
          base_temp_K = (base_temp - 32) * 5/9 + 273.15;  /* convert to degK */
          break;
       }
    }
            
    GNUplot = FALSE;
    prefix_char = noGNUprefix;
    do {

      minval = 1;
      maxval = 4;
      for(;;)
      {
        printf("\nSend data to:\n\t(1) Screen\n\t(2) Printer\n\t");
        printf("(3) Disk file\n\t(4) GNUplot data and command files\n\nChoice?[%d] ", destnum);
        get_input( tty_input, TTYmax );
        if( sscanf(tty_input, "%d", &ianswer) != EOF )
           destnum = ianswer;
        if( destnum < minval )
        { 
           printf("\n*** Error - \"%d\" is too small! Minimum allowed is \"%d\" ***\n%c",
                  destnum, minval, Bell );
           destnum = minval;
        }
        else if( destnum > maxval )
        { 
           printf("\n*** Error - \"%d\" is too large! Maximum allowed is \"%d\" ***\n%c",
                  destnum, maxval, Bell );
           destnum = maxval;
        }
        else break;
      }


#ifdef MSDOS
        if (destnum == 1)
            strcpy(fname, "CON");
        else if (destnum == 2)
            strcpy(fname, "PRN");
#endif
#ifdef UNIX
        if (destnum == 1)
            strcpy(fname, "/dev/tty");
        else if (destnum == 2)
            strcpy(fname, "/dev/lp");
#endif
#ifdef VMS
        if (destnum == 1)
            strcpy(fname, "TT:");
        else if (destnum == 2)
            strcpy(fname, "LP:");
#endif
#ifdef MAC
        if (destnum == 1){
            console_options.title  = "\pRASP Output";
            console_options.txFont = monaco;
            console_options.txSize = 9;   /* Font size in points */
            console_options.ncols  = 80;
            console_options.nrows  = 25;
            console_options.top    =100;
            console_options.left   = 30;
            if( !MsgSent ) {
               MsgSent = TRUE;
               printf("\n\nScrolling of the display may be throttled with the mouse button.\n");
               sleep(1);
            }
         }
        else if (destnum == 2){
            console_options.title  = "\pRASP Print Window";
            console_options.txFont = courier;
            console_options.txSize = 12;   /* Font size in points */
            console_options.nrows  = 25;
            console_options.ncols  = 80;
            console_options.top    = 50;
            console_options.left   = 10;
        }
#endif
        else if (destnum == 3){
            printf("Enter file name: ");
            scanf("%s", fname);
            dump_buffer();
        }
        else{
            GNUplot = TRUE;
            prefix_char = GNUprefix;
            printf("Enter file name PREFIX: ");
            scanf("%s", fnam_prefix);
            dump_buffer();
            strcpy( fname, fnam_prefix );
            strcat( fname, GNUplot_Dat_Ext );
#if GNUcmd
            strcpy( fnam2, fnam_prefix );
            strcat( fnam2, GNUplot_Cmd_Ext );
            if ((stream2 = fopen(fnam2,"w")) == NULL)
                printf("%s cannot be opened.\n",fnam2);
#endif
        }


#ifdef MAC
        if( (destnum == 3) || (destnum == 4) ){
           if ((stream = fopen(fname,"w")) == NULL)
                printf("%s cannot be opened.\n",fname);
         }
        else{
           if ((stream = fopenc()) == NULL)
                printf("CONSOLE cannot be opened.\n");
           else 
                if (destnum == 2) cecho2printer( stream );
        }
#else
        if ((stream = fopen(fname,"w")) == NULL)
             printf("%s cannot be opened.\n",fname);
#endif


    } while (stream == NULL);

    
    wt = rocketwt;
#if GNUcmd
    if ( GNUplot ) {
       fprintf(stream2,"#\n");
       fprintf(stream2,"# To use this file issue the following command to GNUplot:\n");
       fprintf(stream2,"#             load \"%s\" \n", fnam2);
       fprintf(stream2,"#\n");
       fprintf(stream2,"# The following commands (and many others) can be interactively issued\n");
       fprintf(stream2,"# to GNUplot. For GNUplot help try the HELP command under GNUplot.\n");
       fprintf(stream2,"# Other help is available through the newsgroup COMP.GRAPHICS.GNUPLOT.\n");
       fprintf(stream2,"# A GNUplot FAQ is available from RTFM.MIT.EDU [18.70.0.209] in the\n");
       fprintf(stream2,"# directory /pub/usenet/comp.graphics.gnuplot. Check the FAQ for the\n");
       fprintf(stream2,"# locations of sites distributing the sources and executables for\n");
       fprintf(stream2,"# various platforms.\n");
       fprintf(stream2,"#\n");
       fprintf(stream2,"#\n");
       fprintf(stream2,"set data style lines\n");
       fprintf(stream2,"set xlabel \"Seconds\" 0,0 \n");
       fprintf(stream2,"set ylabel 0,0 \n");
       fprintf(stream2,"set grid\n");

       fprintf(stream2,"set title \"Altitude  -  %s\"\n", fname);
       fprintf(stream2,"set ylabel \"Feet\"\n");
       fprintf(stream2,"plot \"%s\" using 1:2 title ''\n", fname);

       fprintf(stream2,"set title \"Velocity  -  %s\"\n", fname);
       fprintf(stream2,"set ylabel \"Feet/Sec\"\n");
       fprintf(stream2,"plot \"%s\" using 1:3 title ''\n", fname);

       fprintf(stream2,"set title \"Acceleration  -  %s\"\n", fname);
       fprintf(stream2,"set ylabel \"Feet/Sec*Sec\"\n");
       fprintf(stream2,"plot \"%s\" using 1:4 title ''\n", fname);

       fprintf(stream2,"set title \"Weight  -  %s\"\n", fname);
       fprintf(stream2,"set ylabel \"Ounces\"\n");
       fprintf(stream2,"plot \"%s\" using 1:5 title ''\n", fname);

       fprintf(stream2,"set title \"Thrust  -  %s\"\n", fname);
       fprintf(stream2,"set ylabel \"Newtons\"\n");
       fprintf(stream2,"plot \"%s\" using 1:6 title ''\n", fname);
       
       fprintf(stream2,"#\n");
       fprintf(stream2,"# To generate a Postscript file of the plots for printing\n");
       fprintf(stream2,"# issue the following commands to GNUplot:\n");
       fprintf(stream2,"#\n");
       fprintf(stream2,"# 1)          show terminal       (make a note of your current terminal type)\n");
       fprintf(stream2,"# 2)          set terminal postscript \n");
       fprintf(stream2,"# 3)          set output \"%s.ps\" \n", fnam_prefix);
       fprintf(stream2,"# 4)          load \"%s\"         (postscript file is produced here)\n", fnam2);
       fprintf(stream2,"# 5)          set output          (reset from step #2)\n");
       fprintf(stream2,"# 6)          set terminal back-to-whatever-it-was-in-step-1 \n");
       fprintf(stream2,"#\n");
    }
#endif

    for (i = 0; i < stagenum; i++) {
      fprintf(stream,"%c\n", prefix_char);
      fprintf(stream,
              "%cStage   Engine            Bare     Launch  Dia.\n",
               prefix_char);
      fprintf(stream,
              "%cNumber  (Quant)Type       Wt(oz)   Wt(oz)  (in)      Cd\n",
              prefix_char);
      fprintf(stream,
           "%c------  ----------------  ------   ------  --------  ----\n",
             prefix_char);
      fprintf(stream,
              "%c   %d    (%d) %-12s   %.2lf     %.2lf     %.2lf    %.2lf\n",
              prefix_char, i+1, stages[i].numeng, e_info[stages[i].engcnum].Mcode,
              stages[i].weight_Oz, wt / OZ2KG, stages[i].maxdia, drag_coff);
      wt -= stages[i].totalw;
      print_engine_header();
      print_engine_info(stages[i].engcnum);
    }
      if (verbose) {
         fprintf(stream,"%c\n", prefix_char);
         fprintf(stream, "%c%4s %11s %10s %12s %8s %8s %10s\n",prefix_char,
                 "Time", "Altitude", "Velocity",    "Acceleration", "Weight", "Thrust", "Drag");
         fprintf(stream, "%c%4s %10s %10s %12s %9s %9s %9s\n", prefix_char,
                 "(Sec)", "(Feet)", "(Feet/Sec)", "(Feet/Sec^2)", "(Ounces)", "(Newtons)" , "(Factor)");
         fprintf(stream, "%c%4s %10s %10s %10s %10s %10s %10s\n", prefix_char,
                 "-----", "------", "---------", "--------", "--------", "--------", "--------");
       } 
       else fprintf(stream, "%c\n", prefix_char);
}
calc() {

    double floor();
    int j, k;                         /* iteration variable */
    int thrust_index = 0;             /* index into engine table */
    int stage_engcnum,                /* engcnum for current stage */
        this_stage = 0;               /* current stage */
    double stage_time = 0.0,          /* elapsed time for current stage */
           burn_time = 0.0,           /* Elapsed time for motor burn */
           stage_wt,                  /* current wt of this_stage */
           stage_diam;                /* max diam at current time */
    double t_rod = 0.0, v_rod, x_rod; /* launch rod info */
    double c;
    double cc, drag_bias;
    double thrust = 0.0;              /* Thrust */
    double vel = 0.0;                 /* Velocity */
    double accel = 0.0;               /* Acceleration */
    double alt = 0.0;                 /* Altitude */
    double print_vel, print_alt, print_accel, print_mass;  /* used for printing */

    double mass = 0.0;
    double t = 0.0;
    double max_accel = 0.0;
    double max_vel = 0.0;
    double old_thrust = 0.0;       /* last engine thrust  from table */
    double old_time = 0.0;         /* last engine thrust time  from table */
    double d, m1;                  /* temp vars */

    int print_index = 0;           /* used for print loop */
    int launched = 0;              /* indicates rocket has lifted off */

    double delta_t = 0.001;        /* Time interval - 1ms */
    double rho;                      /* Constant used for drag calc */
    double air_density();

    stage_wt = stages[this_stage].totalw;
    stage_engcnum = stages[this_stage].engcnum;
    mass = rocketwt;
    


    for (j = 0; j < stagenum; j++) {
        stage_diam = stages[j].maxdia;
        if (j > 0) {
            if (stage_diam < stages[j-1].maxdia)
                stage_diam = stages[j-1].maxdia;
        }
    }
    d = stage_diam * IN2M;

    for(;;) {
        /* Calculate decreasing air density */ 
        rho = air_density (alt,base_temp_K); 
        c = rho * PI * drag_coff * d * d * 0.125;

        t += delta_t;
        stage_time += delta_t;

                  /* handle staging, if needed */

        if ((t > stages[this_stage].end_stage) && (this_stage < stagenum-1)) {
            thrust_index = 0;
            old_thrust = 0.0;
            old_time = 0.0;
            mass = rocketwt -= stage_wt;
            this_stage++;
            stage_wt = stages[this_stage].totalw;
            stage_engcnum = stages[this_stage].engcnum;
            stage_time = burn_time = 0;
            for (j = this_stage; j < stagenum; j++) {
                stage_diam = stages[j].maxdia;
                if (j > this_stage) {
                    if (stage_diam < stages[j-1].maxdia)
                        stage_diam = stages[j-1].maxdia;
                }
            }
            d = stage_diam * IN2M;
            c = rho * PI * drag_coff * d * d * 0.125;
           if (GNUplot)
              fprintf(stream,"\n");
           fprintf(stream,"%cStage %d Ignition.\n", prefix_char, this_stage+1);
        }

                  /* Handle the powered phase of the boost */

        if ((t >= stages[this_stage].start_burn) &&
            (t <= stages[this_stage].end_burn)) {

               burn_time +=  delta_t;  /* add to burn time */

                /* see if we need to use the next thrust point */
          /* All times are relative to burn time for these calculations */

            if (burn_time > e_info[stage_engcnum].thrust[thrust_index]) {
                old_time = e_info[stage_engcnum].thrust[thrust_index];
                thrust_index++;
                old_thrust = e_info[stage_engcnum].thrust[thrust_index];
                thrust_index++;
            }

/*

Logic to smooth transision between thrust points.  Transitions are
linear rather than discontinuous.

*/

            thrust = e_info[stage_engcnum].thrust[thrust_index + 1] - old_thrust;
            thrust *= (burn_time - old_time) / (e_info[stage_engcnum].thrust[thrust_index] - old_time);
            thrust += old_thrust;
            thrust *= stages[this_stage].numeng;
            m1 = (e_info[stage_engcnum].m2 / e_info[stage_engcnum].t2) * delta_t;
            stage_wt -= m1;
            mass = rocketwt -= m1;
          } else thrust = 0;

/*

Crude approximations for MACH 1 Transition.
The drag bias will be applied to the subsonic value.
It will be equal to 1.0 for velocities less than MACH 0.8.
From MACH 0.8 to MACH 1.0 it will increase linearly to 2.0 times
the subsonic value.  From MACH 1.0 to MACH 1.2 it will decrease to
1.5 times the subsonic value.  Above MACH 1.2 it will remain at
a value of 1.5

*/
        if (vel < MACH0_8)
          drag_bias = 1;
        else if (vel < MACH1)
          drag_bias = (1.0 + ((vel - MACH0_8) / (MACH1 - MACH0_8)));
        else if (vel < MACH1_2)
          drag_bias = (2.0 - 0.5*((vel - MACH1) / (MACH1_2 - MACH1)));
        else
          drag_bias = 1.5;

        cc = c * drag_bias;

            /* Simple Newton calculations */

        accel = ((thrust - (cc * vel * vel)) / mass) - G;
        vel = vel + accel * delta_t;
        alt = alt + vel * delta_t;

            /* test for lift-off and apogee */

       if (vel > 0)
           launched = 1;                 /* LIFT-OFF */
       else if (!launched && (vel < 0))
           alt = vel = accel = 0;        /* can't fall off pad! */
       else if (launched && (vel < 0))
           break;                        /* apogee, all done */

        if(alt <= ROD) {
          t_rod = t;
          v_rod = vel * M2FT;
          x_rod = alt * M2FT;
        }

            /* do max evaluations */

        if (accel > max_accel) max_accel = accel;
        if (vel > max_vel) max_vel = vel;

            /* See if we need to print out the current values */

          if (print_index++  == 200) {   /* print every 200 ms */
              print_index = 0;
              print_alt = alt * M2FT;
              print_vel = vel * M2FT;
              print_accel = accel * M2FT;
              print_mass = mass / OZ2KG;
              if (verbose) fprintf(stream,"%4.1lf %10.1lf %10.1lf %10.1lf %10.3lf %10.3lf %10.3lf\n",
                     t + 0.01, print_alt, print_vel, print_accel, print_mass, thrust, drag_bias);
          }
   }
    print_alt = alt * M2FT;
    print_vel = vel * M2FT;
    print_accel = accel * M2FT;
    print_mass = mass / OZ2KG;
    fprintf(stream, "%cMaximum altitude attained was %.1lf feet (%.1lf meters).\n",
            prefix_char, print_alt, alt);
    fprintf(stream, "%cTime to peak altitude = %.2lf seconds.\n",
           prefix_char, t);
    fprintf(stream, "%cMaximum velocity was %.1lf feet/sec.\n", 
           prefix_char, max_vel * M2FT);
    fprintf(stream, "%cMaximum acceleration was %.1lf feet/sec^2.\n",
            prefix_char, max_accel * M2FT);
    fprintf(stream, "%cLaunch rod time = %.2lf,  altitude = %.1lf,  velocity = %.1lf\n",
            prefix_char, t_rod, x_rod, v_rod);
   if( GNUplot ){
      fclose(stream);
#if GNUcmd
      fclose(stream2);
#endif
   }
   else{
#ifdef MAC
         if( destnum != 1  ){
            fclose(stream);
         }
#else    
         fclose(stream);
#endif
   }
}



main(argc, argv)
int argc;
char*argv[];
{
    char ans[3] = "N", canswer[3];
    int i;
    
#ifdef MAC
    console_options.txFont = monaco;
    console_options.txSize = 9;   /* font size in points */
    console_options.procID = 8;
    console_options.nrows  = 25;
    console_options.ncols  = 80;
    console_options.top    = 50;
    console_options.left   = 10;
    console_options.title  = "\pRASP Input";
    console_options.pause_atexit = 0;
    argc = ccommand(&argv);
#endif   

    if (argc > 1) verbose = FALSE;
    
    for ( i=0; i<MAXSTAGE; i++ )    /* initialize defaults for stages */
    {
       stages[i].engcnum    = 0;
       stages[i].numeng     = 1;
       stages[i].start_burn = 0.0;
       stages[i].end_burn   = 0.0;
       stages[i].end_stage  = 0.0;
       stages[i].weight_Oz  = 0.0;
       stages[i].weight_Kg  = 0.0;
       stages[i].totalw     = 0.0;
       stages[i].maxdia     = 0.0;
       stages[i].delay      = 0.0;
       strcpy( stages[i].Mcode, "" );
    }

    printf("\nRASP - Rocket Altitude Simulation Program V%s\n\n", VERSION);

    for (;;)
    {
        Error = FALSE;
        rocketwt = 0;
        nexteng = 0;
        choices();
        if( !Error ) calc();

        if( Error )
           printf("\nDo want to try again (Y/N)? [ %s ]:  ", ans);
        else
           printf("\nDo another one (Y/N)? [ %s ]:  ", ans);

        get_input( tty_input, TTYmax );
        if( sscanf(tty_input, "%s", canswer) != EOF )
           strcpy( ans, canswer );
        if ((strncmp(ans,"y",1) == 0) || (strncmp(ans,"Y",1) == 0)) 
        {
#ifdef MAC
           if( destnum == 1 ) fclose(stream);
#endif            
           continue;
         }
         else exit(0);
     }
}
