7 // Start in a "sane" mode, so unset stuff doesn't freak us out
11 _peaks[0] = _peaks[1] = 1;
15 _widths[i] = 0.01; // half a degree
17 _orient[0] = 1; _orient[1] = 0; _orient[2] = 0;
18 _orient[3] = 0; _orient[4] = 1; _orient[5] = 0;
19 _orient[6] = 0; _orient[7] = 0; _orient[8] = 1;
23 _slatPos = _spoilerPos = _flapPos = 0;
24 _slatDrag = _spoilerDrag = _flapDrag = 1;
31 void Surface::setPosition(float* p)
34 for(i=0; i<3; i++) _pos[i] = p[i];
37 void Surface::getPosition(float* out)
40 for(i=0; i<3; i++) out[i] = _pos[i];
43 void Surface::setChord(float chord)
48 void Surface::setTotalDrag(float c0)
53 float Surface::getTotalDrag()
58 void Surface::setXDrag(float cx)
63 void Surface::setYDrag(float cy)
68 void Surface::setZDrag(float cz)
73 void Surface::setBaseZDrag(float cz0)
78 void Surface::setStallPeak(int i, float peak)
83 void Surface::setStall(int i, float alpha)
88 void Surface::setStallWidth(int i, float width)
93 void Surface::setOrientation(float* o)
100 void Surface::setIncidence(float angle)
105 void Surface::setSlatParams(float stallDelta, float dragPenalty)
107 _slatAlpha = stallDelta;
108 _slatDrag = dragPenalty;
111 void Surface::setFlapParams(float liftAdd, float dragPenalty)
114 _flapDrag = dragPenalty;
117 void Surface::setSpoilerParams(float liftPenalty, float dragPenalty)
119 _spoilerLift = liftPenalty;
120 _spoilerDrag = dragPenalty;
123 void Surface::setFlap(float pos)
128 void Surface::setSlat(float pos)
133 void Surface::setSpoiler(float pos)
138 // Calculate the aerodynamic force given a wind vector v (in the
139 // aircraft's "local" coordinates) and an air density rho. Returns a
140 // torque about the Y axis, too.
141 void Surface::calcForce(float* v, float rho, float* out, float* torque)
143 // Split v into magnitude and direction:
144 float vel = Math::mag3(v);
146 // Handle the blowup condition. Zero velocity means zero force by
150 for(i=0; i<3; i++) out[i] = torque[i] = 0;
154 Math::mul3(1/vel, v, out);
156 // Convert to the surface's coordinates
157 Math::vmul33(_orient, out, out);
159 // "Rotate" by the incidence angle. Assume small angles, so we
160 // need to diddle only the Z component, X is relatively unchanged
161 // by small rotations.
162 out[2] += _incidence * out[0]; // z' = z + incidence * x
164 // Diddle the Z force according to our configuration
165 float stallMul = stallFunc(out);
166 stallMul *= 1 + _spoilerPos * (_spoilerLift - 1);
167 float stallLift = (stallMul - 1) * _cz * out[2];
168 float flaplift = flapLift(out[2]);
170 out[2] *= _cz; // scaling factor
171 out[2] += _cz*_cz0; // zero-alpha lift
175 // Airfoil lift (pre-stall and zero-alpha) torques "up" (negative
176 // torque) around the Y axis, while flap lift pushes down. Both
177 // forces are considered to act at one third chord from the
178 // edge. Convert to local (i.e. airplane) coordiantes and store
181 torque[1] = 0.1667f * _chord * (flaplift - (_cz*_cz0 + stallLift));
183 Math::tmul33(_orient, torque, torque);
185 // The X (drag) force gets diddled for control deflection
186 out[0] = controlDrag(out[2], _cx * out[0]);
188 // Add in any specific Y (side force) coefficient.
191 // Reverse the incidence rotation to get back to surface
193 out[2] -= _incidence * out[0];
195 // Convert back to external coordinates
196 Math::tmul33(_orient, out, out);
198 // Add in the units to make a real force:
199 float scale = 0.5f*rho*vel*vel*_c0;
200 Math::mul3(scale, out, out);
201 Math::mul3(scale, torque, torque);
204 // Returns a multiplier for the "plain" force equations that
205 // approximates an airfoil's lift/stall curve.
206 float Surface::stallFunc(float* v)
208 // Sanity check to treat FPU psychopathology
209 if(v[0] == 0) return 1;
211 float alpha = Math::abs(v[2]/v[0]);
213 // Wacky use of indexing, see setStall*() methods.
214 int fwdBak = v[0] > 0; // set if this is "backward motion"
215 int posNeg = v[2] < 0; // set if the lift is toward -z
216 int i = (fwdBak<<1) | posNeg;
218 float stallAlpha = _stalls[i];
223 stallAlpha += _slatAlpha;
226 if(alpha > stallAlpha+_widths[i])
229 // (note mask: we want to use the "positive" stall angle here)
230 float scale = 0.5f*_peaks[fwdBak]/_stalls[i&2];
233 if(alpha <= stallAlpha)
236 // Inside the stall. Compute a cubic interpolation between the
237 // pre-stall "scale" value and the post-stall unity.
238 float frac = (alpha - stallAlpha) / _widths[i];
239 frac = frac*frac*(3-2*frac);
241 return scale*(1-frac) + frac;
244 // Similar to the above -- interpolates out the flap lift past the
246 float Surface::flapLift(float alpha)
248 float flapLift = _cz * _flapPos * (_flapLift-1);
250 if(alpha < 0) alpha = -alpha;
251 if(alpha < _stalls[0])
253 else if(alpha > _stalls[0] + _widths[0])
256 float frac = (alpha - _stalls[0]) / _widths[0];
257 frac = frac*frac*(3-2*frac);
258 return flapLift * (1-frac) + frac;
261 float Surface::controlDrag(float lift, float drag)
263 // Negative flap deflections don't affect drag until their lift
264 // multiplier exceeds the "camber" (cz0) of the surface. Use a
265 // synthesized "fp" number instead of the actual flap position.
269 fp -= _cz0/(_flapLift-1);
273 // Calculate an "effective" drag -- this is the drag that would
274 // have been produced by an unflapped surface at the same lift.
275 float flapDragAoA = (_flapLift - 1 - _cz0) * _stalls[0];
276 float fd = Math::abs(lift * flapDragAoA * fp);
277 if(drag < 0) fd = -fd;
280 // Now multiply by the various control factors
281 drag *= 1 + fp * (_flapDrag - 1);
282 drag *= 1 + _spoilerPos * (_spoilerDrag - 1);
283 drag *= 1 + _slatPos * (_slatDrag - 1);
288 }; // namespace yasim