mathproc.dll FORTRAN Engineering Subprograms
(Calculating Model Rocket Coefficient of Normal Force & Center of Pressure) ¹
(Includes Downloads)

tir33.for is a collection of FORTRAN subprograms extracted from mathproc.dll to calculate and output Coefficient of Normal Force and Center of Pressure for an aerodynamically stabilized rocket.  The subprograms' passed parameter list is documented within the source code.  The source code file may be combined with mathaux.for.  See Creating FTN95 Project and mathproc Project Properties for setting up the project file and compiling.

  α  : angle of attack, shown as a (Displacement angle) in diagram below;  angle between Flight Direction and rocket's longitudinal axis (referred to as Axis of Symmetry in diagram below)
 cg  : Center of Gravity, rocket's "pivot" point
 cp  : Center of Pressure, point where aerodynamic forces act; calculated by tir33.for subprograms extracted from mathproc.dll for each rocket component;  Technical Information Report TIR-33 refers to each component's relative cp location as ∆x_{component subscript};
x_{component subscript} is the distance from the reference line (usually nose tip) to the component.

  D  : Aerodynamic Drag Force, acts parallel to rocket's Flight Direction; a type of friction acting opposite to the relative motion of an air vehicle flying through the air
  L  : Aerodynamic Lift Force, acts perpendicular to rocket's Flight Direction
  N  : Aerodynamic Normal Force, acts perpendicular to rocket's longitudinal axis (referred to as Axis of Symmetry in diagram below);
when α = 0, N = 0
  A  : Aerodynamic Axial Force, acts parallel to rocket's longitudinal axis

C_D  : Drag Coefficient, dimensionless quantity that is used to quantify all of the complex dependencies of shape, inclination, and flow conditions on air vehicle drag
C_L  : Lift Coefficient, dimensionless quantity describing lift characteristics of a particular shape (or cross-section) of an air vehicle ( C_Lα : Lift Coefficient at angle of attack α)
C_Nα : Coefficient of Normal Force at angle of attack α (sometimes referred to as Normal Force Coeffcient derivative), dimensionless quantity based on shape; calculated by tir33.for subprograms extracted from mathproc.dll for each rocket component

N = CNα ½ ρ V2 α Ar ρ  : Air density V : Airspeed Ar : Reference area indicating size of rocket, usually nose base cross-sectional area; nose base diameter is also used in calculating conical transition CNα When α is small (measured in radians), it is common in Aerodynamics to use α in equations instead of sin (α). Using the simple FORTRAN function at right, can calculate N with the statement:    N = FAREO (CNA, RHO, V, AR) * ALPHA With the appropriate arguments, this function can also be used to calculate D (Do not inlude α), L, etc.

dynamic pressure = ½ ρ V2

L = N cos α - A sin α
D = N sin α + A cos α = C_D ½ ρ V² A_r

D, L, N, and A are vectors;  the tir33.for FORTRAN subprograms calculate the cp and C_Nα.   cp locations and C_Nα must be calculated for each component (using tir33.for) and then summarized (using missle02.pas - see below; requires include files genetyp.inc², nrtyp.inc, and missle.inc).  The components are nose cone, conical shoulder(s), conical boattail(s) (conical shoulders and conical boattails are conical transitions), fins, etc.  For an airplane or rocket to be aerodynamically stable, the cp must lie behind the cg along the longitudinal axis.

3

tir33.for Subprograms
Uses variable bounds one dimensional arrays for the fin data; statically dimensioned two dimensional matrices for the conical transition data, for a maximum of 6 serial stages and 10 conical transitions per serial stage.  Practical design limitations are well below these limits.
Can be called by either FORTRAN or Pascal console programs; can also be called from Windows applications (except for D2R2PRNT (...), PRNTCPFINS (...), and PRNTCPCONTRANS (...)).  To call from a Lazarus application, use missle02.pas procedure CmptMissleCP (...).  Application shown in Lazarus Tab Sheets calls CmptMissleCP (...), which in turns calls CPMRBODY (...) and CPMRFINS (...) in tir33.for and then summarizes the data.

missle.inc, missle02.pas, and tir33.for (including passed parameters) are commented

REAL FUNCTION RLEN (X1, X2, Y1, Y2, Z1, Z2) RESULT (D), REAL FUNCTION MODULUS (A, B, C) RESULT (X) 4 called by SRCONE (...) REAL FUNCTION ATRAPEZOID (A, B, H) RESULT (AR), REAL FUNCTION AELLIPSE (A, B, X, Y) RESULT (AR), REAL FUNCTION SRCYLNDR (R, H) RESULT (S), REAL FUNCTION SRCONE (RB, H) RESULT (S), REAL FUNCTION SRPARABOLIOD (RB, H) RESULT (S), REAL FUNCTION SPTRAPEZOID (A, B, C, D, H, T, BEVEL) RESULT (S), REAL FUNCTION PELLIPSE (A, B) RESULT (P) REAL FUNCTION SPELLIPSE (A, B, X, Y, T, BEVEL) RESULT (S) Calculate planar/surface area5 of various shapes, used in calculating CD and D, called by missle02.dll; not used in calculating CNα and  cp SUBROUTINE D2R2PRNT (A, N, M, IERR) 6 General Purpose Print Subroutine for use in Command Console programs, called by PRNTCPFINS (...) and PRNTCPCONTRANS (...); can be used by both FORTRAN and Pascal command console programs.  Useful for outputting arrays and matrices. Do not call, directly or indirectly, from a Windows program or D command console program.  Calling subprogram passes array/matrix data through matrix A. REAL FUNCTION CPCONICAL (LN) RESULT (XN) Computes the cp location, measured from the nose tip reference line, for a conical shaped nose; called by CPNOSE (...) REAL FUNCTION CPOGIVE (LN) RESULT (XN) Computes the cp location for an ogive shaped nose; called by CPNOSE (...) REAL FUNCTION CPPARABOLIC (LN) RESULT (XN) Computes the cp location for a Parabolic shaped nose; called by CPNOSE (...) SUBROUTINE CPCONTRANS (IERR, L, LCT, D, R1, R2, CNCT, XCT) Computes the cp location and CNα for a single conical transition; called by CPNOSE (...) (for "Honest John" shaped nose); also called by CPMRBODY (...) and PRNTCPCONTRANS (...); FORTRAN allows passing individual subscripted array elements by reference; a conical shoulder produces a positive CNα (increasing stability when aft of the cg), a conical boattail produces a negative CNα (decreasing stability when aft of the cg); each stage may have zero or more conical transitions SUBROUTINE CPNOSE (NSHAPE, IERR, LN, LCB, D, R1, R2, CNN, XN, CNCB, XCB) Computes the cp location and CNα for a variety of nose cone shapes (including "Mercury" capsule 7 and "Honest John" 8 surface-to-surface missile special shapes); determines nose base diameter based on nose shape; called by CPMRBODY (...) SUBROUTINE CPFINS (L, D, R, LROOT, SPAN, LCHORD, M, LTIP, CNFB, XFB, NFINS, IERR) Computes the fin/body cp location and CNα for a single stage; called by PRNTCPFINS (...) and CPMRFINS (...); FORTRAN allows passing individual subscripted array elements by reference SUBROUTINE PRNTCPFINS (MNAME, NSTAGE, D, L, R, LROOT, SPAN, LCHORD, M, LTIP, CNFB, XFB, NFINS, IERR) Can be called by either FORTRAN or Pascal Console Programs, primarily used for testing; calls CPFINS (...) for each stage and outputs the results on the screen SUBROUTINE CPMRBODY (NSTAGE, NCT, NSHAPE, IERR, LN, LCB, D, RN1, RN2, CNN, XN, LNTCT, LCT, R1, R2, CNCT, XCT) The nose, body tubes, and conical transitions comprise the fuselage; Computes the cp location and CNα for the nose by calling CPNOSE (...) and for the conical transitions by calling CPCONTRANS (...) for all conical transitions in every stage (when developing main programs directly/indirectly calling this subroutine, pay extra attention to combined conical transition index and "Honest John" shaped nose processing, if any, on the first stage - certain operations will need to decrement nct[1] after calling CPMRBODY (...); mult-stage w/ "Honest John" nose is still being tested and tweaked; recommend downloading tir33.for periodically for any updates); for "Honest John" 8 idealization, CPMRBODY (...) temporarily adds the conical transition section to the conical transition matrices SUBROUTINE CPMRFINS (NSTAGE, D, L, R, LROOT, SPAN, LCHORD, M, LTIP, CNFB, XFB, NFINS, IERR) Computes the fin/body cp location and CNα for all serial stages by calling CPFINS (...) for each stage SUBROUTINE PRNTCPCONTRANS (MNAME, NSTAGE, NCT, IERR, D, L, LCT, R1, R2, CNCT, XCT) Can be called by either FORTRAN or Pascal Console Programs, primarily used for testing; calls CPCONTRANS (...) for each conical transition for each stage and outputs the results on the screen

Single Stage Rocket Shown Below Stage 1 is the uppermost stage Stage n is the lowermost stage (Indexed from Nose Tip)   d : Reference Diameter cs : Conical Shoulder subscript cb : Conical Boattail subscript   f  : Fin subscript    indicates cp location

1

missle02.pas:

missle02.pas contains model rocket related functions/procedures⁹ compiled into a .dll.  These routines have been successfully called by Object Pascal console programs, Object Pascal Windows programs, and Lazarus programs.  Uses include files genetyp.inc, nrtyp.inc, and missle.inc.

Application shown in Lazarus Tab Sheets calls missle02 procedure CmptMissleCP (...), which in turns calls CPMRBODY (...) and CPMRFINS (...) in tir33.for (mathproc.for) and then summarizes the data via calls to other procedures in missle02.dll: cp = ∑ (CNα[i] * cp[i] ) / ∑ CNα[i]    cp   : overall Center of Pressure CNα[i] : component's Coefficient of Normal             Force  cp[i]  : component's Center of Pressure It is necessary to go through this calculation for the 1st stage launch configuration through the nth stage launch configuration, calculating each stage separately. Warning - the external declaration for CPMRBODY (...) and CPMRFINS (...) refers to mathproc.dll, not tir33.dll; will need to either change mathproc.dll external reference to tir33.dll, or rename tir33.for to mathproc.for.10 procedure missleErrProc (...) may call HWND_APPND_TEXT (...) in hwndio.pas; will need to download and compile into hwndio.dll; see Appending Text to Windows Edit Control for more information.

Use fpcdll.bat to compile; open a Command Prompt (on Win7 - Start, All Programs, Accessories) and type
fpcdll missle02

1. Barrowman, James. TIR-33 Calculating the Center of Pressure
2. Rietman, Edward (1994). Genesis Redux: Experiments Creating Artificial Life. Windcrest (McGraw-Hill). Include file developed from major rewrite of Dawkins8 Biomorphs (Free Pascal dawkins8.dll callable from both Pascal and D programs); contains useful CONSTant and TYPE declarations
3. NASA Rocket Stability
4. http://www.chem.ualberta.ca/~abrown/fortran/section4.html
5. Also requires functions XROTATECW (...), YROTATECW (...), XROTATECCW (...), YROTATECCW (...), and TRIAREA (...) from mathaux.for; called by missle02.dll function FlghtSrfcArea (...) when computing swept back elliptical fin/wing surface area.  Uses ELLIPSE_SEGMENT algorithm from:
Hughes, Gary B., and Mohcine Chraibi (2011).  Calculating Ellipse Overlap Areas
FlghtSrfcArea (...) calculates fin/wing surface area, used in calculating C_D and D
6. Vetterling, William T., Saul A Teukolsky, William H. Press, and Brian P. Flannery (1989). Numerical Recipes: Example Book (Pascal) Revised Edition New York:Press Syndicate of the University of Cambridge.
Extracted from Pascal examples and ported to FORTRAN

7. "Mercury" capsule shaped nose cones must be idealized prior to analysis. CPNOSE (...) idealizes a "Mercury" capsule shaped nose (shown at right, heavy black line) into a conical shape (as shown by dashed line).  The nose tip radius (radius at REFERENCE LINE) and nose base radius is known.  Nose actual length (L - ∆L) is known.  As an alternative to drawing an outline and measuring L and ∆L as given on p. 9 of TIR-33, CPNOSE (...) uses Similar Triangles to calculate both L and ∆L. L = LN + ∆L tan θ = ∆L / R1 = L / R2 = (LN + ∆L) / R2 ∆L = (R1 LN) / (R2 - R1)  θ  : angle between nose tip (line segment 2 3) and dashed line;        angle between nose base (line segment 1 4) and dashed line LN : Actual length of nose (distance from reference line to nose base) R1 : 1st radius from reference line (nose tip radius) R2 : 2nd radius from reference line (nose base radius) CPNOSE (...) calculates the cp for a conical nose of length L and then subtracts ∆L

1

8. "Honest John" shaped nose is broken down into an ogive nose component and a conical boattail component by CPNOSE (...).   L  : Length of ogive section of nose, passed as LN to CPNOSE (...) LCB : Length of conical boattail section of nose, passed as LCB to CPNOSE (...)    R1N : Radius at base of ogive nose section          (same as radius at top of conical boattail section);          1st radius from reference line (nose tip) R2N : Radius at base of conical boattail section;          2nd radius from reference line (nose tip) Refer to p. 9 of TIR-33.

1

9. Procedure CmptMissleCP (...) and related procedures are complete.  missle02.dll also contains additional procedures that are used in other model rocket analysis projects.  As these other procedures are upgraded, a new missle02.pas and/or missle.inc will be uploaded.
10. The same thing applies to functions XROTATECW (...), etc., from download file mathaux.for
To avoid any potential Intellectual Property violations, the entirety of mathproc.for is unavailable for download, since portions were previously published.  Unpublished portions are available for download as tir33.for and mathaux.for (links above).