Rodhern's update 3. #5

Rodhern · 2018-11-10T22:03:58Z

A possible merge as per #4.
I do not use KSP 1.5.1 myself, so I can't even know if this will compile.
The changes are:

Revised interpretations of some stability derivatives, such as the definitions of the 'right' and 'up' components of 'velocity', and inverted signs for Zq, Ze, Xq and Xe to better match the sign convention used for Zu and Xu.
Tweaked stability derivatives solver for level flight AoA.
Reset editor part cache at editor reset to fix the most frequent cause of the stability derivatives changing on craft reload (Values in the stability derivatives panel can change wildly when a craft is reloaded in the editor ferram4/Ferram-Aerospace-Research#177).
Minor updates to localization files.
Cosmetic changes to the stability derivatives GUI window.
Added stability derivatives export button(s) (screenshots available).

Change default AoA search method to one better suited for slow flight scenarios.

…atives' section (minor graphical tweak to GUI). Rename 'siteffG' (situation effective gravity) to 'siteffg' (named 'env.effg' when exported to Octave text file).

dkavolis · 2018-11-11T17:49:05Z

.gitignore

@@ -4,6 +4,7 @@ obj
 *.mdb
 /GameData/FerramAerospaceResearch/Custom*.cfg
 /GameData/FerramAerospaceResearch/FARForceDataUpdate.cfg
+/note.txt


No need to include this in PR

Just need a few minutes to figure out how to do this :-)

I just pushed the 'removed' note.txt to the source branch. I hope this is the Github way that I am supposed to solve the comment.

dkavolis · 2018-11-11T18:48:46Z

FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/InstantConditionSim.cs

-            Vector3d sideways = Vector3.Cross(velocity, liftVector).normalized;
+            Vector3d sideways = -forward * cosAlpha * sinBeta;
+            sideways += right * (cosPhi * cosBeta - sinPhi * sinAlpha * sinBeta); // Rodhern: The same as Vector3.Cross(velocity, liftVector)
+            sideways += up * (cosPhi * sinAlpha * sinBeta + sinPhi * cosBeta);    //  but possibly with better accuracy.


It's better to leave it as it was because Vector3.cross is more concise and better for readability unless the accuracy is much worse

Sure I can change that. Do you feel it is better for readability even when the remark is in the code? I feel there is sort of a nice analogous feel between 'velocity', 'liftVector' and 'sideways' this way. It is up to you though.

We're computing principal axes so cross is easier and shows that explicitly, no one will ever bother computing the last component from scratch (at least on paper) when it is also perpendicular to the other two

dkavolis · 2018-11-11T18:53:31Z

FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/InstantConditionSim.cs

-
-            velocity.Normalize();
+            velocity += right * (sinPhi * sinAlpha * cosBeta + cosPhi * sinBeta); // Rodhern: Changed velocity 'right' and 'up' components
+            velocity += -up * (cosPhi * sinAlpha * cosBeta - sinPhi * sinBeta);   //  to align vector to naive forward direction.


Which order are the rotations applied in? Need to be sure the math is correct before merging

The velocity magnitude is not equal to 1. Perhaps it should be

Vector3d velocity = forward * cosAlpha * cosBeta; velocity += right * (sinPhi * sinAlpha * cosBeta - cosBeta * sinBeta); velocity += -up * (cosPhi * sinAlpha * cosBeta + sinPhi * sinBeta);

http://www.kostasalexis.com/frame-rotations-and-representations.html

[edit: ad velocity magnitude] Are you sure you used the formulas from my code? Did you change the components back to the code from before the merge on purpose?

I wrote a python script to visualize the rotations and yours looks correct. The old implementation not only did not generate unit vectors, they were also not perpendicular and induced additional rotations. I guess this bug went unfixed for so long because large rotations were never used. The simulation should be more accurate with this fix.

dkavolis · 2018-11-11T18:57:06Z

FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/StabilityDerivCalculator.cs

-            pertOutput.Cl *= q * area / mass;
-            pertOutput.Cd *= q * area / mass;
+            pertOutput.Cl *= -q * area / mass; // Rodhern: Replaced 'q' by '-q', so that formulas
+            pertOutput.Cd *= -q * area / mass; //  for Ze and Xe match those for Zu and Xu.


I'd be careful about changing signs, I'm sure ferram had a good reason to set them the way they were. If anything, the negative signs shouldn't be there, the equations follow https://courses.cit.cornell.edu/mae5070/DynamicEquations.pdf (for the most part)

You can see my reason for changing signs here: Kerbal Space Program Forum. It would be great if you can point out where it is that my reasoning fails.

You're probably right, at the very least the signs should be consistent for the same axes. I'm not sure about the need for

Ferram-Aerospace-Research/FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/StabilityDerivCalculator.cs

Lines 235 to 236 in d9665fb

pertOutput.Cl += nominalOutput.Cd;

pertOutput.Cd -= nominalOutput.Cl;

according to the equations.

I have not really paid attention to that. Let us run through some quick back of the envelope calculations and verify that it makes sense to us. Let us look at Zw (i.e. dZ/dw).

The initial state is the aircraft flying horizontally at speed u0. In fact moving at speed u = u0 (forward), w = 0 (downward) and v = 0 (rightward). The angle of attack is some alpha.

Now imagine that we manipulate the direction of travel by rotating it ever so slightly towards the ground. We are in effect increasing alpha, and this will cause two changes that we will monitor. It changes the velocity vector and it changes the downward force. Say the rotation amount is tau radians. As an example of a specific near zero numerical value of tau I will use the symbol 2°.

The velocity deltas are,
dw = sin(tau) * u0 - 0 ≈ 2° * u0,
du = cos(tau) * u0 - u0 ≈ 0,
dv = 0.

The downward force delta (or in our case actually the acceleration delta) is, for some relevant constant k,
dZ = -cos(tau) * k * u0² * CL(AoA= alpha + tau) // the new lift component
+ -sin(tau) * k * u0² * CD(AoA= alpha + tau) // the new drag component
- -cos( 0) * k * u0² * CL(AoA= alpha + 0) // the old downward force
≈ tau * d( -cos(tau) * k * u0² * CL(tau) + -sin(tau) * k * u0² * CD(tau) )
= tau * k * u0² * ((-cos * dCL) + (sin * CL) + (-sin * dCD) + (-cos * CD))
≈ 2° * k * u0² * (-dCL + -CD)

You probably guessed by now that CL(tau) is the coefficient of lift in the given circumstances at an angle of attack equal to alpha plus tau, that CD(tau) is the corresponding coefficient of drag, and that dCL and dCD are the (mathematical) derivatives of CL and CD with respect to the angle of attack parameter (in radians).

Now if q = 1/2 * density * u0² and k = 1/2 * density * area / mass we get
Zw = dZ/dw ≈ 2° * k * u0² * (-dCL + -CD) / (2° * u0)
= -k * u0² * (dCL + CD) / u0
= -(1/2 * density * area / mass) * u0² * (dCL + CD) / u0
= -q * area * (dCL + CD) / (mass * u0)

Now compare to the code:
dCL: pertOutput.Cl = (pertOutput.Cl - nominalOutput.Cl) / (2 * FARMathUtil.deg2rad);
(dCL + CD): pertOutput.Cl += nominalOutput.Cd;
-q * area * (dCL + CD) / (mass * u0): pertOutput.Cl *= -q * area / (mass * u0);

I think it makes sense. Please point out mistakes that I make so that we do not end up confusing other comment readers.

dkavolis · 2018-11-11T19:01:32Z

FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/InstantConditionSim.cs

@@ -158,7 +158,7 @@ public void GetClCdCmSteady(InstantConditionSimInput input, out InstantCondition
                if (!(w && w.part))
                    continue;

-                w.ComputeForceEditor(velocity.normalized, input.machNumber, 2);
+                w.ComputeForceEditor(velocity, input.machNumber, 2);


velocity is no longer normalized here, this will affect force computations

This question is not meant in a harsh way; I am just a bit clumsy writing it succinct.
Can you produce an example where velocity.normalized is (strictly) closer to normalized than the current ('not normalized') velocity vector is?
If so that may also be a good reason to use Vector3.cross for sideways, that you mentioned in an earlier comment.

Looking at both velocity calculations neither seems to be a pure rotation, the magnitudes are not equal to 1. This suggests that there's an error in the computation.

Ad velocity magnitude (spin off from earlier comments): I am well aware that the 'old' formulas do not equal a magnitude of 1. It was fixed in the original code by approximating by '.Normalize()'. I was under the impression that the velocity formula in my 'ksp131' branch did create a magnitude 1 vector. Do you have a convincing set of arguments (alpha, beta, phi) that shows the discrepancy?

You're right, sorry. Had a mistake in my calculations. Still, the floating point errors can accumulate so I would normalize the vector afterwards anyway.

Sure. We can do that to make it easier to read the intention of the formulas.
(for what it is worth I suspect that the direct formulas use double precision floating point calculations and the Normalize method has to deal with single precision and won't improve anything).

dkavolis · 2018-11-11T19:31:48Z

GameData/FerramAerospaceResearch/FAR.version

@@ -10,7 +10,7 @@
    "MAJOR" : 0,
    "MINOR" : 15,
    "PATCH" : 9,
-    "BUILD" : 5
+    "BUILD" : 6


No need to change the version, I will do it when I publish a new release

Roger, will set it back to 5.

dkavolis · 2018-11-11T20:24:52Z

FerramAerospaceResearch/FARMathUtil.cs

+        LblSkipRight:
+            if (f1 > 0) return mfobj.LinearSolution(x0, f0, x1, f1);
+            double x2 = Clamp<double>(x1 + xstepsize, 0d, rightedge);
+            if (x2 == x1) return mfobj.BrentSolve("Reached far right edge.");


Comparison of 2 doubles should be done to a tolerance

dkavolis · 2018-11-11T20:28:54Z

Thanks for your PR. I just have a few questions and suggestions.

You replaced some spaces in localization files with explicit unicode, any reason a simple space is not sufficient? You have also removed some spaces such that values following them will appear without one.
Instead of Loop in export it might be better to use Sweep.
The input values of altitude and Mach number are taken from a file deep in KSP folder, I think it would be better to have a separate tab open where the derivatives are for inputting the values manually (which could be saved in config.xml). I was thinking of having an input list for each, something like in Paraview:

where the user could input values one by one or as a lin/logscale.
The export format is inefficient (outputs too much data in terms of definitions) and will probably be difficult to read with existing tools (I'm thinking Python). What do you think about outputting the results in csv format? Body name, aircraft properties and inertia could be left in the header and the rest in a table format. It would make importing the file that much easier into existing tools.
The exported files are stored in FerramAerospaceResearh/Plugins/PluginData which can be confusing. It might be better to store all the files in FerramAerospaceResearch/StabilityDerivatives or something similar, the default file names should also be changed so they don't overwrite each other, maybe "{CraftName}_{datetime}" such that they are unique.
Running the export at multiple points freezes the game for a short time, this should be offloaded to a different thread.
What are the benefits of your SillySearchMethod over BrentsMethod? The latter is robust, fast and popular root finding algorithm.

Rodhern · 2018-11-11T20:39:20Z

* You replaced some spaces in localization files with explicit unicode, any reason a simple space is not sufficient? You have also removed some spaces such that values following them will appear without one.

In my KSP 1.3.1 version only the explicit unicode spaces show up in the GUI. That is the reason they are made explicit. The ones that are removed are to show clearer that the characters will not appear.

* Instead of `Loop` in export it might be better to use `Sweep`.

That makes sense; thanks.

* The input values of altitude and Mach number are taken from a file deep in KSP folder, I think it would be better to have a separate tab open where the derivatives are for inputting the values manually (which could be saved in config.xml). I was thinking of having an input list for each, something like in Paraview:
  ![image](https://user-images.githubusercontent.com/12998363/48317339-41911280-e5e8-11e8-9b06-b842f01e9c75.png)
  where the user could input values one by one or as a lin/logscale.

Sounds good to me, but I am not able to code that myself.

* The export format is inefficient (outputs too much data in terms of definitions) and will probably be difficult to read with existing tools (I'm thinking Python). What do you think about outputting the results in csv format? Body name, aircraft properties and inertia could be left in the header and the rest in a table format. It would make importing the file that much easier into existing tools.

Sound good too. Again I am not sure I have the skill and/or time to make it so. Feel free to take what can be reused and create some better export functionality.

* The exported files are stored in FerramAerospaceResearh/Plugins/PluginData which can be confusing. It might be better to store all the files in FerramAerospaceResearch/StabilityDerivatives or something similar, the default file names should also be changed so they don't overwrite each other, maybe "{CraftName}_{datetime}" such that they are unique.

:-)

* Running the export at multiple points freezes the game for a short time, this should be offloaded to a different thread.

:-)

* What are the benefits of your `SillySearchMethod` over `BrentsMethod`? The latter is robust, fast and popular root finding algorithm.

The benefit is that BrentsMethod won't (as easily) find the maximum CL, which in turns produces an inability to solve the AoA for slow flight. The SillySearchMethod did better in my examples.
[Edit: You can compile two versions, one with BrentsMethod and one with SillySearchMethod, then see how slow you can go and still get an AoA solution in the SPH with the same plane. It would be interesting to see how the methods compare for your plane designs.]

dkavolis · 2018-11-11T20:47:00Z

Sounds good to me, but I am not able to code that myself.

I will probably give it a try then when I have more time.

Sound good too. Again I am not sure I have the skill and/or time to make it so. Feel free to take what can be reused and create some better export functionality.

This shouldn't be too difficult to implement, once everything is done a contour plot might be a nice additional feature.

The benefit is that BrentsMethod won't (as easily) find the maximum CL, which in turns produces an inability to solve the AoA for slow flight. The SillySearchMethod did better in my examples.

Then a switch based on Mach number could be used.

Rodhern · 2018-11-11T20:52:41Z

This shouldn't be too difficult to implement, once everything is done a contour plot might be a nice additional feature.

Agreed, that would be pretty cool.

Then a switch based on Mach number could be used.

Excellent idea. So it would kind of go: First use Brent, if that fails then go Silly ?

dkavolis · 2018-11-11T20:58:05Z

Excellent idea. So it would kind of go: First use Brent, if that fails then go Silly ?

Simply if mach number is below, say, 0.3, use SillySearchMethod, otherwise use BrentsMethod. The best switch value should be slightly larger than the largest BrentsMethod fails on.

Rodhern · 2018-11-11T21:07:14Z

Simply if mach number is below, say, 0.3, use SillySearchMethod, otherwise use BrentsMethod. The best switch value should be slightly larger than the lowest BrentsMethod fails on.

Sounds fine. My only concern was that the threshold speed might be quite craft dependent.

dkavolis · 2018-11-11T21:09:23Z

You can then add a check for failure but a fixed Mach switch should be faster on average.

Rodhern · 2018-11-11T21:11:30Z

You can then add a check for failure but a fixed Mach switch should be faster on average.

Sounds good. I will see when I get time to look at the code again, and get it fixed, maybe next weekend, maybe later. If you feel like tinkering be my guest and go ahead before then.

dkavolis · 2018-11-11T21:20:07Z

Also, the sweep should use a new CalculateStabilityDerivs method which does not recompute vehicle properties every time.

dkavolis · 2018-11-17T13:23:52Z

Do you mind submitting separate PRs for each of the changes? It would make keeping track of changes easier in case ferram ever returns and some of them can be merged already.

Rodhern · 2018-11-17T14:42:23Z

Do you mind submitting separate PRs for each of the changes? It would make keeping track of changes easier in case ferram ever returns and some of them can be merged already.

I guess I can do that. I will plan to probably do that some time.
(I am already a bit surprised at how Github does PRs; that the branch tips just keep chugging along).

dkavolis · 2018-12-22T09:22:28Z

Do you mind submitting the changes without sweep/export of the stability derivatives (can be a single PR) so I can merge them in?

Rodhern · 2018-12-31T14:15:10Z

Do you mind submitting the changes without sweep/export of the stability derivatives (can be a single PR) so I can merge them in?

That sounds like a plan to me. I will try to find some time in the new year to make an updated pull request to one of your newer branches.

I will do it as a single new untested (because I do not play any of the newest KSP versions and cannot check to see if the merged source compiles) PR, once I have removed / commented out the stability derivatives export buttons in the GUI, and added (some or most of) the code readability enhancing normalizations and tweaks discussed in the comments above.

Rodhern · 2019-01-22T18:15:06Z

Closing this PR; it is to be replaced by a newer PR.

Rodhern added 14 commits August 8, 2017 11:00

Rebase Rodhern's update 1 onto branch 'ksp131'.

808c4ab

Update project references.

9ec4bd2

Add a few explicit format characters.

cdbed43

Merge 'ksp131' onto branch 'ksp130'.

6e1957a

Add ModularFlightIntegrator to .gitignore .

dc4bee4

Merge 'ksp131' onto branch 'master'.

703abdd

Add stability derivatives export button (not functional).

c51b2d3

Add stability derivatives export variables.

eea5af4

Prepare stability derivatives file export (not functional).

8a91640

Enable stability derivatives file export.

fa53a87

Change default AoA search method to one better suited for slow flight scenarios.

Shorten stability derivatives panel shadow in the 'Longitudinal Deriv…

02f4fde

…atives' section (minor graphical tweak to GUI). Rename 'siteffG' (situation effective gravity) to 'siteffg' (named 'env.effg' when exported to Octave text file).

Add stability derivatives export loop button.

8bb48f9

Add file cache for stability derivatives export loop.

b1d6751

Merge 'Rodhern/ksp131' onto branch 'dkavolis/master'.

4bb6acc

dkavolis reviewed Nov 11, 2018

View reviewed changes

Remove 'note.txt' from .gitignore.

8e9c2c4

dkavolis reviewed Nov 11, 2018

View reviewed changes

Reverted update to FAR.version.

508aa93

dkavolis reviewed Nov 11, 2018

View reviewed changes

Rodhern closed this Jan 22, 2019

Rodhern deleted the dkavolis branch January 25, 2019 19:23

	pertOutput.Cl += nominalOutput.Cd;
	pertOutput.Cd -= nominalOutput.Cl;

Rodhern's update 3. #5

Rodhern's update 3. #5

Conversation

Rodhern commented Nov 10, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rodhern Nov 11, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rodhern Nov 11, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rodhern Nov 11, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Rodhern Nov 11, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkavolis commented Nov 11, 2018

Rodhern commented Nov 11, 2018 • edited Loading

dkavolis commented Nov 11, 2018

Rodhern commented Nov 11, 2018 • edited Loading

dkavolis commented Nov 11, 2018 • edited Loading

Rodhern commented Nov 11, 2018

dkavolis commented Nov 11, 2018

Rodhern commented Nov 11, 2018

dkavolis commented Nov 11, 2018 • edited Loading

dkavolis commented Nov 17, 2018

Rodhern commented Nov 17, 2018

dkavolis commented Dec 22, 2018

Rodhern commented Dec 31, 2018

Rodhern commented Jan 22, 2019

Rodhern Nov 11, 2018 •

edited

Loading

Rodhern Nov 11, 2018 •

edited

Loading

Rodhern Nov 11, 2018 •

edited

Loading

Rodhern Nov 11, 2018 •

edited

Loading

Rodhern commented Nov 11, 2018 •

edited

Loading

Rodhern commented Nov 11, 2018 •

edited

Loading

dkavolis commented Nov 11, 2018 •

edited

Loading

dkavolis commented Nov 11, 2018 •

edited

Loading