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


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

		     *** AUTHOR: LARRY CURCIO ***

			  <<< VERSION 4.1 >>>
		     LAST UPDATED 19 JUNE 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 WARRANTY 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

			 ** VERSION 4.1 **
	       SCREEN CLEARING FOR DISPLAYS WITH ANSI CURSOR CONTROL
	       OPTIMAL MASS ROUTINE TUNED AND IMPROVED
	       FORM FACTOR CALCULATION CORRECTED
               DRAG PROFILE READIN ROUTINE PERFORMS CORRECTLY ON ALL MACHINES


			 ** VERSION 4.1 **
	       Minor bug in the reading of thrust curves for very large
	       motors is fixed.

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 ADHERES TO ANSI C STANDARD. 
      ON UNIX, FIND ANSI C SWITCH FOR THE CC COMPILER

      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 COMPILING ON ANSI WITH CURSOR CONTROL, SET ANSI TO 1 - OTHERWISE 0 */
/* EXAMPLES OF ANSI CURSOR CONTROL MASCHINES ARE VTXXX WHERE XXX>=100, AND */
/* IBM PC WITH ANSI.SYS INSTALLED AS A DRIVER. (THOUH IBM OPTION BETTER) */


#define ANSI 1

#if(TURBODOS == 1)
#define printf  cprintf
#include <conio.h>
#define CLS clrscr()

#elif ANSI == 1
#define CLS vmscls()
#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 time 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;
int clear_toggle=0;/* Clears Input Screen Every Other Stage */

/***************************************************************************/
#if (ANSI == 1)
void vmscls()
{

#define ESC 27
char try[13]= {ESC,'[','2','H',ESC,'[','2','J',ESC,'[','2','H','\0'};

printf("%s", try);
return;
}
#endif
/***************************************************************************/
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(char * instring)
/* *** EDIT FOR VALID NUMERIC *** */
{
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(char * cp)

/* *** ACCEPT NUMERIC FROM KEYBOARD, EDIT FOR VALIDITY, AND CONVERT *** */
{
  float res;

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

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

/* *** ACCEPT A LINE FORM THE KEYBOARD *** */
{
  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(char * instring, float x)
/* *** PROMPT USER FOR NUMERIC VALUE AND ACCEPT *** */
/* *** NUMERIC VALUE IS FLOATING POINT, BUT MORE VERSATIE WITH CASTS *** */
{
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] != '\0') 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;
char linein[10];
printf("%s ", legend);

do{
  get_line(1, linein);
  reply = linein[0];
  }while (strlen(linein) != 1);

#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 4.1 >>>");
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( float fin[], float tin[],
		 int nin, float fout[], float tout[])

/* INTERPOLATE COARSE THRUST CURVES */

{
/* 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) + 1;
  /* see comment before for statement */
  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(char filename[])

/* *** CHECK FOR MACH CORRECTION REQUEST *** */
{
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(float temp)
/* *** FIND MACH 1 FROM TEMPERATURE *** */
{
#define MACH1AT0 331.45
return(MACH1AT0*sqrt(temp/273.15));
}
/***************************************************************************/
void machmaker(float cd[], float temp, int nosesw)
{
/*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(int kind_switch, STAGES config)

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

/* 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(int kind_switch, int istage, STAGE *s)

/* *** PROMPTS USED FOR DRAG DATA AND PARSES THE REPLY *** */
/* *** DRAG MAY BE A Cd VALUE, A FORM FACTOR, 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 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))
  {
#if(INUNITS == 'M')
  printf("Drag Coefficient, K (kg*m/sec^2), or File (%s) ",s->dragdollar);
#else
  printf("Drag Coefficient, K (slug*ft/sec^2), or File (%s) ",s->dragdollar);
#endif
  get_line((int) FALSE, s->drag1dollar);
  if(s->drag1dollar[0] != '\0') 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 */

	for (i = 0;
	  ((fscanf(dfp, "%f", &s->cd[i]) != EOF) & (i < DRAGSIZE));
	  i++);
	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] = '\0';

#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(int kind_switch, int istage, STAGE * s)

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

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

float vav;

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

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

ax=PI*s->d*s->d/4.0;

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", vav = 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(ax != 0.0)
    {
    if(s->profile) printf("Trim Factor = %f      ",s->trim);
    else printf("Drag Coefficient = %f      ",s->trim);
			       
#if(INUNITS == 'M')
	  printf("Sea Level K (kg*m/sec^2)*10^4 = %f\n",s->trim*
		  ax*rho*5000.0*s->cd[(int)(vav+0.5)]);
#else
	  printf("Sea Level K (slug*ft/sec^2)*10^4 = %f\n",s->trim*
		  ax*rho*5000.0*s->cd[(int)(vav+0.5)]*2.2*MTOF/32.2);


#endif

    }
  else
    {
    if(s->profile)printf("Trim Factor = %f  ",s->trim);   
      
#if(INUNITS == 'M')
	printf
	  ("Sea Level K (kg*m/sec^2) = %f \n",s->trim*s->cd[(int)(vav+0.5)]);
#else
	printf
	  ("Sea Level K (slug*ft/sec^2) = %f \n",s->trim*s->cd[(int)(vav+0.5)]
	   /*  *2.2*MTOF/32.2 */);
#endif
      }

  }

printf("\n\r");

if(ax > 0.0)
  {
#if(INUNITS == 'M')
  printf("X-Sectional Area (cm^2): %f  ",ax*10000);
  if (!kind_switch)printf("Sea Level K (kg*m/sec^2)*10^4 = %f"
     ,s->trim*s->cd[(int)(vav+0.5)]*ax*rho*5000.0);
#else
  printf("X-Sectional Area (inches^2): %f  ",ax*10000/(INTOCM*INTOCM));
  if (!kind_switch)printf("Sea Level K (slug*ft/sec^2)*10^4 = %f"
     ,s->trim*s->cd[(int)(vav+0.5)]*ax*rho*5000.0*2.2*MTOF/32.2);

#endif
  
printf("\n\r");
  }
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(int kind_switch, STAGES config)

/* *** EXECUTIVE ROUTINE TO CALL STAGE OUTPUT ROUTINE, print_stage *** */
{
int istage;
for (istage=0; istage<NUMSTAGES; istage++)
  {
  print_stage(kind_switch, istage, config.stage[istage]);
  }
return;
}
/***************************************************************************/
void get_stage_data(int istage, int kind_switch, STAGE *s)

/* *** STAGE-SPECIFIC KEYBOARD INPUT ROUTINE *** */
{
int ii, go_on1;

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] != '\0') | (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);
  }

return;
}
/***************************************************************************/
int get_data(STAGES config)

/* *** TAKES WHOLE-ROCKET DATA FROM THE KEYBOARD *** */
/* (STAGE-SPECIFIC DATA TAKEN IN get_stage_data *** */
{
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] != '\0') 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] != '\0')
    {
     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(STAGES config)

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

{
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(double vb, double 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(double vb, double z)
/* z is K/mb */
{
  return(log(z * vb * vb / G + 1) / (z + z));
}
/***************************************************************************/
void boost(STAGE * s)

/* *** BOOST PHASE ROUTINE *** */
{
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;
#if(INUNITS == 'M')
else k = trim / 10000.0;
#else
else k = trim * 32.2 / (10000.0 * 2.2 * MTOF);
#endif
		 
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(int kind_switch, STAGE * s, double delay, int curfew_switch)

/* *** COAST PHASE ROUTINE *** */
/* *** CURFEW IS THE STAGING DELAY TIME, IF ANY *** */
{
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;
#if(INUNITS == 'M')
else k = trim / 10000.0;
#else
else k = trim * 32.2 / (10000.0*2.2*MTOF);
#endif


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(STAGES config, int kind_switch)

/* KNITS BOOST AND COAST PHASES TOGETHER TO MAKE A WHOLE FLIGHT */
/* EFFECTS SIMPLE STAGING DELAYS AND AIRSTART DELAYS */
{
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(STAGES config)
{
double mdd, muu, midm;
double delta, pnt, pntlo, pnthi, maxm;
int iopt, isw, isw1, i;
float totimpulse;
char reply;
 /* ** OPTIMAL MASS ROUTINE ** */
 /* ** Mdd IS M CORRESPONDING TO DOWNWARD SLOPE ** */
 /* ** Muu IS M CORRESPONDING TO UPWARD SLOPE ** */

#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;

 delta = .001;
 iopt = 1;
 isw = isw1 = 1;
 printf("\n\r");
#if(TURBODOS == 1)
clreol();
#endif

 while (isw1)
 {
 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;
 } /* while isw */

  isw1 = 0;


 } /* while isw1 */

if (iopt > MAXITMAS) printf("*** NO CONVERGENCE *** \n\r");
else if((pntlo < 0.0) || (pnthi > 0.0) )
	  {
	  /* Note: Approximate algorithm wih fairly large delta */
	  /* One value exceeds midm, and every little bit helps ! */
	  pntlo = -pntlo;
	  if(pntlo > pnthi) midm -= delta;
	  else midm += delta;
          madjust = midm - maxm; pnt=trajectory(config, 2);
	  }
/* madjust=0.0;  */
return(midm);
}
/***************************************************************************/
    int newton(STAGE * s)

    /* *** NEWTON'S METHOD BACKTRACKING PROCEDURE *** */
    {
    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)));
}
/***************************************************************************/
main()
{

STAGES config;

char reply=EOF, old_reply, reply1;
int go_on = TRUE;
int continue1, i, j;
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;
}

