MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
aerodc81.c File Reference
#include "mbconfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ac/f2c.h"
#include "aerodc81.h"
#include "bisec.h"
#include "c81data.h"
Include dependency graph for aerodc81.c:

Go to the source code of this file.

Functions

static int get_coef (int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach, doublereal *c, doublereal *c0)
 
static doublereal get_dcla (int nm, doublereal *m, doublereal *s, doublereal mach)
 
doublereal c81_data_get_coef (int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach)
 
int c81_aerod2 (doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data)
 
int c81_aerod2_u (doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data, long unsteadyflag)
 

Variables

const outa_t outa_Zero
 

Function Documentation

int c81_aerod2 ( doublereal W,
const vam_t VAM,
doublereal TNG,
outa_t OUTA,
c81_data data 
)

Definition at line 100 of file aerodc81.c.

References c81_data::ad, c81_data::al, outa_t::alpha, c81_data::am, grad::atan2(), vam_t::bc_position, outa_t::cd, vam_t::chord, outa_t::cl, outa_t::clalpha, outa_t::cm, grad::cos(), vam_t::density, grad::fabs(), vam_t::force_position, outa_t::gamma, get_coef(), get_dcla(), M_PI, outa_t::mach, c81_data::md, c81_data::ml, c81_data::mm, c81_data::NAD, c81_data::NAL, c81_data::NAM, c81_data::NMD, c81_data::NML, c81_data::NMM, vam_t::sound_celerity, grad::sqrt(), and c81_data::stall.

101 {
102  /*
103  * velocita' del punto in cui sono calcolate le condizioni al contorno
104  */
105  doublereal v[3];
106  doublereal vp, vp2, vtot;
107  doublereal rho = VAM->density;
108  doublereal cs = VAM->sound_celerity;
109  doublereal chord = VAM->chord;
110 
111  doublereal cl = 0., cl0 = 0., cd = 0., cd0 = 0., cm = 0.;
112  doublereal alpha, gamma, cosgam, mach, q;
113  doublereal dcla;
114 
115  doublereal ca = VAM->force_position;
116  doublereal c34 = VAM->bc_position;
117 
118  const doublereal RAD2DEG = 180.*M_1_PI;
119  const doublereal M_PI_3 = M_PI/3.;
120 
121  enum { V_X = 0, V_Y = 1, V_Z = 2, W_X = 3, W_Y = 4, W_Z = 5 };
122 
123  /*
124  * porta la velocita' al punto di calcolo delle boundary conditions
125  */
126  v[V_X] = W[V_X];
127  v[V_Y] = W[V_Y] + c34*W[W_Z];
128  v[V_Z] = W[V_Z] - c34*W[W_Y];
129 
130  vp2 = v[V_X]*v[V_X] + v[V_Y]*v[V_Y];
131  vp = sqrt(vp2);
132 
133  vtot = sqrt(vp2 + v[V_Z]*v[V_Z]);
134 
135  /*
136  * non considera velocita' al di sotto di 1.e-6
137  * FIXME: rendere parametrico?
138  */
139  if (vp/cs < 1.e-6) {
140  TNG[V_X] = 0.;
141  TNG[V_Y] = 0.;
142  TNG[V_Z] = 0.;
143  TNG[W_X] = 0.;
144  TNG[W_Y] = 0.;
145  TNG[W_Z] = 0.;
146 
147  OUTA->alpha = 0.;
148  OUTA->gamma = 0.;
149  OUTA->mach = 0.;
150  OUTA->cl = 0.;
151  OUTA->cd = 0.;
152  OUTA->cm = 0.;
153  OUTA->clalpha = 0.;
154 
155  return 0;
156  }
157 
158  /*
159  * FIXME: gestire quali angoli indicano cosa
160  *
161  * Idea di base: capire in base ad un criterio di linearita'
162  * a quale angolo di incidenza il profilo stalla, quindi
163  * oltre quell'angolo prendere la correzione per flusso trasverso
164  * in un certo modo, altrimenti in un altro ?!?
165  */
166  alpha = atan2(-v[V_Y], v[V_X]);
167  OUTA->alpha = alpha*RAD2DEG;
168  gamma = atan2(-v[V_Z], fabs(v[V_X])); /* come in COE0 (aerod2.f) */
169  /* gamma = atan2(-v[V_Z], vp); */ /* secondo me (?!?) */
170  OUTA->gamma = gamma*RAD2DEG;
171 
172  if (fabs(gamma) > M_PI_3) {
173  /* tanto ne viene preso il coseno ... */
174  gamma = M_PI_3;
175  }
176 
177  cosgam = cos(gamma);
178  mach = (vtot*sqrt(cosgam))/cs;
179  OUTA->mach = mach;
180 
181  /*
182  * Note: all angles in c81 files MUST be in degrees
183  */
184  get_coef(data->NML, data->ml, data->NAL, data->al,
185  OUTA->alpha, mach, &cl, &cl0);
186  get_coef(data->NMD, data->md, data->NAD, data->ad,
187  OUTA->alpha, mach, &cd, &cd0);
188  get_coef(data->NMM, data->mm, data->NAM, data->am,
189  OUTA->alpha, mach, &cm, NULL);
190 
191  dcla = get_dcla(data->NML, data->ml, data->stall, mach);
192 
193 /*
194  * da COE0 (aerod2.f):
195  *
196  ASLRF = ASLOP0
197  IF(DABS(ALFA).LT.1.D-6) GOTO 10
198  ASLRF = CLIFT/(ALFA*COSGAM)
199  IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0
200  10 CLIFT = ASLRF*ALFA
201  *
202  */
203 
204  /*
205  * in soldoni: se si e' oltre il tratto lineare, prende la
206  * secante con l'angolo corretto per la freccia (fino a 60
207  * gradi) e poi ricalcola il cl con l'angolo vero; in questo
208  * modo il cl viene piu' grande di circa 1/cos(gamma) ma solo
209  * fuori del tratto lineare.
210  */
211  dcla *= RAD2DEG;
212  if (fabs(alpha) > 1.e-6) {
213  doublereal dclatmp = (cl - cl0)/(alpha*cosgam);
214  if (dclatmp < dcla) {
215  dcla = dclatmp;
216  }
217  }
218  cl = cl0 + dcla*alpha;
219 
220  OUTA->cl = cl;
221  OUTA->cd = cd;
222  OUTA->cm = cm;
223  OUTA->clalpha = dcla;
224 
225  q = .5*rho*chord*vp2;
226 
227  TNG[V_X] = -q*(cl*v[V_Y] + cd*v[V_X])/vp;
228  TNG[V_Y] = q*(cl*v[V_X] - cd*v[V_Y])/vp;
229  TNG[V_Z] = -q*cd0*v[V_Z]/vp;
230  TNG[W_X] = 0.;
231  TNG[W_Y] = -ca*TNG[V_Z];
232  TNG[W_Z] = q*chord*cm + ca*TNG[V_Y];
233 
234  /*
235  * Radial drag (TNG[V_Z]) consistent with Harris, JAHS 1970
236  * and with CAMRAD strip theory section forces
237  */
238 
239  return 0;
240 }
doublereal cd
Definition: aerodc81.h:66
doublereal density
Definition: aerodc81.h:108
int NMD
Definition: aerodc81.h:160
#define M_PI
Definition: gradienttest.cc:67
int NAM
Definition: aerodc81.h:166
doublereal sound_celerity
Definition: aerodc81.h:109
int NMM
Definition: aerodc81.h:165
doublereal gamma
Definition: aerodc81.h:63
int NML
Definition: aerodc81.h:146
doublereal bc_position
Definition: aerodc81.h:112
doublereal * mm
Definition: aerodc81.h:167
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
doublereal alpha
Definition: aerodc81.h:62
int NAL
Definition: aerodc81.h:147
static int get_coef(int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach, doublereal *c, doublereal *c0)
Definition: aerodc81.c:637
doublereal mach
Definition: aerodc81.h:64
doublereal force_position
Definition: aerodc81.h:111
doublereal * ad
Definition: aerodc81.h:163
doublereal chord
Definition: aerodc81.h:110
doublereal cl
Definition: aerodc81.h:65
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
doublereal clalpha
Definition: aerodc81.h:69
doublereal * stall
Definition: aerodc81.h:157
doublereal * am
Definition: aerodc81.h:168
GradientExpression< UnaryExpr< FuncCos, Expr > > cos(const GradientExpression< Expr > &u)
Definition: gradient.h:2978
doublereal * al
Definition: aerodc81.h:149
doublereal cm
Definition: aerodc81.h:67
static doublereal get_dcla(int nm, doublereal *m, doublereal *s, doublereal mach)
Definition: aerodc81.c:770
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
doublereal * md
Definition: aerodc81.h:162
doublereal * ml
Definition: aerodc81.h:148
int NAD
Definition: aerodc81.h:161

Here is the call graph for this function:

int c81_aerod2_u ( doublereal W,
const vam_t VAM,
doublereal TNG,
outa_t OUTA,
c81_data data,
long  unsteadyflag 
)

Definition at line 243 of file aerodc81.c.

References c81_data::ad, c81_data::al, outa_t::alf1, outa_t::alf2, outa_t::alpha, c81_data::am, grad::atan2(), vam_t::bc_position, outa_t::cd, vam_t::chord, outa_t::cl, outa_t::clalpha, outa_t::cm, grad::cos(), outa_t::dam, outa_t::dan, outa_t::dcm, outa_t::dcn, vam_t::density, grad::exp(), grad::fabs(), vam_t::force_position, outa_t::gamma, get_coef(), get_dcla(), M_PI, outa_t::mach, c81_data::md, c81_data::ml, c81_data::mm, c81_data::mstall, c81_data::NAD, c81_data::NAL, c81_data::NAM, c81_data::NMD, c81_data::NML, c81_data::NMM, grad::pow(), vam_t::sound_celerity, grad::sqrt(), and c81_data::stall.

Referenced by TheodorsenAeroData::AssJac(), C81AeroData::GetForces(), C81MultipleAeroData::GetForces(), and C81InterpolatedAeroData::GetForces().

245 {
246  /*
247  * velocita' del punto in cui sono calcolate le condizioni al contorno
248  */
249  doublereal v[3];
250  doublereal vp, vp2, vtot;
251  doublereal rho = VAM->density;
252  doublereal cs = VAM->sound_celerity;
253  doublereal chord = VAM->chord;
254 
255  doublereal cl = 0., cl0 = 0., cd = 0., cd0 = 0., cm = 0.;
256  doublereal alpha, gamma, cosgam, mach, q;
257  doublereal dcla;
258 
259  doublereal ca = VAM->force_position;
260  doublereal c34 = VAM->bc_position;
261 
262  const doublereal RAD2DEG = 180.*M_1_PI;
263  const doublereal M_PI_3 = M_PI/3.;
264 
265  enum { V_X = 0, V_Y = 1, V_Z = 2, W_X = 3, W_Y = 4, W_Z = 5 };
266 
267  /*
268  * porta la velocita' al punto di calcolo delle boundary conditions
269  */
270  v[V_X] = W[V_X];
271  v[V_Y] = W[V_Y] + c34*W[W_Z];
272  v[V_Z] = W[V_Z] - c34*W[W_Y];
273 
274  vp2 = v[V_X]*v[V_X] + v[V_Y]*v[V_Y];
275  vp = sqrt(vp2);
276 
277  vtot = sqrt(vp2 + v[V_Z]*v[V_Z]);
278 
279  /*
280  * non considera velocita' al di sotto di 1.e-6
281  * FIXME: rendere parametrico?
282  */
283 
284  if (vp/cs < 1.e-6) {
285  TNG[V_X] = 0.;
286  TNG[V_Y] = 0.;
287  TNG[V_Z] = 0.;
288  TNG[W_X] = 0.;
289  TNG[W_Y] = 0.;
290  TNG[W_Z] = 0.;
291 
292  OUTA->alpha = 0.;
293  OUTA->gamma = 0.;
294  OUTA->mach = 0.;
295  OUTA->cl = 0.;
296  OUTA->cd = 0.;
297  OUTA->cm = 0.;
298  OUTA->clalpha = 0.;
299 
300  return 0;
301  }
302 
303  /*
304  * FIXME: gestire quali angoli indicano cosa
305  *
306  * Idea di base: capire in base ad un criterio di linearita'
307  * a quale angolo di incidenza il profilo stalla, quindi
308  * oltre quell'angolo prendere la correzione per flusso trasverso
309  * in un certo modo, altrimenti in un altro ?!?
310  */
311  alpha = atan2(-v[V_Y], v[V_X]);
312  OUTA->alpha = alpha*RAD2DEG;
313  gamma = atan2(-v[V_Z], fabs(v[V_X])); /* come in COE0 (aerod2.f) */
314  /* gamma = atan2(-v[V_Z], vp); */ /* secondo me (?!?) */
315  OUTA->gamma = gamma*RAD2DEG;
316 
317  if (fabs(gamma) > M_PI_3) {
318  /* tanto ne viene preso il coseno ... */
319  gamma = M_PI_3;
320  }
321 
322  cosgam = cos(gamma);
323  mach = (vtot*sqrt(cosgam))/cs;
324  OUTA->mach = mach;
325 
326  /*
327  * mach cannot be more than .99 (see aerod.f)
328  */
329  if (mach > .99) {
330  mach = .99;
331  }
332 
333  /*
334  * Compute cl, cd, cm based on selected theory
335  */
336  switch (unsteadyflag) {
337  case 0:
338 
339  /*
340  * Note: all angles in c81 files MUST be in degrees
341  */
342  get_coef(data->NML, data->ml, data->NAL, data->al,
343  OUTA->alpha, mach, &cl, &cl0);
344  get_coef(data->NMD, data->md, data->NAD, data->ad,
345  OUTA->alpha, mach, &cd, &cd0);
346  get_coef(data->NMM, data->mm, data->NAM, data->am,
347  OUTA->alpha, mach, &cm, NULL);
348 
349  dcla = get_dcla(data->NML, data->ml, data->stall, mach);
350 
351 /*
352  * da COE0 (aerod2.f):
353  *
354  ASLRF = ASLOP0
355  IF(DABS(ALFA).LT.1.D-6) GOTO 10
356  ASLRF = CLIFT/(ALFA*COSGAM)
357  IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0
358  10 CLIFT = ASLRF*ALFA
359  *
360  */
361 
362  /*
363  * in soldoni: se si e' oltre il tratto lineare, prende la
364  * secante con l'angolo corretto per la freccia (fino a 60
365  * gradi) e poi ricalcola il cl con l'angolo vero; in questo
366  * modo il cl viene piu' grande di circa 1/cos(gamma) ma solo
367  * fuori del tratto lineare.
368  */
369  dcla *= RAD2DEG;
370  if (fabs(alpha) > 1.e-6) {
371  doublereal dclatmp = (cl - cl0)/(alpha*cosgam);
372  if (dclatmp < dcla) {
373  dcla = dclatmp;
374  cl = cl0 + dcla*alpha;
375  }
376  }
377  break;
378 
379  case 1:
380  return -1;
381 
382  case 2: {
383 
384  /*
385  * Constants from unsteady theory
386  * synthetized by Richard L. Bielawa,
387  * 31th A.H.S. Forum Washington D.C.
388  * May 1975
389  */
390  doublereal A, B, A2, B2, ETA, ASN, ASM,
391  SGN, SGM, SGMAX,
392  DAN, DCN, DAM, DCM,
393  S2, alphaN, alphaM, C1,
394  dcma, dclatan, ALF1, ALF2,
395  cn;
396 
397  const doublereal PN[] = {
398  -3.464003e-1,
399  -1.549076e+0,
400  4.306330e+1,
401  -5.397529e+1,
402  5.781402e+0,
403  -3.233003e+1,
404  -2.162257e+1,
405  1.866347e+1,
406  4.198390e+1,
407  3.295461e+2,
408  };
409 
410  const doublereal QN[] = {
411  1.533717e+0,
412  6.977203e+0,
413  1.749010e+3,
414  1.694829e+3,
415  -1.771899e+3,
416  -3.291665e+4,
417  2.969051e+0,
418  -3.632448e+1,
419  -2.268578e+3,
420  6.601995e+3,
421  -9.654208e+3,
422  8.533930e+4,
423  -1.492624e+0,
424  1.163661e+1
425  };
426 
427  const doublereal PM[] = {
428  1.970065e+1,
429  -6.751639e+1,
430  7.265269e+2,
431  4.865945e+4,
432  2.086279e+4,
433  6.024672e+3,
434  1.446334e+2,
435  8.586896e+2,
436  -7.550329e+2,
437  -1.021613e+1,
438  2.247664e+1,
439  };
440 
441  const doublereal QM[] = {
442  -2.322808e+0,
443  -1.322257e+0,
444  -2.633891e+0,
445  -2.180321e-1,
446  4.580014e+0,
447  3.125497e-1,
448  -2.828806e+1,
449  -4.396734e+0,
450  2.565870e+2,
451  -1.204976e+1,
452  -1.157802e+2,
453  8.612138e+0,
454  };
455 
456  enum {
457  U_1 = 0,
458  U_2 = 1,
459  U_3 = 2,
460  U_4 = 3,
461  U_5 = 4,
462  U_6 = 5,
463  U_7 = 6,
464  U_8 = 7,
465  U_9 = 8,
466  U10 = 9,
467  U11 = 10,
468  U12 = 11,
469  U13 = 12,
470  U14 = 13
471  };
472 
473  /*
474  * This is the static stall angle for Mach = 0
475  * (here a symmetric airfoil is assumed; the real
476  * _SIGNED_ static stall should be considered ...)
477  */
478  const doublereal ASN0 = .22689, ASM0 = .22689;
479 
480  ALF1 = OUTA->alf1;
481  ALF2 = OUTA->alf2;
482 
483  A = .5*chord*ALF1/vp;
484  B = .25*chord*chord*ALF2/vp2;
485 
486  ETA = sqrt(pow(A/.048, 2) + pow(B/.016, 2));
487 
488  if (alpha < 0.) {
489  A = -A;
490  B = -B;
491  }
492 
493  if (ETA > 1.) {
494  A /= ETA;
495  B /= ETA;
496  }
497 
498  A2 = A*A;
499  B2 = B*B;
500 
501  ASN = ASN0*(1. - mach);
502  ASM = ASM0*(1. - mach);
503 
504  SGN = fabs(alpha/ASN);
505  SGM = fabs(alpha/ASM);
506  SGMAX = 1.839 - 70.33*fabs(B);
507  if (SGMAX > 1.86) {
508  SGMAX = 1.86;
509  }
510  if (SGN > SGMAX) {
511  SGN = SGMAX;
512  }
513  if (SGM > SGMAX) {
514  SGM = SGMAX;
515  }
516 
517  DAN = (A*(PN[U_1] + PN[U_5]*SGN) + B*(PN[U_2] + PN[U_6]*SGN)
518  + exp(-1072.52*A2)*(A*(PN[U_3] + PN[U_7]*SGN)
519  +A2*(PN[U_9] + PN[U10]*SGN))
520  + exp(-40316.42*B2)*B*(PN[U_4] + PN[U_8]*SGN))*ASN;
521  DCN = A*(QN[U_1] + QN[U_3]*A2 + SGN*(QN[U_7] + QN[U_9]*A2 + QN[U13]*SGN)
522  + B2*(QN[U_5] + QN[U11]*SGN))
523  + B*(QN[U_2] + QN[U_4]*A2
524  + SGN*(QN[U_8] + QN[U10]*A2 + QN[U14]*SGN)
525  + B2*(QN[U_6] + QN[U12]*SGN));
526  DAM = (A*(PM[U_1] + PM[U_3]*A2 + PM[U_5]*B2 + PM[U10]*SGM + PM[U_7]*A)
527  + B*(PM[U_2] + PM[U_4]*B2 + PM[U_6]*A2
528  + PM[U11]*SGM + PM[U_8]*B + PM[U_9]*A))*ASM;
529 
530  S2 = SGM*SGM;
531 
532  DCM = A*(QM[U_2] + QM[U_8]*A + SGM*(QM[U_4] + QM[U10]*A)
533  + S2*(QM[U_6] + QM[U12]*A))
534  + B*(QM[U_1] + QM[U_7]*B + SGM*(QM[U_3] + QM[U_9]*B)
535  + S2*(QM[U_5] + QM[U11]*B));
536 
537  OUTA->dan = DAN*RAD2DEG;
538  OUTA->dam = DAM*RAD2DEG;
539  OUTA->dcn = DCN;
540  OUTA->dcm = DCM;
541 
542  /*
543  * I think I need to apply this contribution
544  * with the sign of alpha, because otherwise
545  * it gets discontinuous as alpha changes sign;
546  * I definitely need the original reference :(
547  */
548  if (alpha < 0.) {
549  DAN = -DAN;
550  DCN = -DCN;
551  DAM = -DAM;
552  DCM = -DCM;
553  }
554 
555  alphaN = (alpha - DAN)*RAD2DEG;
556  get_coef(data->NML, data->ml, data->NAL, data->al,
557  alphaN, mach, &cl, &cl0);
558  get_coef(data->NMD, data->md, data->NAD, data->ad,
559  alphaN, mach, &cd, &cd0);
560 
561  alphaM = (alpha - DAM)*RAD2DEG;
562  get_coef(data->NMM, data->mm, data->NAM, data->am,
563  alphaM, mach, &cm, NULL);
564 
565  dcla = get_dcla(data->NML, data->ml, data->stall, mach);
566  dcma = get_dcla(data->NMM, data->mm, data->mstall, mach);
567 
568  /* note: cl/alpha in 1/deg */
569  dclatan = dcla;
570  if (fabs(alphaN) > 1.e-6) {
571  dclatan = (cl - cl0)/(alphaN*cosgam);
572  }
573  cl = cl0 + dclatan*alphaN;
574 
575  /* back to 1/rad */
576  dcla *= RAD2DEG;
577  dcma *= RAD2DEG;
578 
579  C1 = .9457/sqrt(1. - mach*mach);
580 
581  /*
582  * the unsteady correction is "cn",
583  * so split it in "cl" and "cd"
584  * (note: if "vp" is too small the routine exits
585  * earlier without computing forces)
586  */
587  cn = dcla*DAN + DCN*C1;
588  cl += cn*v[V_X]/vp; /* cos(x) */
589  cd -= cn*v[V_Y]/vp; /* sin(x) */
590  cm += dcma*DAM + DCM*C1;
591 
592  break;
593  }
594  }
595 
596  /*
597  * Save cl, cd, cm for output purposes
598  */
599  OUTA->cl = cl;
600  OUTA->cd = cd;
601  OUTA->cm = cm;
602  OUTA->clalpha = dcla;
603 
604  /*
605  * Local dynamic pressure
606  */
607  q = .5*rho*chord*vp2;
608 
609  /*
610  * airfoil forces and moments in the airfoil frame
611  */
612  TNG[V_X] = -q*(cl*v[V_Y] + cd*v[V_X])/vp;
613  TNG[V_Y] = q*(cl*v[V_X] - cd*v[V_Y])/vp;
614  TNG[V_Z] = -q*cd0*v[V_Z]/vp;
615  TNG[W_X] = 0.;
616  TNG[W_Y] = -ca*TNG[V_Z];
617  TNG[W_Z] = q*chord*cm + ca*TNG[V_Y];
618 
619  /*
620  * Radial drag (TNG[V_Z]) consistent with Harris, JAHS 1970
621  * and with CAMRAD strip theory section forces
622  */
623 
624  return 0;
625 }
GradientExpression< UnaryExpr< FuncExp, Expr > > exp(const GradientExpression< Expr > &u)
Definition: gradient.h:2975
doublereal cd
Definition: aerodc81.h:66
doublereal density
Definition: aerodc81.h:108
int NMD
Definition: aerodc81.h:160
doublereal dcn
Definition: aerodc81.h:74
#define M_PI
Definition: gradienttest.cc:67
int NAM
Definition: aerodc81.h:166
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
doublereal sound_celerity
Definition: aerodc81.h:109
int NMM
Definition: aerodc81.h:165
doublereal dcm
Definition: aerodc81.h:75
doublereal gamma
Definition: aerodc81.h:63
int NML
Definition: aerodc81.h:146
doublereal alf1
Definition: aerodc81.h:70
doublereal bc_position
Definition: aerodc81.h:112
doublereal * mm
Definition: aerodc81.h:167
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
doublereal alpha
Definition: aerodc81.h:62
doublereal alf2
Definition: aerodc81.h:71
int NAL
Definition: aerodc81.h:147
static int get_coef(int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach, doublereal *c, doublereal *c0)
Definition: aerodc81.c:637
doublereal mach
Definition: aerodc81.h:64
doublereal dan
Definition: aerodc81.h:72
doublereal force_position
Definition: aerodc81.h:111
doublereal * ad
Definition: aerodc81.h:163
doublereal chord
Definition: aerodc81.h:110
doublereal cl
Definition: aerodc81.h:65
doublereal * mstall
Definition: aerodc81.h:158
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
doublereal clalpha
Definition: aerodc81.h:69
doublereal * stall
Definition: aerodc81.h:157
doublereal * am
Definition: aerodc81.h:168
doublereal dam
Definition: aerodc81.h:73
GradientExpression< UnaryExpr< FuncCos, Expr > > cos(const GradientExpression< Expr > &u)
Definition: gradient.h:2978
doublereal * al
Definition: aerodc81.h:149
doublereal cm
Definition: aerodc81.h:67
static doublereal get_dcla(int nm, doublereal *m, doublereal *s, doublereal mach)
Definition: aerodc81.c:770
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
double doublereal
Definition: colamd.c:52
doublereal * md
Definition: aerodc81.h:162
doublereal * ml
Definition: aerodc81.h:148
int NAD
Definition: aerodc81.h:161

Here is the call graph for this function:

doublereal c81_data_get_coef ( int  nm,
doublereal m,
int  na,
doublereal a,
doublereal  alpha,
doublereal  mach 
)

Definition at line 90 of file aerodc81.c.

References c, and get_coef().

Referenced by main().

91 {
92  doublereal c;
93 
94  get_coef(nm, m, na, a, alpha, mach, &c, NULL);
95 
96  return c;
97 }
static int get_coef(int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach, doublereal *c, doublereal *c0)
Definition: aerodc81.c:637
static std::stack< cleanup * > c
Definition: cleanup.cc:59
static const doublereal a
Definition: hfluid_.h:289
double doublereal
Definition: colamd.c:52

Here is the call graph for this function:

static int get_coef ( int  nm,
doublereal m,
int  na,
doublereal a,
doublereal  alpha,
doublereal  mach,
doublereal c,
doublereal c0 
)
static

Definition at line 637 of file aerodc81.c.

References bisec_d(), and grad::fabs().

Referenced by c81_aerod2(), c81_aerod2_u(), and c81_data_get_coef().

639 {
640  int im;
641  int ia, ia0 = -1;
642 
643  while (alpha < -180.) {
644  alpha += 360.;
645  }
646 
647  while (alpha >= 180.) {
648  alpha -= 360.;
649  }
650 
651  mach = fabs(mach);
652 
653  /*
654  * im e' l'indice di m in cui si trova
655  * l'approssimazione per difetto di mach
656  */
657  im = bisec_d(m, mach, 0, nm - 1);
658 
659  /*
660  * ia e' l'indice della prima colonna di a in cui si trova
661  * l'approssimazione per difetto di alpha
662  */
663  if (c0 != NULL) {
664  ia0 = bisec_d(a, 0., 0, na - 1);
665  }
666 
667  ia = bisec_d(a, alpha, 0, na - 1);
668 
669  if (im == nm - 1) {
670  if (c0 != NULL) {
671  if (ia0 == na - 1) {
672  *c0 = a[na*(nm + 1) - 1];
673 
674  } else if (ia0 == -1) {
675  *c0 = a[na*nm];
676 
677  } else {
678  doublereal da;
679 
680  ia0++;
681  da = -a[ia0 - 1]/(a[ia0] - a[ia0 - 1]);
682  *c0 = (1. - da)*a[na*nm + ia0 - 1] + da*a[na*nm + ia0];
683  }
684  }
685 
686  if (ia == na - 1) {
687  *c = a[na*(nm + 1) - 1];
688 
689  } else if (ia == -1) {
690  *c = a[na*nm];
691 
692  } else {
693  doublereal da;
694 
695  ia++;
696  da = (alpha - a[ia - 1])/(a[ia] - a[ia - 1]);
697  *c = (1. - da)*a[na*nm + ia - 1] + da*a[na*nm + ia];
698  }
699 
700  } else if (im == -1) {
701  if (c0 != NULL) {
702  if (ia0 == na - 1) {
703  *c0 = a[na*2 - 1];
704 
705  } else if (ia0 == -1) {
706  *c0 = a[na];
707 
708  } else {
709  doublereal da;
710 
711  ia0++;
712  da = -a[ia0 - 1]/(a[ia0] - a[ia0 - 1]);
713  *c0 = (1. - da)*a[na + ia0 - 1] + da*a[na + ia0];
714  }
715  }
716 
717  if (ia == na - 1) {
718  *c = a[na*2 - 1];
719 
720  } else if (ia == -1) {
721  *c = a[na];
722 
723  } else {
724  doublereal da;
725 
726  ia++;
727  da = (alpha - a[ia - 1])/(a[ia] - a[ia - 1]);
728  *c = (1. - da)*a[na + ia - 1] + da*a[na + ia];
729  }
730 
731  } else {
732  doublereal d;
733 
734  im++;
735  d = (mach - m[im - 1])/(m[im] - m[im - 1]);
736 
737  if (c0 != NULL) {
738  if (ia0 == na) {
739  *c0 = (1. - d)*a[na*(im + 1) - 1] + d*a[na*(im + 2) - 1];
740  } else {
741  doublereal a1, a2, da;
742  a1 = (1. - d)*a[na*im + ia0 - 1] + d*a[na*(im + 1) + ia0 - 1];
743  a2 = (1. - d)*a[na*im + ia0] + d*a[na*(im + 1) + ia0];
744  da = -a[ia0 - 1]/(a[ia0] - a[ia0 - 1]);
745  *c0 = (1. - da)*a1 + da*a2;
746  }
747  }
748 
749  if (ia == na - 1) {
750  *c = (1. - d)*a[na*(im + 1) - 1] + d*a[na*(im + 2) - 1];
751 
752  } else if (ia == -1) {
753  *c = (1. - d)*a[na*im] + d*a[na*(im + 1)];
754 
755  } else {
756  doublereal a1, a2, da;
757 
758  ia++;
759  a1 = (1. - d)*a[na*im + ia - 1] + d*a[na*(im + 1) + ia - 1];
760  a2 = (1. - d)*a[na*im + ia] + d*a[na*(im + 1) + ia];
761  da = (alpha - a[ia - 1])/(a[ia] - a[ia - 1]);
762  *c = (1. - da)*a1 + da*a2;
763  }
764  }
765 
766  return 0;
767 }
int bisec_d(doublereal *v, doublereal val, int lb, int ub)
Definition: bisec.cc:57
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
static std::stack< cleanup * > c
Definition: cleanup.cc:59
static const doublereal a
Definition: hfluid_.h:289
double doublereal
Definition: colamd.c:52

Here is the call graph for this function:

static doublereal get_dcla ( int  nm,
doublereal m,
doublereal s,
doublereal  mach 
)
static

Definition at line 770 of file aerodc81.c.

References bisec_d(), and grad::fabs().

Referenced by c81_aerod2(), and c81_aerod2_u().

771 {
772  int im;
773 
774  mach = fabs(mach);
775 
776  /*
777  * im e' l'indice di m in cui si trova l'approssimazione per eccesso
778  * di mach
779  */
780  im = bisec_d(m, mach, 0, nm - 1);
781 
782  if (im == nm - 1) {
783  return s[3*nm - 1];
784 
785  } else if (im == -1) {
786  return s[2*nm];
787 
788  } else {
789  doublereal d;
790 
791  im++;
792  d = (mach - m[im - 1])/(m[im] - m[im - 1]);
793 
794  return (1. - d)*s[2*nm + im - 1] + d*s[2*nm + im];
795  }
796 }
int bisec_d(doublereal *v, doublereal val, int lb, int ub)
Definition: bisec.cc:57
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
double doublereal
Definition: colamd.c:52

Here is the call graph for this function:

Variable Documentation

const outa_t outa_Zero

Definition at line 87 of file aerodc81.c.