/***************************************************************************/ /* <<< PROGRAM DIGITRAK (tm) >>> *** VERTICAL ROCKET TRAJECTORY SIMULATION AND *** *** DIGITAL BACKTRACKING SYSTEM FOR MODEL AND HPR ROCKETS *** *** AUTHOR: LARRY CURCIO *** <<< VERSION 2.8>>> LAST UPDATED 8 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 INSTANCE OF 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 RUNGE-KUTTA DRAG CORRECTION AT EACH 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 2.3 ** CONTINUOUS TIME IMPLEMENTED IN STAGING DELAYS. 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. ** VERSION 2.4 ** 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 2.5 ** OUTPUT FILE FIELDS ARE LABELED. ALSO, OUTPUT FILE UNITS CAN BE SWITCHED BY INUNITS OPTION. PROGRAM IS BULLET-PROOFED. FIXED MINOR INCONSISTENCY IN THRUST CURVE INTERPOLATION. ** VERSION 2.6 ** FURTHER BULLETPROOFING IN CD ENTRY ALSO A PROBLEM IN THE HANDLING OF 3-STAGE ROCKETS WAS FIXED ** VERSION 2.7 ** DIGITRAK WOULD GO INTO A TIGHT LOOP IF BACKTRAKING WERE INVOKED ON THE FIRST PASS. THA'S FIXED. ** VERSION 2.8 ** FIXED BUG IN TREATMENT OF DRAG PROFILES AFTER BACKTRACKING WITHOUT PROFILES. CONSPICUOUS FLAW WHEN PROFILES ARE USED. NO CONSEQUENCE WHEN THEY ARE NOT **** 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\r"); ****** 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 *** */ #define TURBODOS 0 #if(TURBODOS == 1) #include #define printf cprintf #define CLS clrscr() #else #define CLS printf("\n\n\n\n\n\n") #endif /* LAUNCH PAD ALTITUDE IN ****METERS**** */ #define ALT0 0 /* ALTUNITS (Altitude Units) should be 'F' for feet and 'M' for meters */ #define ALTUNITS 'F' /* INUNITS (Input Units) should be 'M' for metric or 'E' for English */ /* English units are inches and ounces. Metri units are cm and grams */ #define INUNITS 'E' /* TEMPUNITS should be 'F' for fahrenheit and 'C' for centigrade */ /* 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 /* SET TEMPUNITS TO 'C' OR 'F' */ #define TEMPUNITS 'F' /* Coast Phase time interval width */ #define DELTAT0 .01 /* 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 /* Number of intervals between file record output */ #define REPNUM 10 /* 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 STRUCTUER */ #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; 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; double hc, vc, amax, amin; double a,k,z,v,vwork,lastv, h, deltat, t, t1, digimp, f, oldf; double mbar, newm; 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 routine here if you like */ #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 2.8 >>>"); 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 500 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) IMCRIT=IMACHPT9; if(IMACHPT9 > DRAGSIZE) return; fwork = cd[IMACHPT9-1]; fwork1=(float) IMACHPT9; #else 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\r",(*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"); } } 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; printf("\n\r"); (*config.stage[BACKTRAK_STAGE]).profile=profile; } #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]); return(kind_switch); } /***************************************************************************/ void init() /* *** INITIALIZES A WHOLE FLIGHT *** */ { v = DELTA; lastv=t=h=0.0; repcnt=0; 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; { double fwork, A1, A2, A3, A4, V1, deltahalf; double m0hold, H; 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; m = 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; newm = m0 - digimp * mconst; mbar = (m+newm)/2.0; if (fwork > mbar * G) {liftoff = TRUE;} if (liftoff) { lastv = v; #if (CONSTANT_RHO) z = k; #else H=h+ALT0; z=k*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14)))); #endif /* *** Runge Kutta Method *** */ deltahalf=deltat*0.5; A1=(oldf - z*CD[(int)(v+0.5)]*v*v)/m -G; V1=v + deltahalf*A1; A2=(fwork - z*CD[(int)(V1+0.5)]*V1*V1)/mbar -G; V1=v+deltahalf*A2; A3=(fwork - z*CD[(int)(V1+0.5)]*V1*V1)/mbar -G; V1=v+deltat*A3; A4=(f - z*CD[(int)(V1+0.5)]*V1*V1)/newm -G; a=(A1 + 2.0*(A2 + A3) + A4)/6.0; v += a*deltat; m=newm; if (v > maxv) maxv = v; if (a > amax) amax = a; if (a < amin) amin = a; /* if(lastv > v) {v = lastv;} */ h = h + (v + lastv) * 0.5 * deltat; 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; { double deltahalf, deltav, V1, A1, A2, A3, A4, H, curfew1; float m0hold, *CD; int go_on; SCALLA; /* ** copy here for faster addressing in loops ** */ CD=(*s).cd; m0hold=m0; m0 +=madjust; deltat=DELTAT0; if(curfew_switch)curfew1=99999e9; else curfew1=delay; deltahalf = DELTAT0/2.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; curfew1=curfew1+t; mb=m0-mp; k=k/mb; repcnt=0; go_on = TRUE; while(go_on) { lastv = v; if (t+deltat > curfew1) { deltat = curfew1-t; go_on = FALSE; deltahalf=deltat/2.0; } #if (!CONSTANT_RHO) H=h+ALT0; z=k*exp(H * (-.00008187+H*(-3.085e-9 + H*(3.975e-14)))); #endif /* *** Runge Kutta Method *** */ A1= z*CD[(int)(v+0.5)]*v*v + G; V1=v - deltahalf*A1; A2= z*CD[(int)(V1+0.5)]*V1*V1 + G; V1=v - deltahalf*A2; A3= z*CD[(int)(V1+0.5)]*V1*V1 + G; V1=v - deltat*A3; A4= z*CD[(int)(V1+0.5)]*V1*V1 + G; a= - (A1 + 2.0*(A2 + A3) + A4)/6.0; if (a < amin) amin = a; deltav = a*deltat; if(-deltav > v) { go_on = FALSE; if(v != 0.0) { h += FM_coast_alt(v, z*CD[(int)(v+0.5)]); t += FM_coast_time(v,z*CD[(int)(v+0.5)]); v=0.0; } }/* if(-deltav) */ else { v += deltav; h += (v + lastv)*deltahalf; 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) */ (*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(); 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 = .001; iopt = 1; isw = 1; printf("\n\r"); 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 = .015; (*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; }