## Copyright (C) 2005, 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{zc}, @var{omega}, @var{Gamma}] = lvm (@var{z}, @var{alpha})
## Return the collocation points @var{zc}, influence matrix @var{omega},
## and vortex strengths @var{Gamma} for a thin aerofoil with complex
## camber-line @var{z} at incidence @var{alpha} (radians).  The complex
## velocities at the collocation points are recoved with @code{w = exp
## (-1i*alpha) + omega * 1i*Gamma;}.
## @end deftypefn

## ### THEORY
##
## A lumped vortex panel has a vortex at the quarter-point and a
## collocation point at the three-quarter-point, at which impermeability
## is enforced (Katz & Plotkin 2001, pp. 263--272; Moran 2003, pp.
## 101--103; McBain 2005).

## ### REFERENCES
##
## J. Katz & A. Plotkin 2001 Low-Speed Aerodynamics. 2nd edn. Cambridge
## University Press.
##
## J. Moran 2003 An Introduction to Theoretical and Computational
## Aerodynamics. Dover.
##
## G. D. McBain 2006 Course notes for Aerodynamics I,
## http://www.aeromech.usyd.edu.au/~mcbain/AerodynamicsI

## Author: G. D. McBain <geordie.mcbain@aeromech.usyd.edu.au>
## Created: 2005-09-08
## Keywords: aerodynamics, panel methods

function [zc, omega, Gamma] = lvm (z, alpha)

  dz = diff (z); 		# panel segments
  n = length (dz);		# number of panels
  lambda = arg (dz)';		# angle of panels to +ve x axis

  zv = z(1:n) + 0.25*dz;	# vortex locations
  zc = (z(1:n) + 0.75*dz).';	# collocation points

  omega = lvm_influence (zc, zv);
  Gamma = - imag (dmult (exp (1i*lambda) * 1i, omega)) \ sin (lambda-alpha);

endfunction

%!demo# NACA 2400
%! [m, p] = naca4pars (2400);
%! [y, x] = naca4meanline (m, p);
%! [zc, omega, Gamma] = lvm (x'+1i*y', alpha = 0);
%! Cl = 2 * sum (Gamma)
%! Cmc4 = 2 * sum (Gamma .* real (0.25-zc)) * cos (alpha)


%!demo# NACA 2400; convergence
%! [m, p] = naca4pars (2400);
%! printf ("%4s  %7s  %8s\n\n", "n", "Cl", "Cmc4");
%! printf ("(%3d  %7.5f  %- 8.5f)\n\n", Inf, 2*pi*2.077*pi/180, -0.05312);
%! for n = 16:16:256
%!   chi = linspace (pi, 0, n)';
%!   x = (1 + cos (chi))/2;
%!   y = naca4meanline (m, p, x);
%!   [zc, omega, Gamma] = lvm (x'+1i*y', alpha = 0);
%!   Cl = 2 * sum (Gamma);
%!   Cmc4 = 2 * sum (Gamma .* real (0.25-zc)) * cos (alpha);
%!   printf ("%4d  %7.5f  %- 8.5f\n", n, Cl, Cmc4);
%! endfor

%!demo# NACA 2400; surface velocity
%! [m, p] = naca4pars (2400);
%! [y, x] = naca4meanline (m, p);
%! [zc, omega, Gamma] = lvm (x'+1i*y', alpha=(alphad=8)*pi/180);
%! w = exp (-1i*alpha) + omega * 1i*Gamma;
%! [ax,ay] = qquiver (real (zc), imag (zc), real (w), -imag (w), 0.2);
%! xlabel ("x/c"); ylabel ("y/c"); axis ("image");
%! title (sprintf ("NACA 2400 at alpha = %3.1f deg.; surface velocity", alphad))

%!demo# NACA 2400
%! [m, p] = naca4pars (2400);
%! [y, x] = naca4meanline (m, p);
%! for alphad=-10:10
%!   [zc, omega, Gamma] = lvm (x'+1i*y', alphad*pi/180);
%!   Cl = 2 * sum (Gamma);
%!   disp ([alphad, Cl]);
%! endfor
