function generate_figures
global model
model = @Land2016;

addpath('./export_fig','./convert','./velocity_data');
set(0,'defaultlinelinewidth',1.5);
set(0,'DefaultAxesFontSize',10);
set(0,'DefaultTextFontSize',10);

humanCai = [0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.170 0.175 0.181 0.188 0.196 0.205 0.215 0.225 0.235 0.247 0.259 0.270 0.280 0.289 0.297 0.304 0.312 0.319 0.328 0.336 0.345 0.353 0.359 0.364 0.368 0.373 0.377 0.382 0.388 0.394 0.401 0.408 0.413 0.416 0.420 0.423 0.427 0.432 0.438 0.444 0.450 0.455 0.460 0.463 0.465 0.468 0.470 0.475 0.480 0.485 0.491 0.496 0.499 0.502 0.504 0.504 0.506 0.510 0.513 0.519 0.524 0.529 0.531 0.532 0.532 0.532 0.533 0.535 0.538 0.543 0.547 0.551 0.554 0.554 0.554 0.554 0.553 0.555 0.559 0.561 0.565 0.568 0.570 0.569 0.568 0.568 0.566 0.567 0.570 0.574 0.577 0.581 0.581 0.581 0.579 0.578 0.578 0.578 0.579 0.582 0.586 0.588 0.588 0.587 0.585 0.583 0.581 0.580 0.583 0.587 0.590 0.593 0.593 0.592 0.590 0.586 0.584 0.585 0.586 0.590 0.592 0.594 0.593 0.593 0.590 0.588 0.586 0.586 0.587 0.589 0.592 0.594 0.593 0.591 0.589 0.586 0.584 0.583 0.585 0.587 0.589 0.591 0.589 0.588 0.585 0.582 0.580 0.579 0.580 0.582 0.585 0.587 0.586 0.584 0.581 0.579 0.578 0.577 0.578 0.580 0.581 0.583 0.581 0.579 0.576 0.572 0.570 0.569 0.572 0.574 0.576 0.578 0.577 0.574 0.572 0.569 0.567 0.566 0.567 0.569 0.570 0.570 0.569 0.567 0.563 0.560 0.558 0.557 0.558 0.559 0.561 0.562 0.563 0.560 0.557 0.554 0.551 0.550 0.550 0.552 0.553 0.555 0.553 0.550 0.546 0.544 0.541 0.541 0.542 0.543 0.545 0.545 0.543 0.541 0.537 0.534 0.531 0.531 0.530 0.532 0.534 0.535 0.533 0.530 0.527 0.523 0.520 0.518 0.518 0.519 0.521 0.521 0.520 0.518 0.515 0.512 0.509 0.508 0.508 0.508 0.510 0.510 0.508 0.505 0.502 0.499 0.497 0.496 0.496 0.496 0.496 0.496 0.496 0.493 0.490 0.486 0.484 0.483 0.483 0.483 0.484 0.483 0.481 0.479 0.476 0.472 0.469 0.467 0.467 0.467 0.468 0.469 0.466 0.464 0.460 0.457 0.454 0.452 0.453 0.453 0.454 0.454 0.453 0.450 0.447 0.444 0.441 0.439 0.439 0.440 0.441 0.441 0.440 0.437 0.433 0.430 0.427 0.426 0.426 0.426 0.426 0.426 0.424 0.421 0.418 0.415 0.413 0.411 0.410 0.411 0.412 0.411 0.409 0.407 0.404 0.401 0.399 0.398 0.397 0.397 0.397 0.397 0.396 0.393 0.390 0.387 0.384 0.383 0.383 0.383 0.384 0.384 0.382 0.380 0.377 0.374 0.372 0.370 0.369 0.370 0.370 0.370 0.368 0.366 0.364 0.361 0.358 0.357 0.357 0.357 0.357 0.356 0.355 0.353 0.350 0.348 0.345 0.344 0.343 0.344 0.344 0.344 0.343 0.341 0.338 0.336 0.333 0.332 0.332 0.332 0.332 0.332 0.332 0.330 0.327 0.325 0.323 0.322 0.321 0.321 0.321 0.321 0.320 0.318 0.316 0.313 0.311 0.310 0.310 0.310 0.311 0.311 0.310 0.308 0.306 0.303 0.301 0.301 0.300 0.300 0.301 0.301 0.300 0.298 0.296 0.294 0.292 0.291 0.291 0.291 0.292 0.292 0.291 0.289 0.287 0.285 0.284 0.283 0.282 0.283 0.283 0.282 0.281 0.280 0.279 0.277 0.275 0.274 0.274 0.274 0.275 0.274 0.273 0.272 0.270 0.269 0.268 0.267 0.267 0.267 0.267 0.267 0.266 0.265 0.263 0.261 0.260 0.260 0.260 0.260 0.260 0.260 0.259 0.258 0.257 0.255 0.253 0.253 0.253 0.253 0.254 0.254 0.253 0.252 0.250 0.248 0.247 0.246 0.246 0.247 0.247 0.247 0.247 0.246 0.244 0.243 0.241 0.241 0.241 0.242 0.242 0.242 0.242 0.240 0.239 0.237 0.236 0.236 0.236 0.236 0.237 0.237 0.236 0.236 0.234 0.233 0.232 0.231 0.232 0.232 0.233 0.232 0.232 0.231 0.229 0.228 0.227 0.226 0.226 0.227 0.228 0.228 0.227 0.226 0.225 0.224 0.223 0.222 0.222 0.223 0.223 0.223 0.223 0.222 0.221 0.220 0.219 0.218 0.219 0.219 0.220 0.220 0.220 0.219 0.218 0.217 0.216 0.215 0.215 0.216 0.216 0.216 0.216 0.215 0.214 0.213 0.212 0.211 0.212 0.212 0.213 0.213 0.213 0.212 0.211 0.209 0.209 0.209 0.209 0.209 0.210 0.210 0.209 0.208 0.207 0.206 0.205 0.205 0.205 0.206 0.206 0.206 0.206 0.206 0.205 0.204 0.203 0.203 0.203 0.204 0.204 0.204 0.204 0.203 0.202 0.201 0.200 0.200 0.200 0.201 0.201 0.201 0.201 0.200 0.199 0.199 0.198 0.197 0.198 0.198 0.199 0.199 0.199 0.198 0.197 0.196 0.195 0.195 0.196 0.196 0.197 0.197 0.197 0.196 0.195 0.194 0.194 0.193 0.194 0.194 0.195 0.195 0.195 0.194 0.193 0.192 0.191 0.191 0.191 0.192 0.192 0.193 0.192 0.192 0.191 0.190 0.190 0.189 0.190 0.190 0.191 0.191 0.191 0.190 0.189 0.189 0.188 0.188 0.188 0.188 0.189 0.189 0.189 0.188 0.187 0.187 0.186 0.186 0.186 0.187 0.187 0.188 0.187 0.187 0.186 0.185 0.185 0.184 0.185 0.185 0.186 0.186 0.186 0.185 0.184 0.184 0.183 0.183 0.183 0.184 0.185 0.185 0.185 0.184 0.183 0.182 0.182 0.181 0.182 0.182 0.183 0.184 0.184 0.183 0.182 0.181 0.181 0.181 0.181 0.181 0.182 0.182 0.182 0.182 0.181 0.180 0.179 0.179 0.180 0.180 0.181 0.181 0.181 0.180 0.180 0.179 0.178 0.178 0.178 0.179 0.180 0.180 0.180 0.179 0.179 0.178 0.177 0.177 0.177 0.178 0.179 0.179 0.179 0.178 0.178 0.177 0.176 0.176 0.176 0.177 0.178 0.178 0.178 0.177 0.176 0.175 0.175 0.175 0.175 0.176 0.176 0.177 0.177 0.177 0.176 0.175 0.175 0.175 0.175 0.176 0.176 0.177 0.176 0.176 0.175 0.174 0.174 0.174 0.174 0.174 0.175 0.175 0.175 0.175 0.174 0.173 0.173 0.173 0.173 0.174 0.174 0.174 0.174 0.174 0.173 0.172 0.172 0.172 0.172 0.172 0.173 0.173 0.173 0.173 0.172 0.171 0.171 0.171 0.171 0.172 0.172 0.173 0.173 0.172 0.171 0.171 0.170 0.170 0.170 0.171 0.172 0.172 0.172 0.172 0.171 0.170 0.170 0.169 0.170 0.171 0.171 0.171 0.171 0.171 0.170 0.169 0.169 0.169 0.169 0.170 0.170 0.171 0.171 0.170 0.170 0.169 0.168 0.168 0.168 0.169 0.170 0.170 0.170 0.170 0.169 0.168 0.168 0.167 0.168 0.169 0.169 0.170 0.170 0.169 0.168 0.168 0.167 0.167 0.168 0.168 0.169 0.169 0.169 0.169 0.168 0.167 0.167 0.167 0.167 0.168 0.168 0.169 0.169 0.168 0.167 0.167 0.166 0.167 0.167 0.168 0.168 0.168 0.168 0.168 0.167 0.166 0.166 0.166 0.166 0.166 0.167 0.168 0.167 0.167 0.166 0.166 0.165 0.165 0.166 0.167 0.167 0.168 0.168 0.167 0.166 0.166 0.165 0.165 0.166 0.166 0.167 0.167 0.167 0.167 0.166 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.167 0.166 0.166 0.165 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.166 0.166 0.165 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.165 0.166 0.167 0.167 0.166 0.165 0.165 0.164 0.164 0.165 0.165 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.164 0.165 0.166 0.166 0.166 ];
t_ca = 0:1000;
[t,y] = ode15s(@(t,y) model(t,y,humanCai,'intact',1,0), 0:1:(10*1000), model() );
tp('Model stats:',y(end-1000:end,1)*480);

lambdas = [0.9 0.95 1 1.05 1.1];
figure(101); 
clf; 
subplot(1,2,2); hold on;
for l = lambdas
  [t,y] = ode15s(@(t,y) model(t,y,humanCai,'intact',l,0), 0:1:(10*1000), model() );
  if l==1; lw=2; else lw=1; end
  plot(t-9e3,y(:,1)*480,'LineWidth',lw);
end

xlabel('Time (ms)');
ylabel('Force (kPa)');
legend( arrayfun(@(l) sprintf('\\lambda = %g',l),lambdas,'UniformOutput',0));
xlim([0 1000]);
set(gca,'XTick',[0 250 500 750 1000],'YTick',[0 25 50 75 100]);

subplot(1,2,1); hold on;
plot(t_ca,humanCai);
xlabel('Time (ms)');
ylabel('Intracellular calcium concentration (\mu{}M)');

ylim([0 0.6]);
set(gca,'XTick',[0 250 500 750 1000],'YTick',[0 0.15 0.3 0.45 0.6]);

myfig2(figure(101), './Fig_fit_ca_twitch.png',[7 2.5]*2,225);


% --- velocity
figure(102); clf;

t={}; tension={}; traces={};
[t{1},tension{1},traces{1}]  = plot_data(0);
[t{2},tension{2},traces{2}]  = plot_data(1);
[t{3},tension{3},traces{3}]  = plot_data(2);

Cai = 30;

[~,y] = ode15s(@(t,y) model(t,y,Cai,'skinned',1,0),0:10000,model());
y00 = y;   

[score, tension_out,rel] = score_par(y00,t,tension);

figure(102); clf;
for i=1:3
 subplot(1,4,i); hold on;
 plot(t{i},tension{i}/30,'Color',[ 0    0.4470    0.7410],'LineWidth',2);
 plot(t{i}(1:end-1),tension_out{i}/30,'Color',[ 0.8500    0.3250    0.0980],'LineWidth',2);
 
 for j=1:length(traces{i})
  plot(t{i},traces{i}{j}/30,'Color',[0.05 0.05 0.05]*(j+4),'LineWidth',1);
 end
 plot(t{i},tension{i}/30,'Color',[ 0    0.4470    0.7410],'LineWidth',1);
 plot(t{i}(1:end-1),tension_out{i}/30,'Color',[ 0.8500    0.3250    0.0980],'LineWidth',2);
end
legend('Mean data','Model fit','Individual traces');
subplot(1,4,1); ylabel('Normalized force');

for i=1:3
 subplot(1,4,i);    
 xlabel('Time (ms)');
 xlim([-75 2000]); set(gca,'XTick',[0 500 1000 1500 2000],'YTick',[0.5 0.75 1 1.25 1.5]);
 ylim([0.5 1.6]);
end

%myfig2(figure(102), './Fig_fit_vel.png',[7 2.5]*2,225);

% --- constant shortening

files = {'constant_shortening/16318001_corrd.abf','constant_shortening/16318002_corrd.abf','constant_shortening/16318003_corrd.abf','constant_shortening/16318004_corrd.abf','constant_shortening/16318005_corrd.abf','constant_shortening/16318006_corrd.abf','constant_shortening/16318007_corrd.abf','constant_shortening/16318008_corrd.abf','constant_shortening/16318009_corrd.abf','constant_shortening/16318010_corrd.abf','constant_shortening/16318011_corrd.abf'};


area = [52*13 52*13 52*13 52*13 52*13 52*13 29*14 29*14 29*14 29*14 29*14];

sl0 = [2.1 2.1 2.06  2.1  2.3*ones(1,11)];
slrest = 1.9;

length_change =  [8.5 8.5 8.5 10.0 11.4 11.4 10.0 8.5 7.1 7.1 8.5]' * [1 1 1 1]; % regular, 1-11
                  
length_dur = [ 100	200	300 0
              100	200	300 0
              150 200 250 300
              150 200 250 300
              150 200 250 300
              100 200 300 0
              100 200 300 0
              100 200 300 0
              100 200 300 0 % 009
              150 200 250 300
              150 200 250 300
              150 200 250 300 % 011
               ];
             
          
t0 = 8310;

tsh = 11000;
trx_sh = 14000;
trx_long  = 19500;

          
x = 0;
lambdas = {};
fnorms = {};
ts = {}; tend={};
lx=[]; fe=[];
for i=1:length(files)
 [d,si,h] = abf2load(files{i});  
 t =  (0:size(d,1)-1)*(si/1000);
 f = 1000 * squeeze(d(:,1,:)) / area(i);
 l = squeeze(d(:,2,:)); 
 
 x0 = x+1;
 fni = [];
 for j=1:size(f,2)
   x=x+1; 
    ts{x} = t;
    [~,ti] = min(abs(t-(t0+length_dur(i,j))));
    [~,tsi] = min(abs(t-t0));
    tend{x} = t(ti);
    fprintf('%d/%d: mean f = %f\n',i,j,mean(f(t < t0-5 & t > t0-40,j)));
    
    frx_sh(x) = mean( f( abs(t-trx_sh)<50 ,j) );
    frx_long(x) = mean( f( abs(t-trx_long)<50 ,j) );
    fsh(x) = mean( f( abs(t-tsh)<50 ,j) );
    fe(x) = mean( f(ti:ti+100,j)); % ~ 20ms
    
    fs(x) = mean( f(tsi-100:tsi,j)); % ~ 20ms
    
    cor_fe(x) = (fe(x) - frx_sh(x)) ./ (fs(x) - frx_long(x));
    dur(x) = length_dur(i,j);     lx(x) = length_change(i,j)/100;
 end
end

figure(102); subplot(1,4,4)
plot(lx./dur,cor_fe,'o', 0.1./(100:10:300), rel(:,1)./rel(:,2), 'x-',[0.2/200 0.2/150 0.2/150 0.2/200 0.2/200],[0.9 0.9 0.5 0.5 0.9],'--'); % [100 200 300]
ylim([0 0.91]);
xlabel('Shortening velocity (/ms)'); ylabel('Relative force at end of shortening'); legend('Experimental data','Skinned Model prediction','Expected force in organ','Location','NorthWest')


ax1 = gca;
ax1_pos = ax1.Position; % position of first axes
ax2 = axes('Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','left',...
    'Color','none');
xlabel('Sarcomere length shortening speed (\mu{}m/s)')
scale = 2.3*1000;

ax2.XLim = ax1.XLim*scale;
ax2.XTick=ax1.XTick*scale
ax2.YTick=[];


myfig2(figure(102), './Fig_fit_vel_cmb.png',[10 2.5]*2,225); 



function [tpt,rt50,rt95,mif,mxf] = tp(tag,fi)
 [mxf,tpt] = max(fi);
 f0 = fi-min(fi);   
 f0(1:tpt) = f0(tpt);
 t50 = find(f0 <= 0.5  * max(f0),1);
 t90 = find(f0 <= 0.1  * max(f0),1);
 t95 = find(f0 <= 0.05 * max(f0),1);
 if isempty(t95); t95=1000; end
 onset = find(diff(fi)>0.01,1); % ?
 fprintf('%s\tTPT: %.1f\tRT50: %.1f\tRT90: %.1f\tRT95: %.1f\tMax %.1f\tMin %.3f\n',tag,tpt-onset,t50-tpt,t90-tpt,t95-tpt,mxf,min(fi));
 rt50 = t50-tpt;
 rt95 = t95-tpt;
 tpt = tpt-onset;
 mif = min(fi);
 

  
  
function [score,tension_out,rel,allrel] = score_par(y00,t,tension)
  global model
  stretch = [0.5 1 2];
  step_dur = 10;
  step_t   = 2000;        
  Cai = 30;
  xb0   = [0.95993   0.97022  0.9748];
  
  ldp  = [ 5.25  6.6 4.5;
          -2.4 -2.4 -2.4];       
  tscale= [0.8498    0.8396    0.8388];
  dif = [0 0 0]; 
       
  for i=1:3
    Tfit = tension{i};
    flen = ones(size(Tfit));
    flen(step_dur:step_t) = stretch(i)/100 + 1; % stretched
    st =  interp1([0 step_dur],[1 flen(6+step_dur)],0:step_dur);
    flen(1:step_dur+1) = st;
    flen(step_t:step_t+step_dur) = st(end:-1:1);

    ti = t{i};
    f = @(tt,y) model(tt,y,Cai,'skinned',interp1(ti,flen,tt), (interp1(ti,flen,tt+1e-3)-interp1(ti,flen,tt))/1e-3 , ldp(:,i)); % TODO: LD
    y0 = y00( find(y00(:,1)/y00(end,1) < xb0(i), 1,'last'), :) % 
    
    [tt,y] = ode15s(f,t{i}(1:end-1),y0,odeset('MaxStep',1)); % reduce xb0
    tatp = cell2mat(cellfun(@(ty) get_ta(f,ty),num2cell([tt y],2),'UniformOutput',0));
    ta = tatp(:,1);
    tp = tatp(:,2);

    tot_tension = (ta + tp)*tscale(i);
    tension_out{i} = tot_tension;
    [mt(i),mi] = max(tension{i});
    mm(i) = max(tot_tension);
    mm0(i) = tot_tension(1);
    dif(i) = norm(tot_tension - tension{i}(1:end-1));

  end
  q=1; i=1;
  goal = [0.2 0.4 0.5];
  for te = 100:10:300
    f2 = @(tt,y) model(tt,y,Cai,'skinned', flen(1) * (1 - (0.1/te)*tt) , -0.1/te , ldp(:,end));
    [tt2,y2] = ode15s(f2,0:te,y0,odeset('MaxStep',1)); % reduce xb0
    tatp2 = cell2mat(cellfun(@(ty) get_ta(f,ty),num2cell([tt2 y2],2),'UniformOutput',0));
    ta2 = tatp2(:,1);
    rel(i,:) = [ta2(end) ta2(1)]; i=i+1;
    if te==100||te==200||te==300
      shs(q) = 50*abs( (ta2(end)/ta2(1)) - goal(q) );
      q=q+1;
    end
  end 
  fprintf('score = %f + %f\n',sum(dif),sum(abs(shs)))
  score = sum(dif) + sum(abs(shs));

  
  
function [t,ta,tp]=get_ta(model,ty)
  [~,~,ta,tp] = model(ty(1),ty(2:end));
  t = [ta tp];
  
