r_0 = 200;
v = 5;
phi = pi/50;
a_1 = 20;
a_2 = 30 ;
a_3 = 40 ;
T_N = 20;
T_S = 25;
T_mean = (T_N+T_S)*0.5;
T_fac = (T_N-T_S)/pi;
beta = 50;
sigma=0.05;

x=[0:5:1000];
y=[0:5:800];
T=zeros(length(y),length(x));
Tran=T;
figure(1);
dt=1;
t=0;
while (1)
    
    for i=1:length(y)
        for j=1:length(x)
            theta=atan2(y(i)-400,x(j)-500);
            r= sqrt((x(j)-500)^2+(y(i)-400)^2);
            xi = r_0+v*t + a_1 * cos(theta) + a_2*sin(2*theta) + a_3*sin(3*theta+phi*t);
            T(i,j) =  T_mean + T_fac*atan((r-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
