%====================================================================================================
%	Forward Model==================================================================================
%====================================================================================================
%========================================================
	function u=heat_forward(uo,K,N,uvel,vvel)
%========================================================	
% Heat equation
% du/dt = a ( d^2u/dx^2 + d^2u/dy^2 )

% 2 D
Q=0;
u=uo;
for t=1:N
      u = heat_forward_step(u,t,K,uvel,vvel);
	u = heat_forward_bc(u,t,Q);
end	
return

%========================================================
	function u = heat_forward_step(u,t,K,uvel,vvel);
%========================================================	
       
      [I,J,T]=size(u);	
	dt=0.01;
	dx=0.05	;
	rhs_diff = zeros(I,J);
	rhs_adv = zeros(I,J);
		
	%for x=2:I-1
	%for y=2:J-1
	x=2:I-1;
	y=2:J-1;
	rhs_diff(x,y) = + K*( u(x+1,y,t) - 2*u(x,y,t) + u(x-1,y,t) )/dx^2 ...
	           + K*( u(x,y+1,t) - 2*u(x,y,t) + u(x,y-1,t) )/dx^2  ;
	rhs_adv(x,y) =  - uvel(x,y).*(u(x+1,y,t) - u(x-1,y,t) )/(2*dx) ...
	           - vvel(x,y).*(u(x,y+1,t) - u(x,y-1,t) )/(2*dx) ;
		     
	u(x,y,t+1)  =  u(x,y,t) + (rhs_diff(x,y)+rhs_adv(x,y))*dt    ;
	
	%end
	%end
return	
	
%========================================================
	function 	u = heat_forward_bc(u,t,Q);
%========================================================	
	% Boundary
      [I,J,T]=size(u);	
	%for x=1:I
	x=1:I;
%		u(x,J,t+1)=Q;
%		u(x,1,t+1)=Q;
		u(x,J,t+1)=u(x,J-1,t+1);
		u(x,1,t+1)=u(x,1+1,t+1);
	%end
	%for y=1:J
	y=1:J;
%		u(I,y,t+1)=Q;
%		u(1,y,t+1)=Q;
		u(I,y,t+1)=u(I-1,y,t+1);
		u(1,y,t+1)=u(1+1,y,t+1);
	%end

return