10 REM ** BSAFE-B PROGRAM: THRUST CURVE FROM CHARGE DIMENSIONS ** 20 REM ** FOR BOWL-SHAPED NOZZLES ONLY ** 30 REM ** ALLOWS TOTAL DEGRADATION OF BORE ** 40 REM ** LARRY CURCIO, JUNE 1988 ** 50 REM ============================================================== 60 DEFDBL Z:REM ** REMOVE IF NOT SUPPORTED BY YOUR SYSTEM ** 70 KEY OFF:REM FOR IBM PC AND CLONES. REMOVE FOR OTHER SYSTEMS 80 NITER=500:REM ** NUMBER OF ITERATIONS: ADJUST TO TASTE ** 90 DIM F(NITER),T(NITER) 100 PI=3.14159265358979#:TWOPI=PI+PI:HALFPI=PI/2! 110 SCRNWDTH=80:REM ** CHARACTER WIDTH OF SCREEN ** 120 FLAT=90:FLATCAP=1:ROUND=2 130 REALBIG=10000:REM ** LARGE NUMBER TO AVOID LEXICALLY-UPWARD GOTO'S ** 140 T$=" ":REM INITIALIZE TITLE 150 DEF FNSQ(A)=A*A:REM ****** INTERSECTION OF LINE WITH CIRCLE BELOW ****** 160 DEF FNQUAD1(A,B,C)=(-B+SQR(B*B-4*A*C))/(A+A) 165 DEF FNTEST1(A,B,C)=B*B-4*A*C:REM ** TESTS FOR REAL ROOTS @ PRECISION ** 167 DEF FNTEST0(A,B,R,M,C)=FNTEST1((M*M+1),2*((C-B)*M-A),FNSQ(C-B)-R*R+A*A) 170 DEF FNINTSCT(A,B,R,M,C)=FNQUAD1((M*M+1),2*((C-B)*M-A),FNSQ(C-B)-R*R+A*A) 180 DEF FNK(A,B,R1,R2)=-A*A+R1*R1-R2*R2-B*B 190 REM ****** INTERSECTION OF TWO CIRCLES, CENTER OF ONE IS ORIGIN ****** 200 DEF FNINT2(A,B,K,R)=FNQUAD1(4*(A*A+B*B),4*A*K,(K*K-4*B*B*R*R)) 205 DEF FNTEST(A,B,K,R)=FNTEST1(4*(A*A+B*B),4*A*K,(K*K-4*B*B*R*R)) 210 DEF FNTS(R,RHO,A1,A2)=TWOPI*R*(R*(SIN(A2)-SIN(A1))+RHO*(A2-A1)) 220 DEF FNTV(R,RHO,A1,A2)=TWOPI*R*(R*R*(SIN(A2)-SIN(A1))/3+R*RHO*(A2-A1)/2) 230 DEF FNVNOZZL(R,RHO,A1,A2,HNOZZLE)=FNTV(R,RHO,A1,A2)+PI*RHO*RHO*HNOZZLE 240 REM ********* EXECUTIVE LOOP ********* 250 FOR BIGLOOP=1 TO REALBIG 260 CLS:REM ** CLEAR SCREEN ** 270 GOSUB 340 280 FOR II=1 TO REALBIG 290 PRINT "ANOTHER CURVE";:INPUT A$ 300 IF A$="n" OR A$="N" THEN 2210 310 IF A$="y" OR A$="Y" THEN 330 320 NEXT II 330 NEXT BIGLOOP:REM ********* END EXECUTIVE LOOP ********* 340 PRINT "TITLE ";:INPUT WORK$:IF WORK$>"" THEN T$=WORK$ 350 CENTER=(SCRNWDTH-LEN(T$))/2:CLS:PRINT TAB(CENTER);T$ 360 PRINT "NAME OF OUTPUT FILE ('NUL,' FOR NONE) ";:INPUT F$ 370 IF F$="" THEN F$="NUL":REM ** 'NUL' IS DOS FOR DUMMY FILE ** 380 OPEN "O",1,F$:REM ** YOUR EQUIVALENT MAY VARY ** 390 PRINT "COMBUSTION CHAMBER RADIUS IN CM (";BIGR;") "; 400 INPUT WORK:IF WORK>0 THEN BIGR=WORK 410 PRINT "RADIUS OF FLAT BOTTOM (0=NONE) (";RFLAT;") "; 420 INPUT WORK$:IF WORK$>"" THEN RFLAT=VAL(WORK$) 430 PRINT "TYPE OF CAP: 1=FLAT, 2=ROUNDED (";ICAP;") "; 440 INPUT WORK:IF WORK>0 THEN ICAP=WORK 450 IF ICAP <> FLATCAP AND ICAP <> ROUND THEN 430 460 PRINT "BASE RADIUS IN CM (";R0BOTTOM;") "; 470 INPUT WORK$:IF WORK$>"" THEN R0BOTTOM=VAL(WORK$) 480 PRINT "TOP RADIUS IN CM (";R0TOP;") "; 490 INPUT WORK$:IF WORK$>"" THEN R0TOP=VAL(WORK$) 500 PRINT "BORE LENGTH IN CM (";H;") "; 510 INPUT WORK$:IF WORK$>"" THEN H=VAL(WORK$) 520 PRINT "LENGTH OF CHARGE IN CM (";LENCHRGE;") "; 530 INPUT WORK:IF WORK>0 THEN LENCHRGE=WORK 540 PRINT "TOTAL IMPULSE IN NEWTON SECONDS (";TOTIMP;") "; 550 INPUT WORK:IF WORK>0 THEN TOTIMP=WORK 560 PRINT "PROPELLANT MASS IN GRAMS (";MPROP;") "; 570 INPUT WORK:IF WORK>0 THEN MPROP=WORK 580 REM ** INITIAL CALCULATIONS ** 590 THETA=0:IF H<>0 THEN THETA=ATN((R0BOTTOM-R0TOP)/H) 600 COSTHETA=COS(THETA):SINTHETA=SIN(THETA):TANTHETA=SINTHETA/COSTHETA 610 M=0:IF TANTHETA<>0 THEN M=-1/TANTHETA 620 PRINT "BURNING TIME IN SECONDS (";TTOT;") "; 630 BIGR1=BIGR-R0TOP:BIGR2=BIGR-R0BOTTOM:TOPA=PI*R0TOP*R0TOP 640 T=0:RBT=0:RBT2=0:RBT0=0:RHOTOP=R0TOP:IF ICAP=FLATCAP THEN 660 650 RBT0=R0TOP/COSTHETA:TOPA=0:RHOTOP=0:BIGR1=BIGR:RBT2=RBT0*RBT0 660 INPUT WORK$:IF WORK$>"" THEN TTOT=VAL(WORK$) 670 REM==================================================================== 680 L=R0BOTTOM-R0TOP:L=SQR(H*H+L*L) 690 DELTAT=TTOT/NITER 700 VE=TOTIMP/MPROP:REM ** EFFECTIVE EXHAUST VELOCITY ** 710 REM ** DEVELOPING LENGTH OF FORWARD PROPELLANT CAP ** 720 SPHEREHT=0:REM ** DEVELOPING HEIGHT OF SPERICAL SEGMENT ** 730 IF ICAP=ROUND AND RBT0<>0 THEN SPHEREHT=RBT0-SQR(RBT0*RBT0-R0TOP*R0TOP) 740 LENCAP=LENCHRGE+RBT0-SPHEREHT-H 750 RB=(SQR(LENCAP*LENCAP+BIGR1*BIGR1)-RBT0)/TTOT 760 VCHAMBER=PI*BIGR*BIGR*LENCHRGE:REM ** VOLUME OF CHAMBER ** 770 REM ** ADJUST VOLUME FOR NOZZLE ** 780 RROUND=BIGR-RFLAT:REM ** ROUND PART RADIUS & HEIGHT OF NOZZLE ** 790 NOZANGLE=HALFPI:IF RROUND=0 THEN 810 800 NOZANGLE=ATN((BIGR-R0BOTTOM)/RROUND) 810 NOZANG1=HALFPI-NOZANGLE:NOZANG10=NOZANG1:REM ** HOLD ORIGINAL ANGLE ** 820 LNOZZLE=BIGR2:IF RROUND>0 THEN LNOZZLE=SQR(RROUND*RROUND+BIGR2*BIGR2) 830 IF NOZANG1>THETA THEN 850 840 BIGR3=BIGR2/SIN(NOZANGLE):GOTO 880:REM** IF NOZZLE DOESN'T SHORTEN WALL ** 850 IF L/COS(NOZANGLE+THETA)0 THEN 1050 1040 MAXL=L-RROUND:GOTO 1055 1050 MAXL=L:IF NOZANG1>THETA THEN MAXL=(BIGR-MAXTR)/SINTHETA 1055 IF MAXL>0 THEN 1060:REM ** CHECK FOR CAP TOUCHDOWN ** 1056 PSI=(RROUND-H)/(RB*TMAXF+RBT0):PSI=ATN(PSI/SQR(1-PSI*PSI)) 1057 MAXA=FNTS(RBT0+RB*TMAXF,RHOTOP,PSI,HALFPI)+TOPA:GOTO 1080 1060 MAXA=FNTS(RBT0+RB*TMAXF,RHOTOP,THETA,HALFPI)+PI*MAXL*(MAXBR+MAXTR)+TOPA 1070 IF NOZANG11.5 THEN 1470 1220 NOZANGLE=HALFPI:NOZANG1=0:IF RBT+R0BOTTOM<=RFLAT THEN 1470 1230 REM ** INTERSECTION OF FRUSTUM CROSS-SECTION WITH BOWL **** 1240 IF BIGR3BIGR THEN X=BIGR:REM ** ELIMINATES EXPLOSIONS FOR PRECISION ** 1300 Y=RROUND-SQR(RROUND*RROUND-FNSQ(X-RFLAT)):X=X-R0BOTTOM 1310 L1=H-Y:IF SINTHETA>0 THEN L1=(X-RTOP+R0BOTTOM)/SINTHETA 1320 IF L1>L THEN 1370 1330 L=L1:IF L<0 THEN 1420 1340 IF SINTHETA=0 THEN 1390 1350 Y=X*M+H+RBT/SINTHETA-M*R0TOP 1360 GOTO 1390 1370 IF TEST>=0 THEN 1376 1372 TEST=FNTEST(RFLAT-R0BOTTOM,RROUND,FNK(RFLAT-R0BOTTOM,RROUND,RROUND,RBT),RBT) 1374 X=RBT+RBT0:IF TEST<0 THEN 1380 1376 X=FNINT2(RFLAT-R0BOTTOM,RROUND,FNK(RFLAT-R0BOTTOM,RROUND,RROUND,RBT),RBT) 1380 IF X<=RBT THEN Y=SQR(RBT*RBT-X*X) 1390 NOZANGLE=HALFPI:IF Y<>0 THEN NOZANGLE=ATN(X/Y) 1400 NOZANG1=HALFPI-NOZANGLE 1410 GOTO 1470 1420 IF TEST>=0 THEN 1430 1423 TEST=FNTEST(RFLAT-RHOTOP,RROUND-H,FNK(RFLAT-RHOTOP,RROUND-H,RROUND,RBT+RBT0),RBT+RBT0) 1426 X=RBT+RBT0:IF TEST<0 THEN 1440 1430 X=FNINT2(RFLAT-RHOTOP,RROUND-H,FNK(RFLAT-RHOTOP,RROUND-H,RROUND,RBT+RBT0),RBT+RBT0) 1440 Y=SQR(FNSQ(RBT+RBT0)-X*X) 1450 L=0:MARK=1.5:PSI=0:IF X<>0 THEN PSI=ATN(Y/X) 1460 A=FNTS(RBT+RBT0,RHOTOP,PSI,HALFPI)+TOPA:GOTO 1720 1470 IF BIGRNOZANG1 THEN A=A+FNTS(RBT,R0BOTTOM,NOZANG1,THETA) 1520 IF MARK<1 THEN MARK=1 1530 GOTO 1720 1540 REM ** TRAILING SPHERE HITS BUT FRUSTUM BASE IS CLEAR: COWLICK ** 1550 IF THETA<=NOZANG10 THEN 1720 1560 THETA1=ATN(SQR(RBT*RBT-BIGR2*BIGR2)/BIGR2) 1570 A=A+FNTS(RBT,R0BOTTOM,THETA1,THETA):MARK=2 1580 GOTO 1720 1590 REM ** FRUSTUM BASE HAS HIT BUT TOP HAS NOT: DESCENDING PORTION ** 1600 IF BIGR=F THEN 1760 1750 DMAXF=F:DTMAXF=T 1760 T=T+DELTAT:RBT=RB*T:RBT2=(RBT+RBT0):RBT2=RBT2*RBT2:ZDIGIMP=ZDIGIMP+F 1770 NEXT I 1780 CLOSE #1:REM USE EQUIVALENT ON YOUR SYSTEM ** 1790 REM ===================================================================== 1800 REM ** PLOT ROUTINE: PC-DEPENDENT ** 1810 REM ** OMIT OR REWRITE FOR OTHER SYSTEMS ** 1820 YSIZE=200:XSIZE=640:MAXY=YSIZE-1:MAXX=XSIZE-1 1830 YSCALE=DMAXF:XSCALE=TTOT 1840 ZDIGIMP=DELTAT*(ZDIGIMP-.5*F(O)-.5*F):REM TRAPAZOIDAL IMPULSE 1850 DIGIMP=ZDIGIMP:REM ** CONVERT TO SINGLE, AVOID 'PRINTUSING' ** 1860 PRINT "DIGITAL MAXIMUM THRUST, TIME= ";DMAXF;" ";DTMAXF 1870 PRINT "DIGITAL IMPULSE = ";DIGIMP 1880 PRINT:PRINT "PRESS TO CONTINUE ";:INPUT"",A$ 1890 FOR II=1 TO REALBIG 1900 CLS 1910 PRINT "MAX X=";TTOT;" MAX Y=";DMAXF:PRINT:PRINT 1920 PRINT "XSCALE=";XSCALE;" YSCALE=";YSCALE:PRINT:PRINT 1930 PRINT"ENTER NEW XSCALE OR PRESS RETURN TO ACCEPT":INPUT GGG:PRINT 1940 IF GGG<>0 THEN XSCALE=GGG 1950 PRINT"ENTER NEW YSCALE OR PRESS RETURN TO ACCEPT":INPUT GGG 1960 IF GGG<>0 THEN YSCALE=GGG 1970 GOSUB 2050 1980 FOR JJ=1 TO REALBIG 1990 PRINT "REPLOT ";:INPUT A$ 2000 IF A$="y" OR A$="Y" THEN 2030 2010 IF A$="n" OR A$="N" THEN 2040 2020 NEXT JJ 2030 NEXT II 2040 RETURN 2050 SCREEN 0:SCREEN 2:REM ** SWITCH TO GRAPHICS WITH FAILSAFE CLEAR ** 2060 LINE (0,0)-(0,MAXY):LINE -(MAXX,MAXY):REM ** DRAW AXES ** 2070 X1=0:Y1=0 2080 REM ** MAKE 10 TIC MARKS ON EACH AXIS ** 2090 FOR I=1 TO 10 2100 PRESET(0,(MAXY-Y1)):PRESET(X1,(MAXY-0)) 2110 X1=X1+XSIZE/10:Y1=Y1+YSIZE/10 2120 NEXT I 2130 PSET(MAXX*T(0)/XSCALE,MAXY-(MAXY*F(0))/YSCALE) 2140 FOR I=1 TO NITER 2150 X=T(I):Y=F(I):REM ** DRAW LINE FROM LAST POINT TO THIS POINT ** 2160 LINE-(MAXX*X/XSCALE,MAXY-(MAXY*Y)/YSCALE) 2170 NEXT I 2180 A$=INKEY$:IF A$="" THEN 2180:REM ** HOLD GRAPH UNTIL KEY PRESSED ** 2190 SCREEN 0:REM ** FLIP BACK TO TEXT SCREEN ** 2200 RETURN 2210 END