10 REM ** BSAFE THRUST CURVE FROM CHARGE DIMENSIONS ** 20 REM ** LARRY CURCIO, JUNE 1988 ** 30 REM ============================================================== 40 DEFDBL Z:REM ** REMOVE IF NOT SUPPORTED BY YOUR SYSTEM ** 50 KEY OFF:REM FOR IBM PC AND CLONES. REMOVE FOR OTHER SYSTEMS 60 NITER=500:REM ** NUMBER OF ITERATIONS: ADJUST TO TASTE ** 70 DIM F(NITER),T(NITER) 80 PI=3.14159265358979#:TWOPI=PI+PI:HALFPI=PI/2! 90 SCRNWDTH=80:REM ** CHARACTER WIDTH OF SCREEN ** 100 FLAT=90:FLATCAP=1:ROUND=2:INZ=90 110 REALBIG=10000:REM ** LARGE NUMBER TO AVOID LEXICALLY-UPWARD GOTO'S ** 120 T$=" ":REM INITIALIZE TITLE 130 DEF FNTS(R,RHO,A1,A2)=TWOPI*R*(R*(SIN(A2)-SIN(A1))+RHO*(A2-A1)) 135 DEF FNQUAD(A,B,C)=(-B+SQR(B*B-4*A*C))/(A+A) 137 DEF FNSQ(A)=A*A 140 REM ********* EXECUTIVE LOOP ********* 150 FOR BIGLOOP=1 TO REALBIG 160 CLS:REM ** CLEAR SCREEN ** 170 GOSUB 240 180 FOR II=1 TO REALBIG 190 PRINT "ANOTHER CURVE";:INPUT A$ 200 IF A$="n" OR A$="N" THEN 1920 210 IF A$="y" OR A$="Y" THEN 230 220 NEXT II 230 NEXT BIGLOOP:REM ********* END EXECUTIVE LOOP ********* 240 PRINT "TITLE ";:INPUT WORK$:IF WORK$>"" THEN T$=WORK$ 250 CENTER=(SCRNWDTH-LEN(T$))/2:CLS:PRINT TAB(CENTER);T$ 260 PRINT "NAME OF OUTPUT FILE ('NUL,' FOR NONE) ";:INPUT F$ 270 IF F$="" THEN F$="NUL":REM ** 'NUL' IS DOS FOR DUMMY FILE ** 280 OPEN "O",1,F$:REM ** YOUR EQUIVALENT MAY VARY ** 290 PRINT "COMBUSTION CHAMBER RADIUS IN CM (";BIGR;") "; 300 INPUT WORK:IF WORK>0 THEN BIGR=WORK 310 PRINT "NOZZLE ANGLE (0=NORMAL TO FRUSTUM) (";INZ;") "; 320 INPUT WORK$:IF WORK$>"" THEN INZ=VAL(WORK$) 330 PRINT "TYPE OF CAP: 1=FLAT, 2=ROUNDED (";ICAP;") "; 340 INPUT WORK:IF WORK>0 THEN ICAP=WORK 350 IF ICAP <> FLATCAP AND ICAP <> ROUND THEN 330 360 PRINT "BASE RADIUS IN CM (";R0BOTTOM;") "; 370 INPUT WORK$:IF WORK$>"" THEN R0BOTTOM=VAL(WORK$) 380 PRINT "TOP RADIUS IN CM (";R0TOP;") "; 390 INPUT WORK$:IF WORK$>"" THEN R0TOP=VAL(WORK$) 400 PRINT "BORE LENGTH IN CM (";H;") "; 410 INPUT WORK$:IF WORK$>"" THEN H=VAL(WORK$) 420 PRINT "LENGTH OF CHARGE IN CM (";LENCHRGE;") "; 430 INPUT WORK:IF WORK>0 THEN LENCHRGE=WORK 440 PRINT "TOTAL IMPULSE IN NEWTON SECONDS (";TOTIMP;") "; 450 INPUT WORK:IF WORK>0 THEN TOTIMP=WORK 460 PRINT "PROPELLANT MASS IN GRAMS (";MPROP;") "; 470 INPUT WORK:IF WORK>0 THEN MPROP=WORK 480 REM ** INITIAL CALCULATIONS ** 490 THETA=0:IF H<>0 THEN THETA=ATN((R0BOTTOM-R0TOP)/H) 500 COSTHETA=COS(THETA):SINTHETA=SIN(THETA) 510 PRINT "BURNING TIME IN SECONDS (";TTOT;") "; 520 BIGR1=BIGR-R0TOP:BIGR2=BIGR-R0BOTTOM:TOPA=PI*R0TOP*R0TOP 530 T=0:RBT=0:RBT2=0:RBT0=0:RHOTOP=R0TOP:IF ICAP=FLATCAP THEN 550 540 RBT0=R0TOP/COSTHETA:TOPA=0:RHOTOP=0:BIGR1=BIGR:RBT2=RBT0*RBT0 550 INPUT WORK$:IF WORK$>"" THEN TTOT=VAL(WORK$) 560 REM====================================================================== 570 L=R0BOTTOM-R0TOP:L=SQR(H*H+L*L) 580 DELTAT=TTOT/NITER 590 VE=TOTIMP/MPROP:REM ** EFFECTIVE EXHAUST VELOCITY ** 600 REM ** DEVELOPING LENGTH OF FORWARD PROPELLANT CAP ** 610 SPHEREHT=0:REM ** DEVELOPING HEIGHT OF SPERICAL SEGMENT ** 620 IF ICAP=ROUND AND RBT0<>0 THEN SPHEREHT=RBT0-SQR(RBT0*RBT0-R0TOP*R0TOP) 630 LENCAP=LENCHRGE+RBT0-SPHEREHT-H 640 RB=(SQR(LENCAP*LENCAP+BIGR1*BIGR1)-RBT0)/TTOT 650 REM MIN DISTANCE FOR CAP TO TOUCH WALL 660 VCHAMBER=PI*BIGR*BIGR*LENCHRGE:REM ** VOLUME OF CHAMBER ** 670 REM ** ADJUST VOLUME FOR NOZZLE ** 680 NOZANGLE=TWOPI*INZ/360:IF NOZANGLE=0 THEN NOZANGLE=HALFPI-THETA 690 DELTAL=0:SINNOZ=0 700 DIFF=HALFPI-NOZANGLE-THETA:TANDIFF=0:IF DIFF>0 THEN TANDIFF=TAN(DIFF) 710 IF TANDIFF<>0 THEN H1=H+(BIGR-R0TOP)*TAN(THETA) 780 NOZANG1=HALFPI-NOZANGLE:SINNOZ=SIN(NOZANGLE):BIGR3=BIGR2/SINNOZ 790 VNOZZLE=0:HNOZZLE=0:MAXL=L 800 IF INZ=FLAT THEN 870 810 IF INZ=0 AND THETA=0 THEN 870 820 HNOZZLE=BIGR2/TAN(NOZANGLE):IF TANDIFF=0 THEN 850 825 TANPHI=TAN(NOZANG1) 830 BIGR3=SIN(THETA+NOZANGLE)*SQR(BIGR2*BIGR2+HNOZZLE*HNOZZLE) 841 REM ** FOR MAXL, SEE DELTAL FORMULA IN MAIN LOOP. NOTE BIGR3=RBT*TMAXF ** 842 MAXL=L-BIGR3*TANDIFF 845 IF MAXL<0 THEN BIGR3=SQR(FNSQ(BIGR-RHOTOP)+FNSQ(HNOZZLE-H))+RHOTOP-R0TOP 850 VNOZZLE=PI*(BIGR*BIGR+R0BOTTOM*R0BOTTOM+R0BOTTOM*BIGR)*HNOZZLE/3! 860 VCHAMBER=VCHAMBER-PI*BIGR*BIGR*HNOZZLE+VNOZZLE:REM ** CORRECTED VOLUME ** 870 REM ** DEVELOP BORE VOLUME ** 880 VBORE=0 890 VBORE=PI*SPHEREHT*SPHEREHT*(3!*RBT0-SPHEREHT)/3!:REM FROM CRC MATH TABLES 900 VBORE=VBORE+PI*H*(R0BOTTOM*R0BOTTOM+R0TOP*R0TOP+R0TOP*R0BOTTOM)/3! 910 VPROP=VCHAMBER-VBORE 920 RHOPROP=MPROP/VPROP:REM ** PROPELLANT DENSITY ** 930 BIGK=VE*RB*RHOPROP:REM ** CONVERSION CONSTANT BETWEEN AREA AND THRUST ** 940 VE=VE*1000!:REM ** CONVERT EXHAUST VELOCITY TO METERS/SECOND ** 950 ZDIGIMP=0:REM ** INITIALIZE DIGITAL IMPULSE (DOUBLE PRECISION) ** 960 TMAXF=BIGR3/RB:REM ** TIME OF MAXIMUM THRUST ** 970 REM ** DEVELOP MAX BOTTOM RADIUS, TOP RADIUS, SURFACE AREA ** 980 MAXBR=BIGR:IF TANDIFF=0 THEN MAXBR=R0BOTTOM+RB*TMAXF*COSTHETA 990 MAXTR=R0TOP+RB*TMAXF*COSTHETA 991 PSI=THETA 992 IF MAXL>0 THEN 1000 993 COSPSI=(BIGR-RHOTOP)/(RB*TMAXF+RBT0) 994 IF ABS(COSPSI-1)<1E-5 THEN COSPSI=1 995 SINPSI=SQR(1-COSPSI*COSPSI) 996 PSI=ATN(SINPSI/COSPSI) 997 MAXL=0 1000 MAXA=FNTS(RBT0+RB*TMAXF,RHOTOP,PSI,HALFPI)+PI*MAXL*(MAXBR+MAXTR)+TOPA 1010 IF INZ<>0 AND TANDIFF=0 THEN MAXA=MAXA+FNTS(BIGR3,R0BOTTOM,NOZANG1,THETA) 1020 ASYMPA=PI*BIGR*BIGR:REM ** ASYMPTOTIC THRUST FROM CHAMBER AREA *** 1030 PRINT:PRINT "MAXIMUM THRUST IN NEWTONS= ";MAXA*BIGK 1040 PRINT "TIME OF MAXIMUM THRUST IN SECONDS= ";TMAXF 1050 PRINT "ASYMPTOTIC THRUST IN NEWTONS= ";ASYMPA*BIGK 1060 PRINT "PROPELLANT DENSITY IN GM/CM^3= ";RHOPROP 1070 PRINT "PROPELLANT VOLUME IN CM^3= ";VPROP 1080 PRINT "BURNING RATE IN CM/SEC= ";RB 1090 PRINT "LENGTH OF FORWARD CAP IN CM= ";LENCAP-RBT0 1100 PRINT "EXHAUST VELOCITY IN M/SEC= ";VE 1102 PRINT#1,BIGR;" ";INZ;" ";ICAP;" ";R0BOTTOM;" ";R0TOP;" "; 1103 PRINT#1,H;" ",LENCHRGE;" BSAFEX" 1105 PRINT#1,TOTIMP,MPROP,TTOT,NITER 1110 REM================= ** MAIN LOOP ** =================================== 1120 FOR I=0 TO NITER 1130 RTOP=R0TOP+RBT*COSTHETA:RBOTTOM=R0BOTTOM+RBT*COSTHETA 1135 IF MARK>2.5 THEN 1169 1140 IF TANDIFF=0 THEN 1170 1141 DELTAL=RBT*TANDIFF 1142 IF DELTAL<=L THEN 1169 1143 H1=H+(R0BOTTOM-RHOTOP)*TANPHI 1144 COSPSI=FNQUAD((TANPHI*TANPHI+1)*RBT2,-2*TANPHI*H1*(RBT+RBT0),H1*H1-RBT2) 1146 SINPSI=SQR(1-COSPSI*COSPSI):TANPSI=SINPSI/COSPSI:PSI=ATN(TANPSI) 1149 RTOP=(RBT+RBT0)*COSPSI+RHOTOP:IF RTOP>BIGR THEN 1370 1152 A=FNTS(RBT+RBT0,RHOTOP,PSI,HALFPI)+TOPA:MARK=2.5 1155 GOTO 1450 1169 RBOTTOM=SQR(DELTAL*DELTAL+RBT*RBT)*SINNOZ+R0BOTTOM 1170 REM ** BEFORE WALLS HIT AT BOTTOM ** 1180 IF BIGR0 THEN 1450 1230 REM ** BEFORE WALLS HIT BY TRAILING SPHERE: ASCENDING PORTION ** 1240 IF BIGR30 THEN A=A+FNTS(RBT,R0BOTTOM,NOZANG1,THETA) 1260 GOTO 1450 1270 REM ** TRAILING SPHERE HITS BUT FRUSTUM BASE IS CLEAR: COWLICK ** 1280 IF INZ=0 THEN 1450:REM THIS SHOULDN'T BE REQUIRED 1290 THETA1=ATN(SQR(RBT*RBT-BIGR2*BIGR2)/BIGR2) 1300 A=A+FNTS(RBT,R0BOTTOM,THETA1,THETA):MARK=2 1310 GOTO 1450 1320 REM ** FRUSTUM BASE HAS HIT BUT TOP HAS NOT: DESCENDING PORTION ** 1330 IF BIGR TO CONTINUE ";:INPUT"",A$ 1600 FOR II=1 TO REALBIG 1610 CLS 1620 PRINT "MAX X=";TTOT;" MAX Y=";MAXA*BIGK:PRINT:PRINT 1630 PRINT "XSCALE=";XSCALE;" YSCALE=";YSCALE:PRINT:PRINT 1640 PRINT"ENTER NEW XSCALE OR PRESS RETURN TO ACCEPT":INPUT GGG:PRINT 1650 IF GGG<>0 THEN XSCALE=GGG 1660 PRINT"ENTER NEW YSCALE OR PRESS RETURN TO ACCEPT":INPUT GGG 1670 IF GGG<>0 THEN YSCALE=GGG 1680 GOSUB 1760 1690 FOR JJ=1 TO REALBIG 1700 PRINT "REPLOT ";:INPUT A$ 1710 IF A$="y" OR A$="Y" THEN 1740 1720 IF A$="n" OR A$="N" THEN 1750 1730 NEXT JJ 1740 NEXT II 1750 RETURN 1760 SCREEN 0:SCREEN 2:REM ** SWITCH TO GRAPHICS WITH FAILSAFE CLEAR ** 1770 LINE (0,0)-(0,MAXY):LINE -(MAXX,MAXY):REM ** DRAW AXES ** 1780 X1=0:Y1=0 1790 REM ** MAKE 10 TIC MARKS ON EACH AXIS ** 1800 FOR I=1 TO 10 1810 PRESET(0,(MAXY-Y1)):PRESET(X1,(MAXY-0)) 1820 X1=X1+XSIZE/10:Y1=Y1+YSIZE/10 1830 NEXT I 1840 PSET(MAXX*T(0)/XSCALE,MAXY-(MAXY*F(0))/YSCALE) 1850 FOR I=1 TO NITER 1860 X=T(I):Y=F(I):REM ** DRAW LINE FROM LAST POINT TO THIS POINT ** 1870 LINE-(MAXX*X/XSCALE,MAXY-(MAXY*Y)/YSCALE) 1880 NEXT I 1890 A$=INKEY$:IF A$="" THEN 1890:REM ** HOLD GRAPH UNTIL KEY PRESSED ** 1900 SCREEN 0:REM ** FLIP BACK TO TEXT SCREEN ** 1910 RETURN 1920 END