5 Surface::Surface( Version * version ) :
8 // Start in a "sane" mode, so unset stuff doesn't freak us out
12 _peaks[0] = _peaks[1] = 1;
16 _widths[i] = 0.01; // half a degree
18 _orient[0] = 1; _orient[1] = 0; _orient[2] = 0;
19 _orient[3] = 0; _orient[4] = 1; _orient[5] = 0;
20 _orient[6] = 0; _orient[7] = 0; _orient[8] = 1;
25 _slatPos = _spoilerPos = _flapPos = 0;
26 _slatDrag = _spoilerDrag = _flapDrag = 1;
29 _flapEffectiveness = 1;
35 void Surface::setPosition(float* p)
38 for(i=0; i<3; i++) _pos[i] = p[i];
41 void Surface::getPosition(float* out)
44 for(i=0; i<3; i++) out[i] = _pos[i];
47 void Surface::setChord(float chord)
52 void Surface::setTotalDrag(float c0)
57 float Surface::getTotalDrag()
62 void Surface::setXDrag(float cx)
67 void Surface::setYDrag(float cy)
72 void Surface::setZDrag(float cz)
77 float Surface::getXDrag()
82 void Surface::setBaseZDrag(float cz0)
87 void Surface::setStallPeak(int i, float peak)
92 void Surface::setStall(int i, float alpha)
97 void Surface::setStallWidth(int i, float width)
102 void Surface::setOrientation(float* o)
109 void Surface::setIncidence(float angle)
114 void Surface::setTwist(float angle)
119 void Surface::setSlatParams(float stallDelta, float dragPenalty)
121 _slatAlpha = stallDelta;
122 _slatDrag = dragPenalty;
125 void Surface::setFlapParams(float liftAdd, float dragPenalty)
128 _flapDrag = dragPenalty;
131 void Surface::setSpoilerParams(float liftPenalty, float dragPenalty)
133 _spoilerLift = liftPenalty;
134 _spoilerDrag = dragPenalty;
137 void Surface::setFlap(float pos)
142 void Surface::setFlapEffectiveness(float effectiveness)
144 _flapEffectiveness = effectiveness;
147 double Surface::getFlapEffectiveness()
149 return _flapEffectiveness;
153 void Surface::setSlat(float pos)
158 void Surface::setSpoiler(float pos)
163 // Calculate the aerodynamic force given a wind vector v (in the
164 // aircraft's "local" coordinates) and an air density rho. Returns a
165 // torque about the Y axis, too.
166 void Surface::calcForce(float* v, float rho, float* out, float* torque)
168 // Split v into magnitude and direction:
169 float vel = Math::mag3(v);
171 // Handle the blowup condition. Zero velocity means zero force by
175 for(i=0; i<3; i++) out[i] = torque[i] = 0;
179 // special case this so the logic below doesn't produce a non-zero
180 // force; should probably have a "no force" flag instead...
181 if(_cx == 0. && _cy == 0. && _cz == 0.) {
182 for(int i=0; i<3; i++) out[i] = torque[i] = 0.;
186 Math::mul3(1/vel, v, out);
188 // Convert to the surface's coordinates
189 Math::vmul33(_orient, out, out);
191 // "Rotate" by the incidence angle. Assume small angles, so we
192 // need to diddle only the Z component, X is relatively unchanged
193 // by small rotations.
194 float incidence = _incidence + _twist;
195 out[2] += incidence * out[0]; // z' = z + incidence * x
197 // Hold onto the local wind vector so we can multiply the induced
200 Math::set3(out, lwind);
202 // Diddle the Z force according to our configuration
203 float stallMul = stallFunc(out);
204 stallMul *= 1 + _spoilerPos * (_spoilerLift - 1);
205 float stallLift = (stallMul - 1) * _cz * out[2];
206 float flaplift = flapLift(out[2]);
208 out[2] *= _cz; // scaling factor
209 out[2] += _cz*_cz0; // zero-alpha lift
213 // Airfoil lift (pre-stall and zero-alpha) torques "up" (negative
214 // torque) around the Y axis, while flap lift pushes down. Both
215 // forces are considered to act at one third chord from the
216 // edge. Convert to local (i.e. airplane) coordiantes and store
219 torque[1] = 0.1667f * _chord * (flaplift - (_cz*_cz0 + stallLift));
221 Math::tmul33(_orient, torque, torque);
223 // The X (drag) force gets diddled for control deflection
224 out[0] = controlDrag(out[2], _cx * out[0]);
226 // Add in any specific Y (side force) coefficient.
229 // Diddle the induced drag
230 Math::mul3(-1*_inducedDrag*out[2]*lwind[2], lwind, lwind);
231 Math::add3(lwind, out, out);
233 // Reverse the incidence rotation to get back to surface
234 // coordinates. Since out[] is now the force vector and is
235 // roughly parallel with Z, the small-angle approximation
236 // must change its X component.
237 if( _version->isVersionOrNewer( Version::YASIM_VERSION_32 )) {
238 out[0] += incidence * out[2];
240 out[2] -= incidence * out[0];
243 // Convert back to external coordinates
244 Math::tmul33(_orient, out, out);
246 // Add in the units to make a real force:
247 float scale = 0.5f*rho*vel*vel*_c0;
248 Math::mul3(scale, out, out);
249 Math::mul3(scale, torque, torque);
255 static const float DEG2RAD = 0.0174532925199;
256 float v[3], force[3], torque[3];
257 float rho = Atmosphere::getStdDensity(0);
264 for(float angle = -90; angle<90; angle += 0.01) {
265 float rad = angle * DEG2RAD;
266 v[0] = spd * -Math::cos(rad);
268 v[2] = spd * Math::sin(rad);
269 calcForce(v, rho, force, torque);
270 float lift = force[2] * Math::cos(rad) + force[0] * Math::sin(rad);
271 //__builtin_printf("%f %f\n", angle, lift);
272 __builtin_printf("%f %f\n", angle, torque[2]);
277 // Returns a multiplier for the "plain" force equations that
278 // approximates an airfoil's lift/stall curve.
279 float Surface::stallFunc(float* v)
281 // Sanity check to treat FPU psychopathology
282 if(v[0] == 0) return 1;
284 float alpha = Math::abs(v[2]/v[0]);
286 // Wacky use of indexing, see setStall*() methods.
287 int fwdBak = v[0] > 0; // set if this is "backward motion"
288 int posNeg = v[2] < 0; // set if the airflow is toward -z
289 int i = (fwdBak<<1) | posNeg;
291 float stallAlpha = _stalls[i];
296 if( _version->isVersionOrNewer( Version::YASIM_VERSION_32 )) {
297 stallAlpha += _slatPos * _slatAlpha;
299 stallAlpha += _slatAlpha;
304 if(alpha > stallAlpha+_widths[i])
307 // (note mask: we want to use the "positive" stall angle here)
308 float scale = 0.5f*_peaks[fwdBak]/_stalls[i&2];
311 if(alpha <= stallAlpha)
314 // Inside the stall. Compute a cubic interpolation between the
315 // pre-stall "scale" value and the post-stall unity.
316 float frac = (alpha - stallAlpha) / _widths[i];
317 frac = frac*frac*(3-2*frac);
319 return scale*(1-frac) + frac;
322 // Similar to the above -- interpolates out the flap lift past the
324 float Surface::flapLift(float alpha)
326 float flapLift = _cz * _flapPos * (_flapLift-1) * _flapEffectiveness;
331 if(alpha < 0) alpha = -alpha;
332 if(alpha < _stalls[0])
334 else if(alpha > _stalls[0] + _widths[0])
337 float frac = (alpha - _stalls[0]) / _widths[0];
338 frac = frac*frac*(3-2*frac);
339 return flapLift * (1-frac);
342 float Surface::controlDrag(float lift, float drag)
344 // Negative flap deflections don't affect drag until their lift
345 // multiplier exceeds the "camber" (cz0) of the surface. Use a
346 // synthesized "fp" number instead of the actual flap position.
350 fp -= _cz0/(_flapLift-1);
354 // Calculate an "effective" drag -- this is the drag that would
355 // have been produced by an unflapped surface at the same lift.
356 float flapDragAoA = (_flapLift - 1 - _cz0) * _stalls[0];
357 float fd = Math::abs(lift * flapDragAoA * fp);
358 if(drag < 0) fd = -fd;
361 // Now multiply by the various control factors
362 drag *= 1 + fp * (_flapDrag - 1);
363 drag *= 1 + _spoilerPos * (_spoilerDrag - 1);
364 drag *= 1 + _slatPos * (_slatDrag - 1);
369 }; // namespace yasim