Python CP Calculator

[blackbrandt]

I wrote this CP Calculator for a class I'm in. Can anyone give it a look and see if you can spot any problems? The numbers I'm getting are coming out similar but not identical to OpenRocket.

With no transitions, I was within 3% of Openrocket. With 1 transition, I was 11% off. I'm curious if this is just a result of the equations I'm using being too simple, or if I have errors in my calculations.

#    _     _      _     _      _     _      _     _      _     _      _     _

[FONT=courier new]#   (c).-.(c)    (c).-.(c)    (c).-.(c)    (c).-.(c)    (c).-.(c)    (c).-.(c)[/FONT]
[FONT=courier new]#    / ._. \      / ._. \      / ._. \      / ._. \      / ._. \      / ._. \[/FONT]
[FONT=courier new]#  __\( Y )/__  __\( Y )/__  __\( Y )/__  __\( Y )/__  __\( Y )/__  __\( Y )/__[/FONT]
[FONT=courier new]# (_.-/'-'\-._)(_.-/'-'\-._)(_.-/'-'\-._)(_.-/'-'\-._)(_.-/'-'\-._)(_.-/'-'\-._)[/FONT]
[FONT=courier new]#    || O ||      || O ||      || O ||      || O ||      || O ||      || O ||[/FONT]
[FONT=courier new]#  _.' `-' '._  _.' `-' '._  _.' `-' '._  _.' `-' '._  _.' `-' '._  _.' `-' '._[/FONT]
[FONT=courier new]# (.-./`-'\.-.)(.-./`-'\.-.)(.-./`-'\.-.)(.-./`-'\.-.)(.-./`-'\.-.)(.-./`-'\.-.)[/FONT]
[FONT=courier new]#  `-'     `-'  `-'     `-'  `-'     `-'  `-'     `-'  `-'     `-'  `-'     `-'[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]#############################################################################[/FONT]
[FONT=courier new]#   Author:  Matt Fletcher                                                  #[/FONT]
[FONT=courier new]#      PID:  None                                                           #[/FONT]
[FONT=courier new]#    Class:  PH412, Spring, 2017                                            #[/FONT]
[FONT=courier new]#  Helpers:  None                                                           #[/FONT]
[FONT=courier new]#                                                                           #[/FONT]
[FONT=courier new]#  Program:  Barrowman Equations Calculator                                 #  [/FONT]
[FONT=courier new]# Due Date:  Feb 13,2017                                                    # [/FONT]
[FONT=courier new]#                                                                           #[/FONT]
[FONT=courier new]# Language:  Python 3.5.2                                                   #[/FONT]
[FONT=courier new]#      IDE:  IDLE                                                           #[/FONT]
[FONT=courier new]#                                                                           #[/FONT]
[FONT=courier new]# Purpose:   Given a user input of dimensions of a rocket, outputs the      #[/FONT]
[FONT=courier new]#            locations of the Center of Pressure relative to the tip of     #[/FONT]
[FONT=courier new]#            the nose cone.                                                 #[/FONT]
[FONT=courier new]#                                                                           #[/FONT]
[FONT=courier new]#Sources:    https://www2.estesrockets.com/pdf/TIR-33_Center_of_Pressure.pdf #[/FONT]
[FONT=courier new]#            https://www.rocketmime.com/rockets/Barrowman.html               #[/FONT]
[FONT=courier new]#   "Bugs":  None                                                           #[/FONT]
[FONT=courier new]#   "Undocumented features": Fin ASCII art not appearing correctly          #[/FONT]
[FONT=courier new]#############################################################################[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]#####################################################[/FONT]
[FONT=courier new]#                Begin Infinite Loop                #[/FONT]
[FONT=courier new]#####################################################[/FONT]
[FONT=courier new]execute_program = True[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]while execute_program:[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]    #####################################################[/FONT]
[FONT=courier new]    #                     Imports                        #[/FONT]
[FONT=courier new]    #####################################################[/FONT]
[FONT=courier new]    from math import sqrt[/FONT]



[FONT=courier new]    #####################################################[/FONT]
[FONT=courier new]    #                  User Inputs                        #[/FONT]
[FONT=courier new]    #####################################################[/FONT]

[FONT=courier new]    #Determining what units will be used[/FONT]

[FONT=courier new]    #Comparator to check if valid units have been selected.[/FONT]
[FONT=courier new]    valid_units = False[/FONT]
[FONT=courier new]    while valid_units == False:[/FONT]
[FONT=courier new]        print ("What units would you like to use?")[/FONT]
[FONT=courier new]        print ("1. Inches")[/FONT]
[FONT=courier new]        print ("2. Centimeters")[/FONT]
[FONT=courier new]        units=str(input("Enter your choice here>>> "))[/FONT]

[FONT=courier new]        #Converting choice to string (to be used later)[/FONT]
[FONT=courier new]        if units=="1":[/FONT]
[FONT=courier new]            units="inches"[/FONT]
[FONT=courier new]            valid_units=True[/FONT]
[FONT=courier new]        elif units=="2":[/FONT]
[FONT=courier new]            units="centimeters"[/FONT]
[FONT=courier new]            valid_units=True[/FONT]
[FONT=courier new]        #If user does not input a valid choice, run through the menu again. [/FONT]
[FONT=courier new]        else:[/FONT]
[FONT=courier new]            print("Invalid choice, please try again.")[/FONT]


[FONT=courier new]
[/FONT]
[FONT=courier new]    #Requesting user to input all dimensions in their chosen dimension.     [/FONT]
[FONT=courier new]    print ("Please enter all dimensions in "+units+".")[/FONT]

[FONT=courier new]    #################################[/FONT]
[FONT=courier new]    #Nose Cone Parameters[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("Nose Cone Calculations")[/FONT]

[FONT=courier new]    #Shape of the Nose Cone[/FONT]
[FONT=courier new]    print ("Nose cone shape options")[/FONT]
[FONT=courier new]    print ("1. Conical")[/FONT]
[FONT=courier new]    print ("2. Ogive")[/FONT]
[FONT=courier new]    print ("3. Parabolic")[/FONT]

[FONT=courier new]    nose_shape=str(input("Enter your choice here>>> "))[/FONT]

[FONT=courier new]    #Length of the nose cone[/FONT]
[FONT=courier new]    print ("How long is your nose cone in "+units+"?")[/FONT]
[FONT=courier new]    nose_length=float(input("Enter the length here>>> "))[/FONT]

[FONT=courier new]    #Diameter at base of the nose cone[/FONT]
[FONT=courier new]    print ("What is the diameter at the base of the nose cone in "+units+"?")[/FONT]
[FONT=courier new]    nose_base_diameter=float(input("Enter the diameter here>>> "))[/FONT]

[FONT=courier new]    ################################[/FONT]
[FONT=courier new]    #Transition Parameters[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("Transition Calculations")[/FONT]
[FONT=courier new]    print ("Enter the number of conical transitions on the rocket.")[/FONT]
[FONT=courier new]    print("This program currently only supports 1 transition")[/FONT]
[FONT=courier new]    print ("If there are no transitions, enter 0.")[/FONT]
[FONT=courier new]    transition_count=int(input("Enter the number here>>> "))[/FONT]

[FONT=courier new]    if transition_count==0:[/FONT]

[FONT=courier new]        print("There are no transitions in this rocket.")[/FONT]

[FONT=courier new]    if transition_count!=0:[/FONT]
[FONT=courier new]        #Get the length of the transition section[/FONT]
[FONT=courier new]        print("Enter the length of the transition")[/FONT]
[FONT=courier new]        transition_length=float(input("Enter here>>> "))[/FONT]

[FONT=courier new]        #Get the forward diameter of the transition[/FONT]
[FONT=courier new]        print("Enter the forward diameter of the transition")[/FONT]
[FONT=courier new]        transition_forward_diameter=float(input("Enter here>>> "))[/FONT]

[FONT=courier new]        #Get the aft diameter of the transition[/FONT]
[FONT=courier new]        print("Enter the aft diameter of the transition")[/FONT]
[FONT=courier new]        transition_aft_diameter=float(input("Enter here>>> "))[/FONT]

[FONT=courier new]        #Get the distance from the nose tip to the front of the transition[/FONT]
[FONT=courier new]        print("Enter the distance from the nose tip to the front of the transition")[/FONT]
[FONT=courier new]        transition_xp=float(input("Enter here>>> "))[/FONT]
[FONT=courier new]    ################################[/FONT]
[FONT=courier new]    #Fin Parameters[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("")[/FONT]
[FONT=courier new]    print("Fin Calculations")[/FONT]
[FONT=courier new]    print("Please note, all fins must be trapezoidal.")[/FONT]
[FONT=courier new]    print("")[/FONT]

[FONT=courier new]    print("Fin Dimensions are shown in the following diagram")[/FONT]

[FONT=courier new]    #Print fin diagram[/FONT]
[FONT=courier new]    print( '''[/FONT]

[FONT=courier new]                        ____[/FONT]
[FONT=courier new]    |*                    Leading Edge[/FONT]
[FONT=courier new]    |  *L                    to[/FONT]
[FONT=courier new]R     |    *E                    Tip[/FONT]
[FONT=courier new]O    |      *A                Parallel[/FONT]
[FONT=courier new]O    |        *D                To[/FONT]
[FONT=courier new]T    |          *                Body[/FONT]
[FONT=courier new]    |            *E                Tube[/FONT]
[FONT=courier new]E     |* M            *D            [/FONT]
[FONT=courier new]D    |  C * I          *G            [/FONT]
[FONT=courier new]G    |     H  * D        *E        ______[/FONT]
[FONT=courier new]E    |        O   *       |[/FONT]
[FONT=courier new]    |           R    *   |[/FONT]
[FONT=courier new]    |           D     |T C[/FONT]
[FONT=courier new]    |                 |I H[/FONT]
[FONT=courier new]    *             |P O[/FONT]
[FONT=courier new]       T*                 |  R[/FONT]
[FONT=courier new]    E    R*           |  D[/FONT]
[FONT=courier new]            D    A*         |[/FONT]
[FONT=courier new]              G   I* |[/FONT]
[FONT=courier new]              E L*[/FONT]

[FONT=courier new]    |_______SemiSpan_____|[/FONT]
[FONT=courier new]    ''')[/FONT]
[FONT=courier new]    #TODO: Insert fin diagram. [/FONT]

[FONT=courier new]    #Get a valid number of fins (3,4, or 6)[/FONT]
[FONT=courier new]    valid_fin_count=False[/FONT]

[FONT=courier new]    #Cue user to input number of fins. [/FONT]
[FONT=courier new]    while valid_fin_count==False:[/FONT]
[FONT=courier new]        fin_count=int(input("Enter the number of fins >>> "))[/FONT]

[FONT=courier new]        #If a bad number of fins is entered, ask again. [/FONT]
[FONT=courier new]        if fin_count!=3 or fin_count!=4 or fin_count!=6:[/FONT]
[FONT=courier new]            print("The Barrowman equations only support rockets with 3, 4, or 6 fins.")[/FONT]
[FONT=courier new]        #If a correct number of fins is entered, exit the loop and continue. [/FONT]
[FONT=courier new]        if fin_count==3 or fin_count==4 or fin_count==6:            [/FONT]
[FONT=courier new]            valid_fin_count=True[/FONT]


[FONT=courier new]    #Cue user to input the fin root chord. [/FONT]
[FONT=courier new]    print("Enter the length of the fin root chord.")[/FONT]
[FONT=courier new]    fin_root=float(input("Enter here>>> "))[/FONT]


[FONT=courier new]    #Cue user to input the fin tip chord. [/FONT]
[FONT=courier new]    print("Enter the length of the fin tip chord.")[/FONT]
[FONT=courier new]    fin_tip=float(input("Enter here>>> "))[/FONT]


[FONT=courier new]    #Cue user to input the fin semi-span. [/FONT]
[FONT=courier new]    print("Enter the length of the fin semi-span.")[/FONT]
[FONT=courier new]    fin_semispan=float(input("Enter here>>> "))[/FONT]


[FONT=courier new]    #Cue user to input the fin mid chord. [/FONT]
[FONT=courier new]    print("Enter the length of the fin mid chord.")[/FONT]
[FONT=courier new]    fin_midchord=float(input("Enter here>>> "))[/FONT]


[FONT=courier new]    #Cue user to input the distance between fin root leading edge  [/FONT]
[FONT=courier new]    #and fin tip leading edge parallel to body[/FONT]
[FONT=courier new]    print("Input the distance between fin root leading edge and fin tip leading edge parallel to body")[/FONT]
[FONT=courier new]    fin_x_r=float(input("Enter here>>> "))[/FONT]

[FONT=courier new]    #Cue user to input the distance between the tip of the nose cone  [/FONT]
[FONT=courier new]    #and the fin root chord leading edge. [/FONT]
[FONT=courier new]    print("Input the distance distance between the tip of the nose cone and the fin root chord leading edge.")[/FONT]
[FONT=courier new]    fin_x_b=float(input("Enter here>>> "))[/FONT]

[FONT=courier new]    #Cue user to input the radius of the body tube at the end of the rocket. [/FONT]
[FONT=courier new]    print("Enter the radius of the body tube at the aft end.")[/FONT]
[FONT=courier new]    tube_end_radius=float(input("Enter here>>> "))[/FONT]



[FONT=courier new]    #####################################################[/FONT]
[FONT=courier new]    #                  Calculations                        #[/FONT]
[FONT=courier new]    #####################################################[/FONT]

[FONT=courier new]    #Nose Cone CP calculation[/FONT]
[FONT=courier new]    #Coefficient of nose is always 2[/FONT]
[FONT=courier new]    c_nose=2.[/FONT]

[FONT=courier new]    #Calculating distance of CP of nose from tip of nose depending on chosen shape of nose cone. [/FONT]
[FONT=courier new]    #If nose is conical[/FONT]
[FONT=courier new]    if nose_shape=="1":[/FONT]
[FONT=courier new]        x_nose=(2./3)*nose_length[/FONT]
[FONT=courier new]    #If nose is ogive[/FONT]
[FONT=courier new]    elif nose_shape=="2":[/FONT]
[FONT=courier new]        x_nose=0.466*nose_length[/FONT]
[FONT=courier new]    #If nose is parabolic.[/FONT]
[FONT=courier new]    elif nose_shape=="3":[/FONT]
[FONT=courier new]        x_nose=0.5*nose_length[/FONT]

[FONT=courier new]    #Transition CP Calculation[/FONT]
[FONT=courier new]    #This section will only run if there is a transition section in the rocket. [/FONT]
[FONT=courier new]    if transition_count!=0:[/FONT]
[FONT=courier new]        #Transition Coefficient[/FONT]
[FONT=courier new]        c_transition=2*((transition_forward_diameter/nose_base_diameter)**2-(transition_aft_diameter/nose_base_diameter)**2)[/FONT]

[FONT=courier new]        #Transition distance from tip of nose cone[/FONT]
[FONT=courier new]        x_transition=transition_xp+transition_length/3.+(1+(1-(transition_forward_diameter/transition_aft_diameter))/(1-(transition_forward_diameter/transition_aft_diameter)**2))[/FONT]

[FONT=courier new]    #If no transitions, set c and x values of transitions equal to zero. [/FONT]
[FONT=courier new]    else:[/FONT]
[FONT=courier new]        c_transition=0[/FONT]
[FONT=courier new]        x_transition=0[/FONT]
[FONT=courier new]
[/FONT]



[FONT=courier new]    #Fin CP Calculation[/FONT]

[FONT=courier new]    #Fin coefficient[/FONT]
[FONT=courier new]    c_fin=(1+(tube_end_radius)*(tube_end_radius+fin_semispan))*((4*fin_count)*(fin_semispan/nose_base_diameter)**2)/(1+sqrt(1+((2*fin_midchord)/(fin_root+fin_tip))**2))[/FONT]

[FONT=courier new]    #Fin distance from nose[/FONT]
[FONT=courier new]    x_fin=fin_x_b+(fin_x_r)/3.*((fin_root+2*fin_tip)/(fin_root+fin_tip))+1/6.*((fin_root+fin_tip)-((fin_root*fin_tip)/(fin_root+fin_tip)))[/FONT]

[FONT=courier new]    #Sum of coefficients[/FONT]
[FONT=courier new]    c_rocket=c_nose+c_transition+c_fin[/FONT]


[FONT=courier new]    #Finding weighted averages[/FONT]
[FONT=courier new]    x_rocket=(c_nose*x_nose+c_transition*x_transition+c_fin*x_fin)/(c_rocket)[/FONT]

[FONT=courier new]    print("The CP is "+str(round(x_rocket,2))+" "+units+" from the nose cone tip.")[/FONT]

[FONT=courier new]    #####################################################[/FONT]
[FONT=courier new]    #                  Loop Check                        #[/FONT]
[FONT=courier new]    #####################################################[/FONT]

[FONT=courier new]    #Ask user if they would like to try another run.[/FONT]

[FONT=courier new]    print("Would you like to simulate another rocket?.")[/FONT]
[FONT=courier new]    print("1. Yes")[/FONT]
[FONT=courier new]    print("2. No")[/FONT]
[FONT=courier new]    if str(input(">>>>")) == "2":[/FONT]
[FONT=courier new]        execute_program = False[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]
[/FONT]

 

[tHoagland]

What are all of lines preceeded by "#"? Comments? Never seen anyone actually use those...

I know that Barrowman equations are common knowledge but I've never looked at the details myself. Can you share the equations/link you're using for reference?

Look at the line for x_transition. It looks like you might be multiplying instead of dividing...

... aft_diameter)) / ( 1-(transition_forward.....

Unless I missed something

 

[blackbrandt]

Yup. Hashtags are comments. I comment my code to death.

Link to the equations is in the header of the code.
It's also here: VV
https://www.rocketmime.com/rockets/Barrowman.html

 

[tHoagland]

I missed the link in the header; but that is the page I looked at too. Did you look at the line where you define x_transition ( can't get a line number from my phone).

 

[dhbarr]

This may be helpful if you're using straight Barrowman:. https://argoshpr.ch/joomla1/articles/pdf/sentinel39-galejs.pdf

 

[Buckeye]

What are all of lines preceeded by "#"? Comments? Never seen anyone actually use those...

Huh? Never seen comments in code?

Looking good, BB.

 

[blackbrandt]

Thanks! I'm starting to assume that it's not a calculation error, but just due to the simplicity of the program.

Also, if anyone's ever bored and needs some reading, the OpenRocket technical documentation is a super fun read: https://openrocket.sourceforge.net/techdoc.pdf

 

[tHoagland]

Except for the error calculating X_transition as pointed out above....

 

[blackbrandt]

Found that, fixed it, and the numbers are still off a bit. Openrocket says 17.8 inches for a rocket with 1 transition, while the python script says 20.4.... ??

 

[Nytrunner]

I'd blame the teddy bear cabal at the top. They're clearly up to no good.

Which source did you use for the equations?

 

[blackbrandt]

The sources are listed in the header.

 

[markkoelsch]

Best guess is that OR has been tweaked beyond Barrowman.

 

[manixFan]

Have you compared it to RockSim values? I've never compared OR to RS but it would give you another validation source.

You might try checking out:

https://apogeerockets.com/education/downloads/Newsletter238.pdf

which is an article on how Rocksim varies from the Barrowman equations. He also links to additional articles.

and https://documentslide.com/documents/60921375-the-descriptive-geometry-of-nose-cone.html

Which is a pretty cool presentation about nosecone shapes. On page 6 the author has different formulas for calculating the CP of nosecones depending on the shape. So it could be that OR is calculating the transition as a modified nosecone using a different formula than the general Barrowman one you are using.

finally https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20010047838.pdf

Is a link to the actual Masters Thesis by Barrowman. He states the equations assume a conical (he calls it a 'sharp point') nosecone. So if OR is using the actual shape of the nosecone to calc CP and you have specified an ogive or something else, that could be it. If so, change the nose cone shape to conical and see if your numbers get closer. But my guess is that OR is also treating the transition as a special case.

But that's all just a best guess without reading the OR tech docs you linked to. Gotta work for a living as much fun as it would be to read through the docs!

Good luck,


Tony