%
%            Helmholtz_opt.m
%
  
  clear all

%
% Load the data and initial the parameters.  Plot the energy with the initial parameters
% and compare to the experimental data.
%

  load Helmholtz.txt;

  P = Helmholtz(:,1);
  psi_data = Helmholtz(:,2);
  
  alpha_1 = -150;
  alpha_11 = 900;
  alpha_111 = 10;
  param_init = [alpha_1; alpha_11; alpha_111];

  psi_init = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6;

  figure(1)
  plot(P,psi_init,'-b',P,psi_data,'rx','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Polarization')
  ylabel('Helmholtz Energy')
  legend('Model','Data','Location','NorthWest')
  
%
% Optimize the parameters using the function energy_fun and Nelder-Mead optimization
%
  
  modelfun = @(param)energy_fun(param,P,psi_data);
  [param_opt,fval] = fminsearch(modelfun,param_init);

  alpha_1 = param_opt(1);
  alpha_11 = param_opt(2);
  alpha_111 = param_opt(3);

  psi = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6; 
  res = psi - psi_data;

%
% Plot the optimized energy and residuals.
%

  figure(2)
  plot(P,psi,'-b',P,psi_data,'rx','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Polarization')
  ylabel('Helmholtz Energy')
  legend('Model','Data','Location','NorthWest')

  figure(3)
  plot(P,res,'rx',P,0*P,'-k','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Polarization')
  ylabel('Helmholtz Energy')