1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/c81data.cc,v 1.52 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 <cerrno>
35 #include <cstdlib>
36 #include <cstdio>
37 #include <cstring>
38
39 #include <iostream>
40 #include <sstream>
41 #include <iomanip>
42 #include <cmath>
43
44 #include "ac/f2c.h"
45 #include "myassert.h"
46
47 extern "C" {
48 #include "aerodc81.h"
49 }
50
51 #include "c81data.h"
52
53 /*
54 * header NML,NAL,NMD,NAD,NMM,NAM A30,6I2
55 * ML(1),....,ML(NML) 7x,9F7.0 eventualmente su piu' righe
56 * AL(1) CL(1,1),....,CL(1,NML) 10F7.0/(7x,9F7.0)
57 * : : : :
58 * AL(NAL)CL(NAL,1),....,CL(NAL,NML) 10F7.0/(7x,9F7.0)
59 * AD(1) CD(1,1),....,CD(1,NMD) 10F7.0/(7x,9F7.0)
60 * : : : :
61 * AD(NAD)CD(NAD,1),....,CD(NAD,NMD) 10F7.0/(7x,9F7.0)
62 * AM(1) CM(1,1),....,CL(1,NMM) 10F7.0/(7x,9F7.0)
63 * : : : :
64 * AM(NAM)CM(NAM,1),....,CL(NAM,NMM) 10F7.0/(7x,9F7.0)
65 */
66
67 static int
68 do_stall(int NM, int NA, doublereal *a, doublereal *stall, const doublereal dcltol);
69
70 static int
get_int(const char * const from,int & i)71 get_int(const char *const from, int &i)
72 {
73 #ifdef HAVE_STRTOL
74 char *endptr = NULL;
75 errno = 0;
76 i = strtol(from, &endptr, 10);
77 int save_errno = errno;
78 if (endptr != NULL && endptr[0] != '\0' && !isspace(endptr[0])) {
79 return -1;
80
81 } else if (save_errno == ERANGE) {
82 silent_cerr("c81data: warning, int "
83 << std::string(from, endptr - from)
84 << " overflows" << std::endl);
85 return -1;
86 }
87 #else /* !HAVE_STRTOL */
88 i = atoi(buf);
89 #endif /* !HAVE_STRTOL */
90 return 0;
91 }
92
93 #if 0 // unused so far
94 static int
95 get_long(const char *const from, long &l)
96 {
97 #ifdef HAVE_STRTOL
98 char *endptr = NULL;
99 errno = 0;
100 l = strtol(from, &endptr, 10);
101 int save_errno = errno;
102 if (endptr != NULL && endptr[0] != '\0' && !isspace(endptr[0])) {
103 return -1;
104
105 } else if (save_errno == ERANGE) {
106 silent_cerr("c81data: warning, int "
107 << std::string(from, endptr - from)
108 << " overflows" << std::endl);
109 return -1;
110 }
111 #else /* !HAVE_STRTOL */
112 l = atol(buf);
113 #endif /* !HAVE_STRTOL */
114 return 0;
115 }
116 #endif
117
118 static int
get_double(const char * const from,doublereal & d)119 get_double(const char *const from, doublereal &d)
120 {
121 #ifdef HAVE_STRTOD
122 char *endptr = NULL;
123 errno = 0;
124 d = strtod(from, &endptr);
125 int save_errno = errno;
126 if (endptr != NULL && endptr[0] != '\0' && !isspace(endptr[0])) {
127 return -1;
128
129 } else if (save_errno == ERANGE) {
130 silent_cerr("c81data: warning, double "
131 << std::string(from, endptr - from)
132 << " overflows" << std::endl);
133 return -1;
134 }
135 #else /* !HAVE_STRTOD */
136 d = atof(buf);
137 #endif /* !HAVE_STRTOD */
138 return 0;
139 }
140
141 static int
get_c81_vec(std::istream & in,doublereal * v,int ncols)142 get_c81_vec(std::istream& in, doublereal* v, int ncols)
143 {
144 char buf[128];
145
146 if (!in || v == NULL || ncols < 1) {
147 return -1;
148 }
149
150 for (int c = 0; c < ncols; c++) {
151 char tmp[8];
152 int i = c%9;
153
154 if (i == 0) {
155 in.getline(buf, sizeof(buf));
156 }
157
158 /* always skip first */
159 memcpy(tmp, &buf[7*(i + 1)], 7);
160 tmp[7] = '\0';
161
162 if (get_double(tmp, v[c])) {
163 return -1;
164 }
165 }
166
167 return 0;
168 }
169
170 static int
check_vec(doublereal * v,int n)171 check_vec(doublereal* v, int n)
172 {
173 for (int i = 1; i < n; i++) {
174 if (v[i] <= v[i - 1]) {
175 return i;
176 }
177 }
178
179 return 0;
180 }
181
182 static int
get_c81_mat(std::istream & in,doublereal * m,int nrows,int ncols)183 get_c81_mat(std::istream& in, doublereal* m, int nrows, int ncols)
184 {
185 char buf[128];
186
187 if (!in || m == NULL || nrows < 1 || ncols < 1) {
188 return -1;
189 }
190
191 for (int r = 0; r < nrows; r++) {
192 int row = 0;
193 int offset = 0;
194
195 for (int c = 0; c < ncols; c++) {
196 char tmp[8];
197 int i = (c - offset)%(9 - offset + 1);
198
199 if (i == 0) {
200 if (++row == 2) {
201 offset = 1;
202 }
203 in.getline(buf, sizeof(buf));
204 }
205
206 /* skip first except first time */
207 memcpy(tmp, &buf[7*(i + offset)], 7);
208 tmp[7] = '\0';
209
210 if (get_double(tmp, m[r + nrows*c])) {
211 return -1;
212 }
213 }
214 }
215
216 return 0;
217 }
218
219 static int
put_c81_row(std::ostream & out,doublereal * v,int dim,int ncols,int first=0)220 put_c81_row(std::ostream& out, doublereal* v, int dim, int ncols, int first = 0)
221 {
222 int start = 0;
223 const int N = 9;
224
225 if (first) {
226 out << std::setw(7) << v[0];
227 start = dim;
228 ncols--;
229 } else {
230 out << std::setw(7) << "";
231 }
232
233 for (int i = 0; i < (ncols-1)/N; i++) {
234 for (int j = 0; j < N; j++) {
235 out << std::setw(7) << v[start+dim*j];
236 }
237 out << std::endl << std::setw(7) << "";
238 start += dim*N;
239 }
240
241 for (int j = 0; j < (ncols-1)%N+1; j++) {
242 out << std::setw(7) << v[start+dim*j];
243 }
244 out << std::endl;
245
246 return 0;
247 }
248
249 static int
put_c81_vec(std::ostream & out,doublereal * v,int nrows)250 put_c81_vec(std::ostream& out, doublereal* v, int nrows)
251 {
252 if (!out || v == NULL || nrows < 1) {
253 return -1;
254 }
255
256 put_c81_row(out, v, 1, nrows);
257
258 return 0;
259 }
260
261 static int
put_c81_mat(std::ostream & out,doublereal * m,int nrows,int ncols)262 put_c81_mat(std::ostream& out, doublereal* m, int nrows, int ncols)
263 {
264 if (!out || m == NULL || nrows < 1 || ncols < 1) {
265 return -1;
266 }
267
268 for (int i = 0; i < nrows; i++) {
269 put_c81_row(out, m+i, nrows, ncols, 1);
270 }
271
272 return 0;
273 }
274
275 extern "C" void
c81_data_destroy(c81_data * data)276 c81_data_destroy(c81_data* data)
277 {
278 delete[] data->ml;
279 delete[] data->al;
280
281 delete[] data->md;
282 delete[] data->ad;
283
284 delete[] data->mm;
285 delete[] data->am;
286
287 delete[] data->stall;
288 delete[] data->mstall;
289 }
290
291 extern "C" int
c81_data_read(std::istream & in,c81_data * data,const doublereal dcltol,int * ff)292 c81_data_read(std::istream& in, c81_data* data, const doublereal dcltol, int *ff)
293 {
294 char buf[BUFSIZ]; // 81 should suffice
295
296 if (ff) {
297 *ff = 0;
298 }
299
300 /* header */
301 in.getline(buf, sizeof(buf));
302 if (strncasecmp(buf, "# FREE FORMAT", STRLENOF("# FREE FORMAT")) == 0) {
303 if (ff) {
304 *ff = 1;
305 }
306
307 return c81_data_read_free_format(in, data, dcltol);
308 }
309
310 buf[42] = '\0';
311 if (get_int(&buf[40], data->NAM)) {
312 return -1;
313 }
314
315 buf[40] = '\0';
316 if (get_int(&buf[38], data->NMM)) {
317 return -1;
318 }
319
320 buf[38] = '\0';
321 if (get_int(&buf[36], data->NAD)) {
322 return -1;
323 }
324
325 buf[36] = '\0';
326 if (get_int(&buf[34], data->NMD)) {
327 return -1;
328 }
329
330 buf[34] = '\0';
331 if (get_int(&buf[32], data->NAL)) {
332 return -1;
333 }
334
335 buf[32] = '\0';
336 if (get_int(&buf[30], data->NML)) {
337 return -1;
338 }
339
340 if (data->NAM <= 0 || data->NMM <= 0
341 || data->NAD <= 0 || data->NMD <= 0
342 || data->NAL <= 0 || data->NML <= 0) {
343 return -1;
344 }
345
346 buf[30] = '\0';
347 strncpy(data->header, buf, 30);
348 data->header[30] = '\0';
349
350 /* lift */
351 data->ml = new doublereal[data->NML];
352 get_c81_vec(in, data->ml, data->NML);
353 if (check_vec(data->ml, data->NML)) {
354 return -1;
355 }
356
357 data->al = new doublereal[(data->NML+1)*data->NAL];
358 get_c81_mat(in, data->al, data->NAL, data->NML+1);
359 if (check_vec(data->al, data->NAL)) {
360 return -1;
361 }
362
363 /* drag */
364 data->md = new doublereal[data->NMD];
365 get_c81_vec(in, data->md, data->NMD);
366 if (check_vec(data->md, data->NMD)) {
367 return -1;
368 }
369
370 data->ad = new doublereal[(data->NMD+1)*data->NAD];
371 get_c81_mat(in, data->ad, data->NAD, data->NMD+1);
372 if (check_vec(data->ad, data->NAD)) {
373 return -1;
374 }
375
376 /* moment */
377 data->mm = new doublereal[data->NMM];
378 get_c81_vec(in, data->mm, data->NMM);
379 if (check_vec(data->mm, data->NMM)) {
380 return -1;
381 }
382
383 data->am = new doublereal[(data->NMM+1)*data->NAM];
384 get_c81_mat(in, data->am, data->NAM, data->NMM+1);
385 if (check_vec(data->am, data->NAM)) {
386 return -1;
387 }
388
389 /* FIXME: maybe this is not the best place */
390 c81_data_do_stall(data, dcltol);
391
392 return 0;
393 }
394
395 extern "C" int
c81_data_merge(unsigned ndata,const c81_data ** data,const doublereal * upper_bounds,doublereal dCsi,doublereal dcltol,c81_data * i_data)396 c81_data_merge(
397 unsigned ndata,
398 const c81_data **data,
399 const doublereal *upper_bounds,
400 doublereal dCsi,
401 doublereal dcltol,
402 c81_data *i_data)
403 {
404 ASSERT(ndata > 0);
405 ASSERT(data != 0);
406 ASSERT(upper_bounds != 0);
407 ASSERT(dCsi >= -1.);
408 ASSERT(dCsi <= 1.);
409 ASSERT(dcltol > 0.);
410 ASSERT(i_data != 0);
411
412 if (dCsi < upper_bounds[0]) {
413 silent_cerr("cannot find C81 data lower bound for point xi=" << dCsi << std::endl);
414 return -1;
415 }
416
417 int to = 0;
418 for (unsigned i = 1; i < ndata; i++) {
419 if (upper_bounds[i] > dCsi) {
420 to = i;
421 break;
422 }
423 }
424 if (to == 0) {
425 silent_cerr("cannot find C81 data upper bound for point xi=" << dCsi << std::endl);
426 return -1;
427 }
428
429 int from = to - 1;
430 if (unsigned(from) == ndata) {
431 silent_cerr("cannot find C81 data upper bound for point xi=" << dCsi << std::endl);
432 return -1;
433 }
434
435 /* we need to interpolate between data[from]
436 * and data[to] */
437
438 /* we only accept homogeneous data sources,
439 * i.e. same Mach and alpha patterns */
440 if (data[from]->NML != data[to]->NML) {
441 silent_cerr("number of Mach points for Cl between airfoils "
442 << from << " (" << data[from]->NML << ") and "
443 << to << " (" << data[to]->NML << ") do not match"
444 << std::endl);
445 return -1;
446 }
447
448 if (data[from]->NAL != data[to]->NAL) {
449 silent_cerr("number of AoA points for Cl between airfoils "
450 << from << " (" << data[from]->NAL << ") and "
451 << to << " (" << data[to]->NAL << ") do not match"
452 << std::endl);
453 return -1;
454 }
455
456 if (data[from]->NMD != data[to]->NMD) {
457 silent_cerr("number of Mach points for Cd between airfoils "
458 << from << " (" << data[from]->NMD << ") and "
459 << to << " (" << data[to]->NMD << ") do not match"
460 << std::endl);
461 return -1;
462 }
463
464 if (data[from]->NAD != data[to]->NAD) {
465 silent_cerr("number of AoA points for Cd between airfoils "
466 << from << " (" << data[from]->NAD << ") and "
467 << to << " (" << data[to]->NAD << ") do not match"
468 << std::endl);
469 return -1;
470 }
471
472 if (data[from]->NMM != data[to]->NMM) {
473 silent_cerr("number of Mach points for Cm between airfoils "
474 << from << " (" << data[from]->NMM << ") and "
475 << to << " (" << data[to]->NMM << ") do not match"
476 << std::endl);
477 return -1;
478 }
479
480 if (data[from]->NAM != data[to]->NAM) {
481 silent_cerr("number of AoA points for Cm between airfoils "
482 << from << " (" << data[from]->NAM << ") and "
483 << to << " (" << data[to]->NAM << ") do not match"
484 << std::endl);
485 return -1;
486 }
487
488 for (int i = 0; i < data[from]->NML; i++) {
489 if (data[from]->ml[i] != data[to]->ml[i]) {
490 silent_cerr("Mach point " << i << "for Cl of airfoils "
491 << from << " (" << data[from]->ml[i] << ") and "
492 << to << " (" << data[to]->ml[i] << ") differs"
493 << std::endl);
494 return -1;
495 }
496 }
497
498 for (int i = 0; i < data[from]->NAL; i++) {
499 if (data[from]->al[i] != data[to]->al[i]) {
500 silent_cerr("AoA point " << i << "for Cl of airfoils "
501 << from << " (" << data[from]->al[i] << ") and "
502 << to << " (" << data[to]->al[i] << ") differs"
503 << std::endl);
504 return -1;
505 }
506 }
507
508 for (int i = 0; i < data[from]->NMD; i++) {
509 if (data[from]->md[i] != data[to]->md[i]) {
510 silent_cerr("Mach point " << i << "for Cd of airfoils "
511 << from << " (" << data[from]->md[i] << ") and "
512 << to << " (" << data[to]->md[i] << ") differs"
513 << std::endl);
514 return -1;
515 }
516 }
517
518 for (int i = 0; i < data[from]->NAD; i++) {
519 if (data[from]->ad[i] != data[to]->ad[i]) {
520 silent_cerr("AoA point " << i << "for Cd of airfoils "
521 << from << " (" << data[from]->ad[i] << ") and "
522 << to << " (" << data[to]->ad[i] << ") differs"
523 << std::endl);
524 return -1;
525 }
526 }
527
528 for (int i = 0; i < data[from]->NMM; i++) {
529 if (data[from]->mm[i] != data[to]->mm[i]) {
530 silent_cerr("Mach point " << i << "for Cm of airfoils "
531 << from << " (" << data[from]->mm[i] << ") and "
532 << to << " (" << data[to]->mm[i] << ") differs"
533 << std::endl);
534 return -1;
535 }
536 }
537
538 for (int i = 0; i < data[from]->NAM; i++) {
539 if (data[from]->am[i] != data[to]->am[i]) {
540 silent_cerr("AoA point " << i << "for Cm of airfoils "
541 << from << " (" << data[from]->am[i] << ") and "
542 << to << " (" << data[to]->am[i] << ") differs"
543 << std::endl);
544 return -1;
545 }
546 }
547
548 snprintf(i_data->header, sizeof(i_data->header),
549 "interpolated (\"%s\"[%e]->\"%s\"[%e]: [%e])",
550 data[from]->header, upper_bounds[from],
551 data[to]->header, upper_bounds[to], dCsi);
552
553 i_data->NML = data[from]->NML;
554 i_data->NAL = data[from]->NAL;
555 i_data->NMD = data[from]->NMD;
556 i_data->NAD = data[from]->NAD;
557 i_data->NMM = data[from]->NMM;
558 i_data->NAM = data[from]->NAM;
559
560 doublereal dw = upper_bounds[to] - upper_bounds[from];
561 doublereal dw_from = (upper_bounds[to] - dCsi)/dw;
562 doublereal dw_to = (dCsi - upper_bounds[from])/dw;
563
564 /* lift */
565 i_data->ml = new doublereal[i_data->NML];
566 for (int i = 0; i < i_data->NML; i++) {
567 i_data->ml[i] = data[from]->ml[i];
568 }
569 i_data->al = new doublereal[(i_data->NML + 1)*i_data->NAL];
570 for (int i = 0; i < i_data->NAL; i++) {
571 i_data->al[i] = data[from]->al[i];
572 }
573 for (int i = i_data->NAL; i < (i_data->NML + 1)*i_data->NAL; i++) {
574 i_data->al[i] = dw_from*data[from]->al[i] + dw_to*data[to]->al[i];
575 }
576
577 /* drag */
578 i_data->md = new doublereal[i_data->NMD];
579 for (int i = 0; i < i_data->NMD; i++) {
580 i_data->md[i] = data[from]->md[i];
581 }
582 i_data->ad = new doublereal[(i_data->NMD + 1)*i_data->NAD];
583 for (int i = 0; i < i_data->NAD; i++) {
584 i_data->ad[i] = data[from]->ad[i];
585 }
586 for (int i = i_data->NAD; i < (i_data->NMD + 1)*i_data->NAD; i++) {
587 i_data->ad[i] = dw_from*data[from]->ad[i] + dw_to*data[to]->ad[i];
588 }
589
590 /* moment */
591 i_data->mm = new doublereal[i_data->NMM];
592 for (int i = 0; i < i_data->NMM; i++) {
593 i_data->mm[i] = data[from]->mm[i];
594 }
595 i_data->am = new doublereal[(i_data->NMM + 1)*i_data->NAM];
596 for (int i = 0; i < i_data->NAM; i++) {
597 i_data->am[i] = data[from]->am[i];
598 }
599 for (int i = i_data->NAM; i < (i_data->NMM + 1)*i_data->NAM; i++) {
600 i_data->am[i] = dw_from*data[from]->am[i] + dw_to*data[to]->am[i];
601 }
602
603 // FIXME: maybe this is not the best place
604 c81_data_do_stall(i_data, dcltol);
605
606 return 0;
607 }
608
609 extern "C" int
c81_data_write(std::ostream & out,c81_data * data)610 c81_data_write(std::ostream& out, c81_data* data)
611 {
612 if (data == 0) {
613 return -1;
614 }
615
616 std::ios::fmtflags tmpflags;
617 tmpflags = out.flags(std::ios::left);
618 out << std::setw(30) << data->header;
619 out.flags(tmpflags);
620 out
621 << std::setw(2) << data->NML
622 << std::setw(2) << data->NAL
623 << std::setw(2) << data->NMD
624 << std::setw(2) << data->NAD
625 << std::setw(2) << data->NMM
626 << std::setw(2) << data->NAM
627 << std::endl;
628
629 put_c81_vec(out, data->ml, data->NML);
630 put_c81_mat(out, data->al, data->NAL, data->NML+1);
631
632 put_c81_vec(out, data->md, data->NMD);
633 put_c81_mat(out, data->ad, data->NAD, data->NMD+1);
634
635 put_c81_vec(out, data->mm, data->NMM);
636 put_c81_mat(out, data->am, data->NAM, data->NMM+1);
637
638 return 0;
639 }
640
641 extern "C" int
read_fc511_row(std::istream & in,doublereal * d,int NC)642 read_fc511_row(std::istream& in, doublereal *d, int NC)
643 {
644 int r;
645
646 for (r = 0; r < NC/8; r++) {
647 char buf[81];
648
649 in.getline(buf, sizeof(buf));
650
651 for (int c = 7; c >= 0; c--) {
652 buf[10*(c + 1)] = '\0';
653 if (get_double(&buf[10*c], d[8*r + c])) {
654 return -1;
655 }
656 }
657 }
658
659 if (NC%8) {
660 char buf[81];
661
662 in.getline(buf, sizeof(buf));
663
664 for (int c = NC%8 - 1; c >= 0; c--) {
665 buf[10*(c + 1)] = '\0';
666 if (get_double(&buf[10*c], d[8*r + c])) {
667 return -1;
668 }
669 }
670 }
671
672 return 0;
673 }
674
675 extern "C" int
read_fc511_mat(std::istream & in,doublereal * d,int NR,int NC)676 read_fc511_mat(std::istream& in, doublereal *d, int NR, int NC)
677 {
678 for (int i = 0; i < NR; i++) {
679 int r;
680
681 for (r = 0; r < NC/8; r++) {
682 char buf[81];
683
684 in.getline(buf, sizeof(buf));
685
686 for (int c = 7; c >= 0; c--) {
687 buf[10*(c + 1)] = '\0';
688 if (get_double(&buf[10*c], d[NR*(8*r + c) + i])) {
689 return -1;
690 }
691 }
692 }
693
694 if (NC%8) {
695 char buf[81];
696
697 in.getline(buf, sizeof(buf));
698
699 for (int c = NC%8 - 1; c >= 0; c--) {
700 buf[10*(c + 1)] = '\0';
701 if (get_double(&buf[10*c], d[NR*(8*r + c) + i])) {
702 return -1;
703 }
704 }
705 }
706 }
707
708 return 0;
709 }
710
711 extern "C" int
read_fc511_block(std::istream & in,int & NA,int & NM,doublereal * & da,doublereal * & dm)712 read_fc511_block(std::istream& in, int &NA, int &NM, doublereal *&da, doublereal *&dm)
713 {
714 char buf[128]; // 81 should suffice; let's make it 128
715
716 in.getline(buf, sizeof(buf));
717
718 buf[10] = '\0';
719 if (get_int(&buf[5], NM)) {
720 return -1;
721 }
722
723 buf[5] = '\0';
724 if (get_int(&buf[0], NA)) {
725 return -1;
726 }
727
728 if (NA <= 0 || NM <= 0) {
729 return -1;
730 }
731
732 dm = new doublereal[NM];
733 doublereal *tmpda = new doublereal[NA*(NM + 1)];
734
735 if (read_fc511_row(in, tmpda, NA)) {
736 return -1;
737 }
738
739 if (read_fc511_row(in, dm, NM)) {
740 return -1;
741 }
742
743 if (read_fc511_mat(in, &tmpda[NA], NA, NM)) {
744 return -1;
745 }
746
747 int pos;
748 in.getline(buf, sizeof(buf));
749 if (get_int(buf, pos)) {
750 return -1;
751 }
752
753 doublereal *posda = new doublereal[2*pos];
754
755 if (read_fc511_row(in, posda, pos)) {
756 return -1;
757 }
758
759 if (read_fc511_row(in, &posda[pos], pos)) {
760 return -1;
761 }
762
763 int neg;
764 in.getline(buf, sizeof(buf));
765 if (get_int(buf, neg)) {
766 return -1;
767 }
768
769 int realNA = (neg + NA + pos);
770 da = new doublereal[realNA*(NM + 1)];
771
772 if (read_fc511_row(in, da, neg)) {
773 return -1;
774 }
775
776 if (read_fc511_row(in, &da[realNA], neg)) {
777 return -1;
778 }
779
780 for (int c = 0; c < 2; c++) {
781 for (int r = 0; r < pos; r++) {
782 da[realNA*c + neg + NA + r] = posda[pos*c + r];
783 }
784 }
785
786 for (int c = 2; c < NM + 1; c++) {
787 for (int r = 0; r < neg; r++) {
788 da[realNA*c + r] = da[realNA + r];
789 }
790 for (int r = 0; r < pos; r++) {
791 da[realNA*c + neg + NA + r] = posda[pos + r];
792 }
793 }
794
795 for (int c = 0; c < (NM + 1); c++) {
796 for (int r = 0; r < NA; r++) {
797 da[realNA*c + neg + r] = tmpda[NA*c + r];
798 }
799 }
800
801 delete[] tmpda;
802 delete[] posda;
803
804 NA = realNA;
805
806 return 0;
807 }
808
809 extern "C" int
c81_data_fc511_read(std::istream & in,c81_data * data,const doublereal dcltol)810 c81_data_fc511_read(std::istream& in, c81_data* data, const doublereal dcltol)
811 {
812 char buf[128]; // 81 should suffice; let's make it 128
813
814 /* header */
815 in.getline(buf, sizeof(buf));
816
817 memcpy(data->header, buf, sizeof(data->header));
818 data->header[STRLENOF(data->header)] = '\0';
819
820 char *p;
821
822 in.getline(buf, sizeof(buf));
823 p = buf;
824
825 while (isspace(p[0])) {
826 p++;
827 }
828
829 if (!p[0] || strncasecmp(p, "cl", 2) != 0 || (p[2] != '\0' && !isspace(p[2]))) {
830 return -1;
831 }
832
833 if (read_fc511_block(in, data->NAL, data->NML, data->al, data->ml)) {
834 return -1;
835 }
836
837 in.getline(buf, sizeof(buf));
838 p = buf;
839
840 while (isspace(p[0])) {
841 p++;
842 }
843
844 if (!p[0] || strncasecmp(p, "cd", 2) != 0 || (p[2] != '\0' && !isspace(p[2]))) {
845 return -1;
846 }
847
848 if (read_fc511_block(in, data->NAD, data->NMD, data->ad, data->md)) {
849 return -1;
850 }
851
852 in.getline(buf, sizeof(buf));
853 p = buf;
854
855 while (isspace(p[0])) {
856 p++;
857 }
858
859 if (!p[0] || strncasecmp(p, "cm", 2) != 0 || (p[2] != '\0' && !isspace(p[2]))) {
860 return -1;
861 }
862
863 if (read_fc511_block(in, data->NAM, data->NMM, data->am, data->mm)) {
864 return -1;
865 }
866
867 /* FIXME: maybe this is not the best place */
868 c81_data_do_stall(data, dcltol);
869
870 return 0;
871 }
872
873 extern "C" int
c81_data_nrel_read(std::istream & in,c81_data * data,const doublereal dcltol)874 c81_data_nrel_read(std::istream& in, c81_data* data, const doublereal dcltol)
875 {
876 char buf[128]; // 81 should suffice; let's make it 128
877
878 // header
879 in.getline(buf, sizeof(buf));
880
881 // keep the first one
882 memcpy(data->header, buf, sizeof(data->header));
883 data->header[STRLENOF(data->header)] = '\0';
884
885 // discard the second one
886 in.getline(buf, sizeof(buf));
887
888 int n;
889 in >> n;
890 if (n != 1) {
891 silent_cerr("NREL format: can only handle one airfoil per file (got " << n << ")" << std::endl);
892 return -1;
893 }
894
895 // discard remaining lines
896 for (int i = 0; i < 12; i++) {
897 // discard
898 in.getline(buf, sizeof(buf));
899 }
900
901 std::streampos pos = in.tellg();
902 for (n = 0; in; n++) {
903 in.getline(buf, sizeof(buf));
904 }
905
906 in.clear();
907 in.seekg(pos);
908 n--;
909
910 if (n < 2) {
911 silent_cerr("insufficient number of data points n=" << n << " in NREL data format (n >= 2 required)" << std::endl);
912 return -1;
913 }
914
915 /* lift */
916 data->NML = 1;
917 data->ml = new doublereal[data->NML];
918 data->ml[0] = 0.;
919 data->NAL = n;
920 data->al = new doublereal[2*data->NAL];
921
922 /* drag */
923 data->NMD = 1;
924 data->md = new doublereal[data->NMD];
925 data->md[0] = 0.;
926 data->NAD = n;
927 data->ad = new doublereal[2*data->NAD];
928
929 /* moment */
930 data->NMM = 1;
931 data->mm = new doublereal[data->NMM];
932 data->mm[0] = 0.;
933 data->NAM = n;
934 data->am = new doublereal[2*data->NAM];
935
936 for (int i = 0; i < n; i++) {
937 doublereal alpha, cl, cd, cm;
938 in >> alpha >> cl >> cd >> cm;
939
940 if (i == 0 && alpha != -180.) {
941 silent_cerr("warning: NREL data format: first alpha=" << alpha << " != -180.)" << std::endl);
942 } else if (i == n - 1 && alpha != 180.) {
943 silent_cerr("warning: NREL data format: last alpha=" << alpha << " != 180.)" << std::endl);
944 }
945
946 data->al[i] = alpha;
947 data->al[data->NAL + i] = cl;
948 data->ad[i] = alpha;
949 data->ad[data->NAD + i] = cd;
950 data->am[i] = alpha;
951 data->am[data->NAM + i] = cm;
952 }
953
954 /* FIXME: maybe this is not the best place */
955 c81_data_do_stall(data, dcltol);
956
957 return 0;
958 }
959
960 extern "C" int
c81_data_read_free_format(std::istream & in,c81_data * data,const doublereal dcltol)961 c81_data_read_free_format(std::istream& in, c81_data* data, const doublereal dcltol)
962 {
963 char buf[128]; // 81 should suffice; let's make it 128
964
965 /* header */
966 in.getline(buf, sizeof(buf));
967 if (strncasecmp(buf, "# FREE FORMAT", STRLENOF("# FREE FORMAT")) == 0) {
968 in.getline(buf, sizeof(buf));
969 }
970
971 char *p = std::strchr(buf, ';');
972 if (p == NULL) {
973 return -1;
974 }
975
976 size_t len = STRLENOF(data->header);
977 if (size_t(p - buf) < len) {
978 len = size_t(p - buf);
979 }
980
981 strncpy(data->header, buf, len);
982 data->header[len] = '\0';
983
984 std::istringstream s(++p);
985
986 s >> data->NML;
987 s >> data->NAL;
988 s >> data->NMD;
989 s >> data->NAD;
990 s >> data->NMM;
991 s >> data->NAM;
992
993 if (data->NAM <= 0 || data->NMM <= 0
994 || data->NAD <= 0 || data->NMD <= 0
995 || data->NAL <= 0 || data->NML <= 0) {
996 return -1;
997 }
998
999 /* lift */
1000 data->ml = new doublereal[data->NML];
1001 for (int c = 0; c < data->NML; c++) {
1002 in >> data->ml[c];
1003 }
1004 if (check_vec(data->ml, data->NML)) {
1005 return -1;
1006 }
1007
1008 data->al = new doublereal[(data->NML + 1)*data->NAL];
1009 for (int r = 0; r < data->NAL; r++) {
1010 /* NOTE: "<=" because the number of columns is data->NML + 1
1011 * for the angle of attack */
1012 for (int c = 0; c <= data->NML; c++) {
1013 in >> data->al[r + data->NAL*c];
1014 }
1015 }
1016 if (check_vec(data->al, data->NAL)) {
1017 return -1;
1018 }
1019
1020 /* drag */
1021 data->md = new doublereal[data->NMD];
1022 for (int c = 0; c < data->NMD; c++) {
1023 in >> data->md[c];
1024 }
1025 if (check_vec(data->md, data->NMD)) {
1026 return -1;
1027 }
1028
1029 data->ad = new doublereal[(data->NMD + 1)*data->NAD];
1030 for (int r = 0; r < data->NAD; r++) {
1031 for (int c = 0; c <= data->NMD; c++) {
1032 in >> data->ad[r + data->NAD*c];
1033 }
1034 }
1035 if (check_vec(data->ad, data->NAD)) {
1036 return -1;
1037 }
1038
1039 /* moment */
1040 data->mm = new doublereal[data->NMM];
1041 for (int c = 0; c < data->NMM; c++) {
1042 in >> data->mm[c];
1043 }
1044 if (check_vec(data->mm, data->NMM)) {
1045 return -1;
1046 }
1047
1048 data->am = new doublereal[(data->NMM + 1)*data->NAM];
1049 for (int r = 0; r < data->NAM; r++) {
1050 for (int c = 0; c <= data->NMM; c++) {
1051 in >> data->am[r + data->NAM*c];
1052 }
1053 }
1054 if (check_vec(data->am, data->NAM)) {
1055 return -1;
1056 }
1057
1058 /* FIXME: maybe this is not the best place */
1059 c81_data_do_stall(data, dcltol);
1060
1061 return 0;
1062 }
1063
1064 extern "C" int
c81_data_write_free_format(std::ostream & out,c81_data * data)1065 c81_data_write_free_format(std::ostream& out, c81_data* data)
1066 {
1067 if (data == 0) {
1068 return -1;
1069 }
1070
1071 out << "# FREE FORMAT" << std::endl;
1072
1073 out << data->header << ";"
1074 << " " << data->NML
1075 << " "<< data->NAL
1076 << " "<< data->NMD
1077 << " "<< data->NAD
1078 << " "<< data->NMM
1079 << " "<< data->NAM
1080 << std::endl;
1081
1082 /* lift */
1083 for (int c = 0; c < data->NML; c++) {
1084 out << data->ml[c] << " ";
1085 }
1086 out << std::endl;
1087 for (int r = 0; r < data->NAL; r++) {
1088 /* NOTE: "<=" because the number of columns is data->NML + 1
1089 * for the angle of attack */
1090 for (int c = 0; c <= data->NML; c++) {
1091 out << data->al[r + data->NAL*c] << " ";
1092 }
1093 out << std::endl;
1094 }
1095
1096 /* drag */
1097 for (int c = 0; c < data->NMD; c++) {
1098 out << data->md[c] << " ";
1099 }
1100 out << std::endl;
1101 for (int r = 0; r < data->NAD; r++) {
1102 for (int c = 0; c <= data->NMD; c++) {
1103 out << data->ad[r + data->NAD*c] << " ";
1104 }
1105 out << std::endl;
1106 }
1107
1108 /* moment */
1109 for (int c = 0; c < data->NMM; c++) {
1110 out << data->mm[c] << " ";
1111 }
1112 out << std::endl;
1113 for (int r = 0; r < data->NAM; r++) {
1114 for (int c = 0; c <= data->NMM; c++) {
1115 out << data->am[r + data->NAM*c] << " ";
1116 }
1117 out << std::endl;
1118 }
1119
1120 return 0;
1121 }
1122
1123 /* data array */
1124 static int c81_ndata = 0;
1125 static c81_data** __c81_pdata = NULL;
1126
1127 extern "C" c81_data*
get_c81_data(long int jpro)1128 get_c81_data(long int jpro)
1129 {
1130 if (__c81_pdata == NULL) {
1131 return NULL;
1132 }
1133 return __c81_pdata[jpro];
1134 }
1135
1136 extern "C" int
c81_data_set(long int jpro,c81_data * data)1137 c81_data_set(long int jpro, c81_data* data)
1138 {
1139 if (__c81_pdata == NULL || jpro >= c81_ndata) {
1140 c81_data** pp = NULL;
1141 int ndata = c81_ndata;
1142
1143 if (jpro >= c81_ndata) {
1144 ndata = jpro+1;
1145 }
1146
1147 pp = new c81_data*[ndata];
1148 if (__c81_pdata != NULL) {
1149 for (int i = 0; i < ndata; i++) {
1150 pp[i] = __c81_pdata[i];
1151 }
1152 delete[] __c81_pdata;
1153 }
1154 __c81_pdata = pp;
1155 c81_ndata = ndata;
1156 } else if (__c81_pdata[jpro] != NULL) {
1157 silent_cerr("profile " << jpro << "already defined" << std::endl);
1158 return -1;
1159 }
1160
1161 __c81_pdata[jpro] = data;
1162
1163 return 0;
1164 }
1165
1166 /*
1167 * sistema i dati di stallo
1168 */
1169 static int
do_stall(int NM,int NA,doublereal * a,doublereal * stall,const doublereal dcltol)1170 do_stall(int NM, int NA, doublereal *a, doublereal *stall, const doublereal dcltol)
1171 {
1172 for (int nm = 0; nm < NM; nm++) {
1173 int start = NA*(nm + 1);
1174 int na = NA/2; // mid point
1175
1176 // look for zero
1177 while (a[na] > 0.) {
1178 na--;
1179 if (na <= 0) {
1180 silent_cerr("do_stall: "
1181 "unable to find negative alpha range "
1182 "for Mach #" << nm + 1 << std::endl);
1183 return -1;
1184 }
1185 }
1186
1187 while (a[na] < 0.) {
1188 na++;
1189 if (na >= NA) {
1190 silent_cerr("do_stall: "
1191 "unable to find positive alpha range "
1192 "for Mach #" << nm + 1 << std::endl);
1193 return -1;
1194 }
1195 }
1196
1197 doublereal a0 = a[na];
1198 doublereal dcl0 = a[start + na];
1199
1200 doublereal dcla0 = (a[start + na + 1] - dcl0)/(a[na + 1] - a0);
1201
1202 /* cerca il punto superiore in cui cessa la linearita' */
1203 stall[nm] = 1.;
1204 stall[NM + nm] = 1.;
1205 stall[2*NM + nm] = 0.;
1206
1207 for (int i = na + 2; i < NA; i++) {
1208 doublereal dcla = (a[start + i] - dcl0)/(a[i] - a0);
1209 if (dcla - dcla0 < -dcla0*dcltol) {
1210
1211 /* alpha di stallo superiore */
1212 stall[nm] = a[i - 1];
1213
1214 /* mette temporaneamente il Cp di stallo */
1215 stall[2*NM + nm] = a[start + i - 1];
1216 break;
1217 }
1218 }
1219
1220 dcla0 = (a[start + na - 1] - dcl0)/(a[na - 1] - a0);
1221
1222 /* cerca il punto inferiore in cui cessa la linearita' */
1223 for (int i = na - 2; i >= 0; i--) {
1224 doublereal dcla = (a[start + i] - dcl0)/(a[i] - a0);
1225 if (dcla - dcla0 < -dcla0*dcltol) {
1226
1227 /* stallo inferiore */
1228 stall[NM + nm] = a[i + 1];
1229
1230 /* sottrae il cl allo stallo inferiore */
1231 stall[2*NM + nm] -= a[start + i + 1];
1232
1233 /* divide per il delta alpha */
1234 stall[2*NM + nm] /= (stall[nm] - stall[NM + nm]);
1235 break;
1236 }
1237 }
1238 }
1239
1240 return 0;
1241 }
1242
1243 int
c81_data_do_stall(c81_data * data,const doublereal dcltol)1244 c81_data_do_stall(c81_data *data, const doublereal dcltol)
1245 {
1246 if (data == NULL || data->NML <= 0 || data->NMM <= 0) {
1247 return -1;
1248 }
1249
1250 data->stall = new doublereal[3*data->NML];
1251 do_stall(data->NML, data->NAL, data->al, data->stall, dcltol);
1252
1253 data->mstall = new doublereal[3*data->NMM];
1254 do_stall(data->NMM, data->NAM, data->am, data->mstall, dcltol);
1255
1256 return 0;
1257 }
1258
1259 static int
flip_one(int NM,int NA,doublereal * v,int ss)1260 flip_one(int NM, int NA, doublereal *v, int ss)
1261 {
1262 int s = -1;
1263 for (int m = 0; m <= NM; m++) {
1264 if (m > 0) {
1265 s = ss;
1266 }
1267
1268 for (int a = 0; a < NA/2; a++) {
1269 doublereal d = v[a];
1270 v[a] = s*v[NA - a - 1];
1271 v[NA - a - 1] = s*d;
1272 }
1273
1274 // NA odd: change sign
1275 if (NA%2) {
1276 v[NA/2] *= s;
1277 }
1278
1279 v += NA;
1280 }
1281
1282 return 0;
1283 }
1284
1285 extern "C" int
c81_data_flip(c81_data * data)1286 c81_data_flip(c81_data *data)
1287 {
1288 int rc;
1289
1290 rc = flip_one(data->NML, data->NAL, data->al, -1);
1291 if (rc) {
1292 return rc;
1293 }
1294
1295 rc = flip_one(data->NMD, data->NAD, data->ad, 1);
1296 if (rc) {
1297 return rc;
1298 }
1299
1300 rc = flip_one(data->NMM, data->NAM, data->am, -1);
1301 return rc;
1302 }
1303
1304