Skip to content

Commit

Permalink
Fixed CenterQuery's torque calculation and removed dead code
Browse files Browse the repository at this point in the history
  • Loading branch information
Benjamin Chung committed Jan 20, 2019
1 parent 753115f commit 32f0dfb
Show file tree
Hide file tree
Showing 5 changed files with 29 additions and 110 deletions.
4 changes: 2 additions & 2 deletions FerramAerospaceResearch/FARAeroComponents/FARVesselAero.cs
Original file line number Diff line number Diff line change
Expand Up @@ -308,7 +308,7 @@ private void CalculateAndApplyVesselAeroProperties()

public void SimulateAeroProperties(out Vector3 aeroForce, out Vector3 aeroTorque, Vector3 velocityWorldVector, double altitude)
{
FARCenterQuery center = new FARCenterQuery();
FARCenterQuery center = new FARCenterQuery(vessel.CoM);

float pressure;
float density;
Expand Down Expand Up @@ -352,7 +352,7 @@ public void SimulateAeroProperties(out Vector3 aeroForce, out Vector3 aeroTorque
}

aeroForce = center.force;
aeroTorque = center.TorqueAt(vessel.CoM);
aeroTorque = center.Torque();
}


Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1519,7 +1519,7 @@ private void CalculateVesselAeroProperties()

public void UpdateSonicDragArea()
{
ferram4.FARCenterQuery center = new ferram4.FARCenterQuery();
ferram4.FARCenterQuery center = new ferram4.FARCenterQuery(Vector3.zero);

Vector3 worldMainAxis = _localToWorldMatrix.MultiplyVector(_vehicleMainAxis);
worldMainAxis.Normalize();
Expand Down
37 changes: 15 additions & 22 deletions FerramAerospaceResearch/FARGUI/FAREditorGUI/EditorAeroCenter.cs
Original file line number Diff line number Diff line change
Expand Up @@ -82,10 +82,6 @@ public void UpdateAeroData(List<FARAeroPartModule> aeroModules, List<FARAeroSect

void UpdateAerodynamicCenter()
{
FARCenterQuery aeroSection, dummy;
aeroSection = new FARCenterQuery();
dummy = new FARCenterQuery();

if((object)EditorLogic.RootPart == null)
return;

Expand All @@ -103,33 +99,33 @@ void UpdateAerodynamicCenter()
vel_fuzz = -0.02f * Vector3.forward;
}

Vector3 vel = (vel_base - vel_fuzz).normalized;

for(int i = 0; i < _currentAeroSections.Count; i++)
{
FARAeroSection section = _currentAeroSections[i];
section.PredictionCalculateAeroForces(1, 0.5f, 100000, 0, 0.005f, vel, aeroSection);
}

FARBaseAerodynamics.PrecomputeGlobalCenterOfLift(aeroSection, dummy, vel, 1);

Vector3 pos = Vector3.zero;//rootPartTrans.position;
float mass = 0;
for(int i = 0; i < EditorLogic.SortedShipList.Count; i++)
for (int i = 0; i < EditorLogic.SortedShipList.Count; i++)
{
Part p = EditorLogic.SortedShipList[i];
float tmpMass = p.mass + p.GetResourceMass();
mass += tmpMass;
pos += p.partTransform.position * tmpMass;
}
pos /= mass;
FARCenterQuery aeroSection = new FARCenterQuery(pos);
FARCenterQuery dummy = new FARCenterQuery(pos);

Vector3 vel = (vel_base - vel_fuzz).normalized;

for (int i = 0; i < _currentAeroSections.Count; i++)
{
FARAeroSection section = _currentAeroSections[i];
section.PredictionCalculateAeroForces(1, 0.5f, 100000, 0, 0.005f, vel, aeroSection);
}

FARBaseAerodynamics.PrecomputeGlobalCenterOfLift(aeroSection, dummy, vel, 1);

Vector3 avgForcePos = Vector3.zero;

Vector3 force0, moment0;
force0 = aeroSection.force;
moment0 = aeroSection.TorqueAt(pos);
avgForcePos += aeroSection.GetPos();
moment0 = aeroSection.Torque();

//aeroSection.force = -aeroSection.force;
//aeroSection.torque = -aeroSection.torque;
Expand All @@ -148,13 +144,10 @@ void UpdateAerodynamicCenter()

Vector3 force1, moment1;
force1 = aeroSection.force;
moment1 = aeroSection.TorqueAt(pos);
avgForcePos += aeroSection.GetPos();
moment1 = aeroSection.Torque();

aeroSection.ClearAll();

avgForcePos *= 0.5f;

Vector3 deltaForce = force1 - force0;
Vector3 deltaMoment = moment1 - moment0;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ public void GetClCdCmSteady(InstantConditionSimInput input, out InstantCondition
MAC += w.GetMAC() * w.S;
b_2 += w.Getb_2() * w.S;
}
FARCenterQuery center = new FARCenterQuery();
FARCenterQuery center = new FARCenterQuery(CoM);
for (int i = 0; i < _currentAeroSections.Count; i++)
{
_currentAeroSections[i].PredictionCalculateAeroForces(2, (float)input.machNumber, 10000, 0, 0.005f, velocity.normalized, center);
Expand All @@ -213,7 +213,7 @@ public void GetClCdCmSteady(InstantConditionSimInput input, out InstantCondition
output.Cy += Vector3d.Dot(centerForce, sideways);
output.Cd += -Vector3d.Dot(centerForce, velocity);

Vector3d centerMoment = -center.TorqueAt(CoM) * 1000;
Vector3d centerMoment = -center.Torque() * 1000;

output.Cm += Vector3d.Dot(centerMoment, sideways);
output.Cn += Vector3d.Dot(centerMoment, liftVector);
Expand Down
92 changes: 9 additions & 83 deletions FerramAerospaceResearch/LEGACYferram4/FARCenterQuery.cs
Original file line number Diff line number Diff line change
Expand Up @@ -56,49 +56,25 @@ public class FARCenterQuery
// Torque needed to compensate if force were applied at origin.
public Vector3d torque = Vector3d.zero;

// Weighted average of force positions used as aid in choosing the
// single center location on the line of physically equivalent ones.
// Reference point about which to calculate torques
public Vector3d pos = Vector3d.zero;
public double amount = 0.0;

/*public Vector3d GetACPosition()
public FARCenterQuery(Vector3 about)
{
Vector3d bVec = new Vector3d();
bVec.x = force.y * torque.z - force.z * torque.y;
bVec.y = force.x * torque.z - force.z * torque.x;
bVec.z = force.x * torque.y - force.y * torque.x;
double approxA01, approxA02, approxA12; //commonly used force directions
approxA01 = -1 / (force.z * force.z);
approxA02 = -1 / (force.y * force.y);
approxA12 = -1 / (force.x * force.x);
Vector3d acPos = new Vector3d();
acPos.x = force.x * force.x / (force.y * force.y * force.z * force.z) * bVec.x + bVec.y * approxA01 + bVec.z * approxA02;
acPos.y = force.y * force.y / (force.x * force.x * force.z * force.z) * bVec.y + bVec.x * approxA01 + bVec.z * approxA12;
acPos.x = force.z * force.z / (force.y * force.y * force.x * force.x) * bVec.z + bVec.y * approxA12 + bVec.x * approxA02;
acPos *= 0.5;
acPos += GetPos();
return acPos;
}*/
this.pos = about;
}

public void ClearAll()
{
force = Vector3d.zero;
torque = Vector3d.zero;
pos = Vector3d.zero;
amount = 0;
}

// Record a force applied at a point
public void AddForce(Vector3d npos, Vector3d nforce)
{
double size = nforce.magnitude;
force += nforce;
torque += Vector3d.Cross(npos, nforce);
pos += npos * size;
amount += size;
torque += Vector3d.Cross(npos - pos, nforce);
}

// Record an abstracted torque or couple; application point is irrelevant
Expand All @@ -107,60 +83,10 @@ public void AddTorque(Vector3d ntorque)
torque += ntorque;
}

// Merge two force sets
public void AddAll(FARCenterQuery q2)
{
force += q2.force;
torque += q2.torque;
pos += q2.pos;
amount += q2.amount;
}

// Returns a center of weight-like average of force positions.
// Unless all forces are strictly parallel it doesn't mean much.
public Vector3d GetPos()
{
return amount > 0 ? pos / amount : Vector3d.zero;
}

public void SetPos(Vector3d npos)
{
pos = npos;
amount = 1;
}

// Compensating torque at different origin.
public Vector3d TorqueAt(Vector3d origin)
{
return torque - Vector3d.Cross(origin, force);
}

// Returns a point that requires minimal residual torque
// (or even 0 if possible) and is closest to origin.
// Any remaining torque is always parallel to force.
public Vector3d GetMinTorquePos(Vector3d origin)
{
double fmag = force.sqrMagnitude;
if (fmag <= 0) return origin;

return origin + Vector3d.Cross(force, TorqueAt(origin)) / fmag;
}

public Vector3d GetMinTorquePos()
{
return GetMinTorquePos(GetPos());
}

// The physics engine limits torque that can be applied to a single
// object. This tries to replicate it based on results of experiments.
// In practice this is probably not necessary for FAR, but since this
// knowledge has been obtained, might as well turn it into code.
public static float TorqueClipFactor(Vector3 torque, Rigidbody body)
// Compute torque.
public Vector3d Torque()
{
Vector3 tq = Quaternion.Inverse(body.rotation * body.inertiaTensorRotation) * torque;
Vector3 tensor = body.inertiaTensor;
float acceleration = new Vector3(tq.x / tensor.x, tq.y / tensor.y, tq.z / tensor.z).magnitude;
return Mathf.Max(1.0f, acceleration * Time.fixedDeltaTime / body.maxAngularVelocity);
return torque;
}
}
}

0 comments on commit 32f0dfb

Please sign in to comment.