ploplot = 0 ; 
%%% toggles various plots set 1 
firstorder  = 1 ; 
%%% to skip  1rst order set 0 
secondorder  = 1 ; 
%% to try second order scheme for the fearless SOFF 
%% if 1 first order must have run before 

if ( 1 == firstorder ) 
%% first order sheme Rouy Tourin discretisation 
%% and Sweeping ordering 

% computational domain 0,1  by 0,1 
xmin = 0; xmax = 1. ;

% Discretization  and data 
N = 41 ; 
Big = 300;
x = linspace(xmin, xmax, N);
y = x;
h = x(2)-x(1);
[X,Y] = meshgrid(x,y) ; 
n = refract(X,Y) ;   %% refraction index only 1 homogeneous for now 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% init sub solution 
phi = Big *ones(N,N) ; 
%% BC outgoing if not something else

% BC Dirichlet 
%(1:N, 1) =  0  ; 
%phi(1:N, N) = 0 ; 
%phi(1, 1:N) = 0 ; 
%phi(N, 1:N)= 0 ; 
%%% 

[phi,exact] = source(X,Y,phi,0.25,0.223,0.7,0.75,h) ; %% enforce 2 source points and compute exact solution 

%% iteration parameters 
conv = Big ;
tol = 1e-14 ; 
niter = 0;

phiold  = phi ; %% iteration variables 


'first order' 

while ( conv > tol)  %% iterations 
niter = niter + 1
%%  Sweep in four directions 

phi = sweep(phiold,n,h,N,[2:N-1],[2:N-1]) ; 
phi = sweep(phi,n,h,N,[2:N-1],[N-1:-1:2]) ; 
phi = sweep(phi,n,h,N,[N-1:-1:2],[2:N-1]) ; 
phi = sweep(phi,n,h,N,[N-1:-1:2],[N-1:-1:2]) ; 

%% test convergence for the iterative prociss
conv  = max(max(abs(phiold-phi)))/max(max(abs(phiold))) 
%% L1 errr
l1err = sum(sum(abs(phi(2:N-1,2:N-1)-exact(2:N-1,2:N-1))))/sum(sum(abs(exact(2:N-1,2:N-1))))
%% iterate
phiold  = phi ;

%% some plots 
if (ploplot == 1) 
figure(1)
contour(x(2:N-1), y(2:N-1), phi(2:N-1,2:N-1), 20,'b-') 
pause 
contour(x(2:N-1), y(2:N-1), exact(2:N-1,2:N-1), 20,'r-') 
pause 
surf(x(2:N-1), y(2:N-1), abs(phi(2:N-1,2:N-1)-exact(2:N-1,2:N-1))) 
colorbar
axis([0 1 0 1 0 0.01])
pause 
end %if 

end %% while 

figure(1)
contour(x(2:N-1), y(2:N-1), phi(2:N-1,2:N-1), 20,'b-') 
pause 
contour(x(2:N-1), y(2:N-1), exact(2:N-1,2:N-1), 20,'r-') 
pause 
surf(x(2:N-1), y(2:N-1), abs(phi(2:N-1,2:N-1)-exact(2:N-1,2:N-1))) 
colorbar
pause 


end %% if first order 

if (secondorder == 1) 


'second order'

%% SOFF : Second Order For the Fearless
%% second order correction based on the compact stencil 
%% with curvature correction 


phi2  = phi; %% variable for  second order scheme
phiold = phi ; %% init with first order 


conv  = Big ; 
tol = 1e-14 ; 
niter = 0;

while ( conv > tol)  %% iterations 
niter = niter + 1

%% sweep3b chooses stencil based on the determinant 
%% sweep3 chooses according 1rst order upwind direction 

[phi2,DD] = sweep3b(phiold,n,h,N,phi,[2:N-1],[2:N-1]) ; 
[phi2,DD] = sweep3b(phi2,n,h,N,phi,[2:N-1],[N-1:-1:2]) ; 
[phi2,DD] = sweep3b(phi2,n,h,N,phi,[N-1:-1:2],[2:N-1]) ; 
[phi2,DD] = sweep3b(phi2,n,h,N,phi,[N-1:-1:2],[N-1:-1:2]) ; 


conv  = max(max(abs(phiold-phi2)))/max(max(abs(phiold))) 
%% just check error here convergence does not make sense  max does not make sense
l1err = sum(sum(abs(phi2(2:N-1,2:N-1)-exact(2:N-1,2:N-1))))/sum(sum(abs(exact(2:N-1,2:N-1))))

phiold = phi2 ; 

if (ploplot == 1 ) 
figure(2) 
surf(x(2:N-1), y(2:N-1), DD(2:N-1,2:N-1))  
axis([0 1 0 1 -0.1 0.1])
pause 
contour(x(2:N-1), y(2:N-1), phi2(2:N-1,2:N-1),20) 
colorbar
pause 
end  %if 

end %% while 

figure(2) 
contour(x(2:N-1), y(2:N-1), phi2(2:N-1,2:N-1), 20,'b-') 
pause 
contour(x(2:N-1), y(2:N-1), exact(2:N-1,2:N-1), 20,'r-') 
pause 
surf(x(2:N-1), y(2:N-1), abs(phi2(2:N-1,2:N-1)-exact(2:N-1,2:N-1))) 
colorbar
pause 

%% difference between 1rst and second 
figure(3) 
surf(x(2:N-1), y(2:N-1), abs(phi2(2:N-1,2:N-1)-phi(2:N-1,2:N-1))) 
colorbar
axis([0 1 0 1 0 0.01])

end  %% if second order