## Copyright (C) 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{w} = cusssp_influence (@var{zc}, @var{z}, @var{d})
##
## Compute the complex velocity influence matrix @var{w} for the column
## of points @var{zc} caused by panels on the segments of row of
## complex lengths @var{d} running from row of points @var{z} with unit
## strength. The matrix @var{w} acts on a column of complex panel
## strengths to produce the complex velocities at the points.
##
## @end deftypefn

## ### BACKGROUND THEORY
##
## The complex velocity at zc of a line source at z is
## Q/(2*pi*(zc-z)).
##
## If we have a sheet from z0 to z1 of such line sources of strength q
## dz (for q constant) the resultant complex velocity field is w =
## (q/(2*pi)) * int_{z0}^{z1} dz / (zc-z).  Now d log (zc - z) = -
## dz / (zc - z) so w
##
##  w = -q/(2*pi) * int_{z=z0}^{z=z1} d log (zc - z)
##
##    = -q/(2*pi) * [ log (zc - z)]_{z=z0}^{z=z1}
##
##    = -q/(2*pi) * [ log (zc - z1) - log (zc - z0)]
##
##    = q/(2*pi) * [ log (zc - z0) - log (zc - z1)] .
##
## Now, the branch cut in the log function must be considered.  In
## Octave, it's taken along the negative real axis.  There are no branch
## cuts in the flow if it's continuous, so we can only allow a branch
## cut along the segment of the panel itself.  This can be arranged if
## the branch cuts of the two logarithms cancel along the ray extending
## from the panel segment in one direction; e.g. from z1 to z0.
##
## Consider the simple case where z0 = 0 and z1 = 1 so that
##
##    = q/(2*pi) * (log (zc') - log (zc' - 1)) ,
##
## and both branch cuts extend in the negative x direction.  Then we
## conformally map this flow field to the general panel from z0 to z1.
## The complex potential is
##
##  W = q/(2*pi) * (zc' * log (zc') - (zc' - 1) * log (zc' - 1))
##  ,
##
## and the linear mapping that takes zc' = 0 and zc' = 1 to zc =
## z0 and zc = z1, respectively, is 
##
##  zc' = (zc - z0) ./ (z1 - z0) .
##
## Thus in the zc plane the complex potential is
##
##  W = q/(2*pi) * ((zc-z0)./(z1-z0) * log ((zc-z0)./(z1-z0)) 
##
##               -  ((zc-z0)./(z1-z0) - 1) * log ((zc-z0)./(z1-z0) - 1))
##
##    = q/(2*pi) * ((zc-z0) * log ((zc-z0)./(z1-z0)) 
##
##               -  (zc-z1) * log ((zc-z1)./(z1-z0)) ) / (z1-z0)
##
## and the complex velocity is w = dW/d(zc)
##
##  w = q/(2*pi*(z1-z0)) * (log ((zc-z0)./(z1-z0))
##
##                        - log ((zc-z1)./(z1-z0)) ) .

## Author: G. D. McBain
## Created: 2006-08-11
## Keywords: aerodynamics, ideal flow, complex velocity, panel methods

function w = cusssp_influence (zc, z, d)

  zc = repmat (zc(:), [1, (length (z))]);
  z = repmat (z, [(length (zc)), 1]);
  d = repmat (d, [(length (zc)), 1]);

  w = (log ((zc - z) ./ d) - log ((zc - z - d) ./ d)) ./ (2*pi*d);

endfunction

%!demo#: uniform vortex panel from 0 to 1+1i
%! wf = inline ("1i * reshape (cusssp_influence (z(:), 0, 1+1i), size (z))");
%! qquiver_w (-2-2i, 2+2i, 40, wf, 0.2);
