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
c81_data_get_coef(int nm,doublereal * m,int na,doublereal * a,doublereal alpha,doublereal mach)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
c81_aerod2(doublereal * W,const vam_t * VAM,doublereal * TNG,outa_t * OUTA,c81_data * data)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
c81_aerod2_u(doublereal * W,const vam_t * VAM,doublereal * TNG,outa_t * OUTA,c81_data * data,long unsteadyflag)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
get_coef(int nm,doublereal * m,int na,doublereal * a,doublereal alpha,doublereal mach,doublereal * c,doublereal * c0)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
get_dcla(int nm,doublereal * m,doublereal * s,doublereal mach)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
get_stall(int nm,doublereal * m,doublereal * s,doublereal mach,doublereal * dcpa,doublereal * dasp,doublereal * dasm)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
832