function [cps,cls,cms]=thickpot(xp,yp,alpha)
%% The function thickpot calculates a potential flow around 2D airfoils.
%% The condition of parallell flow is fulfilled in the middle points
%% of the segments.
%% The algorithm is from a report by L. E. Eriksson,
%% Calculation of two-dimensional potential flow wall interference for
%% multicomponent airfoils in closed low speed windtunnels.
%% FFA AU-1116, 1975.
%% I came across this in a course in aerodynamics at KTH taught by
%% Sven-Olof Ridder and Sven Hedman in 1990.
%% This implementation is made from a description in Hedman's course notes.
%% It takes only one profile although the method can take several.
%% The original formulas are considerable simplified with MATLAB's capability
%% for complex numbers.
%%
%% Copyright (C) 1994 Anders Lennartsson
%% Error Reports in case of UNCHANGED versions to
%% Anders Lennartsson, abel@mech.kth.se
%%
%% Coordinate system:
%% The origin is at the nose and the x-axis points along the line to the
%% trailing edge.
%%
%% Arguments
%% xp    x-coordinates in row vector for points 1, 2,..., N+1
%% yp    y-coordinates in row vector for points 1, 2,..., N+1
%%       (the point N+1 should be identical to point 1)
%% alpha a vector of angle of attacks relative to freestream
%%       obviously a single value will do
%%
%% initialization of variables
twopi=2*pi;
% n is the number of end points, first and last identified separately
n=size(xp,2);
% n1 is the number of segments
n1=n-1;
% make calculation for each alpha
na=max(size(alpha));
for jj=1:na
% first find coordinates relative to freestream, inclined by alpha
xv=xp*cos(alpha(jj))+yp*sin(alpha(jj));
yv=-sin(alpha(jj))*xp+cos(alpha(jj))*yp;
% differentials over segments in v system
dxv=xv(1,2:n)-xv(1,1:n1);
dyv=yv(1,2:n)-yv(1,1:n1);
% angles for local coordinate system in segments
thi=atan2(dyv,dxv);
si=sin(thi);
ci=cos(thi);
% length of segments
di=sqrt(dxv.*dxv+dyv.*dyv);
% midpoints of segments in v system
xmv=(xv(1,2:n)+xv(1,1:n1))/2;
ymv=(yv(1,2:n)+yv(1,1:n1))/2;
% dximv from origin i to midpointvector x
% dyimv from origin i to midpointvector y
% aij   incidence matrix
aij=zeros(n1);
% distances to middle points from first segment origin in vind system
dximv=xmv-xv(1);
dyimv=ymv-yv(1);
% in first segment system
dximi=dximv*ci(1)+dyimv*si(1);
dyimi=-si(1)*dximv+dyimv*ci(1);
% going complex
imagxi=dximi+i*dyimi;
imagveli=(1+(1-imagxi/di(1)).*(log(imagxi)-log(imagxi-di(1))))*i/twopi;
% extracting velocity components
ui=real(imagveli);
vi=-imag(imagveli);
% transforming to vind system
uiv=ui*ci(1)-vi*si(1);
viv=ui*si(1)+vi*ci(1);
% projecting on normal to each segment
vij=-uiv.*si+viv.*ci;
% influence from first segment of gamma1
aij(:,1)=vij';
for ii=2:n1,
  % coordinates from segment ii to middle points in vind system
  dximv=xmv-xv(ii);
  dyimv=ymv-yv(ii);
  % tranforming to ii system
  dximi=dximv*ci(ii)+dyimv*si(ii);
  dyimi=-si(ii)*dximv+dyimv*ci(ii);
  % moving previous ii coordinates to ii-1
  imagxiold=imagxi;
  % going complex with new coordinates
  imagxi=dximi+i*dyimi;
  % induced velocities from segment ii in middle points
  imagveli=(1+(1-imagxi/di(ii)).*(log(imagxi)-log(imagxi-di(ii))))*i/twopi;
  % induced velocitis from segment ii-1 in middle points
  imagveliold=(-1+(imagxiold/di(ii-1)).*(log(imagxiold)-log(imagxiold-di(ii-1))))*i/twopi;
  % adding influences and extracting velocity components
  ui=real(imagveli);
  vi=-imag(imagveli);
  uiold=real(imagveliold);
  viold=-imag(imagveliold);
  % transforming to vind system and summing total influence
  uiv=ui*ci(ii)-vi*si(ii)+uiold*ci(ii-1)-viold*si(ii-1);
  viv=ui*si(ii)+vi*ci(ii)+uiold*si(ii-1)+viold*ci(ii-1);
  % projecting on normal to each segment
  vij=-uiv.*si+viv.*ci;
  aij(:,ii)=vij';
end
% imagxi has complex coordinates for last segment
% calculate influence from -gamma1
ii=n1 %%%%%MODIFICATION FOR MATLAB TO SCILAB TRANSLATOR%%%%%%%%
imagveli=(-1+(imagxiold/di(ii-1)).*(log(imagxiold)-log(imagxiold-di(ii-1))))*i/twopi;
ui=real(imagveli);
vi=-imag(imagveli);
% transforming to vind system and summing total influence
uiv=ui*ci(n1)-vi*si(n1);
viv=ui*si(n1)+vi*ci(n1);
% projecting on normal to each segment
vij=-uiv.*si+viv.*ci;
% subtracting influence from first column of aij
aij(:,1)=aij(:,1)-vij';
guinf=aij\si';
cp=1-guinf.*guinf;
cp=[ cp(1:n1,1) ; cp(1,1) ]';
mcp=(cp(1,2:n)+cp(1,1:n1)).*di/2;
cl=sum(mcp.*ci);
cm=-sum(mcp.*ci.*xmv)-sum(mcp.*si.*ymv);

%%%%%MODIFICATION FOR MATLAB TO SCILAB TRANSLATOR%%%%
cps(jj,1:size(cp,2))=cp; %cps(jj,:)=cp;
cls(jj)=cl;
cms(jj)=cm;
end % for loop in alpha
