amp = 100;
xi_0 = 400;
wavelength = 300;
period = 10;
alpha = 1000;
beta = 100
k=2*pi/wavelength;
omega= 2*pi/period;
T_N = 20;
T_S = 25;
T_mean = (T_N+T_S)*0.5;
T_fac = (T_N-T_S)/pi;
sigma=0.05;

x=[0:5:1000];
y=[0:5:800];
T=zeros(length(y),length(x));
Tran=T;
figure(1);
dt=0.1;
t=0;
while (1)
    
    for i=1:length(y)
        for j=1:length(x)
             xi = xi_0 + amp*exp(-x(j)/alpha)*sin(x(j)*k-omega*t);
            T(i,j) =  T_mean + T_fac*atan((y(i)-xi)/beta);
            Tran(i,j) = gauss_real(T(i,j),sigma);
          
        end
    end
    figure(1);
    pcolor(x,y,T);
    shading('flat');
    title(['Time = ' num2str(t) ' sec']);
    drawnow;
    figure(2);
    pcolor(x,y,Tran);
    shading('flat');
    title(['Time = ' num2str(t) ' sec']);
    drawnow;
    pause(dt);
    t=t+dt;
end
