%
%------------------------------------------------------------------
% FTE_FFT_premodecleaner
%  
% Matlab example to demonstrate the use of several FT_FFT_xxx.m
% functions for an FFT based simulation of a 3-mirror
% triangular cavity used with a Laguerre Guass input beam.
% Based on an OSCAR example by Jerome Degallaix.
%
% Part of the SimTools package
% Andreas Freise 05.01.2010 afreise@googlemail.com 
%-------------------------------------------------------------------
%

clear all;

verbose=0;

% TODO improve saving of txt data files


% Parameters of the input laser beam, q value is taken from Finesse file
Laser=FT_init_gauss_param(1064e-9,1,'w0',362.1e-6,'z',-970.4e-3);
q_start=Laser.q(1);


% Init Grid for FFT 
Grid=FT_init_grid(256,256,0.007,0.007,0,0);

% Init three mirror cavity
s1=FT_init_space('L',0.02);
s2=FT_init_space('L',0.2);
s3=FT_init_space('L',0.2);

alpha1=2.86;
alpha2=0.5*(90-alpha1);
%alpha1=0.0;
%alpha2=0.0;

m_in  = FT_init_mirror('name','m_in' ,'T',500e-6,'L',30e-6, 'diameter',0.0254, 'alpha',alpha2);
m_out = FT_init_mirror('name','m_out','T',500e-6,'L',30e-6, 'diameter',0.0254, 'alpha',alpha2);
m_end = FT_init_mirror('name','m_end','T',100e-6,'L',0.0,   'diameter',0.0254, 'alpha',alpha1, 'Rc', 1);

m_in  = FT_add_mirror_map(m_in,  Grid, Laser.lambda);
m_out = FT_add_mirror_map(m_out, Grid, Laser.lambda);
m_end = FT_add_mirror_map(m_end, Grid, Laser.lambda);

cavity=FT_init_cavity(s1,m_out,s2,m_end,s2,m_in);

% Create input field

field = complex(zeros(Grid.xpoints,Grid.ypoints,'double'));

field = FT_LG_cos_field(Laser,3,3,Grid.xaxis,Grid.yaxis,[0,0,0]);
%field = FT_LG_cos_field(Laser,0,0,Grid.xaxis,Grid.yaxis,[0,0,0]);
field = field/sqrt(FT_power_in_field(field,Grid.xaxis,Grid.yaxis));

% Propagate input field through mode-matching telescope
field = FT_FFT_propagate(field,Grid, Laser.lambda, 0.25, 1);
field = FT_FFT_lens(field,Grid, Laser.lambda, -2);
field = FT_FFT_propagate(field,Grid, Laser.lambda, 0.08, 1);
field = FT_FFT_lens(field,Grid, Laser.lambda, 2);
field = FT_FFT_propagate(field,Grid, Laser.lambda, 0.68, 1);

% propagate 'q' through the same system
Mlens1 = FT_ABCD_lens(-2);
Mlens2 = FT_ABCD_lens(2);
q = q_start+0.25;
q = FT_q_ABCD_q(Mlens1, q, 1, 1);
q = q + 0.08;
q = FT_q_ABCD_q(Mlens2, q, 1, 1);
q = q + 0.68;
% transmit `q' to m_out
q = q + 0.02;


% Pre-compute all roundtrips fields
[req_roundtrips]=FT_FFT_required_cavity_roundtrips(cavity,1e-3);
%roundtrips=100;
roundtrips=req_roundtrips;


dfile = FT_FFT_precompute_roundtrips(cavity,field,Grid,Laser.lambda,roundtrips,0);
dfile = FT_FFT_reorganise_tmp_files(dfile, Grid);

filename=sprintf('precompute_%s.mat',dfile.datestr);
save(filename,'Laser','Grid','dfile','cavity','field','q');

% Do a full cavity scan, using precomputed data from above
fprintf(' ==== Full cavity scan ====\n');

[req_points]=FT_FFT_required_scan_accuracy(cavity);
points=1000;
set_roundtrips=round(0.1*roundtrips);
fprintf(' Estimated requirement: %d roundtrips, %d points for full scan.\n',req_roundtrips,req_points);
fprintf(' Using %d roundtrips, %d points.\n',set_roundtrips, points);
[result.scan,scan_step] = FT_FFT_scan_range(Laser.lambda/2, Laser.lambda, points);
result.power = FT_FFT_cavity_scan(dfile,Grid,Laser.lambda,set_roundtrips,result.scan);
save_xy=[result.scan' result.power'];
save('scan1.txt','-ASCII','save_xy');

if (verbose)
figure
  semilogy(result.scan/Laser.lambda,result.power);
  title('Cavity Scan')
  xlabel('tuning [1/lambda]')
  ylabel('circulating power')
  drawnow;
end
  
[result.max_power,result.max_scan_idx] = max(result.power);
result.operating_point= result.scan(result.max_scan_idx);
fprintf(' -> operating point %f [1/lambda] at max power %f [W] \n',result.operating_point/Laser.lambda, result.max_power); 

% Zooming in 1/2 ...
fprintf(' ==== Zooming in on maximum peak, range = 25 FWHM ====\n');
points=400;
set_roundtrips=round(0.25*roundtrips);
set_range=Laser.lambda/cavity.F*25;
fprintf(' Scanning around %f [1/lambda], range = %g [1/lambda], roundtrips = %d, points = %d.\n', ...
        result.operating_point/Laser.lambda,set_range/Laser.lambda, set_roundtrips, points);
[result.zoom1, zoom_step] = FT_FFT_scan_range(result.operating_point, set_range , points);
result.power_zoom1=FT_FFT_cavity_scan(dfile,Grid,Laser.lambda,set_roundtrips,result.zoom1);
save_xy=[result.zoom1' result.power_zoom1'];
save('zoom1.txt','-ASCII','save_xy');


if (verbose)
figure
  plot(result.zoom1/Laser.lambda,result.power_zoom1);
  title('Cavity Zoom')
  xlabel('tuning [1/lambda]')
  ylabel('circulating power')
  drawnow;
end
  
[result.max_power,result.max_zoom_idx] = max(result.power_zoom1);
result.operating_point= result.zoom1(result.max_zoom_idx);
fprintf(' -> operating point %f [1/lambda] at max power %f [W] \n',result.operating_point/Laser.lambda, result.max_power); 

% Zooming in 2/2 ...
fprintf(' ==== Zooming in on maximum peak, getting down to 1pm accuracy ====\n');
points=round(2*zoom_step/0.5e-12);
set_roundtrips=dfile.roundtrips;
set_range=2*zoom_step;
fprintf(' Scanning around %f [1/lambda], range = %g [1/lambda], roundtrips = %d, points = %d.\n', ...
        result.operating_point/Laser.lambda,set_range/Laser.lambda, set_roundtrips, points);
[result.zoom2, zoom_step] = FT_FFT_scan_range(result.operating_point, set_range, points);
result.power_zoom2 = FT_FFT_cavity_scan(dfile,Grid,Laser.lambda,set_roundtrips,result.zoom2);
save_xy=[result.zoom2' result.power_zoom2'];
save('zoom2.txt','-ASCII','save_xy');

if (verbose)
figure
 plot(result.zoom2/Laser.lambda,result.power_zoom2);
 title('Cavity Zoom')
 xlabel('tuning [1/lambda]')
 ylabel('circulating power')
 drawnow;
end
 
[result.max_power,result.max_zoom_idx] = max(result.power_zoom2);
result.operating_point= result.zoom2(result.max_zoom_idx);
fprintf(' -> operating point %f [1/lambda] at max power %f [W] \n',result.operating_point/Laser.lambda, result.max_power); 

% Computing buildup on operating point

[result] = FT_FFT_results_buildup(cavity, field, Grid, dfile, Laser.lambda, result, 0);
if (verbose)
figure
plot(result.power_buildup);
title('Power Buildup')
xlabel('iteration')
ylabel('circulating power')
drawnow;
end

result.f_in=field;

% Compute remaining fields (trans, refl) of the cavity
result=FT_FFT_results_cavity_fields(result, cavity, Grid, Laser, 1);

filename=sprintf('results_%s.mat',dfile.datestr);
save(filename,'Laser','Grid','dfile','cavity','result');

fprintf('\n === Results \n')
fprintf(' Operating point [1/lambda]: %f \n',result.operating_point/Laser.lambda); 
fprintf(' Operating point [m]:        %g \n',result.operating_point); 
fprintf(' \n');
fprintf(' Circulating power [W]: %f \n', FT_power_in_field(result.f_circ,Grid.xaxis,Grid.yaxis));
fprintf(' Reflected power   [W]: %f \n', FT_power_in_field(result.f_refl,Grid.xaxis,Grid.yaxis));
fprintf(' Transmitted power [W]: %f \n', FT_power_in_field(result.f_trans,Grid.xaxis,Grid.yaxis));

% Do mode decomposition
[gp] = FT_init_gauss_param(Laser.lambda,1,'q',q);
max_mode_order = 9;
[gc,n_coeff]=FT_init_gauss_coefficients('HG', max_mode_order);
fprintf(' Computing mode coefficients for transmitted field...\n');
[gc] = FT_mode_content(gc, gp, result.f_trans,Grid.xaxis,Grid.yaxis);
FT_print_gauss_coefficients(gc,1e-3,0);
fprintf(' Total power in modes %.3g (should be 1)!\n',sum(abs(gc.value(:)).^2));

if (verbose)
figure
clf;    
subplot(2,2,1)
FT_FFT_plot_field(field,1)
title('Input field')
subplot(2,2,2)
FT_FFT_plot_field(result.f_circ,1)
title('Circulating field')
subplot(2,2,3)
FT_FFT_plot_field(result.f_refl,1)
title('Reflected field')
subplot(2,2,4)
FT_FFT_plot_field(result.f_trans,1)
title('Transmitted field')
drawnow;  
end

