/***************************************************************************/ /* <<< PROGRAM DIGITRAK (tm) >>> *** VERTICAL ROCKET TRAJECTORY SIMULATION AND *** *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS *** *** AUTHOR: LARRY CURCIO *** <<< VERSION 4.2 >>> 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.0 ** 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. ** VERSION 4.2 ** Bug in Acceleration as printed in output file 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 #include #include #include /*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 0 #if(TURBODOS == 1) #define printf cprintf #include #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 1001 /* 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.2 >>>"); 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; icomp_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;icd[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; istageprofile; /* 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;icd[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 ( = 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 = (f - z*v*v)/m - G; 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*v*v - 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; istagev=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; iimpulse; } 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; itrim=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;icd[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; }