/***************************************************************************/
/*

		      <<< PROGRAM DIGITRAK (tm) >>>
	   *** VERTICAL ROCKET TRAJECTORY SIMULATION AND ***
       *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS ***

		     *** AUTHOR: LARRY CURCIO ***

			  <<< VERSION 3.6>>>
		     LAST UPDATED 11 MARCH 1994

	       BUG REPORTS, QUESTIONS, SPELLING CORRECTIONS,
		 AND COMMENTS CAN BE SENT TO ME DIRECTLY:

	      FREE DISTRIBUTION AS LONG AS NO PROFIT IS MADE
		  FEE FOR EACH FOR-PROFIT SALE: $25.00
	       PROGRAM PROVIDED AS-IS. NO WARRANRY IS MADE

		   *** INTERNET: LC2B@DELPHI.COM ***

		       *** SNAILMAIL ADDRESS ***
			    P.O. BOX 81124
		      PITTSBURGH, PA 15217-81124

		 *** COPYRIGHT LAWRENCE M. CURCIO ***
		       *** 16 NOVEMBER 1993 ***


			 *** DESCRIPTION ***
VERTICAL ROCKET TRAJECTORY SIMULATION PROGRAM WITH FOLLOWING FEATURES:

1) SUPPORTS UP TO THREE STAGES WITH OR WITHOUT STAGING DELAYS

2) SUPPORTS SIMPLE STAGING DELAYS OR AIR START DELAYS, WHERE SEPARATION
   OCCURS AT LOWER STAGE CUTOFF, BUT UPPER STAGE IGNITION IS DELAYED.

3) USES CLOSED-FORM SOLUTIONS AT EACH DIGITAL INTERVAL

4) DEPRECIATES AIR DENSITY WITH ALTITUDE

5) USES EMPIRICAL DRAG DIVERGENCE MODEL PROPOSED BY BENGEN IN
    MANDELL, CAPORASO AND BENGEN: _TOPICS IN ADVANCED MODEL ROCKETRY_
    (MODEL GOOD FOR QUALITATIVE EFFECTS, AND IMPLICATION OF MACH EFFECTS)

6) SEPARATE TREATMENT OF ROUND AND OPTIMALLY SHARP NOSES IN DRAG DIVERGENCE

7) ADJUSTMENT OF AIR DENSITY AND MACH NUMBER TO TEMPERATURE

8) USES DRAG COEFFICIENTS, SHAPE CONSTANTS, OR DRAG PROFILES

9) USES FINE-INTERVAL OR COARSE-INTERVAL THRUST CURVES. INTERPOLATES
   COARSE-INTERVAL CURVES.

10) TEMPORAL BACKTRACKING(tm) CAPABILITY:
    CALCULATION OF TRAJECTORY STATS AND DRAG COEFFICIENT FROM ASCENT TIME

11) OPTICAL BACKTRACKING(tm) CAPABILITY
    CALCULATION OF TRAJECTORY STATS AND DRAG COEFFICIENT FROM MAXIMUM ALTITUDE

12) OPTIONAL USE OF DRAG PROFILES IN TEMPORAL OR OPTICAL BACKTRACKING(tm)

13) ESTIMATION OF OPTIMAL MASS FOR SINGLE-STAGE ROCKETS OR STAGED ROCKETS
    WITH NO STAGING DELAYS

14) WRITTEN BY A SUPERSTITIOUS AUTHOR, WHO DOESN'T LIKE TO HAVE 13 LISTED
    FEATURES.


			** VERSION 2.0 **
      ADDED RUNGE-KUTTA, AND FIXED SOME TABS BEFORE #if'S TO
      PLEASE VARIOUS UNIX CC COMPILERS.

			** VERSION 2.1 **
      ADDED FEHSKENS-MALEWICKI INTEGRALS FOR LAST TIME INTERVAL
      THIS DISSOCIATED TIMING ACCURACY AND INTERVAL WIDTH, CREATING
      CONTINUOUS TIME. CONTINUOUS TIME ALLOWED HIGHER TOLERANCES FOR
      BACKTRACKING, AND ALSO WIDER COAST INTERVALS. RESULT IS A FASTER,
      MORE ACCURATE BACKTRACKING TOOL.

 			 ** VERSION 2.2 **
       FIXED INSUBSTANTIAL BUGS IN CONTINUOUS TIME ROUTINE AND IN MASS
       INITIALIZATION. IMPLEMENTED ENGLISH INPUT OPTION (#define INUNITS).
       ALLOWED LAUNCHER > SEA LEVEL (#define ALT0). Fixed bug in
       BACKTRK_STAGE routine that would bomb on a single stage rocket
       after an upper stage in a multistage rocket was backtracked. Also,
       header records in the output file are preceeded by '#' to make them
       gnuplot compatible.

			 ** VERSION 3.0 **
	FULL FEHSKENS-MALEWICKI TREATMENT AT EACH INTERVAL. VARIABLE
	COAST PHASE TIME INTERVALS.

			 ** VERSION 3.1 **
	CONTINUOUS TIME EXTENDED TO STAGING DELAYS. THE DEFAULT NUMBER
	OF INTERVALS IN THE BOOST PHASE REDUCED FROM 500 TO 400.

			 ** VERSION 3.2 **
	FIXED BUG IN DRAG PROFILE HANDLING PROCEDURE. IMPLEMENTED COLOR
	FOR TURBO DOS. ALSO IMPLEMENTED TURBODOS SWITH, WHICH SWITCHES
	BETWEEN DOS (#define TURBODOS 1) AND UNIX OR GENERIC
	(#define TURBODOS 0)

			 ** VERSION 3.3 **
	OUTPUT FILE FIELDS ARE LABELED. OUTPUT FILE UNITS ARE DETERMINED
	BY INUNITS OPTION. PROGRAM IS BULLET-PROOFED. FIXED MINOR
	INCONSISTENCY IN THRUST CURVE INTERPOLATION.

			 ** VERSION 3.4 **
	       FURTHER BULLETPROOFING IN CD ENTRY ALSO A PROBLEM
	       IN THE HANDLING OF 3-STAGE ROCKETS WAS FIXED

			 ** VERSION 3.5 **
	       DIGITRAK WOULD GO INTO A TIGHT LOOP IF BACKTRAKING
	       WERE INVOKED ON THE FIRST PASS. THA'S FIXED.

			 ** VERSION 3.6 **
	       FIXED BUG IN TREATMENT OF DRAG PROFILES AFTER
	       BACKTRACKING WITHOUT PROFILES. CONSPICUOUS FLAW
	       WHEN PROFILES ARE USED. NO CONSEQUENCE WHEN THEY ARE NOT



NOTE: BACKTRACKED DRAG COEFFICIENTS ARE *VERY* SENSITIVE TO MOTOR PERFORMANCE
      VARIATION AND MEASUREMENT ERRORS. DRAG COEFFICIENTS *MUST* THEREFORE BE
      BACKTRACKED FROM AVERAGE ALTITUDES FROM SEVERAL FLIGHTS. ALTITUDES
      THEMSELVES MAY BE MEASURED DIRECTLY, OR THROUGH TEMPORAL BACKTRACKING.
      BACKTRACKED ALTITUDES ARE FAR LESS SENSITIVE TO MOTOR PERFORMANCE
      VARIATION THAN BACKTRACKED DRAG COEFFICIENTS ARE.

      MULTI-STAGE VEHICLES SHOULD BE BACKTRACKED BOTTOM-UP, NOT TOP-DOWN.
      MOST 2-STAGE VEHICLE BACKTRACKING IS EXTREMELY INSENSITIVE TO Cd OF
      FIRST STAGE.

		  **** COMPILATION NOTES ****
      SCAN PROGRAM FOR COMPILE OPTION FLAGS (PLURAL). CHANGE #define'S
      ACCORDINGLY. FIRST SECTION ALSO HAS A (ONE) HACKERS' AREA. THESE
      #define'S AREN'T NORMALLY CHANGED, BUT HACKERS WILL BE
      HACKERS. LOOK FOR THE WORD, HACKER.

      THIS VERSION COMPILES IN DOS TURBO C (1 WARNING), WINDOWS TURBO C++
      (SCADS OF WARNINGS FOR NO PROTOTYPES), AND UNIX AT CMU
      (CLEAN COMPILE WHEN #define CLS IS CHANGED TO printf("\n\n\n");

      ****** FOR COMPILERS THAT DON'T SUPPORT #if'S ******
      1) PUT { AT END OF EACH #if LINE
      2) REMOVE # FROM EACH #if (CREATING SIMPLE if)
      3) CHANGE #else TO:   } else {
      4) CHANGE #endif TO }
      5) FOLLOW INSTRUCTIONS AT AREA MARKED  #if EXCEPTIONS
      6) COMPILE
      7) IGNORE MYRIAD WARNINGS FOR NULL if'S, CODES NO EFFECT,
	 AND UNREACHABLE CODES. (COMPILERS THINK THEY KNOW EVERYTHING.)

      THERE IS A FUNCTION CALLED color, IN WHICH YOU CAN PUT
      A SHORT ROUTINE TO SET BACKGROUND COLOR IN YOUR SYSTEM. APOLOGIES
      FOR THIS SPELLING ARE DUE TO THOSE IN THE UK. APOLGIES ARE DUE TO
      EVERYONE FOR THE REST OF THE SPELLING.

      DIGITRAK CAN USE THRUST CURVES FROM SIMLAB (tm). IT CAN ALSO USE
      THRUST CURVES CONVERTED FROM RASP.ENG USING THE PROGRAM, MARKCURVE.C.
      TO USE MARKCURVE, FIRST CHANGE ALL TABS IN RASP.ENG TO SPACES. CALL
      THE NEW FILE ANYTHING YOU LIKE. PUT IT IN A *FRESH DIRECTORY* WITH
      THE MARKCURVE EXECUTABLE. BOWCOO THRUST CURVE FILES RESULT. EACH
      BEGINS WITH THE LETTERS, 'MJ'. THE NEXT LETTER DENOTES THE
      MANUFACTURER, ACCORDING TO RASP TRADITION. THE LAST LETTERS DENOTE
      THE ENGINE TYPE. THE FILE TYPE OF EACH THRUST CURVE IS .ENG.
      DIGITRAK APPENDS THIS TO THE THRUST CURVE NAME THAT YOU SPECIFY
      DURING YOUR INTERACTION WITH IT. THUS -

      TO SPECIFY AN AEROTECH E15, USE THE CODE: MJAE15
      TO SPECIFY AN ESTES E155 (STAND BACK!), USE: MJEE15

      MOST SIMLAB CURVES BEGIN WITH THE PREFIX, NEW

      SIMLAB D1 IS NEWD12

      EXCEPTIONS:
      THE AEROTECH E25 SIMLAB CURVE IS AEROE25. ALSO
      THERE IS AN FLAT, AVERAGE-THRUST C6 CURVE FOR DIAGNOSTICS,
      WHICH IS AVEC6.

      NOTE: AS OF THIS WRITING, THE MJAE25 AND MJAD21 CURVES ARE
      CATALOG VALUES, AND NOT S&T VALUES. USE THE SIMLAB AEROE25
      INSTEAD OF THE MJAE25.

	       ***  16 FEBRUARY 1994 ***
*/
/***************************************************************************/
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <string.h>

/*COMPILE OPTION*/
/* ALL #define's DOWN TO *** END COMPILE OPTIONS *** */

/* IF COMPILING ON TURBO C FOR DOS, SET TURBODOS TO 1 - OTHERWISE 0 */
#define TURBODOS 0

#if(TURBODOS == 1)
#define printf  cprintf
#include <conio.h>
#define CLS clrscr()
#else
#define CLS printf("\n\n\n\n\n\n")
#endif

/* LAUNCH PAD ALTITUDE IN ****METERS**** */
#define ALT0 0;

/* ALTUNITS should be 'F' for feet and 'M' for meters */
#define ALTUNITS 'F'

/* INPUT UNITS ARE 'E' FOR ENGLISH AND 'M' FOR METRIC */
#define INUNITS 'E'

/* SET TO CONSTANT_RHO 1 TO TURN OFF DENSITY/ALTITUDE DEPRECIATION */
/* USEFUL FOR DIAGNOSTIC CHECKS USING CLOSED-FORM COAST PHASE FORMULAS */
#define CONSTANT_RHO 0

/* TEMPUNITS should be 'F' for fahrenheit and 'C' for centigrade */
/* SET TEMPUNITS TO 'C' OR 'F' */
#define TEMPUNITS 'F'

/* SUBSONIC Coast Phase time interval width */
#define DELTAT0 .1
/* X-SONIC AND HIGHER COAST PHASE INTERVAL WIDTH */
#define DELTAT1 .01

/* CLS is the Screen Clear function. Define as clrscr for TURBO C */
/* define as printf("\n\n\n\n\n") for UNIX */



/* Number of intervals between file record output */
#define REPNUM 10

/* *** END COMPILE OPTIONS *** */

/* *** HACKERS' AREA *** */
/* (DOWN TO *** END HACKERS' AREA ***) */
/* *** *HACKERS* MAY WISH TO CHANGE THESE THINGS *** */
/* *** RECOMMEND THEY BE LEFT ALONE *** */

/* MAXIMUM NUMBER OF DRAG ITERATIONS */
#define MAXIT 10

/* MAXIMUM NUMBER OF STAGES */
#define MAXSTAGE 3


/* Size of Drag Profile */
#define DRAGSIZE 1000

/* INTERPOLATION CUTOFF defines num. intervals needing no interpolation */
/* That is, it's the cutoff between long and short thrust curves */
/* It all arises from different standards in the modroc world */
#define INTERPOLATION_CUTOFF 100

/* TOGGLE controls how long thrust curves are read (f, t or t, f) */
#define TOGGLE 0

/* TOGGLE1 controls how short thrust curves are read (f, t or t, f) */
#define TOGGLE1 1

/* Number of documentation lines in drag profile */
#define NUMDOC 5

/* Number of repetitions in Newton's method backtracking procedure */
#define NEWTREP 75

/* Tolerance in Newton's method backtracking procedure */
#define NEWTTOL .000001

/*Precision check. Tells when to stop drag iterations */
/* Also used as initial value of V to avoid singularity */
#define DELTA .000001

/* String Length */
#define SLEN 255

/* Maximum size of Thrust-Time tables */
#define ARRAYSIZE 1001

/* *** END HACKERS' AREA *** */


/* YOU PROBABLY WON'T WANT TO CHANGE THESE THINGS */
/* Missing or null value */
#define MISSING -9.9999E20
/* INITIAL DEFAULT TEMPERATURE IN F AND C */
#define TEMPF 68.0
#define TEMPC 20.0

#define TRUE 1
#define FALSE 0
#define PI 3.14159265358979

/* GRAVITATIONAL CONSTANT IN m/sec^2 */
#define G 9.805

/* ** AIR DENSITY AT 0 DEGREES C at sea level DO NOT CHANGE */
/* SEE #define ALT0 FOR LAUNCHES ABOVE/BELOW SEA LEVEL */
#define RHO 1.2929

/* METERS TO FEET CONVERSION */
#define MTOF 3.28084

/* OUNCES TO GRAMS CONVERSION */
#define OZTOG 28.35

/* INCHES TO CM CONVERSION */
#define INTOCM 2.54

/* POUNDS TO NEWTONS */
#define LBTONEWT 4.45

/* GUTS OF THE STAGE STRUCTURE */
#define DEFSTRING char dddollar[SLEN],dragdollar[SLEN], dragfile[SLEN]; \
char olddragdollar[SLEN], drag1dollar[SLEN];\
char edollar[SLEN];\
char e1dollar[SLEN], engfile[SLEN];\
double m0, m01;\
float d1, d, cd[DRAGSIZE], trim, comp_time, comp_alt, impulse;\
float mp, tb, cd1, a2, tt1[ARRAYSIZE],ff[ARRAYSIZE];\
int numint, profile, signal;\
double digimp, amax, amin, maxv, a,k,v,vc, hc, h,t,curfew, tcoast;\
double airstart, airstarth, airstartv

/*
DIGITRAK WAS ORIGINALLY WRITTEN FOR ONE STAGE. IMPRUDENT USE OF
GLOBAL VARIABLES MADE CONVERSION TO MULTIPLE STAGES HELLISH. DECIDED
TO COPY SCALARS OUT OF STRUCTURES IN SUBROUTINES. DECIDED IT WASN'T SO
BAD SINCE MOST FULL THIRD-GENERATION COMPILERS COPY SCALARS INTO
SUBROUTINES (AND BACK OUT IF NECESSARY) AND USE POINTERS FOR ARRAYS.
PROGRAMS CAN ACTUALLY RUN FASTER THIS WAY - ESPECIALLY WITH VIRTUAL
STORAGE, SEGMENTED STORAGE, OR REGISTER MEMORY OPTIMIZATION.
*/

#define SCALL1(x) x = (*s).x
#define SRET1(x) (*s).x = x

/* strings take longer, but used little and makes language simpler */
#define SCALL2(x) strcpy(x, (*s).x)
#define SRET2(x) strcpy((*s).x,x)

#define SCALLA SCALL1(m01);SCALL1(m0);SCALL1(d1);SCALL1(d);SCALL1(trim);\
SCALL1(comp_time);SCALL1(comp_alt);SCALL1(impulse);SCALL1(mp);SCALL1(tb);\
SCALL1(numint);SCALL1(profile); SCALL1(a); SCALL1(k);SCALL1(v);SCALL1(a2);\
SCALL1(h);SCALL1(t); SCALL1(curfew);SCALL1(cd1);SCALL1(digimp);\
SCALL1(amax); SCALL1(amin); SCALL1(maxv);SCALL1(vc);SCALL1(hc);

#define SCALLB SCALL2(drag1dollar); SCALL2(edollar);\
SCALL2(e1dollar); SCALL2(engfile), SCALL2(olddragdollar),SCALL2(dragfile)

#define SRETA SRET1(m01); SRET1(m0); SRET1(d1); SRET1(d); SRET1(trim);\
SRET1(comp_time); SRET1(comp_alt); SRET1(impulse); SRET1(mp); SRET1(tb);\
SRET1(numint); SRET1(profile), SRET1(curfew);SRET1(digimp);\
SRET1(a); SRET1(k);SRET1(v);SRET1(a2);SRET1(vc);SRET1(hc);\
SRET1(h); SRET1(t); SRET1(cd1); SRET1(amax); SRET1(amin); SRET1(maxv)

#define SCALLC SCALL1(a); SCALL1(k);SCALL1(v); SCALL1(h); SCALL1(t);\
SCALL1(amax); SCALL1(amin); SCALL1(maxv);SCALL1(vc);SCALL1(hc);\
Scall1(digimp)

#define SRETC SRET1(a); SRET1(k);SRET1(v); SRET1(h); SRET1(t);\
SRET1(amax); SRET1(amin); SRET1(maxv);SRET1(vc);SRET1(hc);\
SRET1(digimp);


#define SRETB SRET2(drag1dollar); SRET2(edollar);\
SRET2(e1dollar); SRET2(engfile), SRET2(olddragdollar), SRET2(dragfile)

typedef struct
{
DEFSTRING;
} STAGE;


typedef struct
{
STAGE *stage[MAXSTAGE];
}STAGES;


float fff[INTERPOLATION_CUTOFF], ttt1[INTERPOLATION_CUTOFF];
float comp_time, comp_alt, temp, temp1;
double maxv, mcrit;
float fwork, rho, cd1;
int kind_switch, iswitch, irepcount, profile, nosesw=1, oldnosesw;
int BACKTRAK_STAGE, BS1=1, mach_correct=0;
float slope, standard, compare, last_compare, last_cd;
long int i, istage, inum, NUMSTAGES=1;
int repcnt, liftoff, numint, bigloop;
char ofdollar[SLEN], ddollar[SLEN], dragdollar[SLEN];
char dddollar[SLEN], dragfile[SLEN];
char drag1dollar[SLEN], olddragdollar[SLEN], adollar[SLEN], e1dollar[SLEN];
char engfile[SLEN], edollar[SLEN];
float trim = 1.0e0;
int go_on, iscaff;
double madjust, m0, m01;
float d, d1, a1, a2, a3, a4, ax, k1, tb, mp, mb;
float impulse;
double m, mconst, MACH1, MACH2;
double hc, vc, amax, amin;
double a,k,z,v,vwork, h, deltat, t, t1, digimp, f, oldf;
double mbar, oldm;
double halfdrag, oldhalfdrag, oldt1, curfew;
char old_kind[SLEN] = "P", kind[SLEN], trimdollar[SLEN];
FILE *ofp=NULL, *ifp, *dfp, *fopen();
int clear_toggle=0;/* Clears Input Screen Every Other Stage */
/***************************************************************************/
void color_routine()
/*COMPILE OPTION */
/* *** CAN INSERT ROUTINE HERE TO SET BACKGROUND COLOR *** */
{
/* NOTE
    Insert your own color routine here.
 */

#if (TURBODOS == 1)
  textmode(C80);
  textcolor(WHITE);
  textbackground(BLUE);
#endif
  return;
}
/***************************************************************************/
int editnum(instring)
/* *** EDIT FOR VALID NUMERIC *** */
char *instring;
{
int  numpt=0, numblank=0, numneg=0, notblank=0;

while((*instring !='\0') &&(*instring != '\r'))
  {
  if((*instring < '0') || (*instring >'9'))
    {

    if(*instring == ' ')
      {
      if(notblank) ++numblank;
      else numblank=1;
      }
    else  /* (*instring == ' ') */
    if(*instring == '-')
      {
      ++numneg;
      if(notblank)return(FALSE);
      }
    else if(*instring == '.')++numpt;

    else return(FALSE); /* non-numeric character found */

    notblank=1;
    }
  else notblank=1;
  ++instring;
  }
  if(numblank > 2)return(FALSE);
  if(numpt > 1)   return(FALSE);
  if(numneg > 1)  return(FALSE);
  return(TRUE);
}
/***************************************************************************/
float get_val(cp)

/* *** ACCEPT NUMERIC FROM KEYBOARD, EDIT FOR VALIDITY, AND CONVERT *** */

char *cp;
{
  float res;

  if(!editnum(cp))return(MISSING);
  sscanf(cp, "%f", &res);
  return(res);
}

/***************************************************************************/
void get_line(upcase,cp)

/* *** ACCEPT A LINE FORM THE KEYBOARD *** */

int upcase;
char *cp;
{
  int c;

  for (;;)
    {
    c = getchar();
    if ((c >= 'a') && (c <= 'z') && (upcase))
	c = c + 'A' - 'a';
    if (c == EOF || c == '\n')
      {
      *cp = 0;
      return;
      }
    *cp++ = c;
    }
}

/***************************************************************************/
float prompt(instring, x)
/* *** PROMPT USER FOR NUMERIC VALUE AND ACCEPT *** */
/* *** NUMERIC VALUE IS FLOATING POINT, BUT MORE VERSATIE WITH CASTS *** */
char *instring;
float x;
{
int go_on;
float y;
char response[SLEN];

go_on=TRUE;

while (go_on)
  {
  y=x;
  printf("%s (%f) ",instring, x);
  get_line((int) TRUE, response);
  if( response[0] != NULL) y = get_val(response);
  if (y > MISSING)
#if(TURBODOS == 0)
     {
     go_on=FALSE;
     }
#else
    {
    go_on=FALSE;
    gotoxy(1,wherey()-1);
    clreol();
    printf("%s (%f) \n\r",instring, y);
    }
  else
    {
    gotoxy(1,wherey()-1);
    clreol();
    }
#endif
  }

return(y);
}
/***************************************************************************/
char get_reply(legend)

/* *** ACCEPT SINGLE CHARACTER REPLY AND PUT INTO UPPER CASE *** */

char legend[];
{
char reply;

printf("%s ", legend);
reply = EOF;
while((reply == EOF)|(reply == '\n')) reply = getchar();

(void) getchar();
if ((reply >='a')&(reply <= 'z')) reply = reply+'A'-'a';

#if(TURBODOS == 1)
gotoxy(1,wherey()-1);
printf("%s %c \r", legend, reply);
#endif

return(reply);
}
/***************************************************************************/
void banner()

/* *** DISPLAY COPYRIGHT BANNER *** */

{
char adollar[SLEN];
CLS;
printf("\n\n");
printf(" ****************************************************************");
printf("************\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                       <<< PROGRAM DIGITRAK (tm) >>>   ");
printf("                   *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *            *** VERTICAL ROCKET TRAJECTORY SIMULATION AND ***");
printf("             *\n\r");
printf(" *      *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS");
printf(" ***       *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                           <<< VERSION 3.6 >>>");
printf("                            *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                     *** AUTHOR:  LARRY CURCIO ***");
printf("                        *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                    *** INTERNET: LC2B@DELPHI.COM *** ");
printf("                    *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                        *** SNAILMAIL ADDRESS ***  ");
printf("                       *\n\r");
printf(" *                             P.O. BOX 81124        ");
printf("                       *\n\r");
printf(" *                        PITTSBURGH, PA 15217-81124 ");
printf("                       *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" *                  *** COPYRIGHT LAWRENCE M. CURCIO ***");
printf("                    *\n\r");
printf(" *                        *** 16 NOVEMBER 1993 ***      ");
printf("                    *\n\r");
printf(" *                                                               ");
printf("           *\n\r");
printf(" ****************************************************************");
printf("************");
printf("\n\n");

get_line((int) TRUE, adollar);
return;
}
/***************************************************************************/
int interpolate( fin, tin, nin, fout, tout)

/* INTERPOLATE COARSE THRUST CURVES */

float fin[], tin[], fout[], tout[];
int nin;
{
/* Does linear interpolation in an array of points */
/* NUMBER OF LINEARLY INTERPOLATED POINTS */
#define NUMBER 400

float tb, deltat, slope, intercept, t, fwork, deltat1;
int j, i, m, nout;

tb=tin[nin];
deltat=tb/NUMBER;
nout = NUMBER;
if (deltat > DELTAT0)
  {
  nout = tb/DELTAT0;
  if (nout >= ARRAYSIZE) nout = ARRAYSIZE-1;
  deltat = tb/((float) nout);
  }
t=0.0;
nout = 0;
for(i=1; i<=nin; i++)
  {
  fwork=tin[i]-tin[i-1];
  m = (int)(fwork/deltat + 0.5);
  /* see comment before for staement */
  if (tin[i] == tb) m++;
  deltat1= fwork/((float) m);
  t = tin[i-1];
  slope = (fin[i]-fin[i-1])/fwork;
  intercept=fin[i-1];
  /* note: last point in interval not done except in last interval */
  for(j =1; j< m; j++)
    {
    fout[nout] = slope*deltat1*j + intercept;
    tout[nout] = t;
    t +=deltat1;
    nout++;
    }
  }
/*  nout--;*/
  fout[nout]=fin[nin];
  tout[nout]=tin[nin];
return(nout);
}
/****************************************************************************/
int check_for_correct(filename)

/* *** CHECK FOR MACH CORRECTION REQUEST *** */

char filename[];
{
int cswitch;
char *ipoint;
/* search for '/C' at end of string, indicating mach correction desired */
ipoint = strrchr(filename, '/');
cswitch = FALSE;
if ((ipoint != NULL) && ((strcmp(ipoint, "/C") ==0 ) ||
     (strcmp(ipoint, "/c") == 0)))
  {
  /* Found => mark for mach correction and terminate string for file name */
  cswitch = TRUE;
  *ipoint = '\0';
  }
return(cswitch);
}
/***************************************************************************/
float get_mach1(temp)
/* *** FIND MACH 1 FROM TEMPERATURE *** */
float temp;
{
#define MACH1AT0 331.45
return(MACH1AT0*sqrt(temp/273.15));
}
/***************************************************************************/
void machmaker(cd, temp, nosesw)
int nosesw;
float cd[], temp;
{
/*COMPILE OPTION*/
/*
This routine accounts for drag divergence in the manner described by
Mandell, Caporaso, and Bengen in _TOPICS IN ADVANCED MODEL ROCKETRY_,
first edition. The authors use a critical mach number of 0.9. I have
converted their model to use a critical mach number of 0.8, but I have
given the user a choice of either value.

To choose a critical mach number, #define MCRIT as either 8 (for 0.8) or 9
(for 0.9). Since some c compilers running under UNIX cannot understand a
#if comparison to a decimal constant, I had to use whole numbers instead
of the values themselves.
*/
#define MCRIT 8


/* DON'T CHANGE THIS */
#define SHARP 1

float MACH1, MACHPT9, MACHPT8, fwork, fwork1, machnum;
int IMCRIT, IMACHPT9, IMACHPT8, IMACH1PT05, IMACH1PT2, i, smaller;


MACH1 = get_mach1(temp);
MACHPT9 = 0.9*MACH1;
IMACHPT9 = (int) MACHPT9;
MACHPT8 = 0.8*MACH1;
IMACHPT8 = (int) MACHPT9;

if ((float) IMACHPT9 < MACHPT9) IMACHPT9++;
if ((float) IMACHPT8 < MACHPT8) IMACHPT8++;

IMACH1PT05 = (int) (MACH1*1.05 + .999999);
IMACH1PT2 =  (int) (MACH1*1.2  + .999999);

#if (MCRIT == 9)
  mcrit=MACHPT9;
  IMCRIT=IMACHPT9;
  if(IMACHPT9 > DRAGSIZE) return;
  fwork = cd[IMACHPT9-1];
  fwork1=(float) IMACHPT9;
#else
  mcrit=MACHPT8;
  IMCRIT=IMACHPT8;
  if(IMACHPT8 > DRAGSIZE) return;
  fwork = cd[IMACHPT8-1];
  fwork1=(float) IMACHPT8;
#endif

if(nosesw == SHARP)
  {
      /* Sharp Nosed Vehicles */

  if(IMACH1PT05 < DRAGSIZE) smaller=IMACH1PT05;
  else smaller = DRAGSIZE;

  for (i=IMCRIT; i<smaller; i++)
      {
      machnum=fwork1/MACH1;
      fwork1 +=1.0;
      /* (0.9 =< M =< 1.05)  */
#if(MCRIT == 9)
      cd[i] = fwork*(1.0 + 35.5*(machnum-0.9)*(machnum-0.9));
#else
      cd[i] = fwork*(1.0 + 12.78*(machnum-0.8)*(machnum-0.8));
#endif

      }
  if(smaller < DRAGSIZE) for (i=smaller; i<DRAGSIZE; i++)
      {
      machnum=fwork1/MACH1;
      fwork1 +=1.0;
      /*  1.05 =< M =< 2.0  */
      cd[i]=fwork*(1.27 + 0.53*exp(-5.2*(machnum-1.05)));
      }
  }
else

  {
	   /* Round-Nosed Vehicles */
  if(IMACH1PT2 < DRAGSIZE) smaller=IMACH1PT2;
  else smaller = DRAGSIZE;
for (i=IMCRIT; i<smaller; i++)
      {
      machnum=fwork1/MACH1;
      fwork1 +=1.0;
      /* (0.9 =< M =< 1.2) */
#if(MCRIT == 9)
      cd[i] = fwork * (1.0 + 4.88*exp(1.1*log(machnum-0.9)));
#else
      cd[i] = fwork * (1.0 + 3.556208468*exp(1.1*log(machnum-0.8)));
#endif
      }
  if(smaller < DRAGSIZE) for (i=smaller; i<DRAGSIZE; i++)
      {
      machnum=fwork1/MACH1;
      fwork1 +=1.0;
      /* 1.2 =< M =< 2.0 */
      cd[i] = fwork * (2.0 + 0.3*exp(-5.75*(machnum-1.2)));
      }
  }
  return;
}
/***************************************************************************/
void set_next_values(kind_switch, config)

/* *** SET DEFAULT VALUES FOR NEXT FLIGHT (RUN) *** */

STAGES config;
int kind_switch;
/* Set values for next time around */
{
int i;

  if(kind_switch != 1) (*config.stage[NUMSTAGES-1]).comp_time = t;
  if(kind_switch != 2)
    {
#if(ALTUNITS == 'F')
     (*config.stage[NUMSTAGES-1]).comp_alt = h * MTOF;
#else
     (*config.stage[NUMSTAGES-1]).comp_alt = h;
#endif
for (i=0; i < NUMSTAGES; i++)
  {
    if(madjust != 0.0)
    {
    (*config.stage[i]).m0 += madjust;
#if(INUNITS == 'M')
    (*config.stage[i]).m01=(*config.stage[i]).m0*1000.0;
#else
    (*config.stage[i]).m01=(*config.stage[i]).m0*1000.0/OZTOG;
#endif
    }
  (*config.stage[i]).tcoast=0.0;
  }
madjust=0.0;
return;
    }
}
/***************************************************************************/
int query_cd(kind_switch, istage, s)

/* *** PROMPTS USED FOR DRAG DATA AND PARSES THE REPLY *** */
/* *** DRAG MAY BE A Cd VALUE, A SHAPE CONSTANT, OR A DRAG PROFILE *** */
/* *** A PROFILE MAY BE READ PER SE OR CORRECTED FOR DRAG DIVERGENCE *** */
/* *** NUMBER ENTERED HERE AND IN DIAMETER => Cd VALUE *** */
/* *** NUMBER ENTERED HERE AND NONE IN DIAMETER => SHAPE CONSTANT */
/* *** STRING ENTERED HERE => FILE *** */
/* *** STRING HERE, NO # IN DIAMETER => TAKE DIAMETER FROM PROFILE *** */
/* *** STRING HERE, # IN DIAMETER => TAKE DIAMETER FROM # *** */

int kind_switch, istage;
STAGE *s;
{
int signal;

/* If not the backtracked stage, treat as a simple prediction */
if (istage != BACKTRAK_STAGE)kind_switch=0;
profile=(*s).profile;

signal = (*s).signal;
if ((kind_switch == 0) | (profile))
  {
  printf("Drag Coefficient or File (%s) ",(*s).dragdollar);
  get_line((int) FALSE, (*s).drag1dollar);
  if((*s).drag1dollar[0] != NULL) strcpy((*s).dragdollar,(*s).drag1dollar);
  else signal=(*s).signal;
  if((strcmp((*s).dragdollar, (*s).olddragdollar) != 0)|
     (nosesw !=oldnosesw))
    {
    if(((*s).dragdollar[0] == '.') | (((*s).dragdollar[0] >= '0' ) &
			    ((*s).dragdollar[0] <= '9')))
      {
      (*s).cd1 = get_val((*s).dragdollar);
      if ((*s).cd1 == 0.0) (*s).cd1 = DELTA;
      for (i = 0;i<DRAGSIZE;i++) (*s).cd[i] = (*s).cd1;
      machmaker((*s).cd, temp, nosesw);
      (*s).a2 = 0.0;
      signal=TRUE;
      }
    else if((strcmp((*s).dragdollar, (*s).olddragdollar) != 0)|(oldnosesw != nosesw))
      {
      strcpy((*s).dragfile,(*s).dragdollar);
      mach_correct = check_for_correct((*s).dragfile);
      dfp = fopen((*s).dragfile, "r"); /* clean this up when debugeged */
      if (dfp == NULL) signal = FALSE;
      else
        {
        fgets(dddollar, SLEN, dfp);/* toss first line */
	fscanf(dfp, "%f %f %f %f", &a1, &(*s).a2, &a3, &a4); /* a2 is diameter*/
        for (i=2; i<= NUMDOC; i++) fgets(dddollar, SLEN, dfp); /* toss rest */
	i= -1;
	while ((fscanf(dfp, "%f", &(*s).cd[i]) != EOF) & (++i < DRAGSIZE));
        fclose(dfp);
	signal=TRUE;
        i--;
        fwork = (*s).cd[i];

        while (++i < DRAGSIZE) (*s).cd[i] = fwork;


	if(mach_correct)machmaker((*s).cd, temp, nosesw);
        }
      }
    strcpy((*s).olddragdollar, (*s).dragdollar);
    }
  } /* if kind_switch == 0 or ... */
  else (*s).olddragdollar[0] = NULL;

#if(TURBODOS == 1)
if ((kind_switch == 0) | (profile))
  {
  gotoxy(1,wherey()-1);
  clreol();
  printf("Drag Coefficient or File (%s) \n\r",(*s).dragdollar);
  }
#endif
(*s).signal = signal;
return(signal);
}

/***************************************************************************/
void print_stage(kind_switch, istage, s)

/* *** PRINTS DATA PERTAINING TO A SPECIFIED STAGE *** */

int kind_switch, istage;
STAGE *s;
{

/* #if EXCEPTION */
/*
IN COMPILERS WITHOUT #if'S, WILL HAVE TO CHOOSE THE #defineS THAT CONFORM
TO CHOSEN TEMPERATURE UNITS. COMMENT OUT THE OTHER DEFINITIONS.
*/

#if (TEMPUNITS == 'F')
#define TEMPADJ *9.0/5.0 + 32.0
char tempunits='F';
#else
#define TEMPADJ +0.0
char tempunits='C';
#endif

/* *** END #if EXCEPTION *** */

#if(TURBODOS == 1)
gotoxy(1,wherey());
clreol();
#endif

if(kind_switch == 3)
  {
  CLS;

  if (NUMSTAGES > 1) printf("\nFor Stage %i \n\n\r",istage+1);

  printf(
  "      ** Estimated Optimal Performance at %f Deg %c **\n\n\r",
     (temp-273.15) TEMPADJ, tempunits);
  printf("Optimal Mass of Stage  = %f Grams;  or %f Ounces\n\n\r",
     1000.0*((*s).m0+madjust), 1000.0*((*s).m0+madjust)/OZTOG);
  kind_switch=0;
  }
  else
  {
  if(NUMSTAGES > 1)
    {
    CLS;
    printf("\nFor Stage %i:\n\n\r",istage+1);
    }
  }
/* Only one stage can be backtracked at a time */
if(istage != BACKTRAK_STAGE)kind_switch=0;

if(((*s).curfew > 0.0)||(istage == NUMSTAGES-1))
  printf("Maximum Altitude (m):     %f   (ft):     %f \n\r",(*s).h,
	(*s).h*MTOF);
printf("Cutoff Altitude  (m):     %f   (ft):     %f \n\r",(*s).hc,
	(*s).hc*MTOF);
if((*s).airstart != 0.0)
  printf("Airstart Altitude(m):     %f   (ft):     %f \n\r",(*s).airstarth,
	(*s).airstarth*MTOF);

printf("\n\r");
printf("Maximum Velocity  (m/sec): %f  (ft/sec): %f \n\r", (*s).maxv,
	(*s).maxv*MTOF);
printf("Average Velocity  (m/sec): %f   (ft/sec): %f \n\r", (*s).h/(*s).t,
        (*s).h*MTOF/(*s).t);
if((*s).airstart != 0.0)
  printf("Airstart Velocity (m/sec): %f   (ft/sec): %f \n\r",(*s).airstartv,
	(*s).airstartv*MTOF);
if(((*s).curfew > 0.0)&&(istage != NUMSTAGES-1)&&(NUMSTAGES > 1))
  printf("Staging Velocity  (m/sec): %f   (ft/sec): %f \n\r", (*s).v,
	(*s).v*MTOF);
printf("Maximum Acceleration (Gs): %f   Minimum:  %f \n\r",(*s).amax/G,
	  (*s).amin/G);
printf("\n\r");
printf("Maximum Mach Number (dimensionless) (%f) \n\r",
	  (*s).maxv/get_mach1(temp));
printf("\n\r");
printf("Flight Time  (sec): %f     Coast Time (sec): %f \n\r",
	   (*s).t,(*s).tcoast);

if(kind_switch)
  {
  if((*s).profile) printf("Trim Factor = %f      ",(*s).trim);
  else printf("Drag Coefficient = %f      ",(*s).trim);
  printf("Sea Level K (kg*m/sec^2) = %f\n",(*s).k*(*s).trim*
          (*s).cd[(int)((*s).h/(*s).t + 0.5)]);
  }

printf("\n\r");
ax=PI*(*s).d*(*s).d/4.0;
if(ax > 0.0) printf("Cross-Sectional Area (cm^2): %f \n\r",ax*10000);
printf("Recorded Impulse (n-sec): %f  Digital Impulse(n-sec): %f \n\r",
	 (*s).impulse, (*s).digimp);
printf("\n\r");
if (istage != NUMSTAGES-1) getchar();
return;
}
/***************************************************************************/
void prtout(kind_switch, config)

/* *** EXECUTIVE ROUTINE TO CALL STAGE OUTPUT ROUTINE, print_stage *** */

int kind_switch;
STAGES config;
{
int istage;
for (istage=0; istage<NUMSTAGES; istage++)
  {
  print_stage(kind_switch, istage, config.stage[istage]);
  }
return;
}
/***************************************************************************/
void get_stage_data(istage, kind_switch, s)

/* *** STAGE-SPECIFIC KEYBOARD INPUT ROUTINE *** */

int istage, kind_switch;
STAGE *s;
{
int ii, go_on1;

/* copy all the stuff into the space */

/*
SCALLA;
SCALLB;
*/

profile=(*s).profile;
/* Backtrack only one stage at a time */
if(istage != BACKTRAK_STAGE)kind_switch=0;

if (NUMSTAGES != 1)
  {
  /* TURBO C FLAGS THIS AS A WARNING FOR THE UPDATE OF clear_toggle*/
  /* "I meant to do that." ... PeeWee Herman (RIP) */

  if((clear_toggle = 1-clear_toggle) && (istage !=  0))CLS;
  printf("\nFor Stage %i \n\r",istage+1);
  }

if (istage != 0)(*s).airstart=(double) prompt("Air Start Delay (sec): ",
					 (float) (*s).airstart);
go_on1 = TRUE;
while(go_on1)
   {

#if(INUNITS == 'M')
   (*s).m01=(double) prompt("Mass in Grams ",(float) (*s).m01);
   (*s).m0=(*s).m01/ 1000.0;
#else
   (*s).m01=(double) prompt("Mass in Ounces ",(float) (*s).m01);
   (*s).m0=(*s).m01*OZTOG/ 1000.0;
#endif
   if((*s).m0 > 0.0) go_on1=FALSE;
#if(TURBODOS == 1)
   else gotoxy(1,wherey()-1);
#endif
   }

#if(INUNITS == 'E')
(*s).d1=prompt("Maximum Diameter in Inches ",(*s).d1);
(*s).d=(*s).d1*INTOCM/100.0;
#else
(*s).d1=prompt("Maximum Diameter in CM ",(*s).d1);
(*s).d=(*s).d1/100.0;
#endif

/* The Cd routine is isolated in query_cd() Returns TRUE when all OK */
if ((kind_switch == 0) | (profile))
while(!query_cd(kind_switch, istage, s))
#if(TURBODOS == 1)
gotoxy(1,wherey()-1);
clreol();
#else
;
#endif

/* Trim factor is a factor by which all Cd's in a profile are multiplied. */
/* Trim factor is varied instead of Cd in backtracking with profiles */

if (kind_switch == 0)
  {
  (*s).trim=prompt("Trim Factor ",(*s).trim);

  }/* if kind_switch == 0 */
else
  {
  if(! profile)
    {
    for (i = 0;i<DRAGSIZE;i++) (*s).cd[i] = 1.0;
    machmaker((*s).cd, temp, nosesw);
    (*s).olddragdollar[0]='\0';
    }
  }

go_on1 = TRUE;
while(go_on1)
  {
  printf("Motor (%s) ",(*s).edollar);
  get_line((int) FALSE, (*s).e1dollar);
  /* don't read file again if same as before */
  if(((*s).e1dollar[0] != NULL) | (ifp == NULL))
    {
    strcpy((*s).edollar,(*s).e1dollar);
    strcpy((*s).engfile,(*s).edollar);
    strcat((*s).engfile,".eng");
    (*s).mp=0.0;
    ifp=fopen((*s).engfile,"r");
    if (ifp == NULL)  go_on1 = TRUE;
    else
      {
      go_on1 = FALSE;

      fgets((*s).dddollar, SLEN, ifp);

      fscanf(ifp,"%f %f %f %i", &(*s).impulse, &(*s).mp, &(*s).tb, &(*s).numint);
      (*s).mp=(*s).mp/1000.0;
      if((*s).numint < INTERPOLATION_CUTOFF)
      {
      /* Thrust Curve Too Coarse. Must Interpolate Linearly */
      for(ii=0;ii<= (*s).numint; ii++)
	{
#if (TOGGLE1)
	fscanf(ifp,"%f %f ", &ttt1[ii], &fff[ii]);
#else
	fscanf(ifp,"%f %f %f", &fff[ii], &ttt1[ii], &fwork);
#endif
	}
	(*s).numint=interpolate( fff, ttt1, (*s).numint, (*s).ff, (*s).tt1);
      }
      else
       {
	for(ii=0;ii<= (*s).numint; ii++)
	  {
#if(TOGGLE)
	  fscanf(ifp,"%f %f %f", (&s).tt1[ii], (&s).ff[ii], &fwork);
#else
	  fscanf(ifp,"%f %f %f", &(*s).ff[ii], &(*s).tt1[ii], &fwork);
#endif
	  }
       }
      fclose(ifp);
      }
    }/* if e1dollar[0] != NULL */
    else if(ifp != NULL)go_on1 = FALSE;
    if((*s).mp >= (*s).m0)
      {
      go_on1=TRUE;
      printf("Propellant Mass Exceeds Launch Mass \n\r");
      printf ("Press Any Key To Continue ");
      (void)getchar();
#if(TURBODOS == 1)
      gotoxy(1,wherey());
      clreol();
       gotoxy(1,wherey()-1);
      clreol();
      gotoxy(1,wherey()-1);
      clreol();
#else
      printf("\n\r");
#endif
      }
#if(TURBODOS == 1)
  gotoxy(1,wherey()-1);
  clreol();
  if(!go_on1)printf("Motor (%s) \n\r",(*s).edollar);
#endif

  }

  if(istage != NUMSTAGES-1)
    (*s).curfew = prompt("Staging Delay (sec): ",(*s).curfew);

printf("\n\r");


if((iswitch))
  {
#if(INUNITS == 'M')
  fprintf(ofp,"# M0= %f  Diameter= %f \n",(*s).m0, (*s).d);
#else
  fprintf(ofp,"# M0= %f  Diameter= %f \n",1000*(*s).m0/OZTOG,
	   100*(*s).d/INTOCM);
#endif
  fprintf(ofp,"# Cd= %s  Trim= %f Engine= %s \n",(*s).dragdollar, (*s).trim, (*s).edollar);
#if(INUNITS == 'M')
  fprintf(ofp,"# Tb= %f  Impulse= %f Mp= %f \n",(*s).tb, (*s).impulse, (*s).mp);
#else
  fprintf(ofp,"# Tb= %f  Impulse= %f Mp= %f \n",(*s).tb,
  (*s).impulse/LBTONEWT, 1000*(*s).mp/OZTOG);
#endif
  fprintf(ofp,"# Staging Delay =%f  Airstart Delay= %f\n#\n",(*s).curfew,
							  (*s).airstart);
  }
/*
SRETA;
SRETB;
*/
return;
}
/***************************************************************************/
int get_data(config)

/* *** TAKES WHOLE-ROCKET DATA FROM THE KEYBOARD *** */
/* (STAGE-SPECIFIC DATA TAKEN IN get_stage_data *** */

STAGES config;
{
int kind_switch, go_on, istage;
float OLDNUMSTAGES,ANUMSTAGES;
CLS;

clear_toggle=0;
madjust=0.0;
OLDNUMSTAGES=NUMSTAGES;
go_on = TRUE;
while(go_on)
  {
  NUMSTAGES=OLDNUMSTAGES;
  ANUMSTAGES = prompt("Number of Stages ",(float)NUMSTAGES);
  if((ANUMSTAGES > 0) && (ANUMSTAGES <= MAXSTAGE)) go_on=FALSE;
  if(fabs(ANUMSTAGES) > MAXSTAGE) NUMSTAGES=MAXSTAGE+1;
  else NUMSTAGES = (int) ANUMSTAGES;
  go_on=FALSE;
  if(ANUMSTAGES != (float)NUMSTAGES) go_on=TRUE;
  if((NUMSTAGES > MAXSTAGE)||(NUMSTAGES <1)) go_on=TRUE;
  }

if(NUMSTAGES == 1) CLS;

go_on = TRUE;
while(go_on)
  {
  /* TP or OP gets you backtracking with profiles */
  printf("'T' = Temporal, 'O' = Optical, 'P' = Predict (%s) ", old_kind);

  get_line((int) TRUE, kind);
  if(kind[0] != NULL) strcpy(old_kind, kind);

  go_on = FALSE;
  if(strcmp(old_kind, "P") == 0)
    {
    kind_switch = 0;
    }
  else if(strcmp(old_kind, "T") == 0)
    {
    kind_switch = 1;
    profile = FALSE;
    }
  else if(strcmp(old_kind, "O") == 0)
    {
    kind_switch = 2;
    profile = FALSE;
    }
  else if(strcmp(old_kind, "TP") == 0)
    {
    kind_switch = 1;
    profile = TRUE;
    }
  else if(strcmp(old_kind, "OP") == 0)
    {
    kind_switch = 2;
    profile = TRUE;
    }
  else go_on = TRUE;
#if(TURBODOS == 1)
  CLS;
  if(!go_on)
    printf("'T' = Temporal, 'O' = Optical, 'P' = Predict (%s)\r\n", old_kind);
#endif
  }

iswitch = 0;
if(kind_switch == 0)
  {
  printf("Output File (<CR> = None) ");
  get_line((int) FALSE, ofdollar);
  if(ofdollar[0] != NULL)
    {
     iswitch = 1;
     ofp = fopen(ofdollar, "w");
    }
#if(TURBODOS == 1)
gotoxy(1,wherey()-1);
clreol();
printf("Output File (%s)\n\r", ofdollar);
#endif

  }

else if(NUMSTAGES > 1)
  {
  BACKTRAK_STAGE=9;
  while((BACKTRAK_STAGE < 0 )|| (BACKTRAK_STAGE > NUMSTAGES-1))
    {
    BS1=  (int) prompt("Stage to backtrak", (float) BS1);
    BACKTRAK_STAGE=BS1-1;
    (*config.stage[BACKTRAK_STAGE]).profile=profile;
    }
  }

else
  {
  BACKTRAK_STAGE=0;
  (*config.stage[BACKTRAK_STAGE]).profile=profile;
  printf("\n\r");
  }
#if TEMPUNITS == 'C'
temp1=prompt("Temperature in Deg C ", temp1);
temp = temp1 + 273.15;
#else
temp1=prompt("Temperature in Deg F ", temp1);
temp = (temp1 - 32.0)*5.0/9.0 + 273.15;
#endif

rho = RHO * 273.15/temp;

/* Nose type is used only in supersonic calculations */
oldnosesw=nosesw;
nosesw=5;
while((nosesw !=1)&(nosesw !=2))
  {
  nosesw=oldnosesw;
  nosesw=
  (int)prompt("Nose Type  1=Sharp 2=Round ", (float)nosesw);
  }

  if (kind_switch == 1)
    {
    (*config.stage[NUMSTAGES-1]).comp_time=
	    prompt("Ascent Time ",(*config.stage[NUMSTAGES-1]).comp_time);
    standard = (*config.stage[NUMSTAGES-1]).comp_time;
    }
    else if(kind_switch == 2)
      {
      (*config.stage[NUMSTAGES-1]).comp_alt=
#if(ALTUNITS == 'M')
      prompt("Altitude in Meters ",(*config.stage[NUMSTAGES-1]).comp_alt);
      standard = (*config.stage[NUMSTAGES-1]).comp_alt;
#else
      prompt("Altitude in Feet ",(*config.stage[NUMSTAGES-1]).comp_alt);
      standard = (*config.stage[NUMSTAGES-1]).comp_alt / MTOF;
#endif
      }
if(iswitch)
  {
  fprintf(ofp,"# Number of Stages = %i\n", (int) NUMSTAGES);
#if(INUNITS == 'M')
  fprintf(ofp,"# Units = Metric \n");
  fprintf(ofp,"# Temperature (Deg C) = %f  ",temp-273.15);
#else
  fprintf(ofp, "# Units = English \n");
  fprintf(ofp,"# Temperature (Deg F) = %f  ",(temp-273.15)*9/5 + 32);
#endif
  if(nosesw == 1)fprintf(ofp,"Nose = Sharp\n#\n");
  else fprintf(ofp,"Nose = Round\n#\n");
  }

for (istage=0; istage< NUMSTAGES; istage++)
  get_stage_data(istage, kind_switch, config.stage[istage]);

#if (TURBODOS == 1)
printf("Press Any Key To Continue");
get_line(FALSE,adollar);
gotoxy(1,wherey()-1);
textcolor(RED+BLINK);
clreol();
printf("Calculating");
textcolor(WHITE);
#endif
return(kind_switch);
}
/***************************************************************************/
void init(config)

/* *** INITIALIZES A WHOLE FLIGHT *** */

STAGES config;
{
v = DELTA;
t=h=0.0;
repcnt=0;
m=(*config.stage[0]).m0 + madjust;
oldhalfdrag=maxv=0.0;
liftoff = FALSE;
f=0;
oldf=0.0;
oldt1=0.0;
}
/***************************************************************************/
/* ** FM COAST TIME ** */
/* ** FOR LAST COAST INTERVAL TO MAKE BACKTRACKED TIME CONTINUOUS ** */
/* ** WITHOUT THIS, INTERVALS MUST BE VERY SMALL ** */
double FM_coast_time(vb, z)
double vb, z;
/* z is K/mb */
{
  return (sqrt(1.0 / (G * z)) * atan(vb * sqrt(z / (G))));
}
/***************************************************************************/
/* ** FM COAST ALTITUDE ** */
/* ** FOR LAST COAST INTERVAL TO MAKE BACKTRACKED ALTITUDE PRECISE ** */
double FM_coast_alt(vb, z)
double vb, z;
/* z is K/mb */
{
  return(log(z * vb * vb / G + 1) / (z + z));
}
/***************************************************************************/
void boost(s)

/* *** BOOST PHASE ROUTINE *** */

STAGE *s;
{
register double work;
double fwork, work2, netf, newv;
double sqrtg, sqrtz, sqrtzovg;
double H, m0hold, Gstar, zstar;
float *CD, *FF, *TT;

SCALLA;

/* ** copy here for faster addressing in loops ** */
CD = (*s).cd;
FF=(*s).ff;
TT=(*s).tt1;

maxv=0;

/* For Optimal Mass, adjust m0 */
m0hold=m0;
m0 += madjust;
oldm=m0;

oldt1=0.0;
if(d == 0.0)
  {
  d = a2/100.0;
  d1=a2;
#if(INUNITS == 'E')
  d1 = d1/INTOCM;
#endif
  }

ax=PI * d * d /4.0;
if( ax != 0.0) k = trim * ax * rho / 2.0;
else k = trim / 10000.0;

digimp=0.0;
mconst=mp / impulse;

if (iswitch)
  {
  ++repcnt;
  if (repcnt >= REPNUM)
    {
#if(INUNITS == 'M')
    fprintf(ofp,"%f  %f  0.0  0.0  %f \n",h, v, t);
#else
    fprintf(ofp,"%f  %f  0.0  0.0  %f \n",h*MTOF, v*MTOF, t);
#endif
    repcnt=0;
    }
  }/* if iswitch */


for(inum=0; inum <=numint; inum++)
{
  oldf=f; f = FF[inum]; t1=TT[inum];
  fwork=0.5*(f+oldf);
  deltat = t1-oldt1;
  oldt1=t1;
  t += deltat;
  digimp = digimp + fwork * deltat;
  m = m0 - digimp * mconst;
  mbar = (m+oldm)/2.0;
  oldm=m;
  if (fwork > mbar * G) {liftoff = TRUE;}
  if (liftoff)
    {
#if (CONSTANT_RHO)
    z = k*CD[(int)(v+0.5)];
#else
    H=h+ALT0;
    z=k*CD[(int)(v+0.5)]*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14))));
#endif

    netf=fwork-mbar*G;
    if(netf > 0.0)
      {
        /* STRAIGHT EXTENDED FM BOOST-PHASE INTEGRALS */
        work=cosh(deltat*sqrt(z*netf)/mbar)+
        v*sqrt(z/netf)*sinh(deltat*sqrt(z*netf)/mbar);
        h+=mbar*log(work)/z;

	work=tanh(deltat*sqrt(z*netf)/mbar);
        work2=v*sqrt(z/netf);
	v=sqrt(netf/z)*(work+work2)/(1+work*work2);
      }
    else
      {
       /* COAST PHASE INTEGRALS WITH ADJUSTED G FOR  netf < 0.0 */
       Gstar = -netf/(mbar);
       sqrtg = sqrt(Gstar);
       zstar=z/mbar;
       sqrtz=sqrt(zstar);
       sqrtzovg=sqrtz/sqrtg;

/*     set work to tan(arg1)    */
/*     where arg1=atan(v*sqrtzovg) - deltat*sqrtz*sqrtg; */
/*     using identity for tan(alpha-beta) */
       work=tan(sqrtz*sqrtg*deltat);
       work=(v*sqrtzovg-work)/(1.0+v*sqrtzovg*work);
       newv=work/sqrtzovg;

/*     work is now the cos^2 arg1 */
/*     cos(theta)^2 = 1/[1 +tan(theta)^2]  */
       work=1.0/(work*work+1.0);
       h += 0.5*log((v*v*zstar/Gstar+1.0)*work)/zstar;
       work=0.5*(v+newv);
       v=newv;
       }
    a = (netf-z*work*work)/mbar;
    if (v > maxv) maxv = v;
    if (a > amax) amax = a;
    if (a < amin) amin = a;
    hc=h;
    vc = v;
    if (iswitch)
      {
      ++repcnt;
      if (repcnt >= REPNUM)
	{
#if(INUNITS == 'M')
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h, v, a,
	    trim*CD[(int)(v+0.5)], t);
#else
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h*MTOF, v*MTOF, a*MTOF,
	    trim*CD[(int)(v+0.5)], t);
#endif
	repcnt=0;
	}
      }/* if iswitch */
    }/* if (liftoff);*/
}


    if (iswitch)
      {
      if (repcnt != 0)
	{
#if(INUNITS == 'M')
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h, v, a,
	    trim*CD[(int)(v+0.5)], t);
#else
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h*MTOF, v*MTOF, a*MTOF,
	    trim*CD[(int)(v+0.5)], t);
#endif
	repcnt=0;
	}
      }/* if iswitch */

/* unadjust m0 */
m0=m0hold;
SRETA;
return;
}
/***************************************************************************/
double coast(kind_switch, s, delay, curfew_switch)

/* *** COAST PHASE ROUTINE *** */
/* *** CURFEW IS THE STAGING DELAY TIME, IF ANY *** */

int kind_switch, curfew_switch;
double delay;
STAGE *s;
{
register double work;
double H, newv, sqrtg, sqrtz, sqrtzovg, curfew1;
float m0hold, *CD;
int go_on, repnum,REPNUM1;
SCALLA;

/* ** copy here for faster addressing in loops ** */
CD=(*s).cd;

m0hold=m0;
m0 +=madjust;

/* DECREASE REPNUM TO REPNUM1 FOR WIDER INTERVAL */
/* DELTAT0 IS WIDE INTERVAL FOR LOW ACCELERATION COAST */
/* DELTAT1 IS NARROW INTERVAL FOR HIGH ACCELERATION COAST */

REPNUM1= (int) 0.999999 + REPNUM*(DELTAT1/DELTAT0);


if(curfew_switch)curfew1=99999e9;
else curfew1=delay;

if(d == 0.0)
  {
  d = a2/100.0;
  d1=a2;
#if(INUNITS == 'E')
  d1 = d1/INTOCM;
#endif
  }

ax=PI * d * d /4.0;
if( ax != 0.0) k = trim * ax * rho / 2.0;
else k = trim / 10000.0;

curfew1=curfew1+t;
mb=m0-mp;
k=k/mb;
repcnt=0;
sqrtg=sqrt(G);
go_on = TRUE;
while(go_on)
    {

    if((v < mcrit)||(v > MACH2))
      {
      deltat=DELTAT0;
      repnum=REPNUM1;
      }
    else
      {
      deltat=DELTAT1;  /* SMALLER INTERVAL FOR DIVERGENCE */
      repnum=REPNUM;
      }
      if (t+deltat > curfew1)
	 {
	 deltat = curfew1-t;
	 go_on = FALSE;
	 }
#if (!CONSTANT_RHO)
    H=h+ALT0;
    z=CD[(int)(v+0.5)]*k*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14))));
#else
    z=CD[(int)(v+0.5)]*k;
#endif
    sqrtz=sqrt(z);

    sqrtzovg=sqrtz/sqrtg;
/*    set work to tan(arg1)    */
/*    where arg1=atan(v*sqrtzovg) - deltat*sqrtz*sqrtg; */
/* using identity for tan(alpha-beta) */
    work=tan(sqrtz*sqrtg*deltat);
    work=(v*sqrtzovg-work)/(1.0+v*sqrtzovg*work);


    if(work < 0.0)
      {
       go_on = FALSE;
       if(v != 0.0)
	 {
	  h += FM_coast_alt(v, z);
	  t += FM_coast_time(v, z);
	  v=0.0;
	 }
      }/* if(arg1 < 0.0) */
    else
      {
/*     newv = tan(arg1)*sqrt(g/z)  */
	newv= work/sqrtzovg;
       /* work is now the cos^2 arg1 */
       /* cos(theta)^2 = 1/[1 +tan(theta)^2]  */
       work=1.0/(work*work+1.0);
       h += 0.5*log((v*v*z/G+1.0)*work)/z;
       work=0.5*(v+newv);
       v=newv;

       a= -z*work*work - G;
       if (a < amin) amin = a;
       t += deltat;
       if(iswitch)
        {
	++repcnt;

	if (repcnt >= repnum)
	   {
#if(INUNITS == 'M')
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h, v, a,
	    trim*CD[(int)(v+0.5)], t);
#else
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h*MTOF, v*MTOF, a*MTOF,
	    trim*CD[(int)(v+0.5)], t);
#endif
	   repcnt=0;
	   }
        }/* if iswitch */
      } /* else (v still positive) */
    } /* else (no boost) */

/* FORCE LAST RECORD OUTPUT FOR WIDE INTERVALS */
if(iswitch)
	{
	if (repcnt != 0) /* OMIT REDUNDANT LAST RECORDS */
	   {
#if(INUNITS == 'M')
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h, v, a,
	    trim*CD[(int)(v+0.5)], t);
#else
	fprintf(ofp,"%f  %f  %f  %f  %f \n",h*MTOF, v*MTOF, a*MTOF,
	    trim*CD[(int)(v+0.5)], t);
#endif
	   repcnt=0;
	   }
        }/* if iswitch */

(*s).tcoast=t-(*s).t;
m0=m0hold;
SRETA;

if (kind_switch == 1) return (t);
else return(h);
}
/***************************************************************************/
float trajectory(config, kind_switch)

/* KNITS BOOST AND COAST PHASES TOGETHER TO MAKE A WHOLE FLIGHT */
/* EFFECTS SIMPLE STAGING DELAYS AND AIRSTART DELAYS */

int kind_switch;
STAGES config;
{
int istage;

init(config);

if(iswitch)fprintf(ofp,"#  Alt       Vel       Acc         Cd      T\n");

for(istage=0; istage<NUMSTAGES; istage++)
  {
  (*config.stage[istage]).v=v;
  (*config.stage[istage]).h=h;
  (*config.stage[istage]).t=t;
  (*config.stage[istage]).amax= MISSING; /* BIG NEGATIVE NUMBER */
   (*config.stage[istage]).amin= - MISSING; /* BIG POSITIVE NUMBER */

  if ((*config.stage[istage]).airstart != 0.0)
     {
     (*config.stage[istage]).airstarth = coast(2, config.stage[istage],
		(*config.stage[istage]).airstart, 0);
     (*config.stage[istage]).v=(*config.stage[istage]).airstartv=v;
     (*config.stage[istage]).h=h;
     (*config.stage[istage]).t=t;
     }
  if ((*config.stage[istage]).numint >0) boost(config.stage[istage]);
  else (*config.stage[istage]).amax=0.0;
  if ((*config.stage[istage]).curfew != 0.0)
    {
    (*config.stage[istage]).v=v;
    (*config.stage[istage]).h=h;
    (*config.stage[istage]).t=t;
    (void) coast(kind_switch, config.stage[istage],
	     (*config.stage[istage]).curfew, 0);
    }
  (*config.stage[istage]).v=v;
  (*config.stage[istage]).h=h;
  (*config.stage[istage]).t=t;

  }
return(coast(kind_switch, config.stage[NUMSTAGES-1], 99999e9, 1));
}
/***************************************************************************/
 float optmass(config)
 STAGES config;
{
double mdd, muu, midm;
double delta, pnt, pntlo, pnthi, maxm;
int iopt, isw, i;
float totimpulse;

 /* ** OPTIMAL MASS ROUTINE ** */
 /* ** Mdd IS M CORRESPONDING TO DOWNWARD SLOPE ** */
 /* ** Muu IS M CORRESPONDING TO UPWARD SLOPE ** */

#define DELTAMAS 0.0001
#define MAXITMAS 100


 maxm=(*config.stage[0]).m0;

 totimpulse=0.0;

 for(i=0; i<NUMSTAGES ;i++)
   {
   totimpulse +=(*config.stage[i]).impulse;
   }

 if (tb > .0001) mdd = totimpulse / (2.0 * G * tb);
 else mdd=totimpulse / (2.0 * G * (*config.stage[NUMSTAGES-2]).tb);

 muu=(*config.stage[0]).m0 - (*config.stage[NUMSTAGES-1]).m0
	 +(*config.stage[NUMSTAGES-1]).mp;
/* printf("muu= %f \n\r",muu*1000.0); */

 delta = .0001;
 iopt = 1;
 isw = 1;
 printf("\n\r");
#if(TURBODOS == 1)
clreol();
#endif
 while (isw)
 {

   midm = (mdd + muu) / 2.0;
   printf("trying m= %f Grams;  or  %f Ounces\n\r",
	   midm*1000.0, midm*1000.0/OZTOG);
   madjust = midm - maxm - delta; pntlo=trajectory(config, 2);
   madjust = midm - maxm + delta; pnthi=trajectory(config, 2);
   madjust = midm - maxm; pnt=trajectory(config, 2);
   pntlo=pnt-pntlo;
   pnthi=pnthi-pnt;
   if ( (pntlo < 0.0) & (pnthi < 0.0) ) mdd = midm;
   else if ( (pntlo > 0.0) & (pnthi > 0.0) ) muu = midm;
   else
   {
     isw = 0;
     for(i=0; i < NUMSTAGES; i++)
       {
       if(((*config.stage[i]).airstartv < 0.0)||
	  (((*config.stage[i]).v < 0.0)&& (i < NUMSTAGES -1)))
	  {
	  isw=1;
	  mdd=midm;
	  }
       }
   }

   iopt = iopt + 1;
   if (iopt > MAXITMAS) isw = 0;
   if(fabs(mdd - muu)/midm < DELTAMAS) isw = 0;
 }
if (iopt > MAXITMAS) printf("*** NO CONVERGENCE *** \n\r");
/* madjust=0.0;  */
return(midm);
}
/***************************************************************************/
    int newton(s)

    /* *** NEWTON'S METHOD BACKTRACKING PROCEDURE *** */

    STAGE *s;
    {
    float delta;

      trim=(*s).trim;
      delta = standard - compare;

      if(trim == last_cd) return(FALSE);
      else if (NEWTTOL > fabs(delta) / standard) return(FALSE);
      else if (irepcount++ >= NEWTREP )
       {
       printf(" ** FAILURE TO CONVERGE! ABNORMAL RESULT! **\n\r");
       printf("\n\r");
       return(FALSE);
       }
      else
      {
	if (compare != last_compare)
		slope = (trim - last_cd) / (compare - last_compare);
	else printf("Same Compare \n\r");

	last_cd = trim;
	trim = trim + delta * slope;
	if (trim <= 0) trim = .00015;
	(*s).trim=trim;
	return(TRUE);
      }

 }
/***************************************************************************/
STAGE *newstage()
{
void *malloc();
return ((STAGE *) malloc(sizeof(STAGE)));
}
/***************************************************************************/
void main()
{

STAGES config;

char reply=EOF, old_reply, reply1;
int go_on = TRUE;
int continue1, i;
int kind_switch = 1;
long li;
char *cp;

color_routine();

/* In DOS and WINDOWS this may be only way of getting enough space */
for (i=0; i<MAXSTAGE; i++)
  {
  config.stage[i] = newstage();
  cp = (char *) config.stage[i];
  for(li=0; li < sizeof(STAGE); li++)
    {
    *cp++ = '\0';
    }
    (*config.stage[i]).trim=1.0;
    (*config.stage[i]).signal=FALSE;
  }

MACH1 = (double) get_mach1(temp);
MACH2=MACH1+MACH1;
madjust=0.0;

banner();

#if(TEMPUNITS == 'C')
  temp1 = TEMPC;
#else
  temp1 = TEMPF;
#endif

ifp=dfp=NULL;

while(go_on)
  {
  CLS;
  kind_switch=get_data(config);
  if (kind_switch == 0) /* Simple Predict */
    {
    (void) trajectory(config, 1);
    if(iswitch)fclose(ofp);
    }
  else
    {
    (*config.stage[BACKTRAK_STAGE]).trim = .0001;
    compare = trajectory(config, kind_switch);
    last_cd = (*config.stage[BACKTRAK_STAGE]).trim;
    (*config.stage[BACKTRAK_STAGE]).trim = 2.0;
    irepcount=1;
    continue1 = TRUE;
    while(continue1)
      {
      last_compare = compare;
      compare = trajectory(config, kind_switch);

      continue1 = newton(config.stage[BACKTRAK_STAGE]);
      }
    }
  prtout(kind_switch, config);

  set_next_values(kind_switch, config);

  if(kind_switch != 0)
    {
    if(!(*config.stage[BACKTRAK_STAGE]).profile)
      {
      sprintf((*config.stage[BACKTRAK_STAGE]).dragdollar, "%f", trim);
      (*config.stage[BACKTRAK_STAGE]).cd1=
	  (*config.stage[BACKTRAK_STAGE]).trim;
      (*config.stage[BACKTRAK_STAGE]).trim = 1.0;
      }
    }

  old_reply=reply;
  reply = 'J';
  while ((reply != 'y')&(reply != 'Y')&
	 (reply != 'n')&(reply != 'N'))
    {

    if ((reply == 'o')|(reply == 'O'))
       {
       old_reply='O';
       if((kind_switch != 0) & (!(*config.stage[BACKTRAK_STAGE]).profile))
	 {
	 /* Must fix drag coefficient if just backtracked */
	 for (i = 0;i<DRAGSIZE;i++)
	    (*config.stage[BACKTRAK_STAGE]).cd[i] = trim;
	 machmaker((*config.stage[BACKTRAK_STAGE]).cd, temp, nosesw);

	 cd1=trim;
	 trim=1.0;
	 }
       if((*config.stage[NUMSTAGES-1]).cd[30] < .05)
	   printf("Drag Coefficient Too Small for Optimal Mass Routine\n\r");
       else
       {
       m0 = optmass(config);
       prtout(3, config);
       }
       }
      if(old_reply == 'O')
	{
	reply = get_reply("Another Calculation? (Y/N)");
	if (reply == 'O')reply = EOF;
	}
      else
	reply = get_reply("Another Calculation? (Y/N  or O for Optimal Mass)");
     }
     if (reply == 'N') go_on=FALSE;
     else if ((old_reply == 'O')&&((*config.stage[NUMSTAGES-1]).cd[1] >=.05))
       {
       printf("\n\r");

       reply1=EOF;
       while ((reply1 != 'Y') & (reply1 != 'N'))
	 reply1 = get_reply("Retain Optimal Values For Next Run?");

       if(reply1 == 'Y')
	 {
	 madjust=m0-(*config.stage[0]).m0;
	 set_next_values(0, config);
	 }

       }
  }
return;
}
