BUG: Cld divided by zero

[MateusStano]

Pull request type

• Code base additions (bugfix, features)
• Code maintenance (refactoring, formatting, renaming, tests)
• ReadMe, Docs and GitHub maintenance
• Other (please describe):

Pull request checklist

• Code base additions (for bug fixes / features):

  • Tests for the changes have been added
  • Docs have been reviewed and added / updated if needed
  • Lint (black rocketpy) has passed locally and any fixes were made
  • All tests (pytest --runslow) have passed locally

What is the current behavior?

Currently, RocketPy is only calculating the difference between the roll forcing and roll damping coefficients and passing that as the roll moment. Obviously, this led to some issues with the code and inaccuracies with the results

What is the new behavior?

Now RocketPy calculates the roll moment accordingly. All equations and considerations used are documented in the Roll_Equations.pdf file. The calculations for the fins lift coefficient also changed. A more general formula is now used since it allowed for easy access to some of the coefficients needed.

Does this introduce a breaking change?

• Yes
• No

Other information

Unfortunately, I didn't find any data on angular velocity of any subsonic rockets to validate the results from the moment equations. I ended up resorting to simulating the same rocket on RocketPy and on Open Rocket and comparing the results... and they weren't very promising.
In comparison to open rocket, rocketpy's calculated angular velocity were much higher. However, in Open Rocket's documentation they compare their results to experimental ones:

It seems Open Rocket understates the Roll rate. I think it's coherent that RocketPy's results for the roll rate are bigger.

Still, it would be useful to compare RocketPy's roll calculations results to an experimental launch.

[giovaniceotto]

@MateusStano, just to make sure we are in the same page, this problem is caused by the following approximation, correct?

[giovaniceotto]

To be honest, I am not convinced that the solution given is the best one possible. The freestream speed will most likely be zero only when the roll rate omega3 is also zero. In such a case, there is no moment acting on the rocket. Therefore, I believe the original implementation must be used, but we should check if the roll rate omega3 is different from 0 in the first place before performing the calculations.

Furthermore, if we detect that the angle of attack becomes too high (high roll rate and low freestream speed), perhaps we should raise a warning since the equations being used loose validity.

[MateusStano]

this problem is caused by the following approximation, correct?

Yes, this is causing the problem

Furthermore, if we detect that the angle of attack becomes too high (high roll rate and low freestream speed), perhaps we should raise a warning since the equations being used loose validity.

Think everything you said here makes sense. But one question, when should I consider that the angle of attack is too high? Any way you recommend checking that?

[Gui-FernandesBR]

this problem is caused by the following approximation, correct?

Yes, this is causing the problem

Furthermore, if we detect that the angle of attack becomes too high (high roll rate and low freestream speed), perhaps we should raise a warning since the equations being used loose validity.

Think everything you said here makes sense. But one question, when should I consider that the angle of attack is too high? Any way you recommend checking that?

What if, at the end of the simulation, you check if there's a value higher than Xº on the angle of attack vector and if so you print a warning? That X value could be, let's say, 10 degress? With this solution you would perfom only 1 extra if in the simulation,

[giovaniceotto]

I partially agree with @Gui-FernandesBR. I believe implementing a couple of lines of code in the Flight.postProcess method could be useful to calculate the highest angle of attack along the flight and print this detail in the Maximum Values section.

RocketPy does not usually warn users in case something is bad in the simulation... but according to @Projeto-Jupiter/back-end's OKRs, we should do that from now on. I would raise a warning only if the angle of attack becomes greater than 20 degrees, but this is an arbitrary value, for sure.

Another interesting metric would be for how long the angle of attack remains greater than a specified value, be it 10 or 20 degrees.

Anyways, I trust you to establish a reasonable threshold for such warning @MateusStano.

[MateusStano]

Hey @giovaniceotto @Gui-FernandesBR! Changed quite a bit of the roll implementation and now everything seems to be in order. I changed the description of the PR and added some stuff to it. Just a couple of things:

1. I'm having problems with black automatic verification. Seems to be something to do with access permission and I have no idea how to solve this

2. This PR will definitely cause a conflict with Bug: Issue 115 static margin stability #125. I think this one should be merged first and than some alteration will have to be made on Bug: Issue 115 static margin stability #125

[giovaniceotto]

Explendid work @MateusStano. I believe the new code is, overall, significantly more reliable and trustworthy. We can finally say that RocketPy has a great roll model. Furthermore, I really like the fact that you are starting to add compressibility corrections by introducing Mach number as an argument to some functions.

Here are my main comments:

• I know that this isn't new here, but I am concerned that our code is relying quite heavily on the Function.differentiate method. As we know, its implementation is not ready for production and should be enhanced. Something to add to our development roadmap @Gui-FernandesBR.
• I know we have talked about this, by I really like the new names used for some variables, such as renaming cldata to cl and using clalphaSingleFin.
• Any reason λ and tau are being used in lines 594 - 614, instead of lambda and tau or λ and τ? I am afraid that using special chars such as λ and τ can make the code a little bit harder to maintain, but it does improve readability quite significantly!
• I am worried that using fin = [(0, 0, cpz), cl, rollParameters, "Fins"] together with nose = [(0, 0, cpz), cl, "Nose Cone"] can be a bad choice. For fins, the third element of the list is roll parameters, while for the nose cone, the third element is the name. Furthermore, Rocket.addTail, Rocket.addNose, Rocket.addFins all say that ...storing its parameters as part of the aerodynamicSurfaces list. Its parameters are the axial position along the rocket and its derivative of the coefficient of lift in respect to angle of attack., which do not include the roll parameters. Perhaps we should make everything follow the same convetion, such as [(cpx, cpy, cpz), cl, clfDelta*cantAngleRad, cldOmega, "Surface Name"] or something similar? I am not against moving away from lists and starting to store this info in a custom class as well...

Here are my requests:

• Could the description of cl in lines 421, 479, and 565 of Rocket.py be improved? Saying it contains the lift data is quite vague. Perhaps we can specify the inputs and outputs of the function to make it clearer.
• Could you add a comment to line 585 of Rocket.py (aspect ratio AR definition), describing that this is Barrowmans convetion? The 2 multiplying is quite weird to understand otherwise.
• Could functions beta and finNumCorrection have a short docstring instead of comments on top of them? I believe this is more in line with our formatting practices. Furthermore, function finNumCorrection could benefit from a reference to where we got the numbers [2.37, 2.74, 2.99, 3.24].
• Would removing Aref and d from lines 1259 and 1260 of Flight.py be problematic? I believe the convention the remaining of the Flight.uDot method is to use self.rocket.area and 2 * self.rocket.radius explicitly. Checking the code formatter by black, I believe spacing is not a problem. Therefore, using their explicit form can make the code easier to understand.

Final remarks:

• I have not independently verified the equations to make sure there are no typos. I am not the best person to do so considering I was involved when you coded them.

Make sure you re-request my review when you are done!

[Lucas-KB]

Explendid work @MateusStano. I believe the new code is, overall, significantly more reliable and trustworthy. We can finally say that RocketPy has a great roll model. Furthermore, I really like the fact that you are starting to add compressibility corrections by introducing Mach number as an argument to some functions.

Here are my main comments:

• I know that this isn't new here, but I am concerned that our code is relying quite heavily on the Function.differentiate method. As we know, its implementation is not ready for production and should be enhanced. Something to add to our development roadmap @Gui-FernandesBR.
• I know we have talked about this, by I really like the new names used for some variables, such as renaming cldata to cl and using clalphaSingleFin.
• Any reason λ and tau are being used in lines 594 - 614, instead of lambda and tau or λ and τ? I am afraid that using special chars such as λ and τ can make the code a little bit harder to maintain, but it does improve readability quite significantly!
• I am worried that using fin = [(0, 0, cpz), cl, rollParameters, "Fins"] together with nose = [(0, 0, cpz), cl, "Nose Cone"] can be a bad choice. For fins, the third element of the list is roll parameters, while for the nose cone, the third element is the name. Furthermore, Rocket.addTail, Rocket.addNose, Rocket.addFins all say that ...storing its parameters as part of the aerodynamicSurfaces list. Its parameters are the axial position along the rocket and its derivative of the coefficient of lift in respect to angle of attack., which do not include the roll parameters. Perhaps we should make everything follow the same convetion, such as [(cpx, cpy, cpz), cl, clfDelta*cantAngleRad, cldOmega, "Surface Name"] or something similar? I am not against moving away from lists and starting to store this info in a custom class as well...

Here are my requests:

• Could the description of cl in lines 421, 479, and 565 of Rocket.py be improved? Saying it contains the lift data is quite vague. Perhaps we can specify the inputs and outputs of the function to make it clearer.
• Could you add a comment to line 585 of Rocket.py (aspect ratio AR definition), describing that this is Barrowmans convetion? The 2 multiplying is quite weird to understand otherwise.
• Could functions beta and finNumCorrection have a short docstring instead of comments on top of them? I believe this is more in line with our formatting practices. Furthermore, function finNumCorrection could benefit from a reference to where we got the numbers [2.37, 2.74, 2.99, 3.24].
• Would removing Aref and d from lines 1259 and 1260 of Flight.py be problematic? I believe the convention the remaining of the Flight.uDot method is to use self.rocket.area and 2 * self.rocket.radius explicitly. Checking the code formatter by black, I believe spacing is not a problem. Therefore, using their explicit form can make the code easier to understand.

Final remarks:

• I have not independently verified the equations to make sure there are no typos. I am not the best person to do so considering I was involved when you coded them.

Make sure you re-request my review when you are done!

@giovaniceotto about the list and positioning thing you mentioned, we could, of course implement a class to handle aerodynamic surfaces, but a simple alternative woud be to use a dictionary, at least it is more comprehensive and a lot less confusing than having to remember the indices of each element. And also for being easier to maintain, if two branches want to add something new to the aerodynamic surfaces, it becomes a great mess really fast.

[giovaniceotto]

@giovaniceotto about the list and positioning thing you mentioned, we could, of course implement a class to handle aerodynamic surfaces, but a simple alternative would be to use a dictionary, at least it is more comprehensive and a lot less confusing than having to remember the indices of each element. And also for being easier to maintain, if two branches want to add something new to the aerodynamic surfaces, it becomes a great mess really fast.

Agreed. Another option are named tuples. Here is an example from its documentation:

>>> # Basic example
>>> Point = namedtuple('Point', ['x', 'y'])
>>> p = Point(11, y=22)     # instantiate with positional or keyword arguments
>>> p[0] + p[1]             # indexable like the plain tuple (11, 22)
33
>>> x, y = p                # unpack like a regular tuple
>>> x, y
(11, 22)
>>> p.x + p.y               # fields also accessible by name
33
>>> p                       # readable __repr__ with a name=value style
Point(x=11, y=22)

[Gui-FernandesBR]

Let's go guys, we have a lot to-do here. 1st, many thanks for implementing this, @MateusStano !! Now my comments

1. Stano, you mentioned having issues with black, right? Still having those? If so, could you contact to Giovani directly to solve that? (@giovaniceotto @PatrickSampaioUSP please make sure to include some minutes on our next workshop to explain about black formatting again... Not the first time I see someone facing this)
2. @Lucas-KB I liked the idea. Do you think would be easy to implement it here and right now? If it doesn't request too many time, can you commit that to help Stano?
3. I myself can check the equations. I'm not an expert so if I can't understand the code, it is already too complex for new users.
4. Stano, about the comparison RocketPy vs Openrocket, did you save the numbers? Just to clarify how bigger RocketPy values are
5. Regarding automatic tests, I think we could implement some of them here... Any ideas? I will give coulpe of them. @MateusStano are you able to create them or need to delegate to someone else?
   5.A - Dimensional analysis to verify CL is dimensioless.
   5.B - Verify if simulation with cantangle=0 is giving same results as simulations previous performed (you need to check but probably it is already being done).
6. And well, @MateusStano there are sugestions from giovani that I think are easier to solve so I didn't comment, let me know if you have problem wit them.

Stay hungry guys, we are close to solve this one

[Lucas-KB]

2. @Lucas-KB I liked the idea. Do you think would be easy to implement it here and right now? If it doesn't request too many time, can you commit that to help Stano?

I believe it would. I will try it out.

[Lucas-KB]

Actually I don't know if I can contribute because the branch is on @MateusStano's personal github. I can make the changes, but I am not sure how to integrate them here.

[Gui-FernandesBR]

Actually I don't know if I can contribute because the branch is on @MateusStano's personal github. I can make the changes, but I am not sure how to integrate them here.

did you try this? I think it will not give us conflicts

[Lucas-KB]

Ok, I think it is done @Gui-FernandesBR.

[giovaniceotto]

Here is what is left:

• Flight test with a non-null cant angle. @MateusStano
• Equations review. @Gui-FernandesBR

I'll approve this PR as soon as these are done.

[MateusStano]

are the units of this roll rate really in rad/s?

Yes, rad/s.

furthermore, can you describe here why you believe OpenRocket's code has some inaccuracies?

Okay, I think it's for a few reasons but here is a list of all the differences I could find:

1. OpenRocket doesn't apply the interference factors for roll damping and forcing (equations 26 and 27 on the Roll_Equations.pdf).
2. When calculating the roll damping coefficient they use the Cnalpha0 instead of Cnalpha1. Basically, they used equation 19 instead of 20 from the PDF. I believe this can lead to incorrect results since the values of Cnalpha0 can vary significantly in comparison to those in Cnalpha1. Also, physically speaking Cnalpha0 only accounts for the lift generated in the 2D section of the fin, it is more sensible to take into account the value that represents the three-dimensional fin (Cnalpha1) when calculating the roll. Furthermore, they always use 2π/beta for the value of Cnalpha0, which doesn't take into account airfoils (we currently don't do that either but it's worth mentioning).
3. Stall is considered in the roll forcing coefficient. RocketPy doesn't simulate fin stall (yet) I think, so there is no way to implement that in our code. With that said, Calisto's simulation that I used for comparison shouldn't get into stall areas until the apogee so this shouldn't be causing any differences.
4. They consider that when the angular velocity is smaller than 1 the damping coefficient value is zero. This is probably to avoid the very same singularity that we were correcting in this Pull Request.
5. They calculate the integral (in equations 19 and 20 of the PDF) numerically. This is definitely to take into account free form fins that OpenRocket can simulate. However, this integration might be generating some numerical disparities. But it would require some further investigating to be sure of that.

Items 1 and 2 are the most relevant for the result differences I believe. Removing the interference factors and substituting Cnalpha1 for 2π/beta in the roll damping coefficient equation in RocketPy wields the following results:

These are very similar to OpenRocket's. Also, the little bump you can see right before the 5 seconds mark is coherent with the way we defined beta.

[Gui-FernandesBR]

Okay, it is getting intersting guys... awesome!

I'll review equations comparing with the barrowman theory this weekend.

Meanwhile, I think it is important to estabilish a new test just like Giovani mentioned "What would be needed is a new test with a rocket with cantAngle != 0. There currently exist no tests for this. A particularly important one would be when wind is 0, which was the error which helped us find this bug in the first case."

I mean, it's important we ensure the error doesn't be repeated on future PR, right?
How do you think we could approach it, @MateusStano ?

[Gui-FernandesBR]

Let me correct... Actually I think you already created the test file.
Let's move then.

Correct me if I'm mistaken pls

[MateusStano]

Let me correct... Actually I think you already created the test file.

Yep. It is already done :)

I'll review equations comparing with the barrowman theory this weekend.

I recommend checking Roll Equations.pdf and seeing all that is there is correct

[Gui-FernandesBR]

Question: when you make assert test_flight.allInfo() == None on the "tests/test_flight.py" file, what you want is to verify if the Flight class runned well. In case of any trouble the .allInfo() will return some warning.

Is that what is happening?

[Gui-FernandesBR]

I could not find this definition on Barrowman's equations, do you remember where did you find?

Not saying it's wrong, but I didn't understand

[MateusStano]

The definition of Ymac? I think the 'mac' subtext is from OpenRocket's documentation. I'll change it to Yma which is the form used on Barrowman to make the code more consistent.
But here, it is on equation 3-26 on Barrowman:

[Gui-FernandesBR]

This pull request seems aumost done!! Really awesome @MateusStano !!!

I have some questions but, as I said I'm not an expert on aerodynamics formulation, therefore you might face them as simple question from curiosity, and never think I'm say you're wrong.

[Gui-FernandesBR]

Awesome work by Stano!!!

[Gui-FernandesBR]

@giovaniceotto I believe we are ready to merge already. However, as you requested changes we are asking re-approval here.

Could you see this by the weekend? We are holding another bug on merge queue

[giovaniceotto]

Great work @MateusStano and @Gui-FernandesBR.

In my last review, I noticed what could be a bug. @MateusStano, what would happen in case we create a rocket with fins with cantAngle=0 and used the initialSolution argument to launch it with angular velocity? Would a damping moment be calculated and applied?

[MateusStano]

In my last review, I noticed what could be a bug. @MateusStano, what would happen in case we create a rocket with fins with cantAngle=0 and used the initialSolution argument to launch it with angular velocity? Would a damping moment be calculated and applied?

It wouldn't..... but now it works!!
I just removed a cantAngle != 0 check that isn't needed anymore

[Gui-FernandesBR]

Ready to merge @giovaniceotto ?

[giovaniceotto]

Yes! The only checks failing have to do with formatting.

@MateusStano, will you do the honor of merging?