clear all
close all

% Set this flag to save the figures
savefig = 0;

%% some constants. Don't change these.
YAG=1.064e-6;% IR light
nom_Det=nominalDetuning(); % nominal detuning defined in another file
n_air= 1.00027398; %http://refractiveindex.info/?group=GASES&material=Air, but this is irrelevant and is not used.
offsets=(-2.5:0.02:3)*1e-3; % offset of the telescope, in mm
figidx=0;
detuningTolerance=[-0.5 +0.2]*1e-3;

%% input parameters obtained when the source parameters were measured.
% Change these as you wish.

% Taken from MM Scan 0725_src_doublet,
% it has 10 averages.
Mx = 1.10;  % X-axis M^2 value from the MM scan
My = 1.13;  % Y-axis M^2 value from the MM scan
wSH=[4.48 4.325]*1e-3/2 ; % input beam waist radius [X Y] in mm.
% Mx and My are used here to estimate the Gaussian waist.

MM2sourceLauncher=.134+.216+.124+.284; % distance from the MM to the launcher
% when  the source parameters were measured.

MM2sourceWaist= [10.774, 10.839];
% -1.46 and -3.67 are numbers reported by MM. Negative means closer to the
% source, positive means "into the head".
% Since it's negative and its absolute value is larger than
% MM2sourceLauncher, the waist is "into the launcher".


launcher2waist = MM2sourceWaist+ MM2sourceLauncher;
% If this is positive, the waist is away from the launcher toward the TMS.
% If negative, the waist is away from the TMS "into the launcher".

%% MM parameters for double path measurement.

% Distance from the launcher to the 2ndary, from the MM to the secondary, and
% from the 8" retro mirror to the primary.
% Measure these each time you set the experimental apparatus up.
% 10mm error is OK for the first two numbers.
% 50mm error is more than OK for pri2retro.
launcher2sec=0.134+.216+.115+.204+1.712+.500;  % distance from the 2ndary to the source launcher.
MM2sec = launcher2sec-0.134-0.216+.160; % distance from the 2ndary to the MM.
pri2retro=87e-3+688e-3; % retro 2 outer side of the 2ndary plate , plus 2dary plate to the center of primary.
% this is a measured number (Jul24/2013) but a few cm error doesn't affect
% anything.

waistPos=-launcher2waist + launcher2sec; % distance from the secondary to the input waist.
% Positive means that the waist is upstream of the secondary (i.e. closer to the launcher),
% negative means that it is into the telescope.

% Modemaster double path results
Ndata=9;% Note that this doesn't necessarily agree with the number of scans. You need to look in the .std file to find out how many files were used for calculation.
Mx = 1.23;
My = 1.07;
MMW=[0.652 0.655]*1e-3/2; %radius mean [X Y]
%MMW=[0.724 0.826]*1e-3/2; %radius mean [X Y]
MMWstd=[0.0128 0.0140]*1e-3/2; % radius std [X Y]
%
MMZ= -([-0.856, -0.515] + MM2sec); % waist pos mean [X Y]
%MMZ= -([-0.702, -0.702] + MM2sec); % waist pos mean [X Y]
 % Inside [] are the numbers reported by MM.
 % Note somewhat strange sign convention here, this is so that
 % this parameter is interpreted as the real part of the Q parameter of the
 % returning beam at the secondary.
 % MMZ=zero means that the waist is at the 2ndary.
 % Negative means that the waist is towards MM 
 % (e.g. if MMZ=-MM2sec then the waist is at MM).
 % Positive means that the waist is into the telescope (because, from the
 % Q parameter point of view, the beam travelled by MMZ after the waist).
 
MMZstd=[ 0.0209 0.0208];
%
MMDiv=[2.56 2.22] /2 *1e-3; % divergence angle (half), radian
MMDivstd=[0.046 0.069]/2 *1e-3;% divergence angle (half), std
% actual detuning based on interpolation from MM results
off_interpW=[0 0];
off_interpZ=[0 0];
off_interpDiv=[0 0];

%%
% set up the storage container for return beam  geometry
nWaist= 2;
returnROC=zeros(nWaist, length(offsets));
returnWidth=zeros(nWaist, length(offsets));
returnWaist=zeros(nWaist, length(offsets));
returnZ=zeros(nWaist, length(offsets));
returnDiv=zeros(nWaist, length(offsets));

labels={'X', 'Y'};
for kk = 1:length(wSH);
    for ii=1:length(offsets);
        clear mypath;
        %% Source laser beam
        mypath = beamPath;

        mypath.seedz=0-1e-6; % secondary position
        mypath.seedq = beamq(-waistPos(kk) + 1i*pi*(wSH(kk))^2 / YAG ); % q parameter at the secondary
        % Note that positive waistPos means that the waist is closer to the
        % launcher, negative means that it is into the tele, so we're
        % changing sign here.

        %% telescope system
        Tele2ROC= -0.2; % secondary focal length -0.1m
        zTele2 = 1e-6; %on the secondary is the waist of S-H beam when offsets==0.
        
        Tele1ROC=4; %primary focal length 2m
        % nominally correct pri-sec spacing when offsets=0.
        % This is 2.6mm off of a truly afocal system,
        % i.e. distance = 1.9m + 2.6mm.
        % Though it is different from the actual procedure, we assume that
        % the secondary doesn't move and instead the primary is adjusted.
        dTele2_1=((Tele1ROC+Tele2ROC)/2 +nom_Det+offsets(ii));%*n_air;
        zTele1 = zTele2+dTele2_1;
        
%        ReflectorROC1=-1.548671e+03; % ETM-like mirror including the lens of the substrate
        ReflectorROC1=inf; % flat mirror

        dTele1_Reflector=pri2retro;
        zReflector=zTele1+dTele1_Reflector;
        
        % add a return path
        zTele1_return = zReflector+dTele1_Reflector;
        zTele2_return = zTele1_return + dTele2_1;
        zZeroPosition_return = zTele2_return+1e-6; %to make sure that the zero position is after the secondary.
        %%
        mypath.addComponent(component.curvedMirror(Tele2ROC, zTele2, 'sec1'));
        mypath.addComponent(component.curvedMirror(Tele1ROC, zTele1, 'pri1'));
        mypath.addComponent(component.curvedMirror(ReflectorROC1, zReflector, 'refl'));
        mypath.addComponent(component.curvedMirror(Tele1ROC, zTele1_return, 'pri2'));
        mypath.addComponent(component.curvedMirror(Tele2ROC, zTele2_return, 'sec2'));
        %%
        qout = mypath.qPropagate(zZeroPosition_return);
        mypath.targetq=qout;
        returnROC(kk, ii)=qout.radiusOfCurvature;
        returnWidth(kk, ii)=qout.beamWidth;
        returnWaist(kk, ii)=qout.waistSize;
        returnZ(kk, ii)=qout.waistZ; % positive means that the waist is already past (i.e. into the telescope).
        returnDiv(kk, ii)=qout.divergenceAngle;
        mypath.targetz=zZeroPosition_return;
        qReflector = mypath.qPropagate(zReflector);
        ReflectorP=beamPath;
        ReflectorP.seedWaist(1, 1,  YAG/n_air); % this is just to set lambda
        ReflectorP.seedq=qReflector;
        ReflectorP.seedz=zReflector;
        %%
        if offsets(ii)==0 |offsets(ii)==-2;
            figure
            z=-1:0.01:(zZeroPosition_return+1);
            mypath.plotBeamWidth(z);
            hold on
%            set(gca, 'yscale', 'log')
            mypath.plotComponents(z);
            mypath.plotBeams(z);
            ReflectorP.plotBeams(z);
            title(['w=' num2str(wSH(kk)*1e3) 'mm waist injected, ' num2str(offsets(ii) ) ' mm offset in primary-secondary distance, flat ETM'])
            figidx=figidx+1;
            if savefig
                print('-depsc', ['beamsize_flatETM' num2str( figidx) '.ps'])
            end
        end
    end
    
    % find detuning point close to measured W and Z
    [dummy, idx]=min(abs(returnWaist(kk,:)-MMW(kk)));
    if idx(1)==1
        idx1=1; idx2=2;
    elseif idx(1)==length(offsets);
        idx1=length(offsets)-1; idx2=length(offsets);
    elseif (MMW(kk)-returnWaist(kk,idx(1)))*(MMW(kk)-returnWaist(kk, idx(1)+1))<0
        idx1=idx(1); idx2=idx1+1;
    else
        idx2=idx(1); idx1=idx2-1;
    end
    
    if returnWaist(kk, idx2)==returnWaist(kk, idx1)
        off_interpW(kk)=mean(offsets(idx1:idx2));
    else
        off_interpW(kk)=(offsets(idx1)*(returnWaist(kk, idx2)-MMW(kk)) + offsets(idx2)*(MMW(kk)-returnWaist(kk, idx1)))/(returnWaist(kk, idx2)-returnWaist(kk, idx1));
    end

    [dummy, idx]=min(abs(returnZ(kk,:)-MMZ(kk)));
    if idx(1)==1
        idx1=1; idx2=2;
    elseif idx(1)==length(offsets);
        idx1=length(offsets)-1; idx2=length(offsets);
    elseif (MMZ(kk)-returnZ(kk, idx(1)))*(MMZ(kk)-returnZ(kk, idx(1)+1))<0
        idx1=idx(1); idx2=idx1+1;
    else
        idx2=idx(1); idx1=idx2-1;
    end
    if returnZ(kk, idx2)==returnZ(kk, idx1)
        off_interpZ(kk)=mean(offsets(idx1:idx2));
    else
        off_interpZ(kk)=(offsets(idx1)*(returnZ(kk, idx2)-MMZ(kk)) + offsets(idx2)*(MMZ(kk)-returnZ(kk, idx1)))/(returnZ(kk, idx2)-returnZ(kk, idx1));
    end

    [dummy, idx]=min(abs(returnDiv(kk,:)-MMDiv(kk)));
    if idx(1)==1
        idx1=1; idx2=2;
    elseif idx(1)==length(offsets);
        idx1=length(offsets)-1; idx2=length(offsets);
    elseif (MMDiv(kk)-returnDiv(kk, idx(1)))*(MMDiv(kk)-returnDiv(kk, idx(1)+1))<0
        idx1=idx(1); idx2=idx1+1;
    else
        idx2=idx(1); idx1=idx2-1;
    end
    if returnDiv(kk, idx2)==returnDiv(kk, idx1)
        off_interpDiv(kk)=mean(offsets(idx1:idx2));
    else
        off_interpDiv(kk)=(offsets(idx1)*(returnDiv(kk, idx2)-MMDiv(kk)) + offsets(idx2)*(MMDiv(kk)-returnDiv(kk, idx1)))/(returnDiv(kk, idx2)-returnDiv(kk, idx1));
    end
    
    figure(kk+4)
    subplot 311
    plot(offsets, returnWaist(kk, :), 'r');%, offsets, returnWidth(kk, :));
    hold on
    errorbar(off_interpW(kk), MMW(kk),  MMWstd(kk), 'g');
    plot(off_interpW(kk), MMW(kk), 'g*');    
    plot([1 1]*detuningTolerance(1), [min(returnWaist(kk, :)), max(returnWaist(kk, :))], 'k', ...
        [1 1]*detuningTolerance(2), [min(returnWaist(kk, :)), max(returnWaist(kk, :))], 'k')
    ylim([min(returnWaist(kk, :)), max(returnWaist(kk, :))]);
    legend('Gaussian calc from measured input', ['measured \pm std(' num2str(Ndata) ')']);%, 'width@sec')
    ylabel('return waist radius (m)')
    grid on
    title([labels{kk} ', w=' num2str(wSH(kk)*1e3) 'mm waist injected'])
    subplot 312
    plot(offsets, returnZ(kk, :), 'r');
    hold on
    errorbar(off_interpZ(kk), MMZ(kk),  MMZstd(kk), 'g');
    plot(off_interpZ(kk), MMZ(kk),  'g*');
    plot([1 1]*detuningTolerance(1), [min(returnZ(kk, :)), max(returnZ(kk, :))], 'k', ...
        [1 1]*detuningTolerance(2), [min(returnZ(kk, :)), max(returnZ(kk, :))], 'k')
    ylim([min(returnZ(kk, :)), max(returnZ(kk, :))]);
    ylabel('return waist pos (m)')
    grid on
    subplot 313
    plot(offsets, returnDiv(kk, :), 'r')
    hold on
    errorbar(off_interpDiv(kk), MMDiv(kk),  MMDivstd(kk), 'g');
    plot(off_interpDiv(kk), MMDiv(kk),  'g*');
    plot([1 1]*detuningTolerance(1), [min(returnDiv(kk, :)), max(returnDiv(kk, :))], 'k', ...
        [1 1]*detuningTolerance(2), [min(returnDiv(kk, :)), max(returnDiv(kk, :))], 'k')
    ylim([min(returnDiv(kk, :)), max(returnDiv(kk, :))]);

    grid on
    ylabel('return div angle (rad)')
    xlabel('tele offset (m)');
    %FillPage('tall')
    figidx=figidx+1;
    if savefig
        print('-depsc', ['beamsize_flatETM' num2str( figidx) '.ps'])
    end
    
    %%
end
%%
if savefig
    system(['cat beamsize_flatETM?.ps > beamsize_design_flatETM_all_' date '.ps; ps2pdf beamsize_design_flatETM_all_' date '.ps; rm *ps'])
end
