close all
clear all
clf

xmax = 150;
ymax = 150;

sigma = (xmax/10);
tmax = 100*24;
w = (2*pi)/1000;    % frequency of temperature oscillation in time
A = 0;
Ys = 0;
NY0 = ymax;
Ts = 0;
Ta = zeros(xmax,NY0);
NY = 0;
X = zeros(xmax,NY0);
Y = zeros(xmax,NY0);
Ttop = 0;
Ytop = 2;

count = 1;

for i = 1:xmax
    Ys(i) = 0-exp(-((i-(xmax/2))^2)/(2*sigma)^2);    % height of the ground surface (<=0) as a function of x
end

for i =1:xmax
    NY(i) = (Ytop-Ys(i))*NY0/(Ytop-Ys(xmax/2));      % number of y-points from the ground to the top
    for j = 1:NY(i)
        Y(i,j) = ((Ys(i)-Ytop)/(NY(i)-1))*j+(Ytop-(Ys(i)-Ytop)/(NY(i)-1));     % y gridpoint positions
        X(i,j) = i;                                                            % x gridpoint positions
    end
end

d = 6;   % controls average Temp of ground
dd = 15; % conrols average temp of top
e = 10;  % controls swing of the ground
f = 2;   % controls swing of the top 
valley = 5;

Tmax = -1e10;   % temporary Temp variables for finding minimum temp over t and x (see next for loop)
Tmin = 1e10;

% find the minimum temp over t and x
for i = 1:xmax
    Ts1 = (valley*abs(Ys(i))+d)-((Ys(xmax/2)-Ys(i))+e);
    Tt1 = dd-f;
    Tmin = min(min(Tmin,Ts1),Tt1);
end

% find the maximum over t and x, after subtracting the min to make the new
% min 0
for i = 1:xmax
    Ts1 = (valley*abs(Ys(i))+d)+((Ys(xmax/2)-Ys(i))+e)-Tmin;
    Tt1 = dd+f-Tmin;
    Tmax = max(max(Tmax,Ts1),Tt1);
end

% main time loop
for t = 1:tmax
    r = rand(xmax,ymax);    % random numbers for stochasticity
    Ttop = dd+f*sin(w*t)-Tmin;  % temp at the top at the current time t
    % x loop
    for i = 1:xmax
        Ts(i) = (valley*abs(Ys(i))+d)+((Ys(xmax/2)-Ys(i))+e)*sin(w*t+(pi/2))-Tmin;   % temp of the surface as a function of x and t
        A = (Ts(i) - Ttop)/(Ys(i)-Ytop);  % used for linear interpolation below
        % y loop
        for j = 1:ymax
            Tt = (A*Y(i,j)) + Ttop-A*Ytop;   % linearly interpolate between surface and the top to get temp at the current y
            c = Tt/Tmax;                     % local temperature normalized to be between 0 and 1
            % stochastic code to randomize the fog 
            if c <= 0.1 && r(i,j)>0.1    % very cold: higly probably condensation
                Ta(i,j) = 10;
            end
            if c <= 0.25 && r(i,j)>0.25  % cool: somewhat probable condensation
                Ta(i,j) = 10;
            end
            if c <= 0.5 && r(i,j)>0.4    % warm: rare condensation
                Ta(i,j) = 10;
            end
            if c <= 0.75 && r(i,j)>0.6   % quite warm: very rare condensation
                Ta(i,j) = 10;
            end
            if c <= 0.9 && r(i,j)>0.7    % hot: very, very rare condensation
                Ta(i,j) = 10;
            end
            if c >= 0.25 && r(i,j)>0.8   % cool: rare evaporation
                Ta(i,j) = 0;
            end
            if c >= 0.35 && r(i,j)>0.5   % warm: medium likely evaporation
                Ta(i,j) = 0;
            end
            if c >= 0.65 && r(i,j)>0.3    % hot-ish: likely evaporation
                Ta(i,j) = 0;
            end
            if c >= 0.75 && r(i,j)>0.25   % quite warm: very likely evaporation
                Ta(i,j) = 0;
            end
            if c >= 0.9 && r(i,j)>0.1   % hot: very, very likely evaporation
                Ta(i,j) = 0;
            end
        end
    end
    hour = mod(12+((w*t)/(2*pi))*24,24); % hour modulo 24 hours, starting at 12:00
    Ts(1) = Tmax;
    Ts(xmax) = 0;
    scatter(reshape(X,[xmax*ymax,1]),reshape(Y,[xmax*ymax,1]),reshape(Ta+0.001,[xmax*ymax,1]),'b','filled');
    surface('XData',[1:xmax;1:xmax],'YData',[Ys;Ys],'ZData',zeros(2,xmax),'CData',[Ts;Ts],'FaceColor','none','EdgeColor','flat','Marker','none','LineWidth',5);
    ylabel(colorbar,'Temperature (C)');
    title(['Hour in Day: ',num2str(hour,'%5.2f')]);
    xlabel('Longitudinal Distance')
    ylabel('Altitude')
    hold on
    axis ([0 xmax -1 1])
    hold off
    drawnow;
    pause(0.00000001)
    meanTs = mean(Ts);
    count = count+1;
    if mod(t,10) == 0
        count = count+1;
        saveas(gcf,['Cyclic Daily Fog Aggregation',sprintf('%03d',count),'.png'])
    end
end