MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
aerodc81.c
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aerodc81.c,v 1.54 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33 
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <math.h>
37 #include "ac/f2c.h"
38 #include "aerodc81.h"
39 #include "bisec.h"
40 #include "c81data.h"
41 
42 /*
43 cC************************* AEROD *********************
44 c MODIFICATA NEL CALCOLO TNG
45 C= COMPILER (LINK=IBJ$)
46  SUBROUTINE AEROD2(W, VAM, TNG, OUTA, INST, RSPEED, JPRO)
47 C W: vettore lungo 6; velocita' del corpo aerodinamico nel sistema locale,
48 C riferita ad un punto arbitrario
49 C VAM: vettore lungo 6, contiene dati del problema
50 C TNG: vettore lungo 12, presumubilmente il termine noto
51 C (ora l'ho ridotto a 6)
52 C OUTA: vettore lungo 20, usato in i/o con subroutines di aerod
53 C INST: flag di instazionario (>0)
54 C RSPEED -> MODOMEGA: modulo della velocita' di rotazione
55 C JPRO: tipo di profilo
56 C DEFINIZIONI VETTORE VAM
57 C
58 C VAM(1): densita' dell'aria
59 C VAM(2): celerita' del suono
60 C VAM(3): corda
61 C VAM(4): 1/4 corda
62 C VAM(5): 3/4 corda
63 C VAM(6): svergolamento
64 C
65 C W = [v1 v2 v3 w1 w2 w3 ]
66 C con:
67 C 1 direzione di avanzamento,
68 C 2 direzione normale,
69 C 3 direzione lungo l'apertura,
70 C DM*W da' la velocita' nel punto a 3/4 della corda
71 */
72 
73 static int
74 get_coef(int nm, doublereal* m, int na, doublereal* a,
75  doublereal alpha, doublereal mach,
76  doublereal* c, doublereal* c0);
77 
78 static doublereal
79 get_dcla(int nm, doublereal* m, doublereal* s, doublereal mach);
80 
81 #ifdef USE_GET_STALL
82 static int
83 get_stall(int nm, doublereal* m, doublereal* s, doublereal mach,
84  doublereal *dcpa, doublereal *dasp, doublereal *dasm);
85 #endif /* USE_GET_STALL */
86 
87 const outa_t outa_Zero;
88 
89 doublereal
90 c81_data_get_coef(int nm, doublereal* m, int na, doublereal* a, doublereal alpha, doublereal mach)
91 {
92  doublereal c;
93 
94  get_coef(nm, m, na, a, alpha, mach, &c, NULL);
95 
96  return c;
97 }
98 
99 int
100 c81_aerod2(doublereal* W, const vam_t *VAM, doublereal* TNG, outa_t* OUTA, c81_data* data)
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 }
241 
242 int
243 c81_aerod2_u(doublereal* W, const vam_t *VAM, doublereal* TNG, outa_t* OUTA,
244  c81_data* data, long unsteadyflag)
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 }
626 
627 /*
628  * trova un coefficiente dato l'angolo ed il numero di Mach
629  *
630  * il numero di Mach mach viene cercato nell'array m di lunghezza nm
631  * ed interpolato linearmente; quindi il coefficiente alpha viene cercato
632  * nella prima colonna della matrice a di dimensioni na x nm + 1 ed interpolato
633  * linearmente; infine il coefficiente corrispondente alla combinazione di
634  * mach e alpha viene restituito.
635  */
636 static int
637 get_coef(int nm, doublereal* m, int na, doublereal* a, doublereal alpha, doublereal mach,
638  doublereal* c, doublereal* c0)
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 }
768 
769 static doublereal
770 get_dcla(int nm, doublereal* m, doublereal* s, doublereal mach)
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 }
797 
798 #ifdef USE_GET_STALL
799 static int
800 get_stall(int nm, doublereal* m, doublereal* s, doublereal mach,
801  doublereal *dcpa, doublereal *dasp, doublereal *dasm)
802 {
803  int im;
804 
805  mach = fabs(mach);
806 
807  /*
808  * im e' l'indice di m in cui si trova l'approssimazione per eccesso
809  * di mach
810  */
811  im = bisec_d(m, mach, 0, nm - 1);
812  if (im != nm) {
813  im++;
814  }
815 
816  if (im == nm) {
817  *dcpa = s[3*nm - 1];
818  *dasp = s[nm - 1];
819  *dasm = s[2*nm - 1];
820  } else {
821  doublereal d = (mach - m[im - 1])/(m[im] - m[im - 1]);
822 
823  *dcpa = (1. - d)*s[2*nm + im - 1] + d*s[2*nm + im];
824  *dasp = (1. - d)*s[im - 1] + d*s[im];
825  *dasm = (1. - d)*s[nm + im - 1] + d*s[nm + im];
826  }
827 
828  return 0;
829 }
830 #endif /* USE_GET_STALL */
831 
grad::exp
GradientExpression< UnaryExpr< FuncExp, Expr > > exp(const GradientExpression< Expr > &u)
Definition: gradient.h:2975
outa_t::cd
doublereal cd
Definition: aerodc81.h:66
vam_t::density
doublereal density
Definition: aerodc81.h:108
c81_data::NMD
int NMD
Definition: aerodc81.h:160
outa_t::dcn
doublereal dcn
Definition: aerodc81.h:74
M_PI
#define M_PI
Definition: gradienttest.cc:67
c81_data::NAM
int NAM
Definition: aerodc81.h:166
grad::pow
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
vam_t::sound_celerity
doublereal sound_celerity
Definition: aerodc81.h:109
c81_data::NMM
int NMM
Definition: aerodc81.h:165
bisec_d
int bisec_d(doublereal *v, doublereal val, int lb, int ub)
Definition: bisec.cc:57
outa_t
Definition: aerodc81.h:60
outa_t::dcm
doublereal dcm
Definition: aerodc81.h:75
outa_t::gamma
doublereal gamma
Definition: aerodc81.h:63
c81_data::NML
int NML
Definition: aerodc81.h:146
outa_t::alf1
doublereal alf1
Definition: aerodc81.h:70
vam_t::bc_position
doublereal bc_position
Definition: aerodc81.h:112
c81_data::mm
doublereal * mm
Definition: aerodc81.h:167
c81_aerod2_u
int c81_aerod2_u(doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data, long unsteadyflag)
Definition: aerodc81.c:243
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
c81_data
Definition: aerodc81.h:143
outa_t::alpha
doublereal alpha
Definition: aerodc81.h:62
outa_t::alf2
doublereal alf2
Definition: aerodc81.h:71
c81_data::NAL
int NAL
Definition: aerodc81.h:147
get_coef
static int get_coef(int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach, doublereal *c, doublereal *c0)
Definition: aerodc81.c:637
outa_t::mach
doublereal mach
Definition: aerodc81.h:64
outa_t::dan
doublereal dan
Definition: aerodc81.h:72
vam_t::force_position
doublereal force_position
Definition: aerodc81.h:111
c81_data::ad
doublereal * ad
Definition: aerodc81.h:163
vam_t::chord
doublereal chord
Definition: aerodc81.h:110
outa_t::cl
doublereal cl
Definition: aerodc81.h:65
c81_aerod2
int c81_aerod2(doublereal *W, const vam_t *VAM, doublereal *TNG, outa_t *OUTA, c81_data *data)
Definition: aerodc81.c:100
c81_data::mstall
doublereal * mstall
Definition: aerodc81.h:158
c81_data_get_coef
doublereal c81_data_get_coef(int nm, doublereal *m, int na, doublereal *a, doublereal alpha, doublereal mach)
Definition: aerodc81.c:90
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
outa_t::clalpha
doublereal clalpha
Definition: aerodc81.h:69
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
c81_data::stall
doublereal * stall
Definition: aerodc81.h:157
c81_data::am
doublereal * am
Definition: aerodc81.h:168
outa_t::dam
doublereal dam
Definition: aerodc81.h:73
grad::cos
GradientExpression< UnaryExpr< FuncCos, Expr > > cos(const GradientExpression< Expr > &u)
Definition: gradient.h:2978
c81_data::al
doublereal * al
Definition: aerodc81.h:149
a
static const doublereal a
Definition: hfluid_.h:289
outa_t::cm
doublereal cm
Definition: aerodc81.h:67
get_dcla
static doublereal get_dcla(int nm, doublereal *m, doublereal *s, doublereal mach)
Definition: aerodc81.c:770
grad::atan2
GradientExpression< BinaryExpr< FuncAtan2, LhsExpr, RhsExpr > > atan2(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2962
outa_Zero
const outa_t outa_Zero
Definition: aerodc81.c:87
doublereal
double doublereal
Definition: colamd.c:52
c81_data::md
doublereal * md
Definition: aerodc81.h:162
c81_data::ml
doublereal * ml
Definition: aerodc81.h:148
vam_t
Definition: aerodc81.h:107
c81_data::NAD
int NAD
Definition: aerodc81.h:161