
cc1 = r*dr/dt;
viscD = 2*eta.*(((uf(2:end) - uf(1:end-1))/dr).^2 + (u./r).^2 - 1/3*((uf(2:end).*rf(2:end) - uf(1:end-1).*rf(1:end-1))/dr./r).^2);


AT = spdiags([-lambdaf(2:nF).*rf(2:nF)/dr - max(mFlux(1:end-1),0).*[0;cp(1:end-1)], (cp.*cc1.*rho + [lambdaf(2:nF-1).*rf(2:nF-1)/dr;0] + lambdaf(1:nF-1).*rf(1:nF-1)/dr) + max(mFlux(2:end),0).*cp + max(-mFlux(1:end-1),0).*cp ,  -lambdaf(1:nF-1).*rf(1:nF-1)/dr- max(-mFlux(1:end-1).*cp,0)],-1:1,nC,nC);
%AT = spdiags([-lambdaf(2:end).*rf(2:end)/dr - max(mFlux(1:end-1).*[0;cp(1:end-1)]), (cp.*cc1.*rho + [lambdaf(2:end-1).*rf(2:end-1)/dr;lambdaf(end)*rf(end)/0.5/dr] + lambdaf(1:end-1).*rf(1:end-1)/dr) + max(mFlux(2:end),0).*cp + max(-mFlux(1:end-1),0).*cp,  -lambdaf(1:end-1).*rf(1:end-1)/dr - max(-mFlux(1:end-1).*cp,0)],-1:1,nC,nC); %T = 400 end


%beta = betaFunc(T,p);
beta = -1;

bT = cc1.*cp0.*T0.*rho0 + rho.*T.*r*dr.*((cp - cp0)/dt + u.*(cpf(2:end) - cpf(1:end-1))/dr) - beta.*((p - p0)/dt+ u.*gradp).*r*dr + viscD.*r*dr;
%bT(end) -= max(-mFlux(end),0)*hu;

%bT(end) += lambdaf(end)*rf(end)/0.5/dr*400;

bT += initSource.*r*dr;% + max(0,polyval(EgPol,t))*tanhProfile'.*r*dr + Qc.*r*dr/dt.*rho.*(1 - burnt).*(T>Tadb);
Ty = AT\bT;

TInitRes = res_compute(AT,bT,T);

T = (1 - URFt)*T + URFt*Ty;

TFinalRes = res_compute(AT,bT,T);
if cc == 130
  URFp = 0.8;
  URFt = 0.5;
elseif cc == 80
  URFp = 0.1;
  URFt = 0.1;
  %URFi_u = 0.2;
elseif cc == 50
  URFp = 0.3;
  URFt = 0.3;
  %URFi_u = 0.4;
elseif cc >=170
  uInitRes
  TInitRes
  pInitRes
  keyboard
endif

[cp,lambda,M,rho,eta] = updtProp(T,p);

%T = T*(1 - URFt) + URFt*TFromh(T,h,p);
%T = TFromh(T,h,p);
rhof = interpolate(rho);
cpf = interpolate(cp);

%eta = 0;

%T = h./cp;

%cp = 1000;% + T*.001;
%lambda = 0.02 + T.*0;% + T*.0001;
%M = 29 + T.*0;
