clear all
close all

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

%% some constants. Don't change these.
YAG=1.064e-6;% IR light
HeNe=633e-9;
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=(-0.7:0.1:0.3)*1e-3; % offset of the telescope, in mm
figidx=0;
detuningTolerance=[-0.5 +0.2]*1e-3;

%% telescope parameters.
nWaist= 2;
returnROC=zeros(nWaist, length(offsets)); % storage space for return beam geometry  
returnWidth=zeros(nWaist, length(offsets));
returnWaist=zeros(nWaist, length(offsets));
returnZ=zeros(nWaist, length(offsets));
returnDiv=zeros(nWaist, length(offsets));

%% Set MM input parameter correctly if you want to use it for final adjustment.
retro2sec=0.63627 + 12.7e-3;  % distance from the 2nd to the outside of the Tele Back frame plate
MM2retro= 2.695;              % distance from MM head to the telescope plate


%% input parameters. Change these as you wish.
Mx = 0.98;  % X-axis M^2 value from the MM scan
My = 1.04;  % Y-axis M^2 value from the MM scan
MM2sourceretro = 0.47; % distance from the MM to the last steering mirror (used as a retro mirror)
                        % on the pedestal
% Taken from MM Scan 19Oct2011/ScourceScan1_,
% it has 9 averages and calculated the means.
wSH=[5.658/Mx 5.821/My]*1e-3/2; % input beam waist radius [X Y] in mm.
waistPos= -1*[4.813 3.778] + MM2retro + retro2sec - 2*MM2sourceretro; % zero=at the secondary, minus means upstream of the secondary
% the above two are used as a pair and should have the same number of data.


%% Modemaster double path results
% Taken from 19Oct2011/TeleScan3_
fname = 'TeleScan3';
Mx = 1.25;
My = 1.17;
MMW=[0.532 0.565]*1e-3/2; %radius mean [X Y]
MMWstd=[0.0158 0.0196]*1e-3/2; % radius std [X Y]
%
MMZ= - [-1.315, -1.169]-(MM2retro+retro2sec) ; % waist pos mean [X Y]
 % Zero means that the waist is at the 2ndary
 % Positive mean = waist is downstream of the 2ndary towards MM
MMZstd=[ 0.0237 0.0208];
%
MMDiv=[3.17 2.80] /2 *1e-3; % divergence angle (half), radian
MMDivstd=[0.073     0.074]/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];
%%
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

        %% 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 mibrror including the lens of the substrate
        ReflectorROC1=inf; % flat mirror

        dTele1_Reflector=0.9-offsets(ii); % 0.9m is somewhat arbitrary,
        zReflector=zTele1+dTele1_Reflector;
        
        % add a return path
        zTele1_return = zReflector+dTele1_Reflector;
        zTele2_return = zTele1_return + dTele2_1;
        zZeroPosition_return = zTele2_return+zTele2;
        %%
        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;
        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 - ' fname])
            figidx=figidx+1;
            if savefig
                print('-depsc', ['beamsize_flatETM_' fname '_' num2str( figidx) '.eps'])
            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_interpZ(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
    
    hdl = figure(kk+4)
    subplot 311
    plot(offsets, returnWaist(kk, :), 'r', 'LineWidth', 2);%, offsets, returnWidth(kk, :));
    hold on
    errorbar(off_interpW(kk), MMW(kk),  MMWstd(kk), 'g', 'LineWidth', 2);
    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(10)');%, 'width@sec')
    ylabel('return waist radius (m)')
    grid on
    title(['w=' num2str(wSH(kk)*1e3) 'mm waist injected - ' fname])
    subplot 312
    %returnZ = returnZ+(MM2retro+retro2sec);
    %MMZ = MMZ + (MM2retro+retro2sec);
    plot(offsets, returnZ(kk, :), 'r', 'LineWidth', 2);
    hold on
    errorbar(off_interpZ(kk), MMZ(kk),  MMZstd(kk), 'g', 'LineWidth', 2);
    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', 'LineWidth', 2)
    hold on
    errorbar(off_interpDiv(kk), MMDiv(kk),  MMDivstd(kk), 'g', 'LineWidth', 2);
    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;
    orient(hdl, 'landscape');
    if savefig
        print('-depsc', ['beamsize_flatETM_' fname '_' num2str( figidx) '.eps'])
    end
    
    %%
end
%%

if savefig
    system(['cat beamsize_flatETM_' fname '_?.eps > beamsize_design_flatETM_all_' fname '_' date '.eps; /usr/local/bin/ps2pdf beamsize_design_flatETM_all_' fname '_' date '.eps'])
end
