1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 Author: Jon Berndt, Andreas Gaeb, David Culp
7 ------------- Copyright (C) 2011 Jon S. Berndt (jon@jsbsim.org) -------------
9 This program is free software; you can redistribute it and/or modify it under
10 the terms of the GNU Lesser General Public License as published by the Free Software
11 Foundation; either version 2 of the License, or (at your option) any later
14 This program is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
19 You should have received a copy of the GNU Lesser General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 59 Temple
21 Place - Suite 330, Boston, MA 02111-1307, USA.
23 Further information about the GNU Lesser General Public License can also be found on
24 the world wide web at http://www.gnu.org.
27 --------------------------------------------------------------------------------
30 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
37 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
41 #include "models/FGModel.h"
42 #include "math/FGColumnVector3.h"
43 #include "math/FGMatrix33.h"
44 #include "math/FGTable.h"
46 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
50 #define ID_WINDS "$Id: FGWinds.h,v 1.5 2011/09/07 12:21:45 jberndt Exp $"
52 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
58 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
62 /** Models atmospheric disturbances: winds, gusts, turbulence, downbursts, etc.
64 Various turbulence models are available. They are specified
65 via the property <tt>atmosphere/turb-type</tt>. The following models are
67 - 0: ttNone (turbulence disabled)
70 - 3: ttMilspec (Dryden spectrum)
71 - 4: ttTustin (Dryden spectrum)
73 The Milspec and Tustin models are described in the Yeager report cited below.
74 They both use a Dryden spectrum model whose parameters (scale lengths and intensities)
75 are modelled according to MIL-F-8785C. Parameters are modelled differently
76 for altitudes below 1000ft and above 2000ft, for altitudes in between they
77 are interpolated linearly.
79 The two models differ in the implementation of the transfer functions
80 described in the milspec.
82 To use one of these two models, set <tt>atmosphere/turb-type</tt> to 4 resp. 5,
83 and specify values for <tt>atmosphere/turbulence/milspec/windspeed_at_20ft_AGL-fps<tt>
84 and <tt>atmosphere/turbulence/milspec/severity<tt> (the latter corresponds to
85 the probability of exceedence curves from Fig. 7 of the milspec, allowable
86 range is 0 (disabled) to 7). <tt>atmosphere/psiw-rad</tt> is respected as well;
87 note that you have to specify a positive wind magnitude to prevent psiw from
90 Reference values (cf. figures 7 and 9 from the milspec):
92 <tr><td><b>Intensity</b></td>
93 <td><b><tt>windspeed_at_20ft_AGL-fps</tt></b></td>
94 <td><b><tt>severity</tt></b></td></tr>
96 <td>25 (15 knots)</td>
99 <td>50 (30 knots)</td>
102 <td>75 (45 knots)</td>
106 @see Yeager, Jessie C.: "Implementation and Testing of Turbulence Models for
108 href="http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980028448_1998081596.pdf">
109 pdf</a>), NASA CR-1998-206937, 1998
111 @see MIL-F-8785C: Military Specification: Flying Qualities of Piloted Aircraft
115 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
119 class FGWinds : public FGModel {
126 /** Runs the winds model; called by the Executive
127 Can pass in a value indicating if the executive is directing the simulation to Hold.
128 @param Holding if true, the executive has been directed to hold the sim from
129 advancing time. Some models may ignore this flag, such as the Input
130 model, which may need to be active to listen on a socket for the
131 "Resume" command to be given.
132 @return false if no error */
133 bool Run(bool Holding);
134 bool InitModel(void);
135 enum tType {ttNone, ttStandard, ttCulp, ttMilspec, ttTustin} turbType;
137 // TOTAL WIND access functions (wind + gust + turbulence)
139 /// Retrieves the total wind components in NED frame.
140 virtual const FGColumnVector3& GetTotalWindNED(void) const { return vTotalWindNED; }
142 /// Retrieves a total wind component in NED frame.
143 virtual double GetTotalWindNED(int idx) const {return vTotalWindNED(idx);}
145 // WIND access functions
147 /// Sets the wind components in NED frame.
148 virtual void SetWindNED(double wN, double wE, double wD) { vWindNED(1)=wN; vWindNED(2)=wE; vWindNED(3)=wD;}
150 /// Sets a wind component in NED frame.
151 virtual void SetWindNED(int idx, double wind) { vWindNED(idx)=wind;}
153 /// Retrieves the wind components in NED frame.
154 virtual FGColumnVector3& GetWindNED(void) { return vWindNED; }
156 /// Retrieves a wind component in NED frame.
157 virtual double GetWindNED(int idx) const {return vWindNED(idx);}
159 /** Retrieves the direction that the wind is coming from.
160 The direction is defined as north=0 and increases counterclockwise.
161 The wind heading is returned in radians.*/
162 virtual double GetWindPsi(void) const { return psiw; }
164 /** Sets the direction that the wind is coming from.
165 The direction is defined as north=0 and increases counterclockwise to 2*pi (radians). The
166 vertical component of wind is assumed to be zero - and is forcibly set to zero. This function
167 sets the vWindNED vector components based on the supplied direction. The magnitude of
168 the wind set in the vector is preserved (assuming the vertical component is non-zero).
169 @param dir wind direction in the horizontal plane, in radians.*/
170 virtual void SetWindPsi(double dir);
172 virtual void SetWindspeed(double speed);
174 virtual double GetWindspeed(void) const;
176 // GUST access functions
178 /// Sets a gust component in NED frame.
179 virtual void SetGustNED(int idx, double gust) { vGustNED(idx)=gust;}
181 /// Sets a turbulence component in NED frame.
182 virtual void SetTurbNED(int idx, double turb) { vTurbulenceNED(idx)=turb;}
184 /// Sets the gust components in NED frame.
185 virtual void SetGustNED(double gN, double gE, double gD) { vGustNED(eNorth)=gN; vGustNED(eEast)=gE; vGustNED(eDown)=gD;}
187 /// Retrieves a gust component in NED frame.
188 virtual double GetGustNED(int idx) const {return vGustNED(idx);}
190 /// Retrieves a turbulence component in NED frame.
191 virtual double GetTurbNED(int idx) const {return vTurbulenceNED(idx);}
193 /// Retrieves the gust components in NED frame.
194 virtual FGColumnVector3& GetGustNED(void) {return vGustNED;}
196 /** Turbulence models available: ttNone, ttStandard, ttBerndt, ttCulp, ttMilspec, ttTustin */
197 virtual void SetTurbType(tType tt) {turbType = tt;}
198 virtual tType GetTurbType() const {return turbType;}
200 virtual void SetTurbGain(double tg) {TurbGain = tg;}
201 virtual double GetTurbGain() const {return TurbGain;}
203 virtual void SetTurbRate(double tr) {TurbRate = tr;}
204 virtual double GetTurbRate() const {return TurbRate;}
206 virtual void SetRhythmicity(double r) {Rhythmicity=r;}
207 virtual double GetRhythmicity() const {return Rhythmicity;}
209 virtual double GetTurbPQR(int idx) const {return vTurbPQR(idx);}
210 virtual double GetTurbMagnitude(void) const {return Magnitude;}
211 virtual const FGColumnVector3& GetTurbDirection(void) const {return vDirection;}
212 virtual const FGColumnVector3& GetTurbPQR(void) const {return vTurbPQR;}
214 virtual void SetWindspeed20ft(double ws) { windspeed_at_20ft = ws;}
215 virtual double GetWindspeed20ft() const { return windspeed_at_20ft;}
217 /// allowable range: 0-7, 3=light, 4=moderate, 6=severe turbulence
218 virtual void SetProbabilityOfExceedence( int idx) {probability_of_exceedence_index = idx;}
219 virtual int GetProbabilityOfExceedence() const { return probability_of_exceedence_index;}
221 // Stores data defining a 1 - cosine gust profile that builds up, holds steady
222 // and fades out over specified durations.
223 struct OneMinusCosineProfile {
224 bool Running; ///<- This flag is set true through FGWinds::StartGust().
225 double elapsedTime; ///<- Stores the elapsed time for the ongoing gust.
226 double startupDuration; ///<- Specifies the time it takes for the gust startup transient.
227 double steadyDuration; ///<- Specifies the duration of the steady gust.
228 double endDuration; ///<- Specifies the time it takes for the gust to subsude.
229 OneMinusCosineProfile() ///<- The constructor.
239 enum eGustFrame {gfNone=0, gfBody, gfWind, gfLocal};
241 /// Stores the information about a single one minus cosine gust instance.
242 struct OneMinusCosineGust {
243 FGColumnVector3 vWind; ///<- The input normalized wind vector.
244 FGColumnVector3 vWindTransformed; ///<- The transformed normal vector at the time the gust is started.
245 double magnitude; ///<- The magnitude of the wind vector.
246 eGustFrame gustFrame; ///<- The frame that the wind vector is specified in.
247 struct OneMinusCosineProfile gustProfile; ///<- The gust shape (profile) data for this gust.
248 OneMinusCosineGust() ///<- Constructor.
250 vWind.InitMatrix(0.0);
256 /// Stores information about a specified Up- or Down-burst.
258 double ringLatitude; ///<- The latitude of the downburst run (radians)
259 double ringLongitude; ///<- The longitude of the downburst run (radians)
260 double ringAltitude; ///<- The altitude of the ring (feet).
261 double ringRadius; ///<- The radius of the ring (feet).
262 double ringCoreRadius; ///<- The cross-section "core" radius of the ring (feet).
263 double circulation; ///<- The circulation (gamma) (feet-squared per second).
264 struct OneMinusCosineProfile oneMCosineProfile;
267 // 1 - Cosine gust setters
268 /// Initiates the execution of the gust.
269 virtual void StartGust(bool running) {oneMinusCosineGust.gustProfile.Running = running;}
270 ///Specifies the duration of the startup portion of the gust.
271 virtual void StartupGustDuration(double dur) {oneMinusCosineGust.gustProfile.startupDuration = dur;}
272 ///Specifies the length of time that the gust is at a steady, full strength.
273 virtual void SteadyGustDuration(double dur) {oneMinusCosineGust.gustProfile.steadyDuration = dur;}
274 /// Specifies the length of time it takes for the gust to return to zero velocity.
275 virtual void EndGustDuration(double dur) {oneMinusCosineGust.gustProfile.endDuration = dur;}
276 /// Specifies the magnitude of the gust in feet/second.
277 virtual void GustMagnitude(double mag) {oneMinusCosineGust.magnitude = mag;}
278 /** Specifies the frame that the gust direction vector components are specified in. The
279 body frame is defined with the X direction forward, and the Y direction positive out
280 the right wing. The wind frame is defined with the X axis pointing into the velocity
281 vector, the Z axis perpendicular to the X axis, in the aircraft XZ plane, and the Y
282 axis completing the system. The local axis is a navigational frame with X pointing north,
283 Y pointing east, and Z pointing down. This is a locally vertical, locally horizontal
284 frame, with the XY plane tangent to the geocentric surface. */
285 virtual void GustFrame(eGustFrame gFrame) {oneMinusCosineGust.gustFrame = gFrame;}
286 /// Specifies the X component of velocity in the specified gust frame (ft/sec).
287 virtual void GustXComponent(double x) {oneMinusCosineGust.vWind(eX) = x;}
288 /// Specifies the Y component of velocity in the specified gust frame (ft/sec).
289 virtual void GustYComponent(double y) {oneMinusCosineGust.vWind(eY) = y;}
290 /// Specifies the Z component of velocity in the specified gust frame (ft/sec).
291 virtual void GustZComponent(double z) {oneMinusCosineGust.vWind(eZ) = z;}
305 double MagnitudedAccelDt, MagnitudeAccel, Magnitude;
310 double wind_from_clockwise;
311 double spike, target_time, strength;
312 FGColumnVector3 vDirectiondAccelDt;
313 FGColumnVector3 vDirectionAccel;
314 FGColumnVector3 vDirection;
315 FGColumnVector3 vTurbulenceGrad;
316 FGColumnVector3 vBodyTurbGrad;
317 FGColumnVector3 vTurbPQR;
319 struct OneMinusCosineGust oneMinusCosineGust;
321 // Dryden turbulence model
322 double windspeed_at_20ft; ///< in ft/s
323 int probability_of_exceedence_index; ///< this is bound as the severity property
324 FGTable *POE_Table; ///< probability of exceedence table
327 FGColumnVector3 vTotalWindNED;
328 FGColumnVector3 vWindNED;
329 FGColumnVector3 vGustNED;
330 FGColumnVector3 vCosineGust;
331 FGColumnVector3 vBurstGust;
332 FGColumnVector3 vTurbulenceNED;
334 /// Get T, P and rho for a standard atmosphere at the given altitude.
335 void Turbulence(double h);
338 double CosineGustProfile( double startDuration, double steadyDuration,
339 double endDuration, double elapsedTime);
341 virtual void bind(void);
342 void Debug(int from);
345 } // namespace JSBSim
347 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%