clear all

  m = 2.7;
  k = 8.5;
  c = 0.24;

  pt1 = 1.25;
  pt2 = 2.75;
  omega = pt1:.001:pt2;
  omega_n = length(omega);

  sample_n = 10000;
  
  var_m = 0.002^2;
  var_k = 0.001^2;
  var_c = 0.065^2;

  X = zeros(sample_n,3);
  X(:,1) = m + sqrt(var_m)*randn(sample_n,1);
  X(:,2) = k + sqrt(var_k)*randn(sample_n,1);
  X(:,3) = c + sqrt(var_c)*randn(sample_n,1);

  V = diag([var_m var_k var_c]);

  for i=1:omega_n
    denom1 = sqrt((k-m*omega(i)^2)^2 + (c*omega(i))^2);
    denom2 = denom1^3;
    UF_mean(i) = 1/denom1;
    MC = 1./sqrt((X(:,2)-X(:,1)*omega(i)^2).^2 + (X(:,3)*omega(i)).^2);
    MC_ascend = sort(MC);
   
    dfdm = (k-m*omega(i)^2)*(omega(i)^2)/denom2;
    dfdk = -(k-m*omega(i)^2)/denom2;
    dfdc = -c*(omega(i)^2)/denom2;

    cvec = [dfdm; dfdk; dfdc];
    Sigma_SA(i) = sqrt(cvec'*V*cvec);
    SA99(i) = UF_mean(i) + 2.326*Sigma_SA(i); 

    [mu_UF,sigma_UF] = normfit(MC);
    Sigma_MC0(i) = sigma_UF;
    UF_MC(i) = mu_UF;
    UF_MC2(i) = sum(MC)/sample_n;
    Sigma_MC(i) = sqrt(sum((MC-UF_MC2(i)).^2)/(sample_n-1));
    MC95(i) = MC_ascend(.95*sample_n+1) - UF_MC2(i);
    MC05(i) = UF_MC2(i) - MC_ascend(.05*sample_n+1); 
  end  

  omega0 = 1.60;
  denom1 = sqrt((k-m*omega0^2)^2 + (c*omega0)^2);
  denom2 = denom1^3;
  mu_SA0 = 1/sqrt((k-m*omega0^2)^2 + (c*omega0)^2);
  MC0 = 1./sqrt((X(:,2)-X(:,1)*omega0^2).^2 + (X(:,3)*omega0).^2);
  dfdm = (k-m*omega0^2)*(omega0^2)/denom2;
  dfdk = -(k-m*omega0^2)/denom2;
  dfdc = -c*(omega0^2)/denom2;
  cvec = [dfdm; dfdk; dfdc];
  Sigma_SA0 = sqrt(cvec'*V*cvec);
  x0 = 0.57:.001:0.66;
  SA_fit0 = normpdf(x0,mu_SA0,Sigma_SA0);
  MC_density0 = ksdensity(MC0,x0);
  [bandwidth0,MC_density0,mesh0,cdf0] = kde(MC0,128);
%  [mu0,sigma0] = normfit(MC0); 
%  MC_density0 = normpdf(x0,mu0,sigma0);
  

  omega1 = 1.77;
  denom1 = sqrt((k-m*omega1^2)^2 + (c*omega1)^2);
  denom2 = denom1^3;
  mu_SA1 = 1/sqrt((k-m*omega1^2)^2 + (c*omega1)^2);
  MC1 = 1./sqrt((X(:,2)-X(:,1)*omega1^2).^2 + (X(:,3)*omega1).^2);
  dfdm = (k-m*omega1^2)*(omega1^2)/denom2;
  dfdk = -(k-m*omega1^2)/denom2;
  dfdc = -c*(omega1^2)/denom2;
  cvec = [dfdm; dfdk; dfdc];
  Sigma_SA1 = sqrt(cvec'*V*cvec);
  x1 = 0.0:.005:6;
  SA_fit1 = normpdf(x1,mu_SA1,Sigma_SA1);
  MC_density1 = ksdensity(MC1,x1);

  figure(1)
  plot(omega,UF_MC2,'--b',omega,UF_mean,'-r','linewidth',3)
  axis([pt1 pt2 0 3.0])
  set(gca,'Fontsize',[22]);
  xlabel('\omega_F')
  ylabel('Mean')
  legend('Sampling Mean','Perturbation Mean','Location','NorthEast')
  %print -depsc fig5_4a

  figure(2)
  plot(omega,Sigma_MC,'--b',omega,Sigma_SA,'-r','linewidth',3)
  axis([pt1 pt2 0 1.2])
  set(gca,'Fontsize',[22]);
  xlabel('\omega_F')
  ylabel('Standard Deviation')
  legend('\sigma_s','\sigma_p','Location','NorthEast')
  %print -depsc fig5_4b

  figure(3)
  plot(mesh0,MC_density0,'b--','linewidth',3) 
  hold on
  plot(x0,SA_fit0,'-r','linewidth',3)
  hold off
%  plot(x0,MC_density0,'--b',x0,SA_fit0,'-r','linewidth',3)
  axis([.55 .68 0 50]) 
  set(gca,'Fontsize',[22]);
  xlabel('Response y')
  ylabel('PDF')
  %title('\omega_F = 1.0')
  legend('Sampling Method','Perturbation Method')
  %print -depsc fig5_4c

  figure(4)
  plot(x1,MC_density1,'--b',x1,SA_fit1,'r-','linewidth',3)
  hold off
  set(gca,'Fontsize',[22]);
  xlabel('Response y')
  ylabel('PDF')
  %title('\omega_F = 1.138')
  legend('Sampling Method','Perturbation Method')
%  print -depsc fig5_4d
  
  figure(5)
  errorbar(omega(1:40:1501),UF_MC2(1:40:1501),MC05(1:40:1501),MC95(1:40:1501),'linewidth',3)
  axis([pt1 pt2 0 4.5])
  set(gca,'Fontsize',[22]);
  xlabel('\omega_F')
  ylabel('Response y')
%  title('95% Confidence Intervals')