1% Copyright (C) 2008  VZLU Prague, a.s., Czech Republic
2%
3% Author: Jaroslav Hajek <highegg@gmail.com>
4%
5% This file is part of NLWing2.
6%
7% NLWing2 is free software; you can redistribute it and/or modify
8% it under the terms of the GNU General Public License as published by
9% the Free Software Foundation; either version 3 of the License, or
10% (at your option) any later version.
11%
12% This program is distributed in the hope that it will be useful,
13% but WITHOUT ANY WARRANTY; without even the implied warranty of
14% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15% GNU General Public License for more details.
16%
17% You should have received a copy of the GNU General Public License
18% along with this software; see the file COPYING.  If not, see
19% <http://www.gnu.org/licenses/>.
20%
21
22% -*- texinfo -*-
23% @deftypefn{Function File} {wing =} makewing (acs, pols, ref, panels)
24% Creates the wing structure necessary for further computations.
25% @var{acs} is an N-by-5 array specifying the spanwise geometry description.
26% each row contains @code{[zac xac yac chord twist]}
27% pols is a struct array describing the spanwise wing section data
28% distribution. @code{pols(i).z} is the spanwise coordinate, @code{pols(i).cl}
29% is the lift coefficient on local angle of attack dependence, etc.
30% @var{ref} contains the reference quantities.
31% @var{panels} specifies either an approximate number of panels,
32% or directly the z-coordinates of panel vertices. In the latter case,
33% @var{panels}(1) and @var{panels}(end) should match @var{acs}(1,1) and
34% @var{acs}(end,1), respectively.
35% @end deftypefn
36function wing = makewing (ac, pols, ref, np = 80, zac = [])
37
38  ozac = ac(:,1);
39
40  if (isempty (zac))
41    % distribute points
42    if (ref.sym)
43      zmx = ozac(end);
44      fi = linspace (0, pi/2, np+1).';
45      zac = zmx * sin (fi);
46    else
47      zmxl = ozac(1); zmxu = ozac(end);
48      fi = linspace (0, pi/2, np/2+1)';
49      zac = [zmxl * sin(fi(end:-1:2)); zmxu * sin(fi)];
50    endif
51  endif
52
53  wing.zac = zac;
54  aci = interp1 (ozac, ac(:,2:5), zac, "pchip");
55  wing.xac = aci(:,1);
56  wing.yac = aci(:,2);
57  wing.cac = aci(:,3);
58
59  m2 = @(v) (v(1:end-1)+v(2:end))/2;
60  wing.zc = zc = m2 (zac);
61  wing.ch = m2 (aci(:,3));
62  wing.twc = pi/180 * m2 (aci(:,4));
63
64  zpol = [pols.z];
65  if (! issorted (zpol))
66    [zpol, isrt] = sort (zpol);
67    pols = pols(isrt);
68  endif
69
70  % set jj so that zc(jj(i)-1) <= zpol(i) < zc(jj(i))
71  % which, in turn, means
72  % zpol(i) < zc(jj(i)) <= zc(jj(i+1)-1) <= zpol(i+1)
73  jj = lookup (zc, zpol) + 1;
74  % correct boundary values (to include all of zc)
75  jj(1) = 1;
76  jj(end) = length (zc) + 1;
77
78  wing.pidx = jj;
79  wing.pol = pols;
80  wing.np = np;
81
82  wing.a0 = zeros (np, 1);
83  wing.amax = zeros (np, 1);
84  wing.clmax = zeros (np, 1);
85  wing.cf = zeros (np, 1);
86
87  % FIXME: 3.3.50+ will handle discontinuous interpolation.
88  for i=1:length (jj)-1
89    jl = jj(i); ju = jj(i+1)-1;
90    if (jl <= ju)
91      lcf = (zc(jl:ju) - zpol(i)) / (zpol(i+1) - zpol(i));
92      wing.cf(jl:ju) = lcf;
93      wing.a0(jl:ju) = (1-lcf) * pols(i).a0 + lcf * pols(i+1).a0 + wing.twc(jl:ju);
94      wing.amax(jl:ju) = (1-lcf) * pols(i).amax + lcf * pols(i+1).amax + wing.twc(jl:ju);
95      wing.clmax(jl:ju) = (1-lcf) * pols(i).clmax + lcf * pols(i+1).clmax;
96    endif
97  endfor
98
99  for [val,key] = ref
100    wing.(key) = val;
101  endfor
102  dS = wing.ch .* diff (wing.zac);
103  if (! isfield (wing, 'sym'))
104    wing.sym = true;
105  endif
106  if (! isfield (wing, 'area'))
107    area = wing.area = sum (dS);
108    if (wing.sym)
109      wing.area *= 2;
110    endif
111  else
112    area = wing.area / 2;
113  endif
114  if (! isfield (wing, 'cmac'))
115    wing.cmac = sum (dS .* wing.ch) / area;
116  endif
117  if (! isfield (wing, 'span'))
118    wing.span = wing.zac(end) - wing.zac(1);
119    if (wing.sym)
120      wing.span *= 2;
121    endif
122  endif
123  xac = wing.xac; xac = (xac(1:end-1) + xac(2:end))/2;
124  yac = wing.yac; yac = (yac(1:end-1) + yac(2:end))/2;
125  if (! isfield (wing, 'xmac'))
126    wing.xmac = dot (dS, xac) / area;
127    wing.ymac = dot (dS, yac) / area;
128  endif
129  dxmac = xac - wing.xmac;
130  dymac = yac - wing.ymac;
131  wing.rmac = hypot (dymac, dxmac);
132  wing.amac = atan2 (dymac, dxmac);
133endfunction
134