clear;

load inputs.txt;

cav=inputs(1);
ceq=inputs(2);
L=inputs(3);
dt=inputs(4);
D=inputs(5);
npas=inputs(6);

HV0=inputs(7);
Htest=1;
if HV0==0
    Htest=0;
else
    A=inputs(8);
    if A==0
    Htest=0;
    end
end

Va=0.01166;
Vg=0.04666;
Cav=cav/Va;
Ceq=ceq/Va;
m=50;

C=Cav*ones(1,m);
C(1)=Ceq;
CC(1)=Ceq;
CC(m)=Cav;
dN=zeros(1,m);
temps=0;
N=0;
t=0;
Ntot=0;
dx=L/m;
x=0:dx:(m-1)*dx;
S=100;

for j=1:npas
   for i=1:m-1
   	  dN(i)=-D*S*((C(i)-C(i+1))/dx)*dt;
   end
   
   deltaN=dN(1);
   
   for i=2:m-1
      CC(i)=C(i)+(-dN(i-1)+dN(i))/(S*dx);
   end
   C=CC;
   C(m-1:m)=C(m-1);
   N(j+1)=N(j)+deltaN;
   t=t+dt;
   temps(j+1)=t;
   
    if rem(j,1000)==0
       Computing_step=j
    end
   
end

Vf=(N*Vg)/(S*L);

if Htest==1
    M=HV0+A*sqrt(Vf);
    subplot(1,2,1),semilogx(temps/3600,M);
    xlabel('Ageing time (h)');
    ylabel('Mechanical property (MPa)');
    subplot(1,2,2);
end
    semilogx(temps/3600,Vf);
    xlabel('Ageing time (h)');
    ylabel('Gamma-prime volume fraction');


clear HV0 Htest A C CC Cav Ceq Computing_step D L N Ntot S Va Vg cav ceq dN deltaN dt dx i j m npas t temps x inputs

