%  Solve the heat equation by dumbdumb finite difference on
%  geometry of spectral solution.
%
%	OUTPUT:   rhoI(nx,ny,NT). rhoS(nx,ny,NT) ....
%		  hS(nx,ny,NT) ......
%


	clear rho1
	NT = 600;
        nx = 50;
        ny = 50;
        dt = 0.25;
        D = 0.015;
        dx = 2*pi/nx

        X = 2*pi*[1:nx]/nx;


%	SET UP INITIAL CONDITION FOR THE DENSITY
        [X,Y] = meshgrid(X,X);
	xx = X - mean(mean(X));
	yy = Y - mean(mean(Y));
	C0 = exp(-3*((xx+2).^2 + (yy+2.5).^2));
%	C0 = 0.5*(C0 +  exp(-1*((xx+1).^2 + (yy+1).^2)));

	cdiff = dt/(dx*dx);
        cdiff*D

%	Set Up profiles of rhpo1,rho2,rho3
        rho1(:,:,1) = C0;
        rho1(nx+1,:,1) = rho1(2,:,1);
        rho1(:,ny+1,1) = rho1(:,2,1);
        
%	SET UP DIFFUSIVITY PROFILE
	alpha = peaks(nx+1).*peaks(nx+1);
	alpha = sqrt(abs(peaks(nx+1)));
%	alpha = exp(-0.25*((xx).^2 + (yy).^2));
%        alpha(nx+1,:,1) = alpha(2,:,1);
%        alpha(:,ny+1,1) = alpha(:,2,1);

	scale= max(max(alpha));
        alpha = D*alpha/scale;

%	SET UP WIND
        w(1) = 0.05;
        w(2) = 0.05;

%	EVOLVE THE SOLUTION IN TIME for NT time steps
   for it = 1:NT

        for i = 2:nx
         for j = 2:ny
         apx = 0.5*(alpha(i+1,j) + alpha(i,j));
         amx = 0.5*(alpha(i-1,j) + alpha(i,j));
         apy = 0.5*(alpha(i,j+1) + alpha(i,j));
         amx = 0.5*(alpha(i,j-1) + alpha(i,j));

	% COMPUTE FOR RHO1
	 aupx = apx*(rho1(i+1,j,it)-rho1(i,j,it));
	 aumx = apx*(rho1(i,j,it)-rho1(i-1,j,it));
	 aupy = apy*(rho1(i,j+1,it)-rho1(i,j,it));
	 aumy = apy*(rho1(i,j,it)-rho1(i,j-1,it));
	 diff = aupx-aumx + aupy-aumy;
         rho1(i,j,it+1)=rho1(i,j,it)+cdiff*diff;

	% REPEAT ABOVE COMPUTATION FOR OTHER TWO RHOS!!

	% NOW DO BIOLOGY PART WITH HUMANS
	    dhS = 0;   % -dt*cvh*nbh*hS(i,j,nt)*rhoV(i,j,nt)/rho(i,j,nt);
	    dhI = 0;
	    dhV = 0;
	    dhR = 0;

	  hS(i,j,nt+1) = hS(i,j,nt) + dhS;
	  hI(i,j,nt+1) = hI(i,j,nt) + dhI;
	  hV(i,j,nt+1) = hV(i,j,nt) + dhV;
	  hR(i,j,nt+1) = hR(i,j,nt) + dhR;

 	 end
        end

        rho1(1,:,it+1) = rho1(nx,:,it+1);
        rho1(:,1,it+1) = rho1(:,ny,it+1);
        rho1(nx+1,:,it+1) = rho1(2,:,it+1);
        rho1(:,ny+1,it+1) = rho1(:,2,it+1);


        end