clear all; close all; clc


Filename = '20151124_20151123_TCellOld_8_5_1200_1';
load(Filename);
p = pmin;
p(8) = [388.590678892474];
pinit = p;

%This will return the Vmitin parameter to the young model value 


load Tcelldata

% data for the model
%Set the initial conditions
x0(1) = 0.54;  %IP3
x0(2) = 0.07;  %PLCgamma
x0(3) = 10;    %R -- cannot find this initial condition
x0(4) = 0.05;  %Ca_cyt
x0(5) = 0.1;  %h - 0.02 written in Table 5-1, 0.1 written on page 218 in code
x0(6) = [280.449537337208];   %Ca_ER
x0(7) = 0.1;   %Ca_mit


tdata=1:1200;
tspan = [0:0.1:max(tdata)];

%Specify fit type
ff = fittype('exp2');
%General model Exp2:
%f2(x) = a*exp(b*x) + c*exp(d*x), Tau1 is b, which will be used below

%solve the system
[tsol1,xsol1] = ode23s(@(t,x)AK_f1_20151121_TCellOLD_pub_tofit_8Variables(t,x,p),tspan,x0);
%Find at what point the max occurs
qq=max(xsol1(:,4))
%Take the time step and divide by 10 to get the time in s    
peakT=(find(xsol1(:,4)==qq))/10;
%Specify the decay occurs after the peak time
decay = xsol1(10*peakT:end,4);
% size(decay)
% size(tsol1(10*peakT:end))

[c2,gof2] = fit(tsol1(10*peakT:end),decay,ff,'Robust', 'LAR',...
    'Lower',[0,-Inf,0,-Inf],...
    'Upper',[10,1,10,1]); %,...
tau1=c2.b;

tau1inM = 1./(60*tau1);

% save([Filename '_Peak_OLD' ],'peakT')
% save([Filename '_Max_OLD'],'qq')
% save([Filename '_Tau1_OLD'],'tau1')


save(['PeakTimeTauCombined_VmitinRemoved_' Filename])