From 99f39dc1423adc360a04344689a83d2f788aae33 Mon Sep 17 00:00:00 2001
From: Eugeniy Mikhailov <evgmik@gmail.com>
Date: Wed, 20 Mar 2013 16:53:39 -0400
Subject: Matt's changes to use faster beam propagator

---
 fitness.m                        | 60 ++++++++++++++++++++++++-------------
 fitter_check.m                   | 32 ++++++++------------
 gbeam_propagation.m              |  8 ++---
 gbeam_propagation_froward_only.m |  9 ++----
 mode_match.m                     | 26 +++++-----------
 pick_visualization.m             | 10 +++----
 self_gbeam_propagation.m         | 49 +++++++++++++++---------------
 solution_visualization.m         | 64 +++++-----------------------------------
 8 files changed, 101 insertions(+), 157 deletions(-)

diff --git a/fitness.m b/fitness.m
index 55f9d76..e92879f 100644
--- a/fitness.m
+++ b/fitness.m
@@ -3,25 +3,26 @@ function [Energy, Waist, Penalty] = fitness( q_0, q_final, x_final, optics_posit
     x0 = 0;
     Np=20; % # of pts between start position and second lens
     N_collimated = 10; % # of pts between second lens and third lens to be collimated region
+    Energy = 0;
     
     x1=optics_positions(1);
     x2=optics_positions(2);
     x3=optics_positions(3);
     
-    x_array1=linspace(x0,x2,Np);
-    x_array2=linspace(x2, x3, N_collimated);
-    x = cat(2,x_array1,x_array2);
-    
-    q_f_trial_forward = gbeam_propagation(x,q_0,x0,optics_placer(optics_positions, optics_focal_length));
-    [Waist_trial_forward, Radius_trial_forward] = q2wr(q_f_trial_forward, lambda);
-    q_f_trial_backward = gbeam_propagation(x,q_final,x_final,optics_placer(optics_positions, optics_focal_length));
-    [Waist_trial_backward, Radius_trial_backward] = q2wr(q_f_trial_backward, lambda);
-
-    Energy = 0;
-    Penalty_waist_mismatch = sum(abs((Waist_trial_forward-Waist_trial_backward)./min(Waist_trial_forward, Waist_trial_backward)))/Np;
-    
-    Energy = 1e-2*Penalty_waist_mismatch;
-    
+%     x_array1=linspace(x0,x2,Np);
+%     x_array2=linspace(x2, x3, N_collimated);
+%     x = cat(2,x_array1,x_array2);
+%     
+%     q_f_trial_forward = gbeam_propagation(x,q_0,x0,optics_placer(optics_positions, optics_focal_length));
+%     [Waist_trial_forward, Radius_trial_forward] = q2wr(q_f_trial_forward, lambda);
+%     q_f_trial_backward = gbeam_propagation(x,q_final,x_final,optics_placer(optics_positions, optics_focal_length));
+%     [Waist_trial_backward, Radius_trial_backward] = q2wr(q_f_trial_backward, lambda);
+% 
+%     Energy = 0;
+%     Penalty_waist_mismatch = sum(abs((Waist_trial_forward-Waist_trial_backward)./min(Waist_trial_forward, Waist_trial_backward)))/Np;
+%     
+%     Energy = 1e-2*Penalty_waist_mismatch;
+%     
     % penalty calculation
     % do not put lenses too close to each other and end positions
     lens_size=0.03;
@@ -58,14 +59,33 @@ function [Energy, Waist, Penalty] = fitness( q_0, q_final, x_final, optics_posit
     Energy = Energy + penalty_lenses_outside_optical_path;
     
     % make collimated region between 2nd and 3rd lens   
-    %intialize intermediate points between lenses
-    q_intermediate = q_f_trial_forward((x2<x) & (x<x3));
-    lambda_over_waist_sq = (-imag(1./q_intermediate)); %with numerical factor
     
-    coef = 1e-3;
-    penalty_not_collimated_beam = coef * sum(exp((std(lambda_over_waist_sq)/mean(lambda_over_waist_sq)/abs(optics_positions(2) - optics_positions(3)))));
+    [w0, r0] = q2wr(q_0, lambda);
+    [ w, w_pos ] = self_gbeam_propagation( w0, optics_positions, optics_focal_length, x0, lambda );
+    
+    coef = 1;
+    d_object = ((optics_positions(end) - w_pos(end - 1))^2)^.5;
+    d_lens = optics_positions(end) - optics_positions(end - 1);
+    penalty_not_collimated_beam = coef * exp(-(d_object/d_lens)^2);
+    
     Energy = Energy + penalty_not_collimated_beam;
     
-    Penalty = [ penalty_lenses_too_closely_spaced; penalty_lenses_too_closely_spaced; penalty_not_collimated_beam];
+%     % waist at end matches desired waist
+    coef = 10;
+    waist_desired = 3.709E-5;
+    
+    penalty_waist = coef *((waist_desired-w(end))^2)^.5;
+    
+    Energy = Energy + penalty_waist;
+
+    %intialize intermediate points between lenses
+%     q_intermediate = q_f_trial_forward((x2<x) & (x<x3));
+%     lambda_over_waist_sq = (-imag(1./q_intermediate)); %with numerical factor
+%     
+%     coef = 1e-3;
+%     penalty_not_collimated_beam = coef * sum(exp((std(lambda_over_waist_sq)/mean(lambda_over_waist_sq)/abs(optics_positions(2) - optics_positions(3)))));
+%     Energy = Energy + penalty_not_collimated_beam;
+%     
+   Penalty = [ penalty_lenses_too_closely_spaced; penalty_lenses_outside_optical_path; penalty_not_collimated_beam; penalty_waist];
 end
 
diff --git a/fitter_check.m b/fitter_check.m
index f70ded6..c4c5fbf 100644
--- a/fitter_check.m
+++ b/fitter_check.m
@@ -1,5 +1,6 @@
 %Permute all possible lens combinations out of set of lenses
-lens_set = [.075, .203];
+lens_set = [.075, .203, .05, .03];
+lens_set = [.075,.203];
 lens_permutations = pick(lens_set,3,'or');
 
 %Pre-defined Constants
@@ -13,32 +14,23 @@ rf= 1.0E+100 ;
 xf= Ltot;
 q0=wr2q(w0,r0,lambda);
 qf=wr2q(wf,rf,lambda);
-lens_width = .03;
-show_lens_width = 1;
-show_lens_position = 1;
-display_prop = [show_lens_width, show_lens_position];
 %End list
 
 %Mode match
-[ possible_lens_placement, initial_lens_placement, possible_lens_set, possible_sample_energy] = mode_match( q0, qf, Ltot, lambda, lens_permutations, lens_width );
+[ possible_lens_placement, possible_lens_set, possible_sample_energy] = mode_match( q0, qf, Ltot, lambda, lens_permutations );
 
 %Remove similar solutions
 n_truncate = 3; %number of digits in truncated solution
 [ possible_lens_placement_uniq, possible_lens_placement, possible_sample_energy, possible_lens_set, index ] = remove_similar_soln( possible_sample_energy, possible_lens_placement, possible_lens_set, n_truncate );
 
-n_visualizations = 5; %number of best solutions to visualize
-n_hist = 10000; %number of sample points in histogram
-stability_max = 1; %max stability (y-axis) shown on energy vs. stability graph
+n_visualizations = 2; %number of best solutions to visualize
+pick_visualization( possible_sample_energy, possible_lens_placement_uniq, possible_lens_placement, possible_lens_set, index, n_visualizations, q0, qf, Ltot, lambda );
 
-pick_visualization( possible_sample_energy, possible_lens_placement_uniq, possible_lens_placement, possible_lens_set, index, n_visualizations, q0, qf, Ltot, lambda, lens_width, display_prop );
-[stability] = stability_visualization( possible_lens_placement_uniq, q0, qf, xf, possible_lens_placement, possible_lens_set, lambda, n_visualizations, n_hist, index );
-energy_vs_stability( possible_sample_energy, stability, index, stability_max)
+%Visualize fitness function for fixed f2 and f3
+lens_set = [.075, .075, .203];
+f2= 0.40361319425309 ;
+f3= 0.80361319425309 ;
 
-% %Visualize fitness function for fixed f2 and f3
-% lens_set = [.075, .075, .203];
-% f2= 0.40361319425309 ;
-% f3= 0.80361319425309 ;
-% 
-% fitness_simplified=@(x) fitness(q0, qf, Ltot, [x, f2, f3], lens_set, lambda );
-% figure(6)
-% ezplot(fitness_simplified, [0,Ltot])
\ No newline at end of file
+fitness_simplified=@(x) fitness(q0, qf, Ltot, [x, f2, f3], lens_set, lambda );
+figure(6)
+ezplot(fitness_simplified, [0,Ltot])
\ No newline at end of file
diff --git a/gbeam_propagation.m b/gbeam_propagation.m
index 585048d..d81f70e 100644
--- a/gbeam_propagation.m
+++ b/gbeam_propagation.m
@@ -1,4 +1,4 @@
-function [q] = gbeam_propagation(x_pos, q_in, x_in,  optics_elements)
+function q = gbeam_propagation(x_pos, q_in, x_in,  optics_elements)
 % calculate the 'q' parameter of the Gaussian beam propagating through optical
 % 'optics_elements' array along 'x' axis at points 'x_pos'
 %  takes the gaussian beam with initial q_in parameter at x_in
@@ -8,7 +8,7 @@ function [q] = gbeam_propagation(x_pos, q_in, x_in,  optics_elements)
 
 	if any(x_pos >= x_in)
 		% Forward propagation to the right of x_in
-		[q(x_pos >= x_in)] = gbeam_propagation_froward_only(x_pos(x_pos>=x_in), q_in, x_in,  optics_elements);
+		q(x_pos >= x_in) = gbeam_propagation_froward_only(x_pos(x_pos>=x_in), q_in, x_in,  optics_elements);
 	end
 
 	if any(x_pos < x_in)
@@ -45,10 +45,8 @@ function [q] = gbeam_propagation(x_pos, q_in, x_in,  optics_elements)
 		% final assignment of the backwards propagating beam 
 		% which we need to flip back
 		q(x_pos<x_in) = fliplr(q_backw);
-    end
+	end
 
-    
-    
 end
 
 %!test 
diff --git a/gbeam_propagation_froward_only.m b/gbeam_propagation_froward_only.m
index 2cdc1d8..afc1316 100644
--- a/gbeam_propagation_froward_only.m
+++ b/gbeam_propagation_froward_only.m
@@ -1,4 +1,4 @@
-function [q] = gbeam_propagation_froward_only(x_pos, q_in, x_in,  optics_elements)
+function q = gbeam_propagation_froward_only(x_pos, q_in, x_in,  optics_elements)
 % calculate the 'q' parameter of the Gaussian beam propagating through optical
 % 'optics_elements' only in the positive direction  along 'x' axis at points 'x_pos'
 % takes the gaussian beam with initial q_in parameter at x_in
@@ -8,10 +8,6 @@ function [q] = gbeam_propagation_froward_only(x_pos, q_in, x_in,  optics_element
 %
 % all x_pos must be to the right of x_in
 % x_pos must be monotonic!
-    
-    %Initialize q_lens (q at position of lens)
-    q_lens = zeros(1,size(optics_elements,2));
-
 	if (any(x_pos < x_in))
 		error('all beam positions must be to the right of the x_in');
 	end
@@ -49,7 +45,8 @@ function [q] = gbeam_propagation_froward_only(x_pos, q_in, x_in,  optics_element
     index = (x_last_calc <= x_pos);
     x=x_pos(index);
     q(index)= q_after_free_space(q_last_calc,x-x_last_calc);
-      
+       
+	
 	
 end
 
diff --git a/mode_match.m b/mode_match.m
index dafe209..0a36332 100644
--- a/mode_match.m
+++ b/mode_match.m
@@ -1,4 +1,4 @@
-function [ final_possible_lens_placement, initial_possible_lens_placement, possible_lens_set, possible_sample_energy] = mode_match( q0, qf, Ltot, lambda, lens_permutations, lens_width )
+function [ possible_lens_placement, possible_lens_set, possible_sample_energy] = mode_match( q0, qf, Ltot, lambda, lens_permutations )
 %Shuffles lenses into random positions and stores possible solutions
 %   Shuffles over entire lens permutation array for n_shuffles times.
 %   Afterwards, solutions are sorted by Energy and truncated to n_truncate
@@ -6,10 +6,6 @@ function [ final_possible_lens_placement, initial_possible_lens_placement, possi
 
 n_perms = size(lens_permutations,1);
 n_shuffles=20; %number of random placements of lenses
-initial_MaxFunEvals = 1e8;
-initial_MaxIter = 10;
-final_MaxFunEvals = 1e8;
-final_MaxIter = 100;
 
 %Initialize sample arrays
 N = n_perms * n_shuffles;
@@ -18,33 +14,25 @@ possible_lens_set = zeros(N,3);
 possible_sample_energy = zeros(N,1);
 initial_rand_lens_placement=zeros(N,3);
 
+lens_size = .03; % physical size of the lens
+
 for ip = 1:n_perms
     f3=lens_permutations(ip,3);
     x3=Ltot-f3; % last lense transfer collimated region to focused spot
     for is = 1:n_shuffles
         possible_lens_set((ip-1)*n_shuffles + is,:) = lens_permutations(ip,:);
 
-        initial_rand_lens_placement_tmp = sort(lens_width+(x3-2*lens_width)*rand(1,2));
+        initial_rand_lens_placement_tmp = sort(lens_size+(x3-2*lens_size)*rand(1,2));
         initial_rand_lens_placement((ip-1)*n_shuffles + is,:) = [initial_rand_lens_placement_tmp, x3];
     end
 end
 
-
-
-
 parfor i = 1:N
+      
         fitness_simplified=@(x) fitness(q0, qf, Ltot, x, possible_lens_set(i,:), lambda );
-        [ x_sol, energy ] = find_min(i, initial_rand_lens_placement, fitness_simplified, initial_MaxFunEvals, initial_MaxIter )
-                
-        initial_possible_lens_placement(i,:) =x_sol;
+        [x_sol, energy]=fminsearch(fitness_simplified, initial_rand_lens_placement(i,:), optimset('TolFun',1e-5,'MaxFunEvals',1e8,'MaxIter',200));
         
-end
-
-parfor i = 1:N
-        fitness_simplified=@(x) fitness(q0, qf, Ltot, x, possible_lens_set(i,:), lambda );
-        [ x_sol, energy ] = find_min(i, initial_possible_lens_placement, fitness_simplified, final_MaxFunEvals, final_MaxIter )
-                
-        final_possible_lens_placement(i,:) =x_sol;
+        possible_lens_placement(i,:) =x_sol;
         possible_sample_energy(i) = energy;
 
 end
diff --git a/pick_visualization.m b/pick_visualization.m
index 26bd2a6..602e509 100644
--- a/pick_visualization.m
+++ b/pick_visualization.m
@@ -1,4 +1,4 @@
-function [ ] = pick_visualization( fitness_energy, possible_lens_placement_uniq, possible_lens_placement, possible_lens_set, index, n_visualizations, q0, qf, Ltot, lambda, lens_width, display_prop )
+function [ ] = pick_visualization( fitness_energy, possible_lens_placement_uniq, possible_lens_placement, possible_lens_set, index, n_visualizations, q0, qf, Ltot, lambda )
 %Picks n_visualizations of sets of data and graphs each 
 x0 = 0;
 
@@ -6,11 +6,11 @@ n_possible_lens_placement = min(n_visualizations,size(possible_lens_placement_un
 
 for n_graph = 1:n_possible_lens_placement
     figure(n_graph)
-    [w_final_trial, r_final_trial] = solution_visualization(q0, x0, qf, Ltot, optics_placer(possible_lens_placement(index(n_graph),:), possible_lens_set(index(n_graph),:)), lambda, lens_width, display_prop);
+    [w_final_trial, r_final_trial] = solution_visualization(q0, x0, qf, Ltot, optics_placer(possible_lens_placement(index(n_graph),:), possible_lens_set(index(n_graph),:)), lambda);
     
-    str1=sprintf('\n (red) f_1 = %0.4f, x_1 = %0.4f\n',possible_lens_set(index(n_graph),1),possible_lens_placement(index(n_graph),1));
-    str2=sprintf(' (green) f_2 = %0.4f, x_2 = %0.4f\n',possible_lens_set(index(n_graph),2),possible_lens_placement(index(n_graph),2));
-    str3=sprintf(' (blue) f_3 = %0.4f, x_3 = %0.4f\n',possible_lens_set(index(n_graph),3),possible_lens_placement(index(n_graph),3));
+    str1=sprintf('\n f_1 = %0.4f, x_1 = %0.4f\n',possible_lens_set(index(n_graph),1),possible_lens_placement(index(n_graph),1));
+    str2=sprintf('f_2 = %0.4f, x_2 = %0.4f\n',possible_lens_set(index(n_graph),2),possible_lens_placement(index(n_graph),2));
+    str3=sprintf('f_3 = %0.4f, x_3 = %0.4f\n',possible_lens_set(index(n_graph),3),possible_lens_placement(index(n_graph),3));
     tstr='Solution #';
     str_w = ['w_{final}= ',num2str(w_final_trial)];
     str_r = [', r_{final}= ',num2str(r_final_trial)];
diff --git a/self_gbeam_propagation.m b/self_gbeam_propagation.m
index 807303c..10bcf05 100644
--- a/self_gbeam_propagation.m
+++ b/self_gbeam_propagation.m
@@ -1,4 +1,4 @@
-function [ output_args ] = self_gbeam_propagation( w0, x_lens, f, x0, lambda )
+function [ w, w_pos ] = self_gbeam_propagation( w0, x_lens, f, x0, lambda )
 %SELF_GBEAM_PROPAGATION Summary of this function goes here
 %   Detailed explanation goes here
 
@@ -14,31 +14,28 @@ for i = 1:n_lens
     w_pos(i+1) = x_lens(i) + s_new;
 end
 
-x_pos= x_lens;
-x(1)=w_pos(1);
-waist_prop(1)=w(1);
-
-x_temp=linspace(x0 ,x_pos(1));
-waist_prop=[waist_prop w(1)*sqrt(1+((x_temp-w_pos(1))/(pi*w(1)^2/lambda)).^2)];
-x=[x x_temp];
-
-for i=1:n_lens
-    if i == 3
-        break
-    end
-    x_temp=linspace(x_pos(i),x_pos(i+1));     
-    waist_prop=[waist_prop w(i)*sqrt(1+((x_temp-w_pos(i))/(pi*w(i)^2/lambda)).^2)];  
-    x=[x x_temp];
-end
-
-x_temp=linspace(x_pos(end),x_pos(end)+abs(f(end)));  
-waist_prop=[waist_prop w(end)*sqrt(1+((x_temp-w_pos(end))/(pi*w(end)^2/lambda)).^2)];
-x=[x x_temp];
-
-
-figure(1)
-plot(x,waist_prop,'r',x,-waist_prop,'r'); 
-
+% x_pos= x_lens;
+% x(1)=w_pos(1);
+% waist_prop(1)=w(1);
+% 
+% x_temp=linspace(x0 ,x_pos(1));
+% waist_prop=[waist_prop, w(1)*sqrt(1+((x_temp-w_pos(1))/(pi*w(1)^2/lambda)).^2)];
+% x=[x x_temp];
+% 
+% for i=1:n_lens
+%     x_temp=linspace(x_pos(i),x_pos(i+1));     
+%     waist_prop=[waist_prop, w(i)*sqrt(1+((x_temp-w_pos(i))/(pi*w(i)^2/lambda)).^2)];  
+%     x=[x x_temp];
+% end
+% 
+% x_temp=linspace(x_pos(end),x_pos(end)+abs(f(end)));  
+% waist_prop=[waist_prop, w(end)*sqrt(1+((x_temp-w_pos(end))/(pi*w(end)^2/lambda)).^2)];
+% x=[x x_temp];
+% 
+% 
+% figure(1)
+% plot(x,waist_prop,'r',x,-waist_prop,'r'); 
+% 
 
 
 
diff --git a/solution_visualization.m b/solution_visualization.m
index 2b4158c..2cf8ce7 100644
--- a/solution_visualization.m
+++ b/solution_visualization.m
@@ -1,19 +1,11 @@
-function [waste_at_the_end, radius_at_the_end, waist_at_lens_position] = solution_visualization(q0,x0, qf, xf, optics, lambda, lens_width, display_prop)
+function [waste_at_the_end, radius_at_the_end] = solution_visualization(q0,x0, qf, xf, optics, lambda)
 %Propagates beam with given input parameters
 %   Forward propagation and backward propagation are taken in order to
 %   visualize reasonable solutions.
 
-%Array of lens positions
-n_lens = size(optics,2);
-
-for i = 1:n_lens
-    lens_position(i) = optics{i}.x;
-end
-
 x=linspace(x0,xf,1000); % we will calculate beam profile between x0 and xf
-
 %Forward propagation
-q_forward =gbeam_propagation(x,q0,x0,optics);
+q_forward=gbeam_propagation(x,q0,x0,optics);
 [w_forward,r_forward]=q2wr(q_forward, lambda);
 
 %Backward propagation
@@ -21,53 +13,13 @@ q_backward=gbeam_propagation(x,qf,xf,optics);
 [w_backward,r_backward]=q2wr(q_backward, lambda);
 
 %Plot beam profile
-subplot(2,1,1); plot ( ...
-    x,w_forward, '-r',  ...
-    x,-w_forward, '-r', ...
-    x, w_backward, '-.b', ...
-    x, -w_backward, '-.b')
-legend({'forward propagation', '', 'backward propagation', ''})
-
-%Find q and waist at lens positions
-waist_at_lens_position = zeros(1,n_lens);
-q_lens = zeros(1,n_lens);
-
-for i = 1:n_lens
-    q_lens(i) = gbeam_propagation(lens_position(i),q0,x0,optics);
-end
-
-
-for i = 1:n_lens
-    waist_at_lens_position(i) = q2wr(q_lens(i), lambda);
-end
-
-
-%Plot lenses
-color = ['m' 'g' 'c'];
-
-if display_prop(1) == 1
-    for i = 1:n_lens
-        half_width = lens_width/2;
-        corrected_lens_position = lens_position(i)-half_width;
-        corrected_waist = waist_at_lens_position(i)*2;
-        rectangle('Position', [corrected_lens_position,-waist_at_lens_position(i),lens_width,corrected_waist], 'EdgeColor', color(i), 'LineWidth',1);
-    end
-end
-
-if display_prop(2) == 1
-    for i = 1:n_lens
-        x1 = optics{i}.x;
-        y1 = waist_at_lens_position(i);
-        x2 = x1;
-        y2 = -waist_at_lens_position(i);
-        
-        line([x1;x2], [y1;y2], 'LineWidth', 3, 'Color', color(i));
-    end
-end
-
+plot ( ...
+	x,w_forward, '-r',  ...
+	x,-w_forward, '-r', ...
+	x, w_backward, '-.b', ...
+	x, -w_backward, '-.b')
+legend({'forward propagation', '', 'backward propagation', ''})      
 
 [waste_at_the_end,radius_at_the_end]  = q2wr(q_forward(end), lambda);
 
 
-
-
-- 
cgit v1.2.3