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

%  This code is written by Chiu-Yen Kao

%  Copyright is reserved  !!

global xmin xmax ymin ymax Mx My
global gx gy nx ny nxs nys nxe nye dx dy x y trap_mtr meshx meshy pmeshx pmeshy
global s1 s2 nn sigma1 sigma2 r1 r2 alpha1 alpha2 k1 k2 mu rho1 rho2 beta1 beta2 g
global runTime
global c1 c2



% domain information

nn = 6;
sigma1 = 5;
sigma2 = sigma1;
r1 = 1;
r2 = 0.5;
%r2 = 1;
alpha1 = 5;
alpha2 = alpha1;
k1 = 2;
k2 = k1;
%mu = 5;
mu = 5;

rho1 = 1;
rho2 = rho1;
beta1 = 0.05;
beta2 = beta1;
g = 2;

A1 = (alpha1+sigma1+beta1)/mu;
A2 = (alpha2+sigma2+beta2)/mu;




c1 = 0;
c2 = 0;



xmin =  0.0;
xmax =  A1;
ymin =  0.0;
ymax =  A2;


    

Mx = 64; % number of interior stencils in x-aixs +1 
My = 64; % number of interior stencils in y-aixs +1
gx = 1;  % number of ghostpoint in x-axis
gy = 1;  % number of ghostpoint in y-axis
nx = Mx + 2*gx + 1; % number of stencils in x-axis (including ghost)
ny = My + 2*gy + 1; % number of stencils in y-axis (including ghost)

nxs = gx + 1;
nxe = nx - gx;
nys = gy + 1;
nye = ny - gy;

dx = (xmax - xmin) / Mx;
dy = (ymax - ymin) / My;

if dx~=dy
    error('dx and dy is not equal')
end

xs = xmin - gx*dx; % x coordinate of the first ghost point
ys = ymin - gy*dy; % y coordinate of the first ghost point

x = xs:dx:xmax+gx*dx;
y = ys:dy:ymax+gy*dy;

[meshx,meshy] = meshgrid(x,y);


trap_mtr = zeros(nx,ny);
trap_mtr(nxs,nys) = 1;
trap_mtr(nxe,nys) = 1;
trap_mtr(nxs+2:2:nxe-2,nys) = 2;
trap_mtr(nxs+1:2:nxe-1,nys) = 4;
trap_mtr(:,nye) =  trap_mtr(:,nys);
for i  = nys+2:2:nye-2
    trap_mtr(:,i) =  2*trap_mtr(:,nys);
end
for i = nys+1:2:nye-1
    trap_mtr(:,i) =  4*trap_mtr(:,nys);
end
trap_mtr = trap_mtr/9*dx*dy;

pmeshx = meshx'; 
pmeshy = meshy';
gaussian_sigma = 0.01;

xc = 1;
yc = 1;

%xc = 0.019504;
%yc = 1.012956;

%xc = 0.074608;
%yc = 0.074608;

%xc = 1.012956;
%yc = 0.019504;

pp = zeros(size(meshx));


%for i = nxs : nxe
%    for j = nys : nye
%        x1 = x(i);
%        x2 = y(j);
%        pp(i,j)= exp((-(x1-xc)^2-(x2-yc)^2)/2/gaussian_sigma);       
%    end   
%end
% normalize pp
 %   int_pp = sum(sum(pp.*trap_mtr));
%    pp = pp/int_pp;


pp = ones(nx,ny)/A1/A2;
pp = enforceBC(pp);

u = zeros(nx,ny);v = zeros(nx,ny);



runTime = 0;
tFinal = 1;
titerFinal = tFinal*100;

graph_nbr = 25;
t_graph = floor(titerFinal/graph_nbr);

graph_iter = 0;

zero_mtrx = zeros(nx,ny);

p_all = zeros(nx,ny,graph_nbr);
u_all = zeros(nx,ny,graph_nbr);
v_all = zeros(nx,ny,graph_nbr);

for iter = 1:titerFinal
    
    
    tstep = 0.01;
    T_st =  tstep*(iter-1);
    T_end = tstep*iter
   
   
    [u,v] = pp_to_uv(pp);
    
    
    %pause
    pp = main_advection(pp,zero_mtrx,u,v,T_st,T_end) ;
    pp = max(pp,zero_mtrx);
    
    % focus on smaller window
    
    [i,j] = find(pp>1.e-4);
    
    
     
  
    
    
    % normalize pp
    int_pp = sum(sum(pp.*trap_mtr));
    pp = pp/int_pp;
    
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    [temp,max_ind] = max(pp);
    [temp,max_ind2] = max(temp);
    max_ind1 = max_ind(max_ind2);
    x1cri = pmeshx(max_ind1,max_ind2);
    x2cri = pmeshx(max_ind1,max_ind2);
    test1(iter) = - mu*x1cri+(alpha1*x1cri.^nn./(k1^nn+x1cri.^nn)+sigma1*x1cri/(rho1+x1cri))./(1+x2cri/r2)+beta1;
    test2(iter) = - mu*x2cri+(alpha2*x2cri.^nn./(k2^nn+x2cri.^nn)+sigma2*x2cri/(rho2+x2cri))./(1+x1cri/r1)+beta2;    
    
    xxmax(iter) = xmax;
    yymax(iter) = ymax;
    SS1(iter) = s1;
    SS2(iter) = s2;
    
    if mod(iter,10) == 0
        disp('Percent of total run time = ')
        disp(100*iter/titerFinal)
    end
    
    %if mod(iter,t_graph) == 0
    %if mod(iter,5) == 0
        graph_iter = graph_iter+1;
        u_all(:,:,graph_iter) = u;
        v_all(:,:,graph_iter) = v;
        p_all(:,:,graph_iter) = pp;
        
        figure(1)
        subplot(1,2,1);mesh(pmeshx(nxs:nxe,nys:nye),pmeshy(nxs:nxe,nys:nye),...
            pp(nxs:nxe,nys:nye)); colorbar;title('\phi'); axis tight; axis square; view(-25,75)
        hold on;plot3(pmeshx(max_ind1,max_ind2),pmeshy(max_ind1,max_ind2),pp(max_ind1,max_ind2),'ko')
        interv = round((nxe-nxs)/10); 
        drawnow
        subplot(1,2,2);quiver(pmeshx(nxs:interv:nxe,nys:interv:nye),...
            pmeshy(nxs:interv:nxe,nys:interv:nye),...
            u(nxs:interv:nxe,nys:interv:nye),...
            v(nxs:interv:nxe,nys:interv:nye));axis equal;axis tight;title('velocity');view(0,90)
        figure(50);clf;[c,h] = contour(pmeshx,pmeshy,u,[0 0],'b'); clabel(c);
        hold on;[c,h] = contour(pmeshx,pmeshy,v,[0 0],'r');clabel(c);
        hold on;quiver(pmeshx(nxs:interv:nxe,nys:interv:nye),...
            pmeshy(nxs:interv:nxe,nys:interv:nye),...
            u(nxs:interv:nxe,nys:interv:nye),...
            v(nxs:interv:nxe,nys:interv:nye),'g')
        drawnow
        pause(1)
        
       
        %end
   mass(iter) = sum(sum(pp.*trap_mtr)); 
   mass1(iter) = sum(sum(triu(pp).*trap_mtr)); 
   mass2(iter) = sum(sum(tril(pp).*trap_mtr)); 
   
end
figure(2);plot(tstep:tstep:T_end,SS1,tstep:tstep:T_end,SS2);
legend('S1','S2')
figure(3);plot(tstep:tstep:T_end,mass,tstep:tstep:T_end,mass1,tstep:tstep:T_end,mass2);
legend('mass','mass_1','mass_2');
axis([0 T_end -0.2 1.2])