## Copyright (C) 2004, 2006, 2007 G. D. McBain <geordie_mcbain@yahoo.com.au>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA


## -*- texinfo -*-
## @deftypefn {Function File} {} @var{A} = lline (@var{c})
## Compute the expansion coefficients (per unit radian of absolute
## incidence) of the spanwise lift loading for an untwisted wing with
## chord @var{c} at eccentric angle @code{(1:(r-1))*pi/r}, where @code{r
## = length (@var{c})}.
## @end deftypefn

## ### REFERENCES
##
## J. J. Bertin 2002 Aerodynamics for Engineers, 4th edn, Prentice Hall,
## Upper Saddle River, New Jersey.
##
## A. M. Kuethe & C.-Y. Chow 1998 Foundations of Aerodynamics, 5th edn,
## Wiley, New York.

## Author: G. D. McBain
## Created: 2003
## Keywords: aerodynamics, lifting line

function A = lline (c)

  r = (length (c)) + 1;         # degree of approximation
  n = 1:(r-1);                  # row vector of indices
  S = sin (n'*n*pi/r);          # S(i,j) = sin(i*j*pi/r)
  
  B = ((4 ./ (repmat(2*pi*c, 1, r-1))) + (1 ./ (S(:,1))) * n)  .* S;  

  A = B \ ones (r-1, 1);

endfunction

%!demo # Kuethe & Chow (1998, p. 188), ans.: CL = 4.5273, CD = 1.1378
%! r = 8; AR = 6; c = (ones (r-1, 1)) / AR;
%! A = lline (c);
%! CL = pi * AR * A(1)
%! CD = pi * AR * (1:(r-1)) * A.^2

%!demo # Bertin (2002, pp. 249--253); ans.: CL = 0.4654, CD = 0.00776
%! r = 8;
%! AR = 9;
%! taper = 0.4;
%! theta = (1:(r-1))' * pi/r;
%! c = (2 / (AR * (1+taper))) * (1 - (1-taper)*(abs (cos (theta))));
%! A = ((4 - (-1.2)) * (pi/180)) * lline (c);
%! CL = pi * AR * A(1)
%! CD = pi * AR * (1:(r-1)) * A.^2
