close all
clear all

PD=1.1;
n=32; %number averaged carbon no.
nw=n*PD;
Ts=857; %surface temp [k]
Ta=298;
Mw=450.9; %molar mass [g/mol]
[Tm,mu_l]=n_properties(n,nw,Ts,Mw);
rho_l=674; %liquid density [kg/m3]
rho_s=920; %solid fuel density
rho_sref=970;
lambda_l=0.135; %thermal conductivity [W/(m.K)]
cl=2725.5; %specific heat liquid [J/(kg.K)]
cs=2030;
Lm=170; %heat of melting [kJ/kg]
hm=Lm*10^3+cs*(Tm-Ta); %[J/kg]

at=lambda_l/(cl*rho_s)*log(1+cl*(Ts-Tm)/hm);

K=1.6*10^4;
B=4.7;
CB1=2/(2+1.25*B^0.75);
L=1.06; %grain length [m]

mo=4.32; %oxidizer mass flow rate [kg/s]
mf(1,1)=0; %fuel flow rate [kg/s]
dt=0.05;
tb=48; %burn time [s]
t=(0:dt:tb); %time [s]
o=mo*tb; %oxidizer mass [kg]
n=501; %number of grain length finite elements
z=linspace(0,L,n); %vector of grain length elements
dz=L/500;
R(1,1)=0.055; %initial port radius [m]
G(1,1)=mo/(pi*R(1,1)^2); 
Go(1,1)=G(1,1);
for j=2:length(z)
    R(1,j)=R(1,1);
    G(1,j)=(mo+mf(1,j-1))/(pi*R(1,j)^2); %total mass flux [kg/(m2.s2)]
    Go(1,j)=mo/(pi*R(1,j)^2);
    Rent=K*at^2*rho_s^2*rho_l/(mu_l)*CB1^(-2)*B^(-3)*1.1^(-3)*G(1,j)^(0.4)*z(j)^(0.4);
    fi=1+0.61*Rent^0.4;
    rcl=0.03*(6.5*10^-5)^(0.2)/(rho_s)*1.1*B*CB1*G(1,j)^0.8*z(j)^(-0.2);
    r(1,j)=fi*rcl*10^3;
    mf(1,j)=r(1,j)/1000*rho_s*dz*2*pi*R(1,j)+mf(1,j-1);
    r(1,1)=r(1,2);
end

for i=2:length(t)
    mf(i,1)=0;
    R(i,1)=R(i-1,1)+dt*r(i-1,1)/1000;
    G(i,1)=mo/(pi*R(i,1)^2);
    Go(i,1)=G(i,1);
    for j=2:length(z)
        R(i,j)=R(i-1,j)+dt*r(i-1,j)/1000;
        if R(i,j) < 0.16
            G(i,j)=(mo+mf(i,j-1))/(pi*R(i,j)^2); %mass flux [kg/(m2.s2)]
            Go(i,j)=mo/(pi*R(i,j)^2);
            Rent=K*at^2*rho_s^2*rho_l/(mu_l)*CB1^(-2)*B^(-3)*1.1^(-3)*G(i,j)^(0.4)*z(j)^(0.4);
            fi=1+0.61*Rent^0.4;
            rcl=0.03*(6.5*10^-5)^(0.2)/(rho_s)*1.1*B*CB1*G(i,j)^0.8*z(j)^(-0.2);
            r(i,j)=fi*rcl*10^3;
        else
            G(i,j)=(mo+mf(i,j-1))/(pi*R(i,j)^2); %mass flux [kg/(m2.s2)]
            Go(i,j)=mo/(pi*R(i,j)^2);
            r(i,j)=0;
        end
        mf(i,j)=r(i,j)/1000*rho_s*dz*2*pi*R(i,j)+mf(i,j-1);
        r(i,1)=r(i,2);
    end
end
OF=mo./mf(:,end);
Gav=mean(G);
rav=mean(mean(r))
mfav=mean(mf(:,end));
f=sum(dt*mf(:,end));
OFav=o/f
Goav=mean(mean(Go))

Isp=[2276, 2408, 2525, 2552]; %specific impulse [Ns/kg]
mp=mo+mf(:,end);
F=mp*Isp;
Fav=mean(F)
mpav=mean(mp)

[X,Y] = meshgrid(t,z);
Z=R'*100;
figure(1)
surf(X,Y,Z,'LineStyle','none');
zlim([5 16])
xlabel("t [s]",'Interpreter','latex')
ylabel("z [m]",'Interpreter','latex')
zlabel("R [cm]",'Interpreter','latex')
hold on
spacing1 = 40; 
for i = 1:spacing1:length(X(1,:))
    plot3(X(:,i),Y(:,i),Z(:,i),'-k');
end
spacing2=25;
for i = 1:spacing2:length(Y(:,1))
    plot3(X(i,:),Y(i,:),Z(i,:),'-k');
end
hold off

Z=r';
figure(2)
surf(X,Y,Z,'LineStyle','none');
% zlim([5 16])
xlabel("t [s]",'Interpreter','latex')
ylabel("z [m]",'Interpreter','latex')
zlabel("$\dot{r}$ [mm/s]",'Interpreter','latex')
hold on
spacing1 = 40; 
for i = 1:spacing1:length(X(1,:))
    plot3(X(:,i),Y(:,i),Z(:,i),'-k');
end
spacing2=25;
for i = 1:spacing2:length(Y(:,1))
    plot3(X(i,:),Y(i,:),Z(i,:),'-k');
end
hold off

figure(3)
plot(t,mean(r'))
grid
xlabel("t [s]",'Interpreter','latex')
ylabel("$\dot{r}$ [mm/s]",'Interpreter','latex')

figure(4)
plot(t,mean(R')*100)
grid
xlabel("t [s]",'Interpreter','latex')
ylabel("R [cm]",'Interpreter','latex')

figure(5)
plot(t,mf(:,end))
grid
xlabel("t [s]",'Interpreter','latex')
ylabel("$\dot{m}_f$ [kg/s]",'Interpreter','latex')

figure(6)
plot(t,F/1000)
grid
xlabel("t [s]",'Interpreter','latex')
ylabel("F [kN]",'Interpreter','latex')