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