% Look at convergence of the energy calculation, both for uniform charge density
% and for the conducting rod.

function return_condition = convergence()

xmin = -0.5; % coordinate of left end or rod
xmax = +0.5; % coordinate of right end of rod

Ntrials = 9
maxiter = 2000; % maximum number of iterations
tol = sqrt(eps); % tolerance for convergence

cpu_sec=0; % initialize array

figure(7)
hold off
title('Charge displacement')
xlabel('xinitial')
ylabel('xfinal - xinitial')
plot([xmin,xmax],[0,0],'k-')
hold on
% Now ready to accept an overplot with each trial for fig 7...

figure(13) 
hold off
title('Charge density')
xlabel('x')
ylabel('\lambda (x)')
plot([xmin,xmax],[1,1],'k-')
hold on
% Now ready to accept an overplot with each trial for fig 13...

linestyles = ['k-*';'r-^';'g-o';'m-x';'b-+';'k--';'r--';'g--';'m--';'b--'];
   % ...for the charge density plot.

for i = 1:Ntrials
   dx = 0.5 * 2.0^(1-i)
      % It is very important that xmax-xmin/dx be an integer. If not,
      % the first and last charges would not begin at the ends of the
      % rod. Then the energy relaxation ratio calculated later would 
      % be meaningless.
   x0 = (xmin:dx:xmax)'; % initial charge position array. 
   % Note that sqp uses column vectors when calling @constraints.

   N(i) = numel(x0);  % Number of charges
   Ncharges = N(i)   % Print this to the terminal.
   
   Einitial(i) = energy(x0);
   t0 = cputime();
   [xfinal, Efinal(i), info(i), iter(i)] = sqp(x0, @energy, [], @constraints1,...
       xmin, xmax, maxiter, tol);
   cpu_sec(i) = cputime() - t0;
   elapsed_cputime = sum(cpu_sec)
   return_condition = info(i)
   
   if (return_condition != 101) % 101 signifies normal convergence of SQP.
      break;
   end
   
   % FIGURES UPDATE WITH EACH TIRAL SO THAT IF SQP CRASHES, WE'LL HAVE THE 
   % BEST AVAILABLE PLOTS...

   figure(7)
   hold on
   plot(x0, xfinal-x0, linestyles(i,:))
   print('convergence07_chg_disp.pdf','-dpdf') 
   
   max_charge_displacement(i) = max(abs(xfinal-x0));
   figure(8)
   hold off
   semilogx(N, max_charge_displacement, 'b-o')
   hold on
   title('Maximum absolute charge displacement')
   ylabel('Distance')
   xlabel('Number of charges')
   print('convergence08_max_chg_disp.pdf','-dpdf') 
   hold off

   med_charge_displacement(i) = median(abs(xfinal-x0));
   figure(9)
   hold off
   semilogx(N, med_charge_displacement, 'b-o')
   hold on
   title('Median absolute charge displacement')
   ylabel('Distance')
   xlabel('Number of charges')
   print('convergence09_med_chg_disp.pdf','-dpdf') 
   hold off

   lambda = charge_density(xfinal); % estimate of continuous charge density
   figure(13)
   hold on
   plot(xfinal, lambda, linestyles(i,:))
   print('convergence13_chg_density.pdf','-dpdf')
   
   fraction_enhanced(i) = sum(lambda > 1)/N(i);
   figure(14)
   hold off
   semilogx(N, fraction_enhanced, 'b-o')
   hold on
   title('Fraction of charge in density-enhanced regions (\lambda > 1)')
   xlabel('Number of charges')
   print('convergence14_fraction_enhanced.pdf','-dpdf') 
   hold off

   
   max_lambda(i) = max(lambda);
   figure(15)
   hold off
   semilogx(N, max_lambda, 'b-o')
   hold on
   title('Maximum charge density')
   xlabel('Number of charges')
   print('convergence15_max_chg_dens.pdf','-dpdf')
   hold off

   min_lambda(i) = min(lambda);
   figure(16)
   hold off
   semilogx(N, min_lambda, 'b-o')
   hold on
   title('Minimum charge density')
   xlabel('Number of charges')
   print('convergence16_min_chg_dens.pdf','-dpdf')
   hold off
   
   median_abs_dev(i) = median(abs(lambda-1));
   figure(17)
   hold off
   semilogx(N, median_abs_dev, 'b-o')
   hold on
   title('Median of | \lambda - 1 |')
   xlabel('Number of charges')
   print('convergence17_med_abs_dev.pdf','-dpdf') 
   hold off


   
   figure(1)
   hold off
   loglog(N, iter,'b-o')
   hold on
   xlabel('Number of charges')
   ylabel('Number of iterations')
   title('Computational Effort (iterations)')
   print('convergence01_iterations.pdf','-dpdf')
   hold off

   figure(2)
   hold off
   loglog(N, cpu_sec,'b-o')
   hold on
   xlabel('Number of charges')
   ylabel('CPU time (s)')
   title('Computational Effort (CPU time)')
   print('convergence02_time.pdf','-dpdf')
   hold off

   figure(3)
   hold off
   semilogx(N, Einitial,'b-o')
   hold on
   semilogx(N, Efinal,'r-^')
   xlabel('Number of charges')
   ylabel('Energy (arbitrary units)')
   legend('Initial Energy','Final Energy')
   title('Energy Convergence')
   print('convergence03_energy.pdf','-dpdf')
   hold off

   figure(4)
   hold off
   semilogx(N, Efinal ./ Einitial,'b-o')
   hold on
   xlabel('Number of charges')
   ylabel('Final Energy / Initial Energy')
   title('Energy Relaxation Ratio Convergence')
   print('convergence04_eratio.pdf','-dpdf')
   hold off
   
end



