/***************************************************************************/ /* <<< PROGRAM DIGITRAK (tm) >>> *** VERTICAL ROCKET TRAJECTORY SIMULATION AND *** *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS *** *** AUTHOR: LARRY CURCIO *** <<< VERSION 3.6>>> LAST UPDATED 11 MARCH 1994 BUG REPORTS, QUESTIONS, SPELLING CORRECTIONS, AND COMMENTS CAN BE SENT TO ME DIRECTLY: FREE DISTRIBUTION AS LONG AS NO PROFIT IS MADE FEE FOR EACH FOR-PROFIT SALE: $25.00 PROGRAM PROVIDED AS-IS. NO WARRANRY IS MADE *** INTERNET: LC2B@DELPHI.COM *** *** SNAILMAIL ADDRESS *** P.O. BOX 81124 PITTSBURGH, PA 15217-81124 *** COPYRIGHT LAWRENCE M. CURCIO *** *** 16 NOVEMBER 1993 *** *** DESCRIPTION *** VERTICAL ROCKET TRAJECTORY SIMULATION PROGRAM WITH FOLLOWING FEATURES: 1) SUPPORTS UP TO THREE STAGES WITH OR WITHOUT STAGING DELAYS 2) SUPPORTS SIMPLE STAGING DELAYS OR AIR START DELAYS, WHERE SEPARATION OCCURS AT LOWER STAGE CUTOFF, BUT UPPER STAGE IGNITION IS DELAYED. 3) USES CLOSED-FORM SOLUTIONS AT EACH DIGITAL INTERVAL 4) DEPRECIATES AIR DENSITY WITH ALTITUDE 5) USES EMPIRICAL DRAG DIVERGENCE MODEL PROPOSED BY BENGEN IN MANDELL, CAPORASO AND BENGEN: _TOPICS IN ADVANCED MODEL ROCKETRY_ (MODEL GOOD FOR QUALITATIVE EFFECTS, AND IMPLICATION OF MACH EFFECTS) 6) SEPARATE TREATMENT OF ROUND AND OPTIMALLY SHARP NOSES IN DRAG DIVERGENCE 7) ADJUSTMENT OF AIR DENSITY AND MACH NUMBER TO TEMPERATURE 8) USES DRAG COEFFICIENTS, SHAPE CONSTANTS, OR DRAG PROFILES 9) USES FINE-INTERVAL OR COARSE-INTERVAL THRUST CURVES. INTERPOLATES COARSE-INTERVAL CURVES. 10) TEMPORAL BACKTRACKING(tm) CAPABILITY: CALCULATION OF TRAJECTORY STATS AND DRAG COEFFICIENT FROM ASCENT TIME 11) OPTICAL BACKTRACKING(tm) CAPABILITY CALCULATION OF TRAJECTORY STATS AND DRAG COEFFICIENT FROM MAXIMUM ALTITUDE 12) OPTIONAL USE OF DRAG PROFILES IN TEMPORAL OR OPTICAL BACKTRACKING(tm) 13) ESTIMATION OF OPTIMAL MASS FOR SINGLE-STAGE ROCKETS OR STAGED ROCKETS WITH NO STAGING DELAYS 14) WRITTEN BY A SUPERSTITIOUS AUTHOR, WHO DOESN'T LIKE TO HAVE 13 LISTED FEATURES. ** VERSION 2.0 ** ADDED RUNGE-KUTTA, AND FIXED SOME TABS BEFORE #if'S TO PLEASE VARIOUS UNIX CC COMPILERS. ** VERSION 2.1 ** ADDED FEHSKENS-MALEWICKI INTEGRALS FOR LAST TIME INTERVAL THIS DISSOCIATED TIMING ACCURACY AND INTERVAL WIDTH, CREATING CONTINUOUS TIME. CONTINUOUS TIME ALLOWED HIGHER TOLERANCES FOR BACKTRACKING, AND ALSO WIDER COAST INTERVALS. RESULT IS A FASTER, MORE ACCURATE BACKTRACKING TOOL. ** VERSION 2.2 ** FIXED INSUBSTANTIAL BUGS IN CONTINUOUS TIME ROUTINE AND IN MASS INITIALIZATION. IMPLEMENTED ENGLISH INPUT OPTION (#define INUNITS). ALLOWED LAUNCHER > SEA LEVEL (#define ALT0). Fixed bug in BACKTRK_STAGE routine that would bomb on a single stage rocket after an upper stage in a multistage rocket was backtracked. Also, header records in the output file are preceeded by '#' to make them gnuplot compatible. ** VERSION 3.0 ** FULL FEHSKENS-MALEWICKI TREATMENT AT EACH INTERVAL. VARIABLE COAST PHASE TIME INTERVALS. ** VERSION 3.1 ** CONTINUOUS TIME EXTENDED TO STAGING DELAYS. THE DEFAULT NUMBER OF INTERVALS IN THE BOOST PHASE REDUCED FROM 500 TO 400. ** VERSION 3.2 ** FIXED BUG IN DRAG PROFILE HANDLING PROCEDURE. IMPLEMENTED COLOR FOR TURBO DOS. ALSO IMPLEMENTED TURBODOS SWITH, WHICH SWITCHES BETWEEN DOS (#define TURBODOS 1) AND UNIX OR GENERIC (#define TURBODOS 0) ** VERSION 3.3 ** OUTPUT FILE FIELDS ARE LABELED. OUTPUT FILE UNITS ARE DETERMINED BY INUNITS OPTION. PROGRAM IS BULLET-PROOFED. FIXED MINOR INCONSISTENCY IN THRUST CURVE INTERPOLATION. ** VERSION 3.4 ** FURTHER BULLETPROOFING IN CD ENTRY ALSO A PROBLEM IN THE HANDLING OF 3-STAGE ROCKETS WAS FIXED ** VERSION 3.5 ** DIGITRAK WOULD GO INTO A TIGHT LOOP IF BACKTRAKING WERE INVOKED ON THE FIRST PASS. THA'S FIXED. ** VERSION 3.6 ** FIXED BUG IN TREATMENT OF DRAG PROFILES AFTER BACKTRACKING WITHOUT PROFILES. CONSPICUOUS FLAW WHEN PROFILES ARE USED. NO CONSEQUENCE WHEN THEY ARE NOT NOTE: BACKTRACKED DRAG COEFFICIENTS ARE *VERY* SENSITIVE TO MOTOR PERFORMANCE VARIATION AND MEASUREMENT ERRORS. DRAG COEFFICIENTS *MUST* THEREFORE BE BACKTRACKED FROM AVERAGE ALTITUDES FROM SEVERAL FLIGHTS. ALTITUDES THEMSELVES MAY BE MEASURED DIRECTLY, OR THROUGH TEMPORAL BACKTRACKING. BACKTRACKED ALTITUDES ARE FAR LESS SENSITIVE TO MOTOR PERFORMANCE VARIATION THAN BACKTRACKED DRAG COEFFICIENTS ARE. MULTI-STAGE VEHICLES SHOULD BE BACKTRACKED BOTTOM-UP, NOT TOP-DOWN. MOST 2-STAGE VEHICLE BACKTRACKING IS EXTREMELY INSENSITIVE TO Cd OF FIRST STAGE. **** COMPILATION NOTES **** SCAN PROGRAM FOR COMPILE OPTION FLAGS (PLURAL). CHANGE #define'S ACCORDINGLY. FIRST SECTION ALSO HAS A (ONE) HACKERS' AREA. THESE #define'S AREN'T NORMALLY CHANGED, BUT HACKERS WILL BE HACKERS. LOOK FOR THE WORD, HACKER. THIS VERSION COMPILES IN DOS TURBO C (1 WARNING), WINDOWS TURBO C++ (SCADS OF WARNINGS FOR NO PROTOTYPES), AND UNIX AT CMU (CLEAN COMPILE WHEN #define CLS IS CHANGED TO printf("\n\n\n"); ****** FOR COMPILERS THAT DON'T SUPPORT #if'S ****** 1) PUT { AT END OF EACH #if LINE 2) REMOVE # FROM EACH #if (CREATING SIMPLE if) 3) CHANGE #else TO: } else { 4) CHANGE #endif TO } 5) FOLLOW INSTRUCTIONS AT AREA MARKED #if EXCEPTIONS 6) COMPILE 7) IGNORE MYRIAD WARNINGS FOR NULL if'S, CODES NO EFFECT, AND UNREACHABLE CODES. (COMPILERS THINK THEY KNOW EVERYTHING.) THERE IS A FUNCTION CALLED color, IN WHICH YOU CAN PUT A SHORT ROUTINE TO SET BACKGROUND COLOR IN YOUR SYSTEM. APOLOGIES FOR THIS SPELLING ARE DUE TO THOSE IN THE UK. APOLGIES ARE DUE TO EVERYONE FOR THE REST OF THE SPELLING. DIGITRAK CAN USE THRUST CURVES FROM SIMLAB (tm). IT CAN ALSO USE THRUST CURVES CONVERTED FROM RASP.ENG USING THE PROGRAM, MARKCURVE.C. TO USE MARKCURVE, FIRST CHANGE ALL TABS IN RASP.ENG TO SPACES. CALL THE NEW FILE ANYTHING YOU LIKE. PUT IT IN A *FRESH DIRECTORY* WITH THE MARKCURVE EXECUTABLE. BOWCOO THRUST CURVE FILES RESULT. EACH BEGINS WITH THE LETTERS, 'MJ'. THE NEXT LETTER DENOTES THE MANUFACTURER, ACCORDING TO RASP TRADITION. THE LAST LETTERS DENOTE THE ENGINE TYPE. THE FILE TYPE OF EACH THRUST CURVE IS .ENG. DIGITRAK APPENDS THIS TO THE THRUST CURVE NAME THAT YOU SPECIFY DURING YOUR INTERACTION WITH IT. THUS - TO SPECIFY AN AEROTECH E15, USE THE CODE: MJAE15 TO SPECIFY AN ESTES E155 (STAND BACK!), USE: MJEE15 MOST SIMLAB CURVES BEGIN WITH THE PREFIX, NEW SIMLAB D1 IS NEWD12 EXCEPTIONS: THE AEROTECH E25 SIMLAB CURVE IS AEROE25. ALSO THERE IS AN FLAT, AVERAGE-THRUST C6 CURVE FOR DIAGNOSTICS, WHICH IS AVEC6. NOTE: AS OF THIS WRITING, THE MJAE25 AND MJAD21 CURVES ARE CATALOG VALUES, AND NOT S&T VALUES. USE THE SIMLAB AEROE25 INSTEAD OF THE MJAE25. *** 16 FEBRUARY 1994 *** */ /***************************************************************************/ #include #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(TURBODOS == 1) #define printf cprintf #include #define CLS clrscr() #else #define CLS printf("\n\n\n\n\n\n") #endif /* LAUNCH PAD ALTITUDE IN ****METERS**** */ #define ALT0 0; /* ALTUNITS should be 'F' for feet and 'M' for meters */ #define ALTUNITS 'F' /* INPUT UNITS ARE 'E' FOR ENGLISH AND 'M' FOR METRIC */ #define INUNITS 'E' /* SET TO CONSTANT_RHO 1 TO TURN OFF DENSITY/ALTITUDE DEPRECIATION */ /* USEFUL FOR DIAGNOSTIC CHECKS USING CLOSED-FORM COAST PHASE FORMULAS */ #define CONSTANT_RHO 0 /* TEMPUNITS should be 'F' for fahrenheit and 'C' for centigrade */ /* SET TEMPUNITS TO 'C' OR 'F' */ #define TEMPUNITS 'F' /* SUBSONIC Coast Phase time interval width */ #define DELTAT0 .1 /* X-SONIC AND HIGHER COAST PHASE INTERVAL WIDTH */ #define DELTAT1 .01 /* CLS is the Screen Clear function. Define as clrscr for TURBO C */ /* define as printf("\n\n\n\n\n") for UNIX */ /* Number of intervals between file record output */ #define REPNUM 10 /* *** END COMPILE OPTIONS *** */ /* *** HACKERS' AREA *** */ /* (DOWN TO *** END HACKERS' AREA ***) */ /* *** *HACKERS* MAY WISH TO CHANGE THESE THINGS *** */ /* *** RECOMMEND THEY BE LEFT ALONE *** */ /* MAXIMUM NUMBER OF DRAG ITERATIONS */ #define MAXIT 10 /* MAXIMUM NUMBER OF STAGES */ #define MAXSTAGE 3 /* Size of Drag Profile */ #define DRAGSIZE 1000 /* INTERPOLATION CUTOFF defines num. intervals needing no interpolation */ /* That is, it's the cutoff between long and short thrust curves */ /* It all arises from different standards in the modroc world */ #define INTERPOLATION_CUTOFF 100 /* TOGGLE controls how long thrust curves are read (f, t or t, f) */ #define TOGGLE 0 /* TOGGLE1 controls how short thrust curves are read (f, t or t, f) */ #define TOGGLE1 1 /* Number of documentation lines in drag profile */ #define NUMDOC 5 /* Number of repetitions in Newton's method backtracking procedure */ #define NEWTREP 75 /* Tolerance in Newton's method backtracking procedure */ #define NEWTTOL .000001 /*Precision check. Tells when to stop drag iterations */ /* Also used as initial value of V to avoid singularity */ #define DELTA .000001 /* String Length */ #define SLEN 255 /* Maximum size of Thrust-Time tables */ #define ARRAYSIZE 1001 /* *** END HACKERS' AREA *** */ /* YOU PROBABLY WON'T WANT TO CHANGE THESE THINGS */ /* Missing or null value */ #define MISSING -9.9999E20 /* INITIAL DEFAULT TEMPERATURE IN F AND C */ #define TEMPF 68.0 #define TEMPC 20.0 #define TRUE 1 #define FALSE 0 #define PI 3.14159265358979 /* GRAVITATIONAL CONSTANT IN m/sec^2 */ #define G 9.805 /* ** AIR DENSITY AT 0 DEGREES C at sea level DO NOT CHANGE */ /* SEE #define ALT0 FOR LAUNCHES ABOVE/BELOW SEA LEVEL */ #define RHO 1.2929 /* METERS TO FEET CONVERSION */ #define MTOF 3.28084 /* OUNCES TO GRAMS CONVERSION */ #define OZTOG 28.35 /* INCHES TO CM CONVERSION */ #define INTOCM 2.54 /* POUNDS TO NEWTONS */ #define LBTONEWT 4.45 /* GUTS OF THE STAGE STRUCTURE */ #define DEFSTRING char dddollar[SLEN],dragdollar[SLEN], dragfile[SLEN]; \ char olddragdollar[SLEN], drag1dollar[SLEN];\ char edollar[SLEN];\ char e1dollar[SLEN], engfile[SLEN];\ double m0, m01;\ float d1, d, cd[DRAGSIZE], trim, comp_time, comp_alt, impulse;\ float mp, tb, cd1, a2, tt1[ARRAYSIZE],ff[ARRAYSIZE];\ int numint, profile, signal;\ double digimp, amax, amin, maxv, a,k,v,vc, hc, h,t,curfew, tcoast;\ double airstart, airstarth, airstartv /* DIGITRAK WAS ORIGINALLY WRITTEN FOR ONE STAGE. IMPRUDENT USE OF GLOBAL VARIABLES MADE CONVERSION TO MULTIPLE STAGES HELLISH. DECIDED TO COPY SCALARS OUT OF STRUCTURES IN SUBROUTINES. DECIDED IT WASN'T SO BAD SINCE MOST FULL THIRD-GENERATION COMPILERS COPY SCALARS INTO SUBROUTINES (AND BACK OUT IF NECESSARY) AND USE POINTERS FOR ARRAYS. PROGRAMS CAN ACTUALLY RUN FASTER THIS WAY - ESPECIALLY WITH VIRTUAL STORAGE, SEGMENTED STORAGE, OR REGISTER MEMORY OPTIMIZATION. */ #define SCALL1(x) x = (*s).x #define SRET1(x) (*s).x = x /* strings take longer, but used little and makes language simpler */ #define SCALL2(x) strcpy(x, (*s).x) #define SRET2(x) strcpy((*s).x,x) #define SCALLA SCALL1(m01);SCALL1(m0);SCALL1(d1);SCALL1(d);SCALL1(trim);\ SCALL1(comp_time);SCALL1(comp_alt);SCALL1(impulse);SCALL1(mp);SCALL1(tb);\ SCALL1(numint);SCALL1(profile); SCALL1(a); SCALL1(k);SCALL1(v);SCALL1(a2);\ SCALL1(h);SCALL1(t); SCALL1(curfew);SCALL1(cd1);SCALL1(digimp);\ SCALL1(amax); SCALL1(amin); SCALL1(maxv);SCALL1(vc);SCALL1(hc); #define SCALLB SCALL2(drag1dollar); SCALL2(edollar);\ SCALL2(e1dollar); SCALL2(engfile), SCALL2(olddragdollar),SCALL2(dragfile) #define SRETA SRET1(m01); SRET1(m0); SRET1(d1); SRET1(d); SRET1(trim);\ SRET1(comp_time); SRET1(comp_alt); SRET1(impulse); SRET1(mp); SRET1(tb);\ SRET1(numint); SRET1(profile), SRET1(curfew);SRET1(digimp);\ SRET1(a); SRET1(k);SRET1(v);SRET1(a2);SRET1(vc);SRET1(hc);\ SRET1(h); SRET1(t); SRET1(cd1); SRET1(amax); SRET1(amin); SRET1(maxv) #define SCALLC SCALL1(a); SCALL1(k);SCALL1(v); SCALL1(h); SCALL1(t);\ SCALL1(amax); SCALL1(amin); SCALL1(maxv);SCALL1(vc);SCALL1(hc);\ Scall1(digimp) #define SRETC SRET1(a); SRET1(k);SRET1(v); SRET1(h); SRET1(t);\ SRET1(amax); SRET1(amin); SRET1(maxv);SRET1(vc);SRET1(hc);\ SRET1(digimp); #define SRETB SRET2(drag1dollar); SRET2(edollar);\ SRET2(e1dollar); SRET2(engfile), SRET2(olddragdollar), SRET2(dragfile) typedef struct { DEFSTRING; } STAGE; typedef struct { STAGE *stage[MAXSTAGE]; }STAGES; float fff[INTERPOLATION_CUTOFF], ttt1[INTERPOLATION_CUTOFF]; float comp_time, comp_alt, temp, temp1; double maxv, mcrit; float fwork, rho, cd1; int kind_switch, iswitch, irepcount, profile, nosesw=1, oldnosesw; int BACKTRAK_STAGE, BS1=1, mach_correct=0; float slope, standard, compare, last_compare, last_cd; long int i, istage, inum, NUMSTAGES=1; int repcnt, liftoff, numint, bigloop; char ofdollar[SLEN], ddollar[SLEN], dragdollar[SLEN]; char dddollar[SLEN], dragfile[SLEN]; char drag1dollar[SLEN], olddragdollar[SLEN], adollar[SLEN], e1dollar[SLEN]; char engfile[SLEN], edollar[SLEN]; float trim = 1.0e0; int go_on, iscaff; double madjust, m0, m01; float d, d1, a1, a2, a3, a4, ax, k1, tb, mp, mb; float impulse; double m, mconst, MACH1, MACH2; double hc, vc, amax, amin; double a,k,z,v,vwork, h, deltat, t, t1, digimp, f, oldf; double mbar, oldm; double halfdrag, oldhalfdrag, oldt1, curfew; char old_kind[SLEN] = "P", kind[SLEN], trimdollar[SLEN]; FILE *ofp=NULL, *ifp, *dfp, *fopen(); int clear_toggle=0;/* Clears Input Screen Every Other Stage */ /***************************************************************************/ void color_routine() /*COMPILE OPTION */ /* *** CAN INSERT ROUTINE HERE TO SET BACKGROUND COLOR *** */ { /* NOTE Insert your own color routine here. */ #if (TURBODOS == 1) textmode(C80); textcolor(WHITE); textbackground(BLUE); #endif return; } /***************************************************************************/ int editnum(instring) /* *** EDIT FOR VALID NUMERIC *** */ char *instring; { int numpt=0, numblank=0, numneg=0, notblank=0; while((*instring !='\0') &&(*instring != '\r')) { if((*instring < '0') || (*instring >'9')) { if(*instring == ' ') { if(notblank) ++numblank; else numblank=1; } else /* (*instring == ' ') */ if(*instring == '-') { ++numneg; if(notblank)return(FALSE); } else if(*instring == '.')++numpt; else return(FALSE); /* non-numeric character found */ notblank=1; } else notblank=1; ++instring; } if(numblank > 2)return(FALSE); if(numpt > 1) return(FALSE); if(numneg > 1) return(FALSE); return(TRUE); } /***************************************************************************/ float get_val(cp) /* *** ACCEPT NUMERIC FROM KEYBOARD, EDIT FOR VALIDITY, AND CONVERT *** */ char *cp; { float res; if(!editnum(cp))return(MISSING); sscanf(cp, "%f", &res); return(res); } /***************************************************************************/ void get_line(upcase,cp) /* *** ACCEPT A LINE FORM THE KEYBOARD *** */ int upcase; char *cp; { int c; for (;;) { c = getchar(); if ((c >= 'a') && (c <= 'z') && (upcase)) c = c + 'A' - 'a'; if (c == EOF || c == '\n') { *cp = 0; return; } *cp++ = c; } } /***************************************************************************/ float prompt(instring, x) /* *** PROMPT USER FOR NUMERIC VALUE AND ACCEPT *** */ /* *** NUMERIC VALUE IS FLOATING POINT, BUT MORE VERSATIE WITH CASTS *** */ char *instring; float x; { int go_on; float y; char response[SLEN]; go_on=TRUE; while (go_on) { y=x; printf("%s (%f) ",instring, x); get_line((int) TRUE, response); if( response[0] != NULL) y = get_val(response); if (y > MISSING) #if(TURBODOS == 0) { go_on=FALSE; } #else { go_on=FALSE; gotoxy(1,wherey()-1); clreol(); printf("%s (%f) \n\r",instring, y); } else { gotoxy(1,wherey()-1); clreol(); } #endif } return(y); } /***************************************************************************/ char get_reply(legend) /* *** ACCEPT SINGLE CHARACTER REPLY AND PUT INTO UPPER CASE *** */ char legend[]; { char reply; printf("%s ", legend); reply = EOF; while((reply == EOF)|(reply == '\n')) reply = getchar(); (void) getchar(); if ((reply >='a')&(reply <= 'z')) reply = reply+'A'-'a'; #if(TURBODOS == 1) gotoxy(1,wherey()-1); printf("%s %c \r", legend, reply); #endif return(reply); } /***************************************************************************/ void banner() /* *** DISPLAY COPYRIGHT BANNER *** */ { char adollar[SLEN]; CLS; printf("\n\n"); printf(" ****************************************************************"); printf("************\n\r"); printf(" * "); printf(" *\n\r"); printf(" * <<< PROGRAM DIGITRAK (tm) >>> "); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * *** VERTICAL ROCKET TRAJECTORY SIMULATION AND ***"); printf(" *\n\r"); printf(" * *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS"); printf(" *** *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * <<< VERSION 3.6 >>>"); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * *** AUTHOR: LARRY CURCIO ***"); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * *** INTERNET: LC2B@DELPHI.COM *** "); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * *** SNAILMAIL ADDRESS *** "); printf(" *\n\r"); printf(" * P.O. BOX 81124 "); printf(" *\n\r"); printf(" * PITTSBURGH, PA 15217-81124 "); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" * *** COPYRIGHT LAWRENCE M. CURCIO ***"); printf(" *\n\r"); printf(" * *** 16 NOVEMBER 1993 *** "); printf(" *\n\r"); printf(" * "); printf(" *\n\r"); printf(" ****************************************************************"); printf("************"); printf("\n\n"); get_line((int) TRUE, adollar); return; } /***************************************************************************/ int interpolate( fin, tin, nin, fout, tout) /* INTERPOLATE COARSE THRUST CURVES */ float fin[], tin[], fout[], tout[]; int nin; { /* Does linear interpolation in an array of points */ /* NUMBER OF LINEARLY INTERPOLATED POINTS */ #define NUMBER 400 float tb, deltat, slope, intercept, t, fwork, deltat1; int j, i, m, nout; tb=tin[nin]; deltat=tb/NUMBER; nout = NUMBER; if (deltat > DELTAT0) { nout = tb/DELTAT0; if (nout >= ARRAYSIZE) nout = ARRAYSIZE-1; deltat = tb/((float) nout); } t=0.0; nout = 0; for(i=1; i<=nin; i++) { fwork=tin[i]-tin[i-1]; m = (int)(fwork/deltat + 0.5); /* see comment before for staement */ if (tin[i] == tb) m++; deltat1= fwork/((float) m); t = tin[i-1]; slope = (fin[i]-fin[i-1])/fwork; intercept=fin[i-1]; /* note: last point in interval not done except in last interval */ for(j =1; j< m; j++) { fout[nout] = slope*deltat1*j + intercept; tout[nout] = t; t +=deltat1; nout++; } } /* nout--;*/ fout[nout]=fin[nin]; tout[nout]=tin[nin]; return(nout); } /****************************************************************************/ int check_for_correct(filename) /* *** CHECK FOR MACH CORRECTION REQUEST *** */ char filename[]; { int cswitch; char *ipoint; /* search for '/C' at end of string, indicating mach correction desired */ ipoint = strrchr(filename, '/'); cswitch = FALSE; if ((ipoint != NULL) && ((strcmp(ipoint, "/C") ==0 ) || (strcmp(ipoint, "/c") == 0))) { /* Found => mark for mach correction and terminate string for file name */ cswitch = TRUE; *ipoint = '\0'; } return(cswitch); } /***************************************************************************/ float get_mach1(temp) /* *** FIND MACH 1 FROM TEMPERATURE *** */ float temp; { #define MACH1AT0 331.45 return(MACH1AT0*sqrt(temp/273.15)); } /***************************************************************************/ void machmaker(cd, temp, nosesw) int nosesw; float cd[], temp; { /*COMPILE OPTION*/ /* This routine accounts for drag divergence in the manner described by Mandell, Caporaso, and Bengen in _TOPICS IN ADVANCED MODEL ROCKETRY_, first edition. The authors use a critical mach number of 0.9. I have converted their model to use a critical mach number of 0.8, but I have given the user a choice of either value. To choose a critical mach number, #define MCRIT as either 8 (for 0.8) or 9 (for 0.9). Since some c compilers running under UNIX cannot understand a #if comparison to a decimal constant, I had to use whole numbers instead of the values themselves. */ #define MCRIT 8 /* DON'T CHANGE THIS */ #define SHARP 1 float MACH1, MACHPT9, MACHPT8, fwork, fwork1, machnum; int IMCRIT, IMACHPT9, IMACHPT8, IMACH1PT05, IMACH1PT2, i, smaller; MACH1 = get_mach1(temp); MACHPT9 = 0.9*MACH1; IMACHPT9 = (int) MACHPT9; MACHPT8 = 0.8*MACH1; IMACHPT8 = (int) MACHPT9; if ((float) IMACHPT9 < MACHPT9) IMACHPT9++; if ((float) IMACHPT8 < MACHPT8) IMACHPT8++; IMACH1PT05 = (int) (MACH1*1.05 + .999999); IMACH1PT2 = (int) (MACH1*1.2 + .999999); #if (MCRIT == 9) mcrit=MACHPT9; IMCRIT=IMACHPT9; if(IMACHPT9 > DRAGSIZE) return; fwork = cd[IMACHPT9-1]; fwork1=(float) IMACHPT9; #else mcrit=MACHPT8; IMCRIT=IMACHPT8; if(IMACHPT8 > DRAGSIZE) return; fwork = cd[IMACHPT8-1]; fwork1=(float) IMACHPT8; #endif if(nosesw == SHARP) { /* Sharp Nosed Vehicles */ if(IMACH1PT05 < DRAGSIZE) smaller=IMACH1PT05; else smaller = DRAGSIZE; for (i=IMCRIT; i Cd VALUE *** */ /* *** NUMBER ENTERED HERE AND NONE IN DIAMETER => SHAPE CONSTANT */ /* *** STRING ENTERED HERE => FILE *** */ /* *** STRING HERE, NO # IN DIAMETER => TAKE DIAMETER FROM PROFILE *** */ /* *** STRING HERE, # IN DIAMETER => TAKE DIAMETER FROM # *** */ int kind_switch, istage; STAGE *s; { int signal; /* If not the backtracked stage, treat as a simple prediction */ if (istage != BACKTRAK_STAGE)kind_switch=0; profile=(*s).profile; signal = (*s).signal; if ((kind_switch == 0) | (profile)) { printf("Drag Coefficient or File (%s) ",(*s).dragdollar); get_line((int) FALSE, (*s).drag1dollar); if((*s).drag1dollar[0] != NULL) strcpy((*s).dragdollar,(*s).drag1dollar); else signal=(*s).signal; if((strcmp((*s).dragdollar, (*s).olddragdollar) != 0)| (nosesw !=oldnosesw)) { if(((*s).dragdollar[0] == '.') | (((*s).dragdollar[0] >= '0' ) & ((*s).dragdollar[0] <= '9'))) { (*s).cd1 = get_val((*s).dragdollar); if ((*s).cd1 == 0.0) (*s).cd1 = DELTA; for (i = 0;i 1) printf("\nFor Stage %i \n\n\r",istage+1); printf( " ** Estimated Optimal Performance at %f Deg %c **\n\n\r", (temp-273.15) TEMPADJ, tempunits); printf("Optimal Mass of Stage = %f Grams; or %f Ounces\n\n\r", 1000.0*((*s).m0+madjust), 1000.0*((*s).m0+madjust)/OZTOG); kind_switch=0; } else { if(NUMSTAGES > 1) { CLS; printf("\nFor Stage %i:\n\n\r",istage+1); } } /* Only one stage can be backtracked at a time */ if(istage != BACKTRAK_STAGE)kind_switch=0; if(((*s).curfew > 0.0)||(istage == NUMSTAGES-1)) printf("Maximum Altitude (m): %f (ft): %f \n\r",(*s).h, (*s).h*MTOF); printf("Cutoff Altitude (m): %f (ft): %f \n\r",(*s).hc, (*s).hc*MTOF); if((*s).airstart != 0.0) printf("Airstart Altitude(m): %f (ft): %f \n\r",(*s).airstarth, (*s).airstarth*MTOF); printf("\n\r"); printf("Maximum Velocity (m/sec): %f (ft/sec): %f \n\r", (*s).maxv, (*s).maxv*MTOF); printf("Average Velocity (m/sec): %f (ft/sec): %f \n\r", (*s).h/(*s).t, (*s).h*MTOF/(*s).t); if((*s).airstart != 0.0) printf("Airstart Velocity (m/sec): %f (ft/sec): %f \n\r",(*s).airstartv, (*s).airstartv*MTOF); if(((*s).curfew > 0.0)&&(istage != NUMSTAGES-1)&&(NUMSTAGES > 1)) printf("Staging Velocity (m/sec): %f (ft/sec): %f \n\r", (*s).v, (*s).v*MTOF); printf("Maximum Acceleration (Gs): %f Minimum: %f \n\r",(*s).amax/G, (*s).amin/G); printf("\n\r"); printf("Maximum Mach Number (dimensionless) (%f) \n\r", (*s).maxv/get_mach1(temp)); printf("\n\r"); printf("Flight Time (sec): %f Coast Time (sec): %f \n\r", (*s).t,(*s).tcoast); if(kind_switch) { if((*s).profile) printf("Trim Factor = %f ",(*s).trim); else printf("Drag Coefficient = %f ",(*s).trim); printf("Sea Level K (kg*m/sec^2) = %f\n",(*s).k*(*s).trim* (*s).cd[(int)((*s).h/(*s).t + 0.5)]); } printf("\n\r"); ax=PI*(*s).d*(*s).d/4.0; if(ax > 0.0) printf("Cross-Sectional Area (cm^2): %f \n\r",ax*10000); printf("Recorded Impulse (n-sec): %f Digital Impulse(n-sec): %f \n\r", (*s).impulse, (*s).digimp); printf("\n\r"); if (istage != NUMSTAGES-1) getchar(); return; } /***************************************************************************/ void prtout(kind_switch, config) /* *** EXECUTIVE ROUTINE TO CALL STAGE OUTPUT ROUTINE, print_stage *** */ int kind_switch; STAGES config; { int istage; for (istage=0; istage 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= (*s).m0) { go_on1=TRUE; printf("Propellant Mass Exceeds Launch Mass \n\r"); printf ("Press Any Key To Continue "); (void)getchar(); #if(TURBODOS == 1) gotoxy(1,wherey()); clreol(); gotoxy(1,wherey()-1); clreol(); gotoxy(1,wherey()-1); clreol(); #else printf("\n\r"); #endif } #if(TURBODOS == 1) gotoxy(1,wherey()-1); clreol(); if(!go_on1)printf("Motor (%s) \n\r",(*s).edollar); #endif } if(istage != NUMSTAGES-1) (*s).curfew = prompt("Staging Delay (sec): ",(*s).curfew); printf("\n\r"); if((iswitch)) { #if(INUNITS == 'M') fprintf(ofp,"# M0= %f Diameter= %f \n",(*s).m0, (*s).d); #else fprintf(ofp,"# M0= %f Diameter= %f \n",1000*(*s).m0/OZTOG, 100*(*s).d/INTOCM); #endif fprintf(ofp,"# Cd= %s Trim= %f Engine= %s \n",(*s).dragdollar, (*s).trim, (*s).edollar); #if(INUNITS == 'M') fprintf(ofp,"# Tb= %f Impulse= %f Mp= %f \n",(*s).tb, (*s).impulse, (*s).mp); #else fprintf(ofp,"# Tb= %f Impulse= %f Mp= %f \n",(*s).tb, (*s).impulse/LBTONEWT, 1000*(*s).mp/OZTOG); #endif fprintf(ofp,"# Staging Delay =%f Airstart Delay= %f\n#\n",(*s).curfew, (*s).airstart); } /* SRETA; SRETB; */ return; } /***************************************************************************/ int get_data(config) /* *** TAKES WHOLE-ROCKET DATA FROM THE KEYBOARD *** */ /* (STAGE-SPECIFIC DATA TAKEN IN get_stage_data *** */ STAGES config; { int kind_switch, go_on, istage; float OLDNUMSTAGES,ANUMSTAGES; CLS; clear_toggle=0; madjust=0.0; OLDNUMSTAGES=NUMSTAGES; go_on = TRUE; while(go_on) { NUMSTAGES=OLDNUMSTAGES; ANUMSTAGES = prompt("Number of Stages ",(float)NUMSTAGES); if((ANUMSTAGES > 0) && (ANUMSTAGES <= MAXSTAGE)) go_on=FALSE; if(fabs(ANUMSTAGES) > MAXSTAGE) NUMSTAGES=MAXSTAGE+1; else NUMSTAGES = (int) ANUMSTAGES; go_on=FALSE; if(ANUMSTAGES != (float)NUMSTAGES) go_on=TRUE; if((NUMSTAGES > MAXSTAGE)||(NUMSTAGES <1)) go_on=TRUE; } if(NUMSTAGES == 1) CLS; go_on = TRUE; while(go_on) { /* TP or OP gets you backtracking with profiles */ printf("'T' = Temporal, 'O' = Optical, 'P' = Predict (%s) ", old_kind); get_line((int) TRUE, kind); if(kind[0] != NULL) strcpy(old_kind, kind); go_on = FALSE; if(strcmp(old_kind, "P") == 0) { kind_switch = 0; } else if(strcmp(old_kind, "T") == 0) { kind_switch = 1; profile = FALSE; } else if(strcmp(old_kind, "O") == 0) { kind_switch = 2; profile = FALSE; } else if(strcmp(old_kind, "TP") == 0) { kind_switch = 1; profile = TRUE; } else if(strcmp(old_kind, "OP") == 0) { kind_switch = 2; profile = TRUE; } else go_on = TRUE; #if(TURBODOS == 1) CLS; if(!go_on) printf("'T' = Temporal, 'O' = Optical, 'P' = Predict (%s)\r\n", old_kind); #endif } iswitch = 0; if(kind_switch == 0) { printf("Output File ( = None) "); get_line((int) FALSE, ofdollar); if(ofdollar[0] != NULL) { iswitch = 1; ofp = fopen(ofdollar, "w"); } #if(TURBODOS == 1) gotoxy(1,wherey()-1); clreol(); printf("Output File (%s)\n\r", ofdollar); #endif } else if(NUMSTAGES > 1) { BACKTRAK_STAGE=9; while((BACKTRAK_STAGE < 0 )|| (BACKTRAK_STAGE > NUMSTAGES-1)) { BS1= (int) prompt("Stage to backtrak", (float) BS1); BACKTRAK_STAGE=BS1-1; (*config.stage[BACKTRAK_STAGE]).profile=profile; } } else { BACKTRAK_STAGE=0; (*config.stage[BACKTRAK_STAGE]).profile=profile; printf("\n\r"); } #if TEMPUNITS == 'C' temp1=prompt("Temperature in Deg C ", temp1); temp = temp1 + 273.15; #else temp1=prompt("Temperature in Deg F ", temp1); temp = (temp1 - 32.0)*5.0/9.0 + 273.15; #endif rho = RHO * 273.15/temp; /* Nose type is used only in supersonic calculations */ oldnosesw=nosesw; nosesw=5; while((nosesw !=1)&(nosesw !=2)) { nosesw=oldnosesw; nosesw= (int)prompt("Nose Type 1=Sharp 2=Round ", (float)nosesw); } if (kind_switch == 1) { (*config.stage[NUMSTAGES-1]).comp_time= prompt("Ascent Time ",(*config.stage[NUMSTAGES-1]).comp_time); standard = (*config.stage[NUMSTAGES-1]).comp_time; } else if(kind_switch == 2) { (*config.stage[NUMSTAGES-1]).comp_alt= #if(ALTUNITS == 'M') prompt("Altitude in Meters ",(*config.stage[NUMSTAGES-1]).comp_alt); standard = (*config.stage[NUMSTAGES-1]).comp_alt; #else prompt("Altitude in Feet ",(*config.stage[NUMSTAGES-1]).comp_alt); standard = (*config.stage[NUMSTAGES-1]).comp_alt / MTOF; #endif } if(iswitch) { fprintf(ofp,"# Number of Stages = %i\n", (int) NUMSTAGES); #if(INUNITS == 'M') fprintf(ofp,"# Units = Metric \n"); fprintf(ofp,"# Temperature (Deg C) = %f ",temp-273.15); #else fprintf(ofp, "# Units = English \n"); fprintf(ofp,"# Temperature (Deg F) = %f ",(temp-273.15)*9/5 + 32); #endif if(nosesw == 1)fprintf(ofp,"Nose = Sharp\n#\n"); else fprintf(ofp,"Nose = Round\n#\n"); } for (istage=0; istage< NUMSTAGES; istage++) get_stage_data(istage, kind_switch, config.stage[istage]); #if (TURBODOS == 1) printf("Press Any Key To Continue"); get_line(FALSE,adollar); gotoxy(1,wherey()-1); textcolor(RED+BLINK); clreol(); printf("Calculating"); textcolor(WHITE); #endif return(kind_switch); } /***************************************************************************/ void init(config) /* *** INITIALIZES A WHOLE FLIGHT *** */ STAGES config; { v = DELTA; t=h=0.0; repcnt=0; m=(*config.stage[0]).m0 + madjust; oldhalfdrag=maxv=0.0; liftoff = FALSE; f=0; oldf=0.0; oldt1=0.0; } /***************************************************************************/ /* ** FM COAST TIME ** */ /* ** FOR LAST COAST INTERVAL TO MAKE BACKTRACKED TIME CONTINUOUS ** */ /* ** WITHOUT THIS, INTERVALS MUST BE VERY SMALL ** */ double FM_coast_time(vb, z) double vb, z; /* z is K/mb */ { return (sqrt(1.0 / (G * z)) * atan(vb * sqrt(z / (G)))); } /***************************************************************************/ /* ** FM COAST ALTITUDE ** */ /* ** FOR LAST COAST INTERVAL TO MAKE BACKTRACKED ALTITUDE PRECISE ** */ double FM_coast_alt(vb, z) double vb, z; /* z is K/mb */ { return(log(z * vb * vb / G + 1) / (z + z)); } /***************************************************************************/ void boost(s) /* *** BOOST PHASE ROUTINE *** */ STAGE *s; { register double work; double fwork, work2, netf, newv; double sqrtg, sqrtz, sqrtzovg; double H, m0hold, Gstar, zstar; float *CD, *FF, *TT; SCALLA; /* ** copy here for faster addressing in loops ** */ CD = (*s).cd; FF=(*s).ff; TT=(*s).tt1; maxv=0; /* For Optimal Mass, adjust m0 */ m0hold=m0; m0 += madjust; oldm=m0; oldt1=0.0; if(d == 0.0) { d = a2/100.0; d1=a2; #if(INUNITS == 'E') d1 = d1/INTOCM; #endif } ax=PI * d * d /4.0; if( ax != 0.0) k = trim * ax * rho / 2.0; else k = trim / 10000.0; digimp=0.0; mconst=mp / impulse; if (iswitch) { ++repcnt; if (repcnt >= REPNUM) { #if(INUNITS == 'M') fprintf(ofp,"%f %f 0.0 0.0 %f \n",h, v, t); #else fprintf(ofp,"%f %f 0.0 0.0 %f \n",h*MTOF, v*MTOF, t); #endif repcnt=0; } }/* if iswitch */ for(inum=0; inum <=numint; inum++) { oldf=f; f = FF[inum]; t1=TT[inum]; fwork=0.5*(f+oldf); deltat = t1-oldt1; oldt1=t1; t += deltat; digimp = digimp + fwork * deltat; m = m0 - digimp * mconst; mbar = (m+oldm)/2.0; oldm=m; if (fwork > mbar * G) {liftoff = TRUE;} if (liftoff) { #if (CONSTANT_RHO) z = k*CD[(int)(v+0.5)]; #else H=h+ALT0; z=k*CD[(int)(v+0.5)]*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14)))); #endif netf=fwork-mbar*G; if(netf > 0.0) { /* STRAIGHT EXTENDED FM BOOST-PHASE INTEGRALS */ work=cosh(deltat*sqrt(z*netf)/mbar)+ v*sqrt(z/netf)*sinh(deltat*sqrt(z*netf)/mbar); h+=mbar*log(work)/z; work=tanh(deltat*sqrt(z*netf)/mbar); work2=v*sqrt(z/netf); v=sqrt(netf/z)*(work+work2)/(1+work*work2); } else { /* COAST PHASE INTEGRALS WITH ADJUSTED G FOR netf < 0.0 */ Gstar = -netf/(mbar); sqrtg = sqrt(Gstar); zstar=z/mbar; sqrtz=sqrt(zstar); sqrtzovg=sqrtz/sqrtg; /* set work to tan(arg1) */ /* where arg1=atan(v*sqrtzovg) - deltat*sqrtz*sqrtg; */ /* using identity for tan(alpha-beta) */ work=tan(sqrtz*sqrtg*deltat); work=(v*sqrtzovg-work)/(1.0+v*sqrtzovg*work); newv=work/sqrtzovg; /* work is now the cos^2 arg1 */ /* cos(theta)^2 = 1/[1 +tan(theta)^2] */ work=1.0/(work*work+1.0); h += 0.5*log((v*v*zstar/Gstar+1.0)*work)/zstar; work=0.5*(v+newv); v=newv; } a = (netf-z*work*work)/mbar; if (v > maxv) maxv = v; if (a > amax) amax = a; if (a < amin) amin = a; hc=h; vc = v; if (iswitch) { ++repcnt; if (repcnt >= REPNUM) { #if(INUNITS == 'M') fprintf(ofp,"%f %f %f %f %f \n",h, v, a, trim*CD[(int)(v+0.5)], t); #else fprintf(ofp,"%f %f %f %f %f \n",h*MTOF, v*MTOF, a*MTOF, trim*CD[(int)(v+0.5)], t); #endif repcnt=0; } }/* if iswitch */ }/* if (liftoff);*/ } if (iswitch) { if (repcnt != 0) { #if(INUNITS == 'M') fprintf(ofp,"%f %f %f %f %f \n",h, v, a, trim*CD[(int)(v+0.5)], t); #else fprintf(ofp,"%f %f %f %f %f \n",h*MTOF, v*MTOF, a*MTOF, trim*CD[(int)(v+0.5)], t); #endif repcnt=0; } }/* if iswitch */ /* unadjust m0 */ m0=m0hold; SRETA; return; } /***************************************************************************/ double coast(kind_switch, s, delay, curfew_switch) /* *** COAST PHASE ROUTINE *** */ /* *** CURFEW IS THE STAGING DELAY TIME, IF ANY *** */ int kind_switch, curfew_switch; double delay; STAGE *s; { register double work; double H, newv, sqrtg, sqrtz, sqrtzovg, curfew1; float m0hold, *CD; int go_on, repnum,REPNUM1; SCALLA; /* ** copy here for faster addressing in loops ** */ CD=(*s).cd; m0hold=m0; m0 +=madjust; /* DECREASE REPNUM TO REPNUM1 FOR WIDER INTERVAL */ /* DELTAT0 IS WIDE INTERVAL FOR LOW ACCELERATION COAST */ /* DELTAT1 IS NARROW INTERVAL FOR HIGH ACCELERATION COAST */ REPNUM1= (int) 0.999999 + REPNUM*(DELTAT1/DELTAT0); if(curfew_switch)curfew1=99999e9; else curfew1=delay; if(d == 0.0) { d = a2/100.0; d1=a2; #if(INUNITS == 'E') d1 = d1/INTOCM; #endif } ax=PI * d * d /4.0; if( ax != 0.0) k = trim * ax * rho / 2.0; else k = trim / 10000.0; curfew1=curfew1+t; mb=m0-mp; k=k/mb; repcnt=0; sqrtg=sqrt(G); go_on = TRUE; while(go_on) { if((v < mcrit)||(v > MACH2)) { deltat=DELTAT0; repnum=REPNUM1; } else { deltat=DELTAT1; /* SMALLER INTERVAL FOR DIVERGENCE */ repnum=REPNUM; } if (t+deltat > curfew1) { deltat = curfew1-t; go_on = FALSE; } #if (!CONSTANT_RHO) H=h+ALT0; z=CD[(int)(v+0.5)]*k*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14)))); #else z=CD[(int)(v+0.5)]*k; #endif sqrtz=sqrt(z); sqrtzovg=sqrtz/sqrtg; /* set work to tan(arg1) */ /* where arg1=atan(v*sqrtzovg) - deltat*sqrtz*sqrtg; */ /* using identity for tan(alpha-beta) */ work=tan(sqrtz*sqrtg*deltat); work=(v*sqrtzovg-work)/(1.0+v*sqrtzovg*work); if(work < 0.0) { go_on = FALSE; if(v != 0.0) { h += FM_coast_alt(v, z); t += FM_coast_time(v, z); v=0.0; } }/* if(arg1 < 0.0) */ else { /* newv = tan(arg1)*sqrt(g/z) */ newv= work/sqrtzovg; /* work is now the cos^2 arg1 */ /* cos(theta)^2 = 1/[1 +tan(theta)^2] */ work=1.0/(work*work+1.0); h += 0.5*log((v*v*z/G+1.0)*work)/z; work=0.5*(v+newv); v=newv; a= -z*work*work - G; if (a < amin) amin = a; t += deltat; if(iswitch) { ++repcnt; if (repcnt >= repnum) { #if(INUNITS == 'M') fprintf(ofp,"%f %f %f %f %f \n",h, v, a, trim*CD[(int)(v+0.5)], t); #else fprintf(ofp,"%f %f %f %f %f \n",h*MTOF, v*MTOF, a*MTOF, trim*CD[(int)(v+0.5)], t); #endif repcnt=0; } }/* if iswitch */ } /* else (v still positive) */ } /* else (no boost) */ /* FORCE LAST RECORD OUTPUT FOR WIDE INTERVALS */ if(iswitch) { if (repcnt != 0) /* OMIT REDUNDANT LAST RECORDS */ { #if(INUNITS == 'M') fprintf(ofp,"%f %f %f %f %f \n",h, v, a, trim*CD[(int)(v+0.5)], t); #else fprintf(ofp,"%f %f %f %f %f \n",h*MTOF, v*MTOF, a*MTOF, trim*CD[(int)(v+0.5)], t); #endif repcnt=0; } }/* if iswitch */ (*s).tcoast=t-(*s).t; m0=m0hold; SRETA; if (kind_switch == 1) return (t); else return(h); } /***************************************************************************/ float trajectory(config, kind_switch) /* KNITS BOOST AND COAST PHASES TOGETHER TO MAKE A WHOLE FLIGHT */ /* EFFECTS SIMPLE STAGING DELAYS AND AIRSTART DELAYS */ int kind_switch; STAGES config; { int istage; init(config); if(iswitch)fprintf(ofp,"# Alt Vel Acc Cd T\n"); for(istage=0; istage0) boost(config.stage[istage]); else (*config.stage[istage]).amax=0.0; if ((*config.stage[istage]).curfew != 0.0) { (*config.stage[istage]).v=v; (*config.stage[istage]).h=h; (*config.stage[istage]).t=t; (void) coast(kind_switch, config.stage[istage], (*config.stage[istage]).curfew, 0); } (*config.stage[istage]).v=v; (*config.stage[istage]).h=h; (*config.stage[istage]).t=t; } return(coast(kind_switch, config.stage[NUMSTAGES-1], 99999e9, 1)); } /***************************************************************************/ float optmass(config) STAGES config; { double mdd, muu, midm; double delta, pnt, pntlo, pnthi, maxm; int iopt, isw, i; float totimpulse; /* ** OPTIMAL MASS ROUTINE ** */ /* ** Mdd IS M CORRESPONDING TO DOWNWARD SLOPE ** */ /* ** Muu IS M CORRESPONDING TO UPWARD SLOPE ** */ #define DELTAMAS 0.0001 #define MAXITMAS 100 maxm=(*config.stage[0]).m0; totimpulse=0.0; for(i=0; i .0001) mdd = totimpulse / (2.0 * G * tb); else mdd=totimpulse / (2.0 * G * (*config.stage[NUMSTAGES-2]).tb); muu=(*config.stage[0]).m0 - (*config.stage[NUMSTAGES-1]).m0 +(*config.stage[NUMSTAGES-1]).mp; /* printf("muu= %f \n\r",muu*1000.0); */ delta = .0001; iopt = 1; isw = 1; printf("\n\r"); #if(TURBODOS == 1) clreol(); #endif while (isw) { midm = (mdd + muu) / 2.0; printf("trying m= %f Grams; or %f Ounces\n\r", midm*1000.0, midm*1000.0/OZTOG); madjust = midm - maxm - delta; pntlo=trajectory(config, 2); madjust = midm - maxm + delta; pnthi=trajectory(config, 2); madjust = midm - maxm; pnt=trajectory(config, 2); pntlo=pnt-pntlo; pnthi=pnthi-pnt; if ( (pntlo < 0.0) & (pnthi < 0.0) ) mdd = midm; else if ( (pntlo > 0.0) & (pnthi > 0.0) ) muu = midm; else { isw = 0; for(i=0; i < NUMSTAGES; i++) { if(((*config.stage[i]).airstartv < 0.0)|| (((*config.stage[i]).v < 0.0)&& (i < NUMSTAGES -1))) { isw=1; mdd=midm; } } } iopt = iopt + 1; if (iopt > MAXITMAS) isw = 0; if(fabs(mdd - muu)/midm < DELTAMAS) isw = 0; } if (iopt > MAXITMAS) printf("*** NO CONVERGENCE *** \n\r"); /* madjust=0.0; */ return(midm); } /***************************************************************************/ int newton(s) /* *** NEWTON'S METHOD BACKTRACKING PROCEDURE *** */ STAGE *s; { float delta; trim=(*s).trim; delta = standard - compare; if(trim == last_cd) return(FALSE); else if (NEWTTOL > fabs(delta) / standard) return(FALSE); else if (irepcount++ >= NEWTREP ) { printf(" ** FAILURE TO CONVERGE! ABNORMAL RESULT! **\n\r"); printf("\n\r"); return(FALSE); } else { if (compare != last_compare) slope = (trim - last_cd) / (compare - last_compare); else printf("Same Compare \n\r"); last_cd = trim; trim = trim + delta * slope; if (trim <= 0) trim = .00015; (*s).trim=trim; return(TRUE); } } /***************************************************************************/ STAGE *newstage() { void *malloc(); return ((STAGE *) malloc(sizeof(STAGE))); } /***************************************************************************/ void main() { STAGES config; char reply=EOF, old_reply, reply1; int go_on = TRUE; int continue1, i; int kind_switch = 1; long li; char *cp; color_routine(); /* In DOS and WINDOWS this may be only way of getting enough space */ for (i=0; i=.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; }