MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
c81data.cc
Go to the documentation of this file.
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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*
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
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
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
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
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
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 
get_c81_vec
static int get_c81_vec(std::istream &in, doublereal *v, int ncols)
Definition: c81data.cc:142
c81_data::NMD
int NMD
Definition: aerodc81.h:160
read_fc511_row
int read_fc511_row(std::istream &in, doublereal *d, int NC)
Definition: c81data.cc:642
c81_data::NAM
int NAM
Definition: aerodc81.h:166
read_fc511_mat
int read_fc511_mat(std::istream &in, doublereal *d, int NR, int NC)
Definition: c81data.cc:676
c81_data_do_stall
int c81_data_do_stall(c81_data *data, const doublereal dcltol)
Definition: c81data.cc:1244
c81_data_destroy
void c81_data_destroy(c81_data *data)
Definition: c81data.cc:276
c81_data::NMM
int NMM
Definition: aerodc81.h:165
get_double
static int get_double(const char *const from, doublereal &d)
Definition: c81data.cc:119
c81_data_write_free_format
int c81_data_write_free_format(std::ostream &out, c81_data *data)
Definition: c81data.cc:1065
c81_data_flip
int c81_data_flip(c81_data *data)
Definition: c81data.cc:1286
c81_data::NML
int NML
Definition: aerodc81.h:146
c81_data::mm
doublereal * mm
Definition: aerodc81.h:167
get_int
static int get_int(const char *const from, int &i)
Definition: c81data.cc:71
c81_data::header
char header[31]
Definition: aerodc81.h:144
c81_data
Definition: aerodc81.h:143
c81_data_read_free_format
int c81_data_read_free_format(std::istream &in, c81_data *data, const doublereal dcltol)
Definition: c81data.cc:961
c81_data_read
int c81_data_read(std::istream &in, c81_data *data, const doublereal dcltol, int *ff)
Definition: c81data.cc:292
c81_data::NAL
int NAL
Definition: aerodc81.h:147
put_c81_row
static int put_c81_row(std::ostream &out, doublereal *v, int dim, int ncols, int first=0)
Definition: c81data.cc:220
put_c81_vec
static int put_c81_vec(std::ostream &out, doublereal *v, int nrows)
Definition: c81data.cc:250
c81_data::ad
doublereal * ad
Definition: aerodc81.h:163
c81_data_merge
int c81_data_merge(unsigned ndata, const c81_data **data, const doublereal *upper_bounds, doublereal dCsi, doublereal dcltol, c81_data *i_data)
Definition: c81data.cc:396
c81_data::mstall
doublereal * mstall
Definition: aerodc81.h:158
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
get_c81_mat
static int get_c81_mat(std::istream &in, doublereal *m, int nrows, int ncols)
Definition: c81data.cc:183
flip_one
static int flip_one(int NM, int NA, doublereal *v, int ss)
Definition: c81data.cc:1260
for
for(int i=0;i< 3;++i) for(int j=0
c81_ndata
static int c81_ndata
Definition: c81data.cc:1124
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
c81_data_fc511_read
int c81_data_fc511_read(std::istream &in, c81_data *data, const doublereal dcltol)
Definition: c81data.cc:810
put_c81_mat
static int put_c81_mat(std::ostream &out, doublereal *m, int nrows, int ncols)
Definition: c81data.cc:262
STRLENOF
#define STRLENOF(s)
Definition: mbdyn.h:166
get_c81_data
c81_data * get_c81_data(long int jpro)
Definition: c81data.cc:1128
c81_data::stall
doublereal * stall
Definition: aerodc81.h:157
c81_data::am
doublereal * am
Definition: aerodc81.h:168
c81_data_nrel_read
int c81_data_nrel_read(std::istream &in, c81_data *data, const doublereal dcltol)
Definition: c81data.cc:874
check_vec
static int check_vec(doublereal *v, int n)
Definition: c81data.cc:171
c81_data::al
doublereal * al
Definition: aerodc81.h:149
read_fc511_block
int read_fc511_block(std::istream &in, int &NA, int &NM, doublereal *&da, doublereal *&dm)
Definition: c81data.cc:712
a
static const doublereal a
Definition: hfluid_.h:289
__c81_pdata
static c81_data ** __c81_pdata
Definition: c81data.cc:1125
buf
static doublereal buf[BUFSIZE]
Definition: discctrl.cc:333
c81_data_set
int c81_data_set(long int jpro, c81_data *data)
Definition: c81data.cc:1137
doublereal
double doublereal
Definition: colamd.c:52
c81_data_write
int c81_data_write(std::ostream &out, c81_data *data)
Definition: c81data.cc:610
c81_data::md
doublereal * md
Definition: aerodc81.h:162
c81_data::ml
doublereal * ml
Definition: aerodc81.h:148
c81_data::NAD
int NAD
Definition: aerodc81.h:161
do_stall
static int do_stall(int NM, int NA, doublereal *a, doublereal *stall, const doublereal dcltol)
Definition: c81data.cc:1170