clear classes
close all
clear all 
clc

set(0,'DefaultLineLineWidth',2)
set(0,'DefaultAxesFontSize',12)
set(0,'DefaultAxesXGrid','on')
set(0,'DefaultAxesYGrid','on')
set(0,'DefaultAxesZGrid','on')
legFontSize = 10;

%%
printFigs = 1;

gitRoot = '/Users/kissel/ligo/gitcommon/';
svnRoot = '/Users/kissel/ligo/svncommon/SeiSVN/seismic/';
resultsDir = [svnRoot 'Common/SPI/Results/FiberCollimators/'];
toolsDir = [svnRoot 'Common/MatlabTools/'];
jammtDir = [gitRoot 'jammt/'];

fcsn = 'S0272503';
ftsn = 'S3228003';
measDate = {'2026-02-03'...
            '2026-02-03'};

author = ', J. Kissel';
measNote = {'OzOptics';...
            'AxcelPhotonics'};

figTag = [measDate{end} '_spi_fc_' fcsn '_ft_' ftsn '_7p5mPatchCord'];

%% Settings and raw measurements
w0d = 1.05e-3; % Desired waist size, T2400413
z0d = 0; % "The best, most collimated beam has its waist at the collimator lens"
         % and I set up the beam profiler distances with the collimator
         % lens at z = 0 (+/- 5 mm).
lambda = 1064e-9;

% Measured Data (Use D4sig, mean values "as is," i.e. as beam *diameter* in [um])
zPosition_m = [0.508 0.991 1.499 2.007 3.251 4.496 5.41];

xRadius_m =   {[1778  1535  1400  1435  2153  3221  4107]*1e-6/2;... % OzOptics
               [1785  1547  1408  1444  2150  3223  4089]*1e-6/2};   % AxcelPhotonics
yRadius_m =   {[1787  1552  1454  1501  2193  3242  4072]*1e-6/2;... % pre-bake
               [1792  1556  1455  1505  2196  3254  4074]*1e-6/2};   % post-ClassB-bake

for iMeas = 1:size(xRadius_m,1)
    outName{iMeas} = [jammtDir figTag '_' measNote{iMeas} 'Laser.txt'];
    out{iMeas} = [transpose(zPosition_m * 1e2) transpose(xRadius_m{iMeas})*1e6 transpose(yRadius_m{iMeas})*1e6];

    dlmwrite(outName{iMeas},out{iMeas},'\t')
end

% In pre-bake data, collimator flat parallel to optical table.
% In post-bake data, collimator flat is perpendicular to the table, per
% design (see Figure 5 of T2400413).
% So, I flip-flip the 
%     - x (NanoScan Axis 1, Parallel to the Table Surface) and
%     - y (NanoScan Axis 2, Perpendicular to the Table Surface)
% for the pre-bake data set, as the NanoScan axes orientation w.r.t. the
% table was the same for both setups.
xRadius_m_temp = xRadius_m{1};
yRadius_m_temp = yRadius_m{1};
xRadius_m{1} = yRadius_m_temp;
yRadius_m{1} = xRadius_m_temp;

% Set Initial Guess:
measRange_m = [0 20];   % essentially, the desired xlim of the plot
gRadius_m = 1800e-6/2; % Guess based on the data
zGuess_m = 1.5;        % Guess based on the data

%% Fit the data
for iMeas = 1:length(xRadius_m)
    model(iMeas).goo = beamPath;

    % Creat a position vector for the model
    zModel_m = measRange_m(1):.01:measRange_m(2);

    % Initialize the fitting object
    model(iMeas).goo.seedWaist(gRadius_m,zGuess_m);

    % Fit the data
    model(iMeas).goox = model(iMeas).goo.fitBeamWidth(zPosition_m,xRadius_m{iMeas});
    model(iMeas).gooy = model(iMeas).goo.fitBeamWidth(zPosition_m,yRadius_m{iMeas});

    % Assume only one waist within the measurement range, and find its
    % properties
    [model(iMeas).z0x_m, model(iMeas).q0x] = model(iMeas).goox.getWaists(measRange_m(1),measRange_m(2));
    [model(iMeas).z0y_m, model(iMeas).q0y] = model(iMeas).gooy.getWaists(measRange_m(1),measRange_m(2));
    model(iMeas).w0x_m = model(iMeas).q0x.waistSize;
    model(iMeas).w0y_m = model(iMeas).q0y.waistSize;
    model(iMeas).zRx_m = pi*(model(iMeas).w0x_m.^2)/lambda;
    model(iMeas).zRy_m = pi*(model(iMeas).w0y_m.^2)/lambda;
    model(iMeas).wzzRx_m = model(iMeas).w0x_m * sqrt(2);  % w0 * sqrt(1 + (zR / zR)^2)
    model(iMeas).wzzRy_m = model(iMeas).w0y_m * sqrt(2);
    model(iMeas).thetaDx = lambda/(pi*model(iMeas).w0x_m);
    model(iMeas).thetaDy = lambda/(pi*model(iMeas).w0y_m);
    model(iMeas).wffx = (measRange_m(end) - model(iMeas).z0x_m)*tan(model(iMeas).thetaDx);
    model(iMeas).wffy = (measRange_m(end) - model(iMeas).z0y_m)*tan(model(iMeas).thetaDy);
end

%%
figNum = 1; 

xLimits = {[0 6];
           measRange_m};
yLimits = {[-0.25e-3 2.5e-3];
           [-0.25e-3 7.0e-3]};
figExt = {'Near';
          'Far'};

%%
for iMeas = 1:length(xRadius_m)
    for iZoom = 1:2
        figure(figNum)
        figNames{figNum} = [resultsDir figTag '_beamprofile_fit_' measNote{iMeas} '_' figExt{iZoom} 'Field.pdf'];
        figNum = figNum + 1;
        clf
        hold on
        plot(zPosition_m,xRadius_m{iMeas},'bo')
        model(iMeas).goox.plotBeamWidth(zModel_m,'b')
        pw0x = plot(model(iMeas).z0x_m,model(iMeas).w0x_m,'bx');
        pwZx = plot(model(iMeas).z0x_m+model(iMeas).zRx_m,model(iMeas).wzzRx_m,'bs');
        divx = plot([model(iMeas).z0x_m measRange_m(end)],[0 model(iMeas).wffx],'b--');
        plot(zPosition_m,yRadius_m{iMeas},'ro')
        model(iMeas).gooy.plotBeamWidth(zModel_m,'r')
        pw0y = plot(model(iMeas).z0y_m,model(iMeas).w0y_m,'rx');
        pwZy = plot(model(iMeas).z0y_m+model(iMeas).zRy_m,model(iMeas).wzzRy_m,'rs');
        divy = plot([model(iMeas).z0y_m measRange_m(end)],[0 model(iMeas).wffy],'r--');
        hold off
        set(pw0x,'MarkerSize',14)
        set(pw0y,'MarkerSize',14)
        set(pwZx,'MarkerSize',16)
        set(pwZy,'MarkerSize',16)
        leg = legend('Measured Data, X',...
            'Fit, X',...
            ['       (w0x, z0x) = (' num2str(model(iMeas).w0x_m*1e3) ',' num2str(model(iMeas).z0x_m) ') [mm,m]'],...
            'Fit Waist Radius, w0x',...
            'Beam Radius at zRx (+ z0x); w(zR+z0x) = sqrt(2)*w0x',...
            'Divergence, Beam Radius at end = (end - z0x)*atan(\lambda / [\pi * w0x])',...
            'Measured Data, Y',...
            'Fit, Y',...
            ['       (w0y, z0y) = (' num2str(model(iMeas).w0y_m*1e3) ',' num2str(model(iMeas).z0y_m) ') [mm,m]'],...
            'Fit Waist Radius, w0y',...
            'Beam Radius at zRy (+ z0y); w(zR+z0y) = sqrt(2)*w0y',...
            'Divergence, Beam Radius at end = (end - z0y)*atan(\lambda / [\pi * w0y])',...
            'Location','NorthWest');
        xlabel('Distance From Collimator [m]')
        ylabel('Beam Waist [m]')
        set(gca,'XMinorGrid','on','YMinorGrid','on')
        set(leg,'FontSize',legFontSize)
        title({[fcsn ' SuK Fiber Collimator Characterization'];...
            [measDate{iMeas} ' (' measNote{iMeas} '); ' figExt{iZoom} ' Field']})
        xlim(xLimits{iZoom})
        ylim(yLimits{iZoom})

    end
end

%%
iZoom = 1;
legFontSize = 8;

figure(figNum)
figNames{figNum} = [resultsDir figTag '_beamprofile_fit_' figExt{iZoom} 'Field_comparison.pdf'];
figNum = figNum + 1;

subplot(211)
hold on
plot(zPosition_m,xRadius_m{1},'mo')
model(1).goox.plotBeamWidth(zModel_m,'m')
model(1).pw0x = plot(model(1).z0x_m,model(1).w0x_m,'mx');
plot(zPosition_m,xRadius_m{2},'ro')
model(2).goox.plotBeamWidth(zModel_m,'r')
model(2).pw0x = plot(model(2).z0x_m,model(2).w0x_m,'rx');
hold off
set(model(1).pw0x,'MarkerSize',14)
set(model(2).pw0x,'MarkerSize',14)
leg = legend([measNote{1} ' Data, X'],...
             [measNote{1} ' Fit, X'],...
             ['       ' measNote{1} ' (w0x, z0x) = (' num2str(model(1).w0x_m*1e3) ',' num2str(model(1).z0x_m) ') [mm,m]'],...
             [measNote{1} ' Fit Waist Radius, w0x'],...
             [measNote{2} ' Data, X'],...
             [measNote{2} ' Fit, X'],...
             ['       ' measNote{2} ' (w0x, z0x) = (' num2str(model(2).w0x_m*1e3) ',' num2str(model(2).z0x_m) ') [mm,m]'],...
             [measNote{2} ' Fit Waist Radius, w0x'],...
             'Location','NorthWest');
xlabel('Distance From Collimator [m]')
ylabel('Beam Radius [m]')
set(gca,'XMinorGrid','on','YMinorGrid','on')
set(leg,'FontSize',legFontSize)
title({[fcsn ' SuK Fiber Collimator Characterization; X Radii (p-pol) ' figExt{iZoom} ' Field'];...
    [measNote{1} ' (' measDate{1} ') vs. ' measNote{2} '(' measDate{2} '); ' figExt{iZoom} ' Field']})
xlim(xLimits{iZoom})
ylim(yLimits{iZoom})

subplot(212)
hold on
plot(zPosition_m,yRadius_m{1},'co')
model(1).gooy.plotBeamWidth(zModel_m,'c')
model(1).pw0y = plot(model(1).z0y_m,model(1).w0y_m,'cx');
plot(zPosition_m,yRadius_m{2},'bo')
model(2).gooy.plotBeamWidth(zModel_m,'b')
model(2).pw0y = plot(model(2).z0y_m,model(2).w0y_m,'bx');
hold off
set(model(1).pw0y,'MarkerSize',14)
set(model(2).pw0y,'MarkerSize',14)
leg = legend([measNote{1} ' Data, Y'],...
             [measNote{1} ' Fit, Y'],...
             ['       ' measNote{1} ' (w0y, z0y) = (' num2str(model(1).w0y_m*1e3) ',' num2str(model(1).z0y_m) ') [mm,m]'],...
             [measNote{1} ' Fit Waist Radius, w0y'],...
             [measNote{2} ' Data, Y'],...
             [measNote{2} ' Fit, Y'],...
             ['       ' measNote{2} ' (w0y, z0y) = (' num2str(model(2).w0y_m*1e3) ',' num2str(model(2).z0y_m) ') [mm,m]'],...
             [measNote{2} ' Fit Waist Radius, w0y'],...
             'Location','NorthWest');
xlabel('Distance From Collimator [m]')
ylabel('Beam Radius [m]')
set(gca,'XMinorGrid','on','YMinorGrid','on')
set(leg,'FontSize',legFontSize)
title({[fcsn ' SuK Fiber Collimator Characterization; Y Radii (s-pol) ' figExt{iZoom} ' Field'];...
    [measNote{1} ' (' measDate{1} ') vs. ' measNote{2} ' (' measDate{2} '); ' figExt{iZoom} ' Field']})
xlim(xLimits{iZoom})
ylim(yLimits{iZoom})

%%
if printFigs
    addpath(toolsDir)
    mergeCommand = [svnRoot 'Common/Misc/pdfmerge ' resultsDir figTag '_beamprofile_fit.pdf'];

    for iFig = [5 1:4]
        figure(iFig)
        FillPage('w')
        IDfig(author)
        saveas(gcf,figNames{iFig})

        mergeCommand = [mergeCommand ' ' figNames{iFig}];
    end
end
