clear;
clf;


n=550;
t=0;
s=1e5;
phi=101.325; %kPa
Em=18.02/28.966;
EmH=22/28.966;

S=(5e3)^2;%m^2
dt=0.01;
V=zeros(n,1);
V(1)=s^3; % m^3
VH=zeros(n,1);
VH(1)=s^3; % m^3
Rn=350; %W/m2
c=4.19; %kJ/kgC
cp= 1.006; %kJ/kgC
Ts=zeros(n,1);
Ts(1)=20; %C
Ta=zeros(n,1);
Ta(1)=22;%C
es=zeros(n,1);
es(1)=-273.15+exp(77.345+(0.0057*(Ts(1)+273.15))-(7235./(Ts(1)+273.15)))./((Ts(1)+273.15).^8);
ea=zeros(n,1);
ea(1)=-273.15+exp(77.345+(0.0057*(Ta(1)+273.15))-(7235./(Ta(1)+273.15)))./((Ta(1)+273.15).^8);
lhv=2261;%kJ/kg
N=zeros(n,1);
N(1)=(n^3)*c.*Ts(1); %kJ
B=zeros(1,n);
B(1)= (cp*phi*(Ts(1)-Ta(1)))./(Em*lhv*(es(1)-ea(1)));
BH=zeros(1,n);
BH(1)= (cp*phi*(Ts(1)-Ta(1)))./(EmH*lhv*(es(1)-ea(1)));
Tb=zeros(n,1);
Tb(1)=20;
k=zeros(n,1);
k(1)=(Rn+N(1))/(lhv*(1+B(1))+c*(Ts(1)-Tb(1)));
kH=zeros(n,1);
kH(1)=(Rn+N(1))/(lhv*(1+BH(1))+c*(Ts(1)-Tb(1)));
theta=(0);
A=15;
g=0;
for i=1:n
    
    Ta(i+1)=1.1*A*sin(theta)+A;
    Ts(i+1)=0.8*A*sin(theta-0.2)+A;
    
    
    Tb(i+1)=0;
    es(i+1)=-273.15+exp(77.345+(0.0057*(Ts(i+1)+273.15))-(7235./(Ts(i+1)+273.15)))./((Ts(i+1)+273.15).^8);
    ea(i+1)=-273.15+exp(77.345+(0.0057*(Ta(i+1)+273.15))-(7235./(Ta(i+1)+273.15)))./((Ta(i+1)+273.15).^8);
    N(i+1)=(n^3)*c.*Ts(i+1);
    B(i+1)= (cp*phi*(Ts(i+1)-Ta(i+1)))./(Em*lhv*(es(i)-ea(i)));
    BH(i+1)= (cp*phi*(Ts(i+1)-Ta(i+1)))./(EmH*lhv*(es(i)-ea(i)));
    k(i+1)=(Rn+N(i))/(lhv*(1+B(i))+c*(Ts(i+1)-Tb(i+1)));
    kH(i+1)=(Rn+N(i))/(lhv*(1+B(i))+c*(Ts(i+1)-Tb(i+1)));
    V(i+1)=(-k(i)*S*dt)+V(i);
    if V(i+1)<1e13
        V(i+1)=V(i)+12e13;
    end
    VH(i+1)=(-kH(i)*S*dt)+VH(i);
    if VH(i+1)<1e13
        VH(i+1)=VH(i+1)+12e13;
    end
    t=t+1.25*dt;
    
    plot(t,0.6e-13*V(i),'bo-','LineWidth',3)
    plot(t,Ts(i),'ro-','LineWidth',4)
    plot(t,Ta(i),'mo-','LineWidth',4)
    plot(t,0.6e-13*VH(i),'co-','LineWidth',2)
    title('Cumulative Evaporation From Standing Water','FontSize',16)
    xlabel('Time(days)','FontSize',14)
    ylabel('T(C) and Water Volume (6*10^2 m^3)','FontSize',14)
    drawnow
    hold all
    axis ([0 7 -5 60])
    %    pause(.1)
    theta=theta+((1/40)*pi);
    if theta>2*pi
        theta=0;
    end
    g=g+1;
    saveas(gcf,['final_project2_',sprintf('%03d',g),'.png'])
end