Skip to content

Commit

Permalink
Adding capability to aim external forces in inertial frame (JSBSim-Te…
Browse files Browse the repository at this point in the history 
…am#962)
  • Loading branch information
@jonsberndt
jonsberndt authored Sep 24, 2023
1 parent c9b12c4 commit ae4652a
Show file tree
Hide file tree
Showing 8 changed files with 231 additions and 12 deletions.
67 changes: 60 additions & 7 deletions aircraft/ball/ball.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@
<ixy unit="SLUG*FT2"> -0 </ixy>
<ixz unit="SLUG*FT2"> -0 </ixz>
<iyz unit="SLUG*FT2"> -0 </iyz>
<emptywt unit="LBS"> 40 </emptywt>
<emptywt unit="LBS"> 20000 </emptywt>
<location name="CG" unit="IN">
<x> 0 </x>
<y> 0 </y>
Expand Down Expand Up @@ -100,14 +100,15 @@
<external_reactions>
<!-- "Declare" the reefing term -->
<property>fcs/parachute_reef_pos_norm</property>

<property>propulsion/rocket_thrust</property>

<force name="parachute" frame="WIND">
<function>
<product>
<property>aero/qbar-psf</property>
<property>fcs/parachute_reef_pos_norm</property>
<value> 1.0 </value> <!-- Full drag coefficient -->
<value> 20.0 </value> <!-- Full parachute area -->
<value> 10000.0 </value> <!-- Full parachute area -->
</product>
</function>
<!-- The location below is in structural frame (x positive
Expand All @@ -121,19 +122,71 @@
<!-- The direction describes a unit vector. In this case, since
the selected frame is the WIND frame, the "-1" x component
describes a direction exactly opposite of the direction
into the wind vector. That is, the direction specified below
into the WIND vector. That is, the direction specified below
is the direction that the drag force acts in. -->
<direction>
<x>-1</x>
<x> -1 </x>
<y> 0 </y>
<z> 0 </z>
</direction>
</force>

<force name="rocket" frame="INERTIAL">
<function>
<property> propulsion/rocket_thrust </property>
</function>
<location unit="FT">
<x>0</x>
<y>0</y>
<z>0</z>
</direction>
</location>
<!-- The direction that the thrust acts in is defined below in the Pointing
channel of the flight_control section, through the creation of a unit
pointing vector. That is, the direction specified below
is the direction that the drag force acts in. -->
</force>

</external_reactions>

<propulsion/>

<flight_control name="FGFCS"/>
<flight_control name="FGFCS">

<channel name="Pointing">

<fcs_function name="propulsion/tvc_inertial_x">
<function>
<quotient>
<property> velocities/eci-x-fps </property>
<property> velocities/eci-velocity-mag-fps </property>
</quotient>
</function>
<output> external_reactions/rocket/x </output>
</fcs_function>

<fcs_function name="propulsion/tvc_inertial_y">
<function>
<quotient>
<property> velocities/eci-y-fps </property>
<property> velocities/eci-velocity-mag-fps </property>
</quotient>
</function>
<output> external_reactions/rocket/y </output>
</fcs_function>

<fcs_function name="propulsion/tvc_inertial_z">
<function>
<quotient>
<property> velocities/eci-z-fps </property>
<property> velocities/eci-velocity-mag-fps </property>
</quotient>
</function>
<output> external_reactions/rocket/z </output>
</fcs_function>

</channel>

</flight_control>

<aerodynamics>
<axis name="DRAG">
Expand Down
29 changes: 29 additions & 0 deletions aircraft/ball/iss_orbit.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
<?xml version="1.0"?>
<initialize name="reset01" version="2.0">
<!--
This file sets up the spacecraft in a hig inclination elliptical orbit.
2023-09-17T01:31:30.000 -4084.032825462220 -5429.695419910110 -0.556775951837 3.81325195203505 -2.85811815725326 6.00306592847078
-->

<position frame="ECI" unit="KM">
<x> -4084.03 </x>
<y> -5429.69 </y>
<z> -0.5567 </z>
</position>

<orientation unit="DEG" frame="LOCAL">
<yaw> 90.0 </yaw>
</orientation>

<velocity unit="KM/SEC" frame="ECI">
<x> 3.813252 </x>
<y> -2.858118 </y>
<z> 6.003066 </z>
</velocity>

<attitude_rate unit="DEG/SEC" frame="ECI">
<y> 0.0 </y>
<z> 0.0 </z>
</attitude_rate>

</initialize>
120 changes: 120 additions & 0 deletions scripts/ball_orbit_phase.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet type="text/xsl" href="http://jsbsim.sf.net/JSBSimScript.xsl"?>
<runscript xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:noNamespaceSchemaLocation="http://jsbsim.sf.net/JSBSimScript.xsd"
name="Orbital Phasing Test">
<use aircraft="ball" initialize="iss_orbit"/>
<run start="0.0" end="11000" dt="0.005">
<description>
Integrators

0: No integrator (Freeze)
1: Rectangular Euler
2: Trapezoidal
3: Adams Bashforth 2
4: Adams Bashforth 3
5: Adams Bashforth 4

</description>
<property value="3"> simulation/integrator/rate/rotational </property>
<property value="3"> simulation/integrator/rate/translational </property>
<property value="1"> simulation/integrator/position/rotational </property>
<property value="4"> simulation/integrator/position/translational </property>
<property value="1"> simulation/gravity-model </property>
<property value="1"> simulation/output[0]/log_rate_hz </property>

<property> simulation/notify-time-trigger </property>

<event name="Time Notify" persistent="true">
<description>Output message at 10 minute intervals</description>
<notify>
<property caption="Height AGL (ft): "> position/h-agl-ft </property>
<property caption="Altitude (Geod, ft): "> position/geod-alt-ft </property>
<property caption="Latitude (Geod, deg): "> position/lat-geod-deg </property>
<property caption="Vehicle Radius (ft): "> position/radius-to-vehicle-ft </property>
<property caption="Inertial Vel Mag (ft/s): "> velocities/eci-velocity-mag-fps </property>
<property caption="Body X Velocity (ft/s): "> velocities/u-fps </property>
<property caption="Simulation Frame: "> simulation/frame </property>
<property caption="Density: "> atmosphere/rho-slugs_ft3 </property>
<property caption="Angular momentum (ft^2/s)"> orbital/specific-angular-momentum-ft2_sec </property>
<property caption="Inclination (deg) "> orbital/inclination-deg </property>
<property caption="Right ascension (deg) "> orbital/right-ascension-deg </property>
<property caption="Argument of perigee (deg)"> orbital/argument-of-perigee-deg </property>
<property caption="Period (s) "> orbital/period-sec </property>
<property caption="Eccentricity "> orbital/eccentricity </property>
<property caption="Apoapsis radius (ft) "> orbital/apoapsis-radius-ft </property>
<property caption="Periapsis radius (ft) "> orbital/periapsis-radius-ft </property>
<property caption="True anomaly (deg) "> orbital/true-anomaly-deg </property>
</notify>
<condition> simulation/sim-time-sec >= simulation/notify-time-trigger </condition>
<set name="simulation/notify-time-trigger" value="600" type="FG_DELTA"/>
</event>

<event name="Undock burn">
<description> Lower perigee </description>
<condition>
simulation/sim-time-sec gt 10
</condition>
<set name="propulsion/rocket_thrust" value="2000" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

<event name="End Undock burn">
<description> Stop lowering perigee </description>
<condition>
simulation/sim-time-sec gt 12
</condition>
<set name="propulsion/rocket_thrust" value="0" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

<event name="Orbit Lowering Burn">
<description> Lower orbit below station (Thrust is in negative direction here)</description>
<condition>
simulation/sim-time-sec gt 2803
</condition>
<set name="propulsion/rocket_thrust" value="-2000" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

<event name="End Orbit Lowering Burn">
<description> Stop lowering orbit below station </description>
<condition>
simulation/sim-time-sec gt 2808
</condition>
<set name="propulsion/rocket_thrust" value="0" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

<event name="Orbit Circularizing Burn">
<description> Circularize orbit below station (Thrust is in negative direction here)</description>
<condition>
simulation/sim-time-sec gt 5591
</condition>
<set name="propulsion/rocket_thrust" value="-2000" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

<event name="End Orbit Circularizing Burn">
<description> Stop circularizing orbit below station </description>
<condition>
simulation/sim-time-sec gt 5598
</condition>
<set name="propulsion/rocket_thrust" value="0" action="FG_RAMP" tc="1.0"/>
<notify>
<property>simulation/sim-time-sec</property>
</notify>
</event>

</run>
</runscript>
5 changes: 5 additions & 0 deletions src/models/FGExternalForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,8 @@ FGParameter* FGExternalForce::bind(Element *el, const string& magName,
ttype = tLocalBody;
} else if (sFrame == "WIND") {
ttype = tWindBody;
} else if (sFrame == "INERTIAL") {
ttype = tInertialBody;
} else {
cerr << el->ReadFrom()
<< "Invalid frame specified for external " << el->GetName() << ", \""
Expand Down Expand Up @@ -239,6 +241,9 @@ void FGExternalForce::Debug(int from)
case tWindBody:
cout << "WIND";
break;
case tInertialBody:
cout << "INERTIAL";
break;
default:
cout << "ERROR/UNKNOWN";
}
Expand Down
1 change: 1 addition & 0 deletions src/models/FGExternalForce.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,7 @@ CLASS DOCUMENTATION
completes the right handed system. This is modified from a normal
wind frame definition, which is rotated about the Y axis 180 degrees
from this WIND frame.
- INERTIAL frame refers to the ECI frame.
Much of the substance of this class is located in the FGForce base class,
from which this class is derived.
Expand Down
2 changes: 2 additions & 0 deletions src/models/propulsion/FGForce.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,8 @@ const FGMatrix33& FGForce::Transform(void) const
return fdmex->GetAuxiliary()->GetTw2b();
case tLocalBody:
return fdmex->GetPropagate()->GetTl2b();
case tInertialBody:
return fdmex->GetPropagate()->GetTi2b();
case tCustom:
case tNone:
return mT;
Expand Down
2 changes: 1 addition & 1 deletion src/models/propulsion/FGForce.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ class JSBSIM_API FGForce : public FGJSBBase
/// Destructor
virtual ~FGForce();

enum TransformType { tNone, tWindBody, tLocalBody, tCustom };
enum TransformType { tNone, tWindBody, tLocalBody, tInertialBody, tCustom };

virtual const FGColumnVector3& GetBodyForces(void);

Expand Down
17 changes: 13 additions & 4 deletions tests/TestExternalReactions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,12 @@

from JSBSim_utils import JSBSimTestCase, CreateFDM, RunTest, CopyAircraftDef

def getParachuteArea(tree):
parachute_area = 1.0
for value in tree.getroot().findall('external_reactions/force/function/product/value'):
parachute_area *= float(value.text)
return parachute_area

class TestExternalReactions(JSBSimTestCase):
def getLeverArm(self, fdm, name):
lax = (fdm['external_reactions/'+name+'/location-x-in']
Expand Down Expand Up @@ -52,9 +58,12 @@ def test_wind_frame(self):
self.assertAlmostEqual(fdm['external_reactions/parachute/y'], 0.0)
self.assertAlmostEqual(fdm['external_reactions/parachute/z'], 0.0)

tree, _, _ = CopyAircraftDef(script_path, self.sandbox)
parachute_area = getParachuteArea(tree)

while fdm.run():
Tw2b = fdm.get_auxiliary().get_Tw2b()
mag = fdm['aero/qbar-psf'] * fdm['fcs/parachute_reef_pos_norm']*20.0
mag = fdm['aero/qbar-psf'] * fdm['fcs/parachute_reef_pos_norm']*parachute_area
f = Tw2b * np.mat([-1.0, 0.0, 0.0]).T * mag
self.assertAlmostEqual(fdm['forces/fbx-external-lbs'], f[0, 0])
self.assertAlmostEqual(fdm['forces/fby-external-lbs'], f[1, 0])
Expand Down Expand Up @@ -199,8 +208,7 @@ def test_body_frame(self):
def test_moment(self):
script_path = self.sandbox.path_to_jsbsim_file('scripts',
'ball_chute.xml')
tree, aircraft_name, _ = CopyAircraftDef(script_path,
self.sandbox)
tree, aircraft_name, _ = CopyAircraftDef(script_path, self.sandbox)
extReact_element = tree.getroot().find('external_reactions')
moment_element = et.SubElement(extReact_element, 'moment')
moment_element.attrib['name'] = 'parachute'
Expand Down Expand Up @@ -228,10 +236,11 @@ def test_moment(self):
self.assertAlmostEqual(fdm['external_reactions/parachute/n'], mDir[2])

fdm['external_reactions/parachute/magnitude-lbsft'] = -3.5
parachute_area = getParachuteArea(tree)

while fdm.run():
Tw2b = fdm.get_auxiliary().get_Tw2b()
mag = fdm['aero/qbar-psf'] * fdm['fcs/parachute_reef_pos_norm']*20.0
mag = fdm['aero/qbar-psf'] * fdm['fcs/parachute_reef_pos_norm']*parachute_area
f = Tw2b * np.mat([-1.0, 0.0, 0.0]).T * mag
self.assertAlmostEqual(fdm['forces/fbx-external-lbs'], f[0, 0])
self.assertAlmostEqual(fdm['forces/fby-external-lbs'], f[1, 0])
Expand Down

0 comments on commit ae4652a

Please sign in to comment.