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/4
    Ta(i+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)));
    k(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);
    x=1:100;
    plot(x,0.6e-13*V(i),'bo-','LineWidth',3)

    drawnow
    hold all
    axis ([1 n/6 50 100])
   

    theta=theta+((1/40)*pi);
    if theta>2*pi
        theta=0;
    end
    
end
hold all
g=0;
for i=1:n/12
    Ta(i+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)));
    k(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);
    x=1:100;
    plot(x,0.6e-13*V(i),'wo-','LineWidth',3)
     
 
 
  r=0.25*i;
    xx=0:0.001:2*r;
   
    a=r;
    b=85;
    
    y=sqrt(r.^2 -(xx-a).^2)+b;
    yy=-sqrt(r.^2 -(xx-a).^2)+b;
    
    plot(xx,y,'yo','LineWidth',5)
    hold on
    plot(xx,yy,'yo','LineWidth',5)
    axis square
    pause(0.1)
    

    

    drawnow
    hold all
    axis ([1 n/6 50 100])

    %        pause(.1)
    theta=theta+((1/40)*pi);
    if theta>2*pi
        theta=0;
    end
    saveas(gcf,['final_project_cartoon',sprintf('%03d',g),'.png'])
    g=g+1;
end
