%   
%                  spring_sensitivity.m
%

  clear all


%
%  Compute the sensitivities of the sping model (8.2) using the sensitivity equations
%  (8.14) and analytic solutions (8.6).  Here N=p=2. 
%

%
% Set parameters and initial conditions
%

  K = 9;
  C = 0.05;
  params = [K C]; 

  z0 = 2;
  z0_dot = -C;
  z_k0 = 0;
  z_k0_dot = 0;
  z_c0 = 0;
  z_c0_dot = -1;
  

  tf = 6;
  dt = 0.1;
  t_vals = 0:dt:tf;
  Y0 = [z0; z0_dot; z_k0; z_k0_dot; z_c0; z_c0_dot];

%
% Integrate the coupled system using ode45.m.  There are now N(p+1) coupled
% differential equations comprised of the state and sensitivity equations
%

  ode_options = odeset('RelTol',1e-8);
  [t,Y] = ode45(@spring_rhs,t_vals,Y0,ode_options,params);

%
% Extract the states and sensitivities.
%

  z = Y(:,1);  
  z_dot = Y(:,2);
  z_k = Y(:,3);
  z_k_dot = Y(:,4);
  z_c = Y(:,5);
  z_c_dot = Y(:,6);

%
% Compute the analytic sensitivities. 
%

  f0 = sqrt(4*K - C^2);
  f1 = sqrt(K - C^2/4)*t_vals;
  f2 = exp(-C*t_vals/2);
  f3 = (-2*t_vals/f0).*sin(f1);
  f4 = (C*t_vals/f0).*sin(f1) - t_vals.*cos(f1);
  z_k_anal = f2.*f3;
  z_c_anal = f2.*f4; 
 
%
% Plot z_k(t) and z_c(t).
%

  figure(1)
  plot(t,z_k,'-b',t,z_k_anal,'--r',t,0*t,'k-','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Time')
  ylabel('z_K')
  legend('Sensitivity Eq','Analytic','Location','North')
  %print -depsc z_K

  figure(2)
  plot(t,z_c,'-b',t,z_c_anal,'--r',t,0*t,'k-','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Time')
  ylabel('z_C')
  legend('Sensitivity Eq','Analytic','Location','North')
  %print -depsc z_C