Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Rodhern's update 3. #5

Closed
wants to merge 16 commits into from
Closed

Rodhern's update 3. #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
808c4ab
Rebase Rodhern's update 1 onto branch 'ksp131'.
Rodhern Aug 8, 2017
9ec4bd2
Update project references.
Rodhern Aug 12, 2017
cdbed43
Add a few explicit format characters.
Rodhern Aug 14, 2017
6e1957a
Merge 'ksp131' onto branch 'ksp130'.
Rodhern Aug 25, 2017
dc4bee4
Add ModularFlightIntegrator to .gitignore .
Rodhern Nov 7, 2017
703abdd
Merge 'ksp131' onto branch 'master'.
Rodhern Apr 8, 2018
c51b2d3
Add stability derivatives export button (not functional).
Rodhern Apr 13, 2018
eea5af4
Add stability derivatives export variables.
Rodhern Apr 28, 2018
8a91640
Prepare stability derivatives file export (not functional).
Rodhern Apr 30, 2018
fa53a87
Enable stability derivatives file export.
Rodhern May 2, 2018
02f4fde
Shorten stability derivatives panel shadow in the 'Longitudinal Deriv…
Rodhern May 3, 2018
8bb48f9
Add stability derivatives export loop button.
Rodhern May 6, 2018
b1d6751
Add file cache for stability derivatives export loop.
Rodhern May 9, 2018
4bb6acc
Merge 'Rodhern/ksp131' onto branch 'dkavolis/master'.
Rodhern Nov 10, 2018
8e9c2c4
Remove 'note.txt' from .gitignore.
Rodhern Nov 11, 2018
508aa93
Reverted update to FAR.version.
Rodhern Nov 11, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions FerramAerospaceResearch/FARGUI/FAREditorGUI/EditorGUI.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,8 @@ void OnDestroy()
#region EditorEvents
private void ResetEditorEvent(ShipConstruct construct)
{
FARAeroUtil.ResetEditorParts();

if (EditorLogic.RootPart != null)
{
List<Part> partsList = EditorLogic.SortedShipList;
Expand Down
12 changes: 6 additions & 6 deletions FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/InstantConditionSim.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,15 +141,15 @@ public void GetClCdCmSteady(InstantConditionSimInput input, out InstantCondition
AngVel += (sinPhi * alphaDot - cosAlpha * cosPhi * betaDot) * up;

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

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.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

@Rodhern Rodhern Nov 11, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[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?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


//this is negative wrt the ground
Vector3d liftVector = -forward * sinAlpha + right * sinPhi * cosAlpha - up * cosPhi * cosAlpha;

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.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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



for (int i = 0; i < _wingAerodynamicModel.Count; i++)
Expand All @@ -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);
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

@Rodhern Rodhern Nov 11, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

@Rodhern Rodhern Nov 11, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

@Rodhern Rodhern Nov 11, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).


if (clear)
w.EditorClClear(reset_stall);
Expand Down
49 changes: 42 additions & 7 deletions FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/StabilityDerivCalculator.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,33 @@ public StabilityDerivCalculator(InstantConditionSim instantConditionSim)
_instantCondition = instantConditionSim;
}

public StabilityDerivOutput CalculateStabilityDerivs(double u0, double q, double machNumber, double alpha, double beta, double phi, int flapSetting, bool spoilers, CelestialBody body, double alt)
public StabilityDerivExportOutput CalculateStabilityDerivs(CelestialBody body, double alt, double machNumber, int flapSetting, bool spoilers)
{
// Rodhern: Only return a result if one is found.

if (body.GetPressure(alt) > 0)
{
StabilityDerivExportOutput result = CalculateStabilityDerivs(body, alt, machNumber, flapSetting, spoilers, 0, 0, 0);
if (result.outputvals.stableAoAState == "")
return result;
else
return null;
}
else
return null;
}

public StabilityDerivExportOutput CalculateStabilityDerivs(CelestialBody body, double alt, double machNumber, int flapSetting, bool spoilers, double alpha, double beta, double phi)
{
double pressure = body.GetPressure(alt);
double temperature = body.GetTemperature(alt);
double density = body.GetDensity(pressure, temperature);
double sspeed = body.GetSpeedOfSound(pressure, density);
double u0 = sspeed * machNumber;
double q = u0 * u0 * density * 0.5f;

StabilityDerivOutput stabDerivOutput = new StabilityDerivOutput();
StabilityDerivExportVariables stabDerivExport = new StabilityDerivExportVariables();
stabDerivOutput.nominalVelocity = u0;
stabDerivOutput.altitude = alt;
stabDerivOutput.body = body;
Expand Down Expand Up @@ -210,7 +234,7 @@ public StabilityDerivOutput CalculateStabilityDerivs(double u0, double q, double
//Longitudinal Mess
_instantCondition.SetState(machNumber, neededCl, CoM, 0, input.flaps, input.spoilers);

alpha = FARMathUtil.BrentsMethod(_instantCondition.FunctionIterateForAlpha, -30d, 30d, 0.001, 500);
alpha = FARMathUtil.SillySearchMethod(_instantCondition.FunctionIterateForAlpha);
input.alpha = alpha;
nominalOutput = _instantCondition.iterationOutput;
//alpha_str = (alpha * Mathf.PI / 180).ToString();
Expand Down Expand Up @@ -275,8 +299,8 @@ public StabilityDerivOutput CalculateStabilityDerivs(double u0, double q, double
pertOutput.Cd = (pertOutput.Cd - nominalOutput.Cd) / 0.05;
pertOutput.Cm = (pertOutput.Cm - nominalOutput.Cm) / 0.05;

pertOutput.Cl *= q * area * MAC / (2 * u0 * mass);
pertOutput.Cd *= q * area * MAC / (2 * u0 * mass);
pertOutput.Cl *= -q * area * MAC / (2 * u0 * mass); // Rodhern: Replaced 'q' by '-q', so that formulas
pertOutput.Cd *= -q * area * MAC / (2 * u0 * mass); // for Zq and Xq match those for Zu and Xu.
pertOutput.Cm *= q * area * MAC * MAC / (2 * u0 * Iy);

stabDerivOutput.stabDerivs[9] = pertOutput.Cl; //Zq
Expand All @@ -292,8 +316,8 @@ public StabilityDerivOutput CalculateStabilityDerivs(double u0, double q, double
pertOutput.Cd = (pertOutput.Cd - nominalOutput.Cd) / 0.1;
pertOutput.Cm = (pertOutput.Cm - nominalOutput.Cm) / 0.1;

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.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

pertOutput.Cm *= q * area * MAC / Iy;

stabDerivOutput.stabDerivs[12] = pertOutput.Cl; //Ze
Expand Down Expand Up @@ -357,7 +381,18 @@ public StabilityDerivOutput CalculateStabilityDerivs(double u0, double q, double
stabDerivOutput.stabDerivs[23] = pertOutput.Cn; //Nr
stabDerivOutput.stabDerivs[22] = pertOutput.C_roll; //Lr

return stabDerivOutput;
// Assign values to export variables
stabDerivExport.craftmass = mass;
stabDerivExport.envpressure = pressure;
stabDerivExport.envtemperature = temperature;
stabDerivExport.envdensity = density;
stabDerivExport.envsoundspeed = sspeed;
stabDerivExport.envg = _instantCondition.CalculateAccelerationDueToGravity(body, alt);
stabDerivExport.sitmach = machNumber;
stabDerivExport.sitdynpres = q;
stabDerivExport.siteffg = effectiveG;

return new StabilityDerivExportOutput(stabDerivOutput, stabDerivExport);
}

}
Expand Down
32 changes: 32 additions & 0 deletions FerramAerospaceResearch/FARGUI/FAREditorGUI/Simulation/StabilityDerivExportOutput.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
using System;
using System.Collections.Generic;

namespace FerramAerospaceResearch.FARGUI.FAREditorGUI.Simulation
{
class StabilityDerivExportVariables
{
public double craftmass;

public double envpressure; // in kPa
public double envtemperature; // in Kelvin
public double envdensity;
public double envsoundspeed;
public double envg;

public double sitmach;
public double sitdynpres;
public double siteffg; // local gravity corrected for speed
}

class StabilityDerivExportOutput
{
public StabilityDerivOutput outputvals;
public StabilityDerivExportVariables exportvals;

public StabilityDerivExportOutput(StabilityDerivOutput outputvalues, StabilityDerivExportVariables exportvalues)
{
this.outputvals = outputvalues;
this.exportvals = exportvalues;
}
}
}
Loading