% Copyright (c) 2015, Lennon O'Naraigh
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
% 
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer 
% in the documentation and/or other materials provided with the distribution.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 
% OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
% OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function [z1,phi1,z2,phi2,z,phi,eeta,eeta1,eeta2]=get_streamfunction0(vec)

[n_height,~,~,~,N_1,N_21]=fix_all_parameters();

nn=N_1+N_21+2;
vv=vec(1:nn);
ww=vec(nn+(1:nn));


%**************************************************************************
%  Define streamfunction at collocation gridpoints

n1=50;
kk1=1:(n1-3);
eta1=cos((pi/(n1-4))*(kk1-1));

n2=50;
kk2=1:(n2-3);
eta2=cos((pi/(n2-4))*(kk2-1));

eta1=flipdim(eta1,2);
z1 =(1/2)*(eta1-1);

eta2=flipdim(eta2,2);
z2=(n_height/2)*(eta2+1);

%**************************************************************************

phi1=0*z1;
eeta1=0*z1;

for i=1:length(z1)
    
    sum_val=0;
    rum_val=0;
    
    for j=1:N_1+1
        sum_val=sum_val+vv(j)*T(eta1(i),j-1);
        rum_val=rum_val+ww(j)*T(eta1(i),j-1);
    end
    
    phi1(i)=sum_val;
    eeta1(i)=rum_val;
end

phi2=0*z2;
eeta2=0*z2;

for i=1:length(z2)
    
    sum_val=0;
    rum_val=0;
    
    for j=1:N_21+1
        sum_val=sum_val+vv(j+N_1+1)*T(eta2(i),j-1);
        rum_val=rum_val+ww(j+N_1+1)*T(eta2(i),j-1);
    end
    
    phi2(i)=sum_val;
    eeta2(i)=rum_val;
end

z=[z1,z2];
phi=[phi1,phi2];
eeta=[eeta1,eeta2];



end


%**************************************************************************
%**************************************************************************
%**************************************************************************
%**************************************************************************

function t = T(z,k)

t=cos(k*acos(z));

end

