clear all

%==========Parameters of Products==========
m1=0;         %mass [kg]
Tp=273.15+3;  %initial temperature of product [k]

if Tp>273.15

Cp1=4187;     %Specific heat capacity of water[J/kgK]
u1=0;         %Latent heat of coolant[J/kg] 
    
else
    
Cp1=2108;     %Specific heat capacity of ice[J/kgK]
u1=335000;    %Latent heat of coolant[J/kg]   

end
%==========Parameters of Coolants==========
m2=0.80617;%mass of coolant[kg]
Cp2=2108;  %Specific heat capacity of coolant[J/kgK] (Assume only ice)
u2=335000; %Latent heat of coolant[J/kg](Assume only ice)
Cp2b=4187; %Specific heat capacity of coolant after melting[J/kgK] (Water)

%==========Parameters of environment==========
Te=273.15+20;  %External Temperature [K]
Ti=273.15-20;  %Initital Inside Temperature [K]

%==========Parameters of Packaging layers==========
th0 =0; %[m]total thickness of air spaces where both bounding sourfaces are covered with aluminium foil
th1=0;  %[m]total thickness of those air spaces where only one bounding surface is covered with aluminum foil

k1=0.01;  %thermal conductivity 1st material [W/mK] Insulation Layer
th2=0.1; %thickness 1st material [m]0.02745
R1=1/k1;   %thermal resistance 1st material [mK/W]

k2=0.0439;  %thermal conductivity 2nd material [W/mK] Carboard Box
th3=0.00278; %thickness 2nd material [m]
R2=1/k2;   %thermal resistance 2nd material [mK/W]

k3=210;    %thermal conductivity 3rd material [W/mK] AI Foil
th4=0;     %thickness 3rd material [m]
R3=1/k3;   %thermal resistance 3rd material [mK/W]

th5=0;     %thickness of 1st non-foiled air space [m]
th6=0;     %thickness of 2nd non-foiled air space [m]

%==========Parameters of Packaging==========
%Inside dimensions
l1=0.256-2*0.00278-2*0.02745; %length [m]
w1=0.205-2*0.00278-2*0.02745; %width  [m]
h1=0.212-2*0.00278-2*0.02745; %height [m]
Ain=2*(h1*l1)+2*(h1*w1)+2*(l1*w1);
Ein=0.83;   %Emissivity of inside surface 0.83

%Outside dimensions
l2=l1+th2+th3;   %length [m]
w2=w1+th2+th3;   %width  [m]
h2=h1+th2+th3;   %height [m]
Aout=2*(h2*l2)+2*(h2*w2)+2*(l2*w2);
Eout=0.83;  %Emissivity of outside surface0.83

%==========Product dimensions Coolants in this case==========
lp=0.163;   %length [m]
wp=0.091;   %width  [m]
hp=0.032*2; %height [m]
Aprod=2*(hp*lp)+2*(hp*wp)+2*(lp*wp);
Eprod=0.97; %Emissivity of product gel packs cover with polypropylene

%==========Surface Area calculation according to S.J, Choi' model==========
A2=wp*lp;
A1=Aout-A2;
A3=Ain-A2;
A4=Aprod-A2;

%==========Packaging heat capacity(not sure should involve)==========
%rhop= 10.62; %density of insulating materials [kg/m^3]
%v3= h2*l2*w2-h1*l1*w1; %volume of the insulating material [m^3]
%m3=v3*rhop; %mass of insulating materials [kg]
%Cp3= 1700; %heat capacity  of the insulating materials [J/kgK]



%==========Calculations according to S.J, Choi' model==========
%Thermal resistance of wall [m^2.K/W]
R=0.039*th0+0.037*th1+R1*th2+R2*th3+R3*th4; %+((0.217*th5)/(5.918+th5))+((0.217*th6)/(5.918+th6)); 

%Surface heat transfer coefficient [W/m^2.K]
hout=1.778+5.198*Eout;
h3=(3.557*(A4/(A3+A4))+(4.6/((A3/A4)*(1/Ein)+(1/Eprod)-1)));

%Heat Penetration Rate [W/K]
HPR=1/((1/(hout*Aout))+(1/(1/((1/(h3*A3))+(R/(((A1*A3)^0.5))))+(A2/(R))))); 

%Total thermal resistance within the system [K/W]
RQ=1/HPR; 

%==========amount of energy for inside temp go from Ti to 0==========
 a=m1*Cp1*(273.15-Tp)+m2*Cp2*(273.15-Ti);%+m3*Cp3*(Ti-273.15);
%==========amount of energy with latent heat included========== 
 b=(m1*Cp1*(273.15-Tp)+m2*Cp2*(273.15-Ti)+m2*u2+m1*u1);%+m3*Cp3*(Ti-Te)

%==========Mathmatical Model Running==========
 
for t=1:1:700000
   dT(1)=Te-Ti; %Temperature difference at start (initial temp difference)
    i(t)=dT(t)/RQ; %Heat flow across a packaging at time t [W]
    Q(t)=sum(i); %Total heat flow across a packaging [W]
    
 if m2>0 % when there is coolant inside the box
        
     if   dT(t)>(Te-273.15) && Q(t)<a %stage 1 when temp difference less than 0,
          
          T(t)=(Ti-Te)*exp(-t/(RQ*(m1*Cp1+m2*Cp2)))+Te;
          
     elseif Q(t)>b %stage 3 when phase change completely finish
        
          T(t)=(273.15-Te)*exp(-(t-y)/(RQ*(m1*Cp1+m2*Cp2b)))+Te;
      
     else
          T(t)=273.15; %stage 2 when only latent heat, no temp different
    
               if Q(t)~=b %when Q(t) closest to b, specify a time y
        
                  y=[t];  
      
               end    
     end
 else
 T(t)=(Ti-Te)*exp(-(t)/(RQ*(m1*Cp1)))+Te; %when without coolant
   
 end
 dT(t+1)=Te-T(t);
   
               if T(t)>273.15+5
                  return;
               end 
 
end

plot(T-273.15)