clear all
close all
clear classes
addpath alm-master

%NPO_a_21_Nov2024.txt

z_offset = 0.383 + -0.102;
zScan_NPRO = [0:10:100]*1e-3 + z_offset; %this is equivalent to [0,10,20,.....,100];
xWidth_NPRO = [1628.7, 1663.4, 1701.8, 1739.2, 1777.3, 1822.1, 1855.3, 1904.1, 1926.3, 1958.7, 2008.8] * 1e-6/2;
yWidth_NPRO = [1739.5, 1775.0, 1805.1, 1848.5, 1878.2, 1924.4, 1967.9, 2007.8, 2024.7, 2064.9, 2097.5] * 1e-6/2;

%%%% beam profile measurements for NPRO 1661  (found in
%%%% labutils/beam_scans/Data/NPRO_toAMP1_e_27Oct2024.txt) 
% I think this text file (NPRO_26Oct2024.txt) is the profile from the
% leakage of M01
%to make this fit match the waist location from alog 80895 I need to add an
%offset, probably the numbers in text file are rail tick marks
% z_offset = 0.383 + -0.102;
% zScan_NPRO = [0:10:100]*1e-3 + z_offset; %this is equivalent to [0,10,20,.....,100];
% xWidth_NPRO = [1745.9, 1784.4, 1819.8, 1873.0, 1921.9, 1941.0, 2014.3, 2044.9, 2084.5, 2135.2, 2178.2] * 1e-6/2;
% yWidth_NPRO = [1805.2, 1847.2, 1882.4, 1945.7, 1966.7, 2031.8, 2084.5, 2118.2, 2160.6, 2230.2, 2267.7] * 1e-6/2;

% reading alog 80895 this seems to be measurements taken after the
% telescope
zScan_after = [400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500]*1e-3 ; 
xWidth_after = [518.6, 541.4, 576.9, 620.7, 646.8, 686.7, 726.0, 765.2, 804.9, 850.7, 885.7] * 1e-6/2;
yWidth_after = [492.6, 519.5, 548.6, 581.8, 605.2, 637.3, 672.9, 718.8, 735.3, 779.5, 829.0] * 1e-6/2;

%initialize beamPath object that you can use for fitting these measurements
goo_NPRO = beamPath;
% your initial gues of what the waist might be, 
%1st argument is waist diameter w0, second is waist location in your coordinate system,
%third argument is wavelength default is 1064e-9 meters
goo_NPRO.seedWaist(100e-6,-0.1);
goox_NPRO = goo_NPRO.fitBeamWidth(zScan_NPRO,xWidth_NPRO);
gooy_NPRO = goo_NPRO.fitBeamWidth(zScan_NPRO,yWidth_NPRO);

%get some information to display
[locx_NPRO , qx_NPRO]= goox_NPRO.getWaists;
radiusx_NPRO = qx_NPRO.beamWidth;
[locy_NPRO , qy_NPRO]= gooy_NPRO.getWaists;
radiusy_NPRO = qy_NPRO.beamWidth;

%make a fugure to show this data and the fit
zdomain = -0.5:.01:0.5;%linspace(min(zScan),max(zScan),100);

title_string = 'NPRO 1661: horiz waist of ' ...
    + compose("%.1f", radiusx_NPRO*1e6) + 'um at '  + compose("%.3f", locx_NPRO) + ...
    'm vert waist of '+ compose("%.1f", radiusy_NPRO*1e6) + 'um at ' + compose("%.3f", locy_NPRO);
figure(1)
clf
hold on
plot(zScan_NPRO,xWidth_NPRO,'bo')
goox_NPRO.plotBeamWidth(zdomain,'b')
goox_NPRO.plotWaists(zdomain,'b')
plot(zScan_NPRO,yWidth_NPRO,'ro')
gooy_NPRO.plotBeamWidth(zdomain,'r')
gooy_NPRO.plotWaists(zdomain,'r')
xlabel('Distance [m]')
ylabel('Beam Size [m]')
title(title_string)
grid on
hold off
f=gcf;
saveas(gcf, 'beamfitting.png')
%%% now add lenses and a target waist, using alog screenshot of Jammt
%add lenses to both vertical and horizontal path
%                  lens syntax:     (focal length,z position,string label)
for Path = [goox_NPRO, gooy_NPRO]
    Path.addComponent(component.lens(.222,.102,'L1'));
    Path.addComponent(component.lens(-0.334,1.144,'L2'));
    Path.addComponent(component.lens(0.222,1.278,'L3'));
    Path.targetWaist(165e-6,1.794);  % goal of what you want the waist to be

end

%check the overlap
modematch = sqrt(goox_NPRO.targetOverlap*gooy_NPRO.targetOverlap);
disp(['modematching = ',num2str(modematch)])



%now, can we do a fitting to get better overlap?  take an average of the
%waist location and size as a starting point for fitting
avg_waist_loc = (locx_NPRO + locy_NPRO)/2;
avg_waist_radius = (radiusx_NPRO + radiusy_NPRO)/2;
newPath = goox_NPRO.duplicate;
newPath.seedWaist(avg_waist_radius, avg_waist_loc)

%%%
%slide only lens 2 + 3 to leave wasit in EOM the same if we can 
newPath = newPath.optimizePath('L2',[1 1.7],'L3',[1 1.7]);
%now check how this solution works for the actaul X and Y path, using fit waists from
%above
newXPath = newPath.duplicate;
newXPath.seedWaist(radiusx_NPRO, locx_NPRO);
newYPath = newPath.duplicate;
newYPath.seedWaist(radiusy_NPRO, locy_NPRO);
new_modematch = sqrt(newXPath.targetOverlap*newYPath.targetOverlap);
disp(['modematching with new component locations= ',num2str(new_modematch)])

disp('original componenet locations')
goox_NPRO.components
disp('new component locations')
newXPath.components

%change plotting domain 
zdomain = 0:0.01:2;
figure(2)
subplot(2,1,1)
hold on % right now all the plot commands act like the matlab plot command and will overwrite the existing
        % figure unless you turn hold on

% The plot commands actually plots two traces, the top and bottom of the beam.
% The output of the plot commands returns the plot handle of the top so when we make the legend we don't 
% have to put a label on the top and bottom of the beam.
yplot = newYPath.plotBeamWidth(zdomain,'b');

xplot = newXPath.plotBeamWidth(zdomain,'r');
newPath.plotComponents(zdomain,0,'r*');
axis tight
newYPath.plotWaists(zdomain, 'b')
newXPath.plotWaists(zdomain, 'r')
legend([yplot, xplot],'vertical','horizontal') % if we didn't use handles we would need to do 
                              %    legend('Y top','Y bottom','X top','X bottom') or something.
ylabel('Beam width (m)')
grid on
hold off

subplot(2,1,2)
hold on
newYPath.plotGouyPhase(zdomain,'wrap','b');
newXPath.plotGouyPhase(zdomain,'wrap','r');
newPath.plotComponents(zdomain,0,'r*');
axis tight
grid on
hold off

ylabel('Gouy Phase (degrees)')
xlabel('axial distance (m)')
saveas(gcf, 'lens_locations.png')