% muscle load velocity curves

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

% Hill data
g = [2.6,8.05,16.1,27,40.6,51.8,66]';
v = [4.05,2.96,1.84,1.17,0.6,0.31,0]';
 
 %fit the data to (p+a)v=b(p_0-p)

A=[v,g,-ones(7,1)];
R = -g.*v;

B=A'*A;
r=A'*R;
X = B\r;
a = X(1)
b = X(2)
p0 = X(3)/b

p = [0:0.01:1]*p0;

V = b*(p0-p)./(p+a);
figure(1)
plot(g,v,'o',g,v,'.',p,V,'linewidth',2)
axis([0 70 0 5])
xlabel('Load, g','fontsize',20)
ylabel('Velocity of shortening, cm/sec','fontsize',20)

if(cmykflg==1)
  print('../../figs_c/chapt_9/fig18-7_redone','-deps','-cmyk')
 end 

 % now make some plots for the Huxley model
 
 xm = [-3:0.01:0];
 xp = [0:0.01:1];
 G2 = 3.919;
 F1 = 13/16;
 
 %specify phi/V
 Vlist = [0.0001;0.6;2;4];
 for j = 1:4
     Vp = Vlist(j);
 
 
 nI = F1*(1-exp(-1/Vp))*exp(xm*G2/(2*Vp));
 nII = F1*(1-exp((xp.^2-1)/Vp));
 
 
 figure(2)
 plot([xm,xp],[nI, nII],[0,0],[0,1],'--','linewidth',2)
 xlabel('x/h','fontsize',20)
 ylabel('n(x)','fontsize',20)
 set(gca,'XTICK',-3:0.5:1)
 hold on
 end
 text(0.05,0.85,'V=0','fontsize',20)
 text(0.05,0.2,'V=V_{max}','fontsize',20)
 text(0.05,0.34,'V=0.5V_{max}','fontsize',20)
 text(0.05,0.68,'V=0.15V_{max}','fontsize',20)
 hold off
 if(cmykflg==1)
  print('../../figs_c/chapt_9/fig18-13_redone','-deps','-cmyk')
 end 
 
 % calculate the load velocity curve
 Vh = [0:0.01:4];
 F = 1-Vh.*(1-exp(-1./Vh)).*(1+Vh/(2*G2^2));
 
 figure(3)
 plot(F,Vh/4, g/p0,v/5,'o',p/p0,V/5,'--','linewidth',2)
 axis([0 1 0 1])
 ylabel('V/V_{max}','fontsize',20)
 xlabel('Load/Load_{max}','fontsize',20)

 if(cmykflg==1)
  print('../../figs_c/chapt_9/fig18-14_redone','-deps','-cmyk')
 end 