%
%  Project 2
%
%                     Problem 1   
%

  clear all

%
%  Load the aluminum data from the file final_al_data.txt and construct the x datapoints.
%

  xdata = [10 14 18 22 26 30 34 38 42 46 50 54 58 62 66];
  u_amb = 21.29;

% Input dimensions and material constants for the aluminum rod.

  a = 0.95;       % cm
  b = 0.95;       % cm
  L = 70.0;       % cm
  k = 2.37;       % W/cm C
  h = 0.00191;
  Phi = -18.41;
  n = 15;         % Number of measurements

  dk = 1e-8;      % Stepsize in k
  dh = 1e-10;     % Stepsize in h    
  dPhi = 1e-4;    % Stepsize in Phi
  k_p = k + dk;   % Perturbed k
  h_p = h + dh;   % Perturbed h
  Phi_p = Phi + dPhi; % Perturbed Phi

%
%  Compute constants and analytic solution to the steady state heat equation.
%

  gamma = sqrt(2*(a+b)*h/(a*b*k));
  f1 = exp(gamma*L)*(h + k*gamma);
  f2 = exp(-gamma*L)*(h - k*gamma);
  f3 = f1/(f2 + f1);
  c1 = -Phi*f3/(k*gamma);
  c2 = Phi/(k*gamma) + c1;
  uvals_data = c1*exp(-gamma*xdata) + c2*exp(gamma*xdata) + u_amb;

%
% Compute constants and analytic and finite-difference sensitivities with respect to h.
%

  gamma_h = (1/(2*h))*gamma;
  f1_h = exp(gamma*L)*(gamma_h*L*(h+k*gamma) + 1 + k*gamma_h);
  f2_h = exp(-gamma*L)*(-gamma_h*L*(h-k*gamma) + 1 - k*gamma_h);
  den2 = (f1+f2)^2;
  f3_h = (f1_h*(f1+f2) - f1*(f1_h+f2_h))/den2;
  f4 = Phi/(k*gamma*gamma);
  c1_h = f4*gamma_h*f3 - (Phi/(k*gamma))*f3_h;
  c2_h = -f4*gamma_h + c1_h;
  uvals_h_data = c1_h*exp(-gamma*xdata) + c2_h*exp(gamma*xdata) + gamma_h*xdata.*(-c1*exp(-gamma*xdata) + c2*exp(gamma*xdata));
  gamma_hp = sqrt(2*(a+b)*h_p/(a*b*k));
  f1_hp = exp(gamma_hp*L)*(h_p + k*gamma_hp);
  f2_hp = exp(-gamma_hp*L)*(h_p - k*gamma_hp);
  f3_hp = f1_hp/(f2_hp + f1_hp);
  c1_hp = -Phi*f3_hp/(k*gamma_hp);
  c2_hp = Phi/(k*gamma_hp) + c1_hp;
  uvals_data_hp = c1_hp*exp(-gamma_hp*xdata) + c2_hp*exp(gamma_hp*xdata) +u_amb;
  h_fd = (1/dh)*(uvals_data_hp - uvals_data);

%
% Compute constants and analytic and finite-difference sensitivities with respect to Phi.
%

  c1_Phi = -(1/(k*gamma))*f3;
  c2_Phi = (1/(k*gamma)) + c1_Phi;
  uvals_Phi_data = c1_Phi*exp(-gamma*xdata) + c2_Phi*exp(gamma*xdata);
  c1_Phip = -Phi_p*f3/(k*gamma);
  c2_Phip = Phi_p/(k*gamma) + c1_Phip;
  uvals_data_Phip = c1_Phip*exp(-gamma*xdata) + c2_Phip*exp(gamma*xdata) + u_amb;
  Phi_fd = (1/dPhi)*(uvals_data_Phip - uvals_data);

%
% Compute constants and finite-difference sensitivities with respect to k.
%

  gamma_k = sqrt(2*(a+b)*h/(a*b*k_p));
  f1_hk = exp(gamma_k*L)*(h + k_p*gamma_k);
  f2_hk = exp(-gamma_k*L)*(h - k_p*gamma_k);
  f3_hk = f1_hk/(f2_hk + f1_hk);  
  c1_hk = -Phi*f3_hk/(k_p*gamma_k);
  c2_hk = Phi/(k_p*gamma_k) + c1_hk;
  uvals_data_kp = c1_hk*exp(-gamma_k*xdata) + c2_hk*exp(gamma_k*xdata) + u_amb;
  k_fd = (1/dk)*(uvals_data_kp - uvals_data);


%
% Construct the 3x3 sensitivity matrix sens_mat3 and matrix V.  Compute the rank of V.
%

  sens_mat3 = [uvals_Phi_data; uvals_h_data; uvals_data_kp];
  V3 = sens_mat3*sens_mat3';
  rank(V3)

%
% Construct the 2x2 analytic and finite difference sensitivity matrices.
%
  sens_mat = [uvals_Phi_data; uvals_h_data];
  sens_mat_fd = [Phi_fd; h_fd];

%
% Construct the measurement covariance sigma2 and the matrices
% V and V_fd constructed using the analytic and finite difference
% sensitivity relations.
%
  format short e

  V = inv(sens_mat*sens_mat')
  V_fd = inv(sens_mat_fd*sens_mat_fd')

%
% Plot the sensitivities with respect to Phi, h and k.
%

  figure(1)
  plot(xdata,uvals_Phi_data,'o',xdata,Phi_fd,'x','linewidth',3)
  axis([0 70 -5 0])
  set(gca,'Fontsize',[20]);
  xlabel('Distance (cm)')
  ylabel('Sensitivity in \Phi')
  legend('Analytic','Finite Difference','Location','Southeast')

  figure(2)
  plot(xdata,uvals_h_data,'o',xdata,h_fd,'x','linewidth',3)
  axis([0 70 -3.5e4 -0.5e4])
  set(gca,'Fontsize',[20]);
  xlabel('Distance (cm)')
  ylabel('Sensitivity in h')
  legend('Analytic','Finite Difference','Location','Southeast')

  figure(3)
  plot(xdata,uvals_data_kp,'x','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('Distance (cm)')
  ylabel('Sensitivity in k')