function [Fn,Fn0,err,nFnDeln,tharr,curvi,slipar] = ...
            DensSwelVelLowRiJac(zeti,Ric,ipt,ee,Fn0,...
                Tolr,Maxiter,flplot,resol)

% Solves collocation equations for the flow above a Rankine-half body (RHB)
% swollen in response to an increase in the density of source fluid. Uses
% the collocation method of Saffman and Yuen. Outputs arrays of curvilinear
% coordinate and slip on interface, as well as relative errors. Explicitly
% uses the Jacobian to hasten fsolve convergence.
%
% Michel Louge and Cian Carroll 3/12/2012.
%
%	input	-	zeti    =   1-rho/rho'
%                           where   rho  = ambient fluid density
%                                   rho' = source fluid density > rho
%           -   Ric     =   Richardson number 2 zet g H'/rho' U'^2
%                           where   g    = gravitational acceleration
%                                   H'   = tail height
%                                   U'   = tail speed
%           -   ipt    	=   # of collocation pts on swollen interface
%           -   ee      =   exponent to skew theta array toward low angles
%                           set to 1 for no bias; any higher exponent skews
%           -   Fn0     =   column vector first guess for Fn of size 2*ipt
%                           if scalar > 0, guess for all Fn0 components
%                           if scalar < 0, Fn0 found by linearization
%           -   Tolr    =   tolerance
%           -   Maxiter =   maximum number of iterations
%           -   flplot  =   if =1, plot results; if not, do not
%           -   resol   =   increased resolution for plot
%                               = 1 => same resolution
%                               = 2 => twice resolution
%                               = 3 => three times, etc...
%
%	output	-	Fn      =   column vector of size 2*ipt of coefficient Fn
%                           defining the streamfunction in ambient fluid
%           -   Fn0     =   column vector first guess for Fn
%           -   err     =   column vector of errors
%                           for Bernoulli eqs: relative to stagnation pres
%                           for streamline eqs: dim(-) distance of
%                           erroneous streamline to interface streamline
%           -   nFnDeln =   sum(n*Fn*delta^n) (should be zero)
%           -   tharr   =   array of polar angle for plots
%           -   curvi   =   array of curvilinear coordinates relative to b'
%           -   slipar  =   array of dim(-) interface slip
%
%   global variables
%   
%           -   bet   	=   gb/U^2
%                           where   b = gamma/(2 pi U) (m)
%                                   U = uniform upstream flow speed (m/s)
%                                   gamma = source strength (m^2/s)
%           -   delt  	=   swelling coefficient    with U'= delt U
%                                                   and  b'= b/delt
%           -   thei    =   col vec of polar angle of interf pts (size ipt)
%                           (does not include theta = 0 and theta = pi)
%           -   sini    =   col vector sin(thei)
%           -   cosi    =   col vector cos(thei)
%           -   usip2   =   col vector = half the norm squared of speed at
%                           collocation pts on interface in source fluid
%           -   ri      =   col vec of rad coord of interf pts (size ipt)
%           -   ndeln   =   row vec of  n delt^n
%           -   nsini   =   (i,n) mat of  [n sin(n thei)/ri^(n+1)]
%                           where   i = twice the number of collocation pts
%                                       (i.e. two eqs per colloc pt)
%                                   n = number of terms in the series sol
%           -   ncosi   =   (i,n)  mat of [n cos(n thei)/ri^(n+1)]
%           -   ideln   =   (i,n) matr of n*delt^n rows repeated ipt times
%           -   nsinr   =   (i,n)  matrix of [sin(n thei)/ri^n]
%           -   zet     =   zeta = 1-rho/rho'
%           -   onen    =   1x2ipt row vector of ones
%
% Examples, see "MatlabLogCollocation.txt" 
% cd ~/Documents/Michel_files/Papers/CianThesis/CSC_thesis/matlab

global bet delt thei sini cosi usip2 ri ...
        ndeln nsini ncosi ideln nsinr zet onen

%
%   internal variables
%
%           -   narr  	=   column vect of indices of Fn unknowns
%           -   nn      =   index along n for loop
%           -   ii      =   index along i for loop
%           -   rin     =   (i,n) matrix of ri^n
%           -   rin1    =   (i,n) matrix of ri^(n+1)
%           -   narr    =   row vector of indices (1;2;...;n;...) of Fn
%           -   Fn0     =   column vector first guess for Fn of size 2*ipt
%           -   options =   options for fsolve algorithm
%           -   optnew  =   options for fsolve algorithm
%           -   sqti    =   col vector at each colloc pt, see below
%           -   num     =   # of points for table of slip vs angles
%           -   span    =   array of polar angle to calculate curvil coord
%           -   rarr    =   matrix (num,2*ipt) of 1/radii^(n+1)
%           -   nsinar  =   matrix (num,2*ipt)
%           -   ncosar  =   matrix (num,2*ipt)
%           -   sinrat  =   col vect to avoid blowup of sin(t)/t at t=0
%           -   resid   =   col vec of residuals of 2*ipt eqs (size 2*ipt)
%           -   flagFN0 =   flag denoting how Fn0 is known
%
%   variables to linearize the problem (only works
%   at low delta; all statements involving these variables are commented out)
%           -   Ain     =   (i,n) matrix in linearized Bernoulli eqs
%           -   sinn    =   (i,n) matrix, see below
%           -   sinn1   =   (i,n) matrix, see below
%           -   Ci      =   col vect, see below
%           -   Din     =   (i,n) matrix, see below
%           -   ADnn    =   (n,n) matrix, see below
%           -   Ci0     =   (n,1) col vect, see below
%           -   ndelnFn =   (n,1) col vector of n delt^n Fn

%

% exit the program if parameters are out of correct range:

if ipt < 2
    error('# points < 2')
end

if zeti < 0 || zeti == 0
    error('zeta <=0')
end

if Ric < 0 || Ric == 0
    error('Richardson <=0')
end

if floor(resol) ~= resol || resol == 0 || resol < 0
    error('resolution must be a positive integer')
end

% set global variables to pass to functions

zet     =   zeti                ; % initialize zeta as a global variable

bet     =   Ric/(1+Ric)^(3/2)*(1-zet)^(3/2)/2/pi/zet	; % calculate beta
delt 	=   sqrt((1-zet)/(1+Ric))                     	; % calculate delta

% for efficient computations, calculate arrays computed in the function
% called by fsolve and pass to function as global variables:

% column vectors on interface collocation points and corresponding
% quantities at those points; bias toward low angles with exponent ee:

thei    =   (1:1:ipt)'.^ee.*pi./(ipt+1)^ee   	; % col vec of polar angles
sini    =   sin(thei)                         	; % col vec of sin(pol ang)
cosi    =   cos(thei)                         	; % col vec of cos(pol ang)
sqti    =   sqrt(1+(sini./thei).^2- ...
                    2.*sini.*cosi./thei)        ; % col vec of sqrt(Bi)

usip2   =   delt^2/2.*sqti.^2                   ; % col vec of (1/2)U'^2
ri      =   thei./sini./delt                  	; % col vec of radial coord

% calculate (i,n) matrices for efficient computation: 
% for each equation of index i or ipt+i, quantities calculated at 
% interface collocation points i:

narr    =   (1:1:2*ipt)         ; % row vect of indices of Fn unknowns

ndeln   =   (narr.*delt.^narr)	; % row vec  of n*delt^n
onen    =   ones(1,2*ipt)       ; % row vector of ones

rin     =   zeros(ipt,2*ipt)    ; % preallocation of rin for efficiency
rin1    =   zeros(ipt,2*ipt)    ; % preallocation of rin1 for efficiency
sinn    =   zeros(ipt,2*ipt)    ; % preallocation of sinn for efficiency
sinn1   =   zeros(ipt,2*ipt)    ; % preallocation of sinn1 for efficiency
    for nn=1:2*ipt
rin(:,nn)   =   ri.^nn          ; % build matrix (i,n) of ri^n
rin1(:,nn) 	=   ri.^(nn+1)     	; % build matrix (i,n) of ri^(n+1)
sinn(:,nn) =   (sin(thei)./...
                  thei).^nn     ; % mat (i,n) of (sin(thei)/thei))^n
sinn1(:,nn) =   (sin(thei)./...
                  thei).^(nn+1) ; % mat (i,n) of (sin(thei)/thei))^(n+1)
    end
    
nsini   =   (ones(ipt,1)*narr).*...
            sin(thei*narr)./...
        	rin1                ; %  (i,n) matrix [n sin(n thei)/ri^(n+1)]

ncosi   =   (ones(ipt,1)*narr).*...
            cos(thei*narr)./...
        	rin1                ; %  (i,n) matrix [n cos(n thei)/ri^(n+1)]
        
ideln   =   thei*ndeln          ; %  (i,n) matrix of thei*n*delt^n

nsinr   =	sin(thei*narr)./rin	; % (i,n) matrix [sin(n thei)/ri^n]

if size(Fn0,1) == 1
    
    if Fn0 > 0 || Fn0 == 0 % if Fn0 is scalar >= 0, guess for Fn0 vect
        
flagFN0 =   1   ; % uniform guess for Fn0
Fn0     =   Fn0.*ones(2*ipt,1)	; % first guess Fn0 = Fn0*(1;1;...;1)

    else       % if Fn0 is scalar < 0, linearize the pb to find Fn0 vector
        
%
% find a first estimate of Fn from linearized equations:
% Do not use this approach if delta is too large: because Fn ~ ri^n and
% 1/delta < ri^n < (N/2/delta)^N, Fn can be very large at large n...
%
flagFN0 =   2   ; % Fn0 found by linearization

Ain     =   ((sini./thei-cosi)*ones(1,2*ipt)).*...
             ( sinn1.*cos(thei*narr) - (sini./thei)*ones(1,2*ipt) ) + ...
               (sini*ones(1,2*ipt)).*...
                    sinn1.*sin(thei*narr)	; % (i,n) matrix
Ci      =   (1/2)*Ric/(1+Ric).*...
                (thei./pi - sqti.^2)    	; % (1,i) column vector
Din     =   (ones(ipt,1)*(1./narr)).*...
                sinn.*sin(thei*narr) - ...
                    (thei*ones(1,2*ipt))  	; % (i,n) matrix
ADnn    =   [delt.*Ain;Din]     ; % build (n,n) matrix to invert
Ci0     =   [Ci;zeros(ipt,1)]   ; % build (1,n) col vector
ndelnFn =   ADnn\Ci0            ; % invert matrix to find n delt^n Fn
Fn0     =   ndelnFn./narr'./(delt.^(narr')) ; % linear estimate of sol

    end
        
else
    
flagFN0 =   3   ; % Fn0 input

end
    
% call the non-linear solver:

%options =   optimset('Display','iter')  ; % Option to display output

maxiterations   =   Maxiter ; % maximum number of iterations
maxfunctions    =   Maxiter ; % max numb of function evaluations
%maxiterations   =   50000 ; % maximum number of iterations
%maxfunctions    =   50000 ; % max numb of function evaluations
%Tolerance       =   1e-3  ; % termination toler. on funct (ie, on eqn = 0)
%ToleranceX      =   1e-3  ; % termination tolerance on Fn
Tolerance       =   Tolr  ; % termination toler. on funct (ie, on eqn = 0)
ToleranceX      =   Tolr  ; % termination tolerance on Fn

options = optimset('fsolve')	; %options for optimization toolbox
         
optnew  =   optimset(options,'MaxIter',maxiterations,'MaxFunEvals',...
             maxfunctions,'TolFun',Tolerance, 'TolX',ToleranceX,...
             'Jacobian','on') ;

%             'PrecondBandWidth',inf,'LargeScale','off');
%             'Algorithm','trust-region-dogleg',...
%             'Algorithm','trust-region-reflective',...
%             'Algorithm','levenberg-marquardt',...

[Fn] = fsolve(@colloc,Fn0,optnew)   ; % Calls solver

[resid] = colloc(Fn)                ; % find residuals
%
% Residuals in Bernoulli eqs are (p-p')/rho' U^2
%
% Residuals in streamline eqs are approx dim(-) distance
% from actual streamline times [delta sin(thei)^2/sqti]
%

% evaluate solution's accuracy

% calculates sum(n*Fn*delta^n) -- it should vanish
nFnDeln =   (narr.*delt.^narr)*Fn   ;

% errors in Bernoulli Eqs. evaluated relative to stagnation pressure:
err(1:ipt,1)        = resid(1:ipt).*2/(1-zet)                   ;

% errors in streamline Eqs. measured as distance of erroneous streamline
% from the interface streamline:
err(ipt+1:2*ipt,1)  = - resid(ipt+1:2*ipt)./delt.*sqti./sini.^2	;

% Calculate curvilinear coordinate on interface for plotting solution:

span    =   linspace(0,thei(end),resol*ipt+1)	;
num     =   size(span,2)    ; % number of points in integration span

Y0      =   span(1)         ; % initial polar angle

[tharr,curvi]	= ode45(@curvider,span,Y0) ; % ODE solver calling du2dx

% build (num,2*ipt) matrices to calculate slip velocity:

sinrat  =   [1;sin(tharr(2:end))./tharr(2:end)] ; % avoids blowup sin(t)/t

rarr    =   zeros(num,2*ipt)    ; % preallocation of rarr for efficiency
    for nn=1:2*ipt
rarr(:,nn) 	=   (delt*sinrat).^(nn+1)	; % build matrix of 1/r^{n+1)
    end
    
nsinar	=   (ones(num,1)*narr).*...
            sin(tharr*narr).*...
        	rarr                ; % (num,2*ipt) matrx [n sin(n th)/r^(n+1)]

ncosar	=   (ones(num,1)*narr).*...
            cos(tharr*narr).*...
        	rarr                ; % (num,2*ipt) matrx [n cos(n th)/r^(n+1)]

slipar  =   sqrt( (sinrat*(1-delt*nFnDeln)-cos(tharr)+ncosar*Fn).^2  +...
                  (sin(tharr)+nsinar*Fn).^2       ) -...
                    delt.*sqrt(sinrat.^2-2.*sinrat.*cos(tharr)+1)	;               
            
%x   =   [1;tharr(2:end)./tan(tharr(2:end))]/delt	; % distance if needed
        
if flplot == 1

    scrsz   = get(0,'ScreenSize')	; % find screen size of current machine
figure('Position',...
 [1 scrsz(4) scrsz(3) scrsz(4)/3])  ; % reserves screen space for figure

subplot(1,5,1)	% plot slip velocity vs curvilinear coordinate 
plot(curvi,slipar,'o',max(curvi)*[.5 1.5],...
            (1-delt)*[1 1],'--') % to plot a linear plot instead
% show the asymptote in the tail
%semilogx(curvi(2:end),slipar(2:end),'o')
%title('...')
xlabel('curvilinear coord (-)')
ylabel('slip velocity (-)')
axis([0 max(curvi)*1.5 0 1.1*max(slipar)])

subplot(1,5,2)	% plot  Fn's vs n to check order of magnitude
loglog((1:1:2*ipt),abs(Fn),'o')
xlabel('log(n)')
ylabel('log(Fn) ')

subplot(1,5,3)	% plot relative errors vs x for Bernoulli eqs
plot(thei./delt./sini.*cosi,err(1:ipt),'o')
xlabel('x (-)')
ylabel('error in Bernoulli eqs relative to stag pressure')

subplot(1,5,4)	% plot errors vs x for streamline eqs
plot(thei./delt./sini.*cosi,err(ipt+1:2*ipt),'o')
xlabel('x (-)')
ylabel('distance error in streamline eqs relative to stag distance')

subplot(1,5,5)	% subplot for displaying parameters
axis off

text(0.1,.0,['zeta = ',num2str(zet)])
text(0.1,.05,['Ri = ',num2str(Ric)])
text(0.1,.1,['delta = ',num2str(delt)])
text(0.1,.15,['beta = ',num2str(bet)])
text(0.1,.3,['I = ',num2str(ipt)])
text(0.1,.35,['N = ',num2str(2*ipt)])
text(0.1,.4,['plot resolution = ',num2str(resol)])

text(0.1,.55,['bias exponent = ',num2str(ee)])
text(0.1,.65,['n*Fn*Delta^n = ',num2str(nFnDeln),' should vanish'])
text(0.1,.7,['maximum relative equation error = ',...
        num2str(100*max(abs(err))),' %'])

if flagFN0 == 1 % uniform guess for Fn0 vect
    
text(0.1,.8,['uniform guess for Fn0 = ',num2str(Fn0(1))])

else       
        
if flagFN0 == 2 % linearize the pb to find Fn0 vector
    
    text(0.1,.8,['stdev absolute error = ',...
        num2str(std((Fn-Fn0)./(Fn+Fn0).*200)),'%'])
    text(0.1,.85,['mean relat error 2|Fn-Fn0|/|Fn+Fn0| = ',...
        num2str(mean(abs((Fn-Fn0)./(Fn+Fn0).*200))),'%'])
    text(0.1,.9,'Fn0 found by linearizing the problem')

else % Fn0 input in command line
        
    text(0.1,.8,'Fn0 input into the problem')

end

end

    text(0.1,1.0,['tolerance = ',num2str(Tolr),'Jacobian is on'])
    text(0.1,1.05,['Maximum # iterations = ',num2str(Maxiter)])

end

% --------------------------------------------------------------
function [resid,jac] = colloc(Fn) % called by fsolve with  Fn

global bet delt thei sini cosi usip2 ri ...
        ndeln nsini ncosi ideln nsinr zet onen

%	input	-	Fn      =   col vec of unknown  Fn (size 2*ipt)
%
%	output	-	resid   =   col vec of residuals of 2*ipt eqs (size 2*ipt)
%           -   jac     =   2ipt X 2ipt Jacob mat \partialresid/\partial Fn
%
%   global variables    see list above
%
%   internal variables
%
%           -   ucos    =   partial column vector in ui^2
%           -   usin    =   partial column vector in ui^2
%           -   usi2    =   col vector = half the norm squared of speed at
%                           collocation pts on interface in ambient fluid
%           -   resber	=   col vec of resid of Bernoulli eqs (size ipt)
%                               resber represents (p-p')/rho' U^2
%           -   respsi	=   col vec of resid of streamfunc eqs (size ipt)
%           -   jacB    =   Jacobian of Bernoulli equations
%           -   jacS    =   Jacobian of streamline eqs
%

ucos    =   (1./delt-ndeln*Fn)./ri - cosi + ncosi*Fn    ; % col vec
usin    =   sini + nsini*Fn                             ; % col vec

usi2    =   (1/2).*(ucos.^2 + usin.^2)	; % col vec of (1/2)U^2

% Calculate residuals of ipt Bernoulli eqs first (col vec of size ipt):

resber  =   usip2 + bet.*thei./delt -...
            (1-zet).*(usi2 + bet.*thei./delt )      ; % col vec of resid
        
% now, calculate residuals of ipt streamfunction eqs (col vec of size ipt)

respsi  =   (nsinr - ideln)*Fn                      ; % col vec of resids

% concatenate resber and respsi into output col vec of size 2*ipt:
resid   =   [resber;respsi]                         ;

% Use statement below to enforce explicitly n delta^n Fn = 0 as last eq:
% resid   =   [resber;respsi(1:end-1);ndeln*Fn]	;

% Now calculate the Jacobian:

jacB    =   -(1-zet).*(  (ucos*onen).*...
                           (ncosi-(sini./thei)*ndeln*delt) + ...
                         (usin*onen).*...
                            nsini	)	; % Bernoulli Jacobian
jacS    =   nsinr - ideln            	; % Streamfunction Jacobian

% concatenate jacB and jacS into output col vec of size 2*ipt:
jac     =   [jacB;jacS]                 	;

% --------------------------------------------------------------
function Yprime = curvider(th,~)

% Calculates the derivative of dim(-) curvilinear distance 
% with respect to polar angle
%
% MYL 3/1/12
%
%
%	input	-	th		=	polar angle
%	output	-	Yprime	=	derivative of dim(-) curv dist WRT polar angle
%
%   global variables    none
%   

% global  none

% 	other internal variables    none
%

if th < 0.1 % to avoid blow up at vanishing polar angle, expand in series:
Yprime	= 1+2*th^2/9+14*th^4/405                ;
else
Yprime	= sqrt(th^2+(1-th/tan(th))^2)/sin(th)   ; % dim(-) derivative
end
            
% --------------------------------------- %

