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