              
%clear all

recov = 0;
debug = 0;
dbf = 15;
itmax = 200;
tol = 1e-3;
icri = 0;

dg = 2e-3;
pInit = 101325;
Tu = 298;
gamma = 1.37;
Cbd = 200;
Ebd = (4.3 + pInit/1e5/Tu*(136 + 324*dg*1000))^2/Cbd^2/dg/1000;
Ti = Tu*(1 + (60000/Tu - 1)/gamma);
di = 2*sqrt((gamma - 1)/gamma*Ebd/pInit/dg/(1 - Tu/Ti)/pi);
ri = 0.5*di;
Tadb = 1892;
Qc = 1934670;
Le = 1.7;
Su0 = 0.203;

dr = 30e-6;
rmin = 0;%dr;
rmax = 0.24e-2;
dt = 2e-11;%0.2e-9;
nC = round((rmax - rmin)/dr);
nF = nC + 1;
t_end = 100e-6;
				%rk = zeros(ceil(t_end/dt)+1,1);
rk(1) = ri;
rp = rk;
rp2 = rp;
rpCor = rp2;
Energy = 54.8;
switch Energy
  case 2
    EgPol = [-1.9231e+11   3.4208e+07  -4.0267e+04   1.2831e+01];
  case 11.5
    EgPol = [-1.7604e+10   2.6177e+07  -4.0369e+04   2.7137e+01];
  case 54.8
    EgPol = [-2.9110e+08   1.7790e+06  -9.9874e+03   2.9136e+01];
  case 12.4
    EgPol = [-9.6989e+09   1.8887e+07  -3.9422e+04   2.8194e+01];
endswitch

%Egdot = max(0,polyval(EgPol,t));
%%electEnergy
eta_e = 0.5;

rf = (rmin:dr:rmax)';
r = 0.5*(rf(1:end-1)+rf(2:end));
%rf = rf(2:end);

global dT = 1;
if (exist("hfTFile.m"))
  load "hfTFile.m"
else
  dP = [10000*ones(50,1);50000*ones(10,1);100000*ones(40,1);200000*ones(100,1)];
  T=(273.16:dT:60000)';
  global p_i = zeros(1,rows(dP));
  global hfTP = zeros(rows(T),rows(dP));
  global betaP = zeros(rows(T),rows(dP));
  for i = 1:rows(dP)
    p_i(i) = 40000 + sum(dP(1:i));
    p = T.*0 + p_i(i);
    [cp,lambda,M,rho,eta] = updtProp(T,p);    
    hfT = [0;cumsum(cp(1:end-1))]*dT; hfT = hfT(1:1/dT:end);hfTP(:,i) = hfT;
    beta = (log(rho)(2:end) - log(rho)(1:end-1))./(log(T)(2:end) - log(T)(1:end-1));
    beta = [beta(1);beta];
    betaP(:,i) = beta;
    if 0
      hfT = [0;cumsum(cp(1:end-1))]*dT;
      drhohfT = (rho(2:end).*hfT(2:end)-rho(1:end-1).*hfT(1:end-1))/dT;
      drhohfT =[drhohfT;drhohfT(end)];
      drhofT = (rho(2:end) - rho(1:end-1))/dT;
      drhofT =[drhohfT;drhohfT(end)];
      hfT = hfT(1:1/dT:end);
      drhohfT = drhohfT(1:1/dT:end);
      drhofT = drhofT(1:1/dT:end);
      hfTP(:,i) = hfT;
      drhohfTP(:,i) = drhohfT;
      drhofTP(:,i) = drhofT;
      rhofTP(:,i) = rho;
    endif
  endfor
  if 0
    drhohfPT = (rhofTP(:,2:end).*hfTP(:,2:end) - rhofTP(:,1:end-1).*hfTP(:,1:end-1))./dP(1:end-1)';
    drhohfPT = [drhohfPT(:,1),drhohfPT];
  endif
save -binary "hfTFile.m" p_i dT hfTP betaP
endif


T = Tu;
p = pInit;
[cp,lambda,M,rho,eta] = updtProp(T,p);
hu = hFromT(T,p);
rhou = rho;
etau = eta;
Mu = M;
lambdau = lambda;

%Initial cond
				%T = Tu*ones(nC,1);
				%T(1:round(ri/dr)) = Ti;
%T = Ti + (Tu - Ti)*(0.5+0.5*tanh((50*(r'-ri)/rmax)));%plot(r,T);keyboard
tanhProfile = (0.5-0.5*tanh((50*(r'-ri)/rmax)));
%tanhProfile /= sum(0.5-0.5*tanh((50*(r'-ri)/rmax)));
tanhProfile /= (0.5*ri^2);

T = Tu + (Ti - Tu)*exp(-pi/4*(r/ri).^6);
T.*=0; T+= Tu;
T0 = T;
p = T.*0 + p;
%Y = 0 + 4.29e-2*(0.5+0.5*tanh((50*(r'-ri)/rmax)));
burnt = (T>=Tadb);
[cp,lambda,M,rho,eta] = updtProp(T,p);
h = hFromT(T,p);

h0 = h;

%Solve diffusion equation from post-breakdown to T = 3 x Tadb
uf0 = ones(nF,1)*1e-15*0;
u = uf0(2:end);
uf = uf0;
ufavg0 = uf0;
ufavg = ufavg0;
u0 = u;
cp0 = cp;

pf = ones(nF,1).*pInit;
p = pf(2:end);
pp = p*0;
p0 = p;
pf = [p(1);0.5*(p(1:end-1)+p(2:end));pInit];%p(end-1) + (p(end) - p(end-1))*1.5];
gradp = (pf(2:end) - pf(1:end-1))/dr;
gradpf = [0;(p(2:end)-p(1:end-1))/dr;(pInit-p(end))/.5/dr];
gradpfavg = [gradp(1);0.5*(gradp(1:end-1)+gradp(2:end));gradp(end)];

      %uf0 = 0.9*uf0 + 0.1*uf0y;
      %plot(uf0)
      %
      %(qf - [0;qf(1:end-1)] + uFlux.*uf0 - [0;uFlux(1:end-1).*uf0(1:end-1)]);
      %keyboard
rho0 = rho;
rhof = interpolate(rho);
etaf = interpolate(eta);
Mf = interpolate(M);
rhof0 = rhof;
lambdaf = [lambda(1);0.5*(lambda(1:end-1)+lambda(2:end));lambdau];
M0 = M;
uf = uf0;
mFlux = rf.*uf.*rhof;
cpf = [cp(1);0.5*(cp(1:end-1)+cp(2:end));cp(end)];
t = 0;
i = 0;

subplot (3, 1, 1);
hold off
plot(0)
			
subplot (3, 1, 2)
hold off
plot(0)
			
subplot (3, 1, 3)
hold off
plot(0)

uDel = 10 + u*0;
pDel = uDel;
hDel = pDel;

if recov
  load -binary "recover.m" t h u uf ufavg T rho p M i
  h0 = h;
  u0 = u;
  uf0 = uf;
  T0 = T;
  rho0 = rho;
  p0 = p;
  M0 = M;
endif
pInitRes = 0; uInitRes = 0; %rho = p/8314*Mu/Tu; rho0 = rho; p = rho*8314./M.*T; p0 = p;
%h(1:10) *= 10; h0 = h; uInitRes = 0; pInitRes = 0; dt = 1e-7;
while t<t_end
  t += dt;
  i++;
  cc = 0;
  initSource = 243000000*exp(-t/150e-13).*(1-exp(-t/30e-13))*tanhProfile';
  initSource *= 0.5;
  %initSource = max(0,10000000*(1-t*1e8))*tanhProfile';
  URFu = 0.8;
  URFi_u = 1;
  URFp = 0.3;
  URFt = 0.5;%0.2;
  out = 0;
  do

    if i >= 100000
      uInitRes
      TInitRes
      pInitRes
      uy(1)
      uf(1)
      hy(1)
      mFlux(1)
      ppy(1)
      keyboard
    endif

    u_1 = u;
    p_1 = p;
    h_1 = h;


    %uEqn

    %eta = 0;
    momentumEqStaggered
%energyEquation


%TInitRes
%keyboard

%    continuityEquation
  %  pressureCorrEq
    
    %enthalpyEquation
%keyboard
    TEquation
    %pEqOF

    %pCorrEqn
    
    %p = p0.*M0./T0.*T./M;
    %keyboard
if 0
    Au = spdiags(ones(nC,1),0,nC,nC) - spdiags(ones(nC-1,1),-1,nC,nC);
    bu = -(rho - rho0)/dt*dr.*r;
    mFlux(2:nF) = Au\bu; %mFlux = min(mFlux,.02);
    u = 0.5*(mFlux(1:end-1)+mFlux(2:end))./rho./r; uf = mFlux./rhof./[1;rf(2:end)];
endif
    
    if 0
    uDel_1 = uDel;
    uDel = u - u_1;
    pDel_1 = pDel;
    pDel = p - p_1;
    hDel_1 = hDel;
    hDel = h - h_1;
    uLamb = norm(uDel)/norm(uDel_1 +1e-10);
    pLamb = norm(pDel)/norm(pDel_1 +1e-10);
    hLamb = norm(hDel)/norm(hDel_1 +1e-10);

    uErr = norm(uDel)/(uLamb-1);
    pErr = norm(pDel)/(pLamb-1);
    hErr = norm(hDel)/(hLamb-1);
    Err = max(abs([uErr,pErr,hErr]));
    endif

    Err = max([uInitRes;pInitRes;TInitRes]);


    %if uInitRes < 1e-4 && max(abs(pp./p)) < 1e-2 && TInitRes < 1e-5
if Err < tol % || max(abs(pp))<1
  out = 1;
endif
    
    cc++;
        %keyboard
  until ((cc >= itmax) || (out == 1)) % || max(abs(pp))<1e-3
%keyboard
  
%plot (r,T)
  %  pause(0.5)
 % keyboard
%uInitRes
%TInitRes
%pInitRes

%  (rho-rho0)/dt.*r*dr + mFlux(2:end) - mFlux(1:end-1);
%  (rho.*h - rho0.*h0)

  if i == 100
    dt = 5e-10;
  elseif i == 150
    dt = 5e-9;
  elseif i == 200
    dt = 1e-8;
  endif

  if cc == itmax
    disp("Falla de convergencia")
    keyboard
  endif
  if 0
  if i > 100
    cri = max([max(abs(T-T0)./T0),max(abs(p-p0)./p0),max(abs(u-u0)./(u0 + 10))]);
    %cri = max(abs(T-T0)./T0);
    if cri > .3 && i > icri + 5
      icri = i;
      %keyboard
    dt /= 1.5;
    elseif cri < .05 && i > icri + 5
      icri = i;
      %keyboard
    dt *= 1.1;
  endif
  endif
  dt = min(max(dt,1e-12),20e-6);
  endif
if 0
  if i == 50 || i ==100 || i == 150  || i == 300 || i == 450 || i == 600 || i == 750 || i == 900
    dt *= 5;
  endif
  if min(p) < pInit*.7
    dt = 1e-9;
  endif
endif
  
				%if p(1) < 295000 && t>100*dt
%    keyboard
 % endif
  
  h0 = h;
  u0 = u;
  uf0 = uf;
  ufavg0 = ufavg;
  T0 = T;
  rho0 = rho;
  p0 = p;
  M0 = M;
  cp0 = cp;

  save -binary "recover.m" t h u uf ufavg T rho p M i

  
  if (mod(i,1) == 0) && 1
    subplot (3, 1, 1);
    plot (r,T);
    title('T')  
    subplot (3, 1, 2)
    plot (r,u);
    title('u')  
    subplot (3, 1, 3)
    plot (r,p);
    title('p')  
    pause(0.5)
  endif

%  keyboard
   if (mod(i,1) == 0)&&0
    subplot (3, 1, 1);
    plot (r(1:20),T(1:20));
    title('T')  
    subplot (3, 1, 2)
    hold on
    plot (rf(1:20),uf(1:20));
    title('u')  
    subplot (3, 1, 3)
    plot (r(1:20),p(1:20));
    title('p')  
    pause(0.5)
  endif
  
  burnt = max(burnt,T>=Tadb); 
  if 0
    Tfil = sum(T>Tadb);
    rk(i+1) = r(Tfil) + (Tadb - T(Tfil))/(T(Tfil+1)-T(Tfil))*(r(Tfil+1) - r(Tfil));

    sk = (rk(i+1) - rk(i))/dt;
    if (max(T) > 3000)
      eta_e = 0.36 + 0.14*sk^3/(sk^3 + 500);
    else
      eta_e = 0.08 + 0.22*sk^3/(sk^3 + 700);  
    endif
				%D = 
				%Y = Le*rho*cp/lambda
    RR = 2*1.68e-3/rpCor(i);
    RRR = 1 + 1/RR + 1/RR^2 +2/3/RR^3;
    
    rp(i+1) = rp(i) + dt*(sk + 6.3*Su0*max(0.7,rp(i)/(rp(i)+2*1.68e-3)));
    rp2(i+1) = rp2(i) + dt*(sk + 6.3*Su0*max(0.9,rp(i)/(rp(i)+2*1.68e-3)));
    rpCor(i+1) = rpCor(i) + dt*(max(sk, 6.3*Su0/(RRR)));
  endif

endwhile

if 0
subplot (3, 2, 1);
plot (T(1:10));
title('T')  
subplot (3, 2, 2)
plot (u(1:10));
title('u')  
subplot (3, 2, 3)
plot (mFlux(1:10));
title('mFlux')  
subplot (3, 2, 4)
plot (p(1:10));
title('p')  
subplot (3, 2, 5)
plot (pp(1:10));
title('pp')  
subplot (3, 2, 6)
plot (gradp(1:10));
title('gradp')  
%endif

    subplot (1, 1, 1)
    plot (r,u);
    hold on
    plot (rf,uf);
    hold off
endif

