clear all
close all
%%
schnupp=0.05;
green=532e-9;
red=1064e-9;
nR= 1.44963;
nG= 1.46071;
armLen=4000;
figNum=1;
%% Mirrors we'll be potentially using.
EY.ROC=[2245, 2312];
EY.names={'realEY', 'pilotEY'};
EY.Tred=[5e-6, 0.015];
EY.Tgrn=[0.05, 0.0103];
IY.ROC=[2307];
IY.names={'pilotIY'};
IY.Tred=0.0138;
IY.Tgrn=0.0095;
IY.TARred=1-(nR-1)^2/(nR+1)^2; % no AR
IY.TARgrn=1-(nG-1)^2/(nG+1)^2; % no AR
EX.ROC=[2245];
EX.names={'realEX'};
EX.Tred=5E-6;
EX.Tgrn=0.05;
IX.ROC=[1934, 2250.02];
IX.Tred=[0.014, 0.014];
IX.Tgrn=[0.01, 0.01];
IX.TARred=[1, 1-(nR-1)^2/(nR+1)^2]; %the former is with and the latter is without AR
IX.TARgrn=[1, 1-(nG-1)^2/(nG+1)^2];
IX.names={'realIX', 'pilotIX'};

%% nominal ROC of the mirrors.
Ereal.ROC=2245;
Ireal.ROC=1934;

%% First calculate the nominal mode, which is a mean of
% nominal X and nominal Y, which are identical modes except that
% one is shifted by Schnupp asymmetry.
RealArmRoundTripMode=beamPath();
RealArmRoundTripMode.addComponent(component.curvedMirror(Ireal.ROC, 0, 'realITM'));
RealArmRoundTripMode.addComponent(component.curvedMirror(Ereal.ROC, armLen, 'realETM'));
qGreenRealArm=RealArmRoundTripMode.eigenMode(-1e-6, 2*armLen-1e-6, green);
qRedRealArm=RealArmRoundTripMode.eigenMode(-1e-6, 2*armLen-1e-6, red);
% These q parameters are the ones of the beam coming from inside
% the arm just before hitting the ITM.

realArmGreenMode=beamPath();
realArmGreenMode.seedq=qGreenRealArm;
realArmGreenMode.seedz=-1e-6 - schnupp/2; 
% mean nominal mode should be just the X nominal mode shifted by -shnupp/2, or the Y nominal mode shifted by +schnupp/2, they are the same thing. 
f=Ireal.ROC/(nG-1);
realArmGreenMode.addComponent(component.lens(f, -schnupp/2, 'realMeanIX'));

realArmRedMode=beamPath();
realArmRedMode.seedq=qRedRealArm;
realArmRedMode.seedz=-1e-6 -schnupp/2;
f=Ireal.ROC/(nR-1);
realArmRedMode.addComponent(component.lens(f, -schnupp/2, 'realMeanIX'));
        
        
%% X and Y arm with some reality

for ii=1:length(EX.ROC)
    for jj=1:length(IX.ROC)
        idx=jj+(ii-1)*length(IX.ROC);
        
        R=sqrt((1-IX.Tgrn(jj))*(1-EX.Tgrn(ii)));
        XArmGreenFinesse(idx) = pi*sqrt(R)/(1-R);
        XArmGreenTrans(idx) = IX.TARgrn(jj)*IX.Tgrn(jj)*EX.Tgrn(ii)/(1-R)^2;

        R=sqrt((1-IX.Tred(jj))*(1-EX.Tred(ii)));
        XArmRedFinesse(idx) = pi*sqrt(R)/(1-R);
        XArmRedTrans(idx) = IX.TARred(jj)*IX.Tred(jj)*EX.Tred(ii)/(1-R)^2;

        Xname{idx} = [IX.names{jj}, '-', EX.names{ii}];
        XArmRoundTripMode(idx)=beamPath();
        XArmRoundTripMode(idx).addComponent(component.curvedMirror(IX.ROC(jj), 0, IX.names{jj}));
        XArmRoundTripMode(idx).addComponent(component.curvedMirror(EX.ROC(ii), armLen, EX.names{ii}));
        qGreenX(idx)=XArmRoundTripMode(idx).eigenMode(-1e-6, 2*armLen-1e-6, green);
        qRedX(idx)=XArmRoundTripMode(idx).eigenMode(-1e-6, 2*armLen-1e-6, red);
        % These q parameters are the ones of the beam coming from inside
        % the arm just before hitting the ITM.

        XArmGreenMode(idx)=beamPath();
        XArmGreenMode(idx).seedq=qGreenX(idx);
        XArmGreenMode(idx).seedz=-1e-6;
        f=IX.ROC(jj)/(nG-1);
        XArmGreenMode(idx).addComponent(component.lens(f, 0, IX.names{jj}));
        
        XArmRedMode(idx)=beamPath();
        XArmRedMode(idx).seedq=qRedX(idx);
        XArmRedMode(idx).seedz=-1e-6;
        f=IX.ROC(jj)/(nR-1);
        XArmRedMode(idx).addComponent(component.lens(f, 0, IX.names{jj}));
    end
end

for ii=1:length(EY.ROC)
    for jj=1:length(IY.ROC)
        idx=jj+(ii-1)*length(IY.ROC);

        R=sqrt((1-IY.Tgrn(jj))*(1-EY.Tgrn(ii)));
        YArmGreenFinesse(idx) = pi*sqrt(R)/(1-R);
        YArmGreenTrans(idx) = IY.TARgrn(jj)*IY.Tgrn(jj)*EY.Tgrn(ii)/(1-R)^2;

        R=sqrt((1-IY.Tred(jj))*(1-EY.Tred(ii)));
        YArmRedFinesse(idx) = pi*sqrt(R)/(1-R);
        YArmRedTrans(idx) = IY.TARred(jj)*IY.Tred(jj)*EY.Tred(ii)/(1-R)^2;

        Yname{idx} = [IY.names{jj}, '-', EY.names{ii}];
        YArmRoundTripMode(idx)=beamPath();
        YArmRoundTripMode(idx).addComponent(component.curvedMirror(EY.ROC(ii), armLen, EY.names{ii}));
        YArmRoundTripMode(idx).addComponent(component.curvedMirror(IY.ROC(jj), 0, IY.names{jj}));
        qGreenY(idx)=YArmRoundTripMode(idx).eigenMode(-1e-6, 2*armLen-1e-6, green);
        qRedY(idx)=YArmRoundTripMode(idx).eigenMode(-1e-6, 2*armLen-1e-6, red);
        % These q parameters are the ones of the beam coming from inside
        % the arm just before hitting the ITM. Note that Schnupp assymmetry
        % is added later.

        YArmGreenMode(idx)=beamPath();
        YArmGreenMode(idx).seedq=qGreenY(idx);
        YArmGreenMode(idx).seedz=-1e-6-schnupp;
        f=IY.ROC(jj)/(nG-1);
        YArmGreenMode(idx).addComponent(component.lens(f, -schnupp, IX.names{jj}));
        
        YArmRedMode(idx)=beamPath();
        YArmRedMode(idx).seedq=qRedY(idx);
        YArmRedMode(idx).seedz=-1e-6-schnupp;
        f=IY.ROC(jj)/(nR-1);
        YArmRedMode(idx).addComponent(component.lens(f, -schnupp, IX.names{jj}));
        
    end
end

%%
z=-4000:40:2000;
zz=0:40:2000;


for iX=1:length(XArmGreenMode)
    for jY=1:length(YArmGreenMode)
        realq=realArmGreenMode.qPropagate(10);
        XArmGreenMode(iX).targetq=realq;
        XArmGreenMode(iX).targetz=10;
        X_real_coupling = XArmGreenMode(iX).targetOverlap;

        YArmGreenMode(jY).targetq=realq;
        YArmGreenMode(jY).targetz=10;
        Y_real_coupling = YArmGreenMode(jY).targetOverlap;

        xq=XArmGreenMode(iX).qPropagate(10);
        YArmGreenMode(jY).targetq=xq;
        YArmGreenMode(jY).targetz=10;
        XY_coupling = YArmGreenMode(jY).targetOverlap;
        
        figure
        XArmGreenMode(iX).plotBeamWidth(z, 'b');
        YArmGreenMode(jY).plotBeamWidth(z, 'r');
        realArmGreenMode.plotBeamWidth(zz, 'g');
        title(['Green, ' Xname{iX} '-' Yname{jY}]);
        text(-3700, 0.055, ['X: X-nominal coupling ' num2str(X_real_coupling) ', X-Y coupling ' num2str(XY_coupling) ', F=' num2str(XArmGreenFinesse(iX)), ', T=' num2str(XArmGreenTrans(iX))]);
        text(-3700, 0.05, ['Y: Y-nominal coupling ' num2str(Y_real_coupling)  ', F=' num2str(YArmGreenFinesse(jY)), ', T=' num2str(YArmGreenTrans(jY))]);
        text(-3700, 0.045, 'nominal');
        line([-4000, -3800], [0.0557, 0.0557], 'Color', 'b');
        line([-4000, -3800], [0.0507, 0.0507], 'Color', 'r');
        line([-4000, -3800], [0.0457, 0.0457], 'Color', 'g');
        grid on;
        print('-depsc', ['mode' num2str(figNum) '.eps']);
        figNum=figNum+1;

        realq=realArmRedMode.qPropagate(10);
        XArmRedMode(iX).targetq=realq;
        XArmRedMode(iX).targetz=10;
        X_real_coupling = XArmRedMode(iX).targetOverlap;

        YArmRedMode(jY).targetq=realq;
        YArmRedMode(jY).targetz=10;
        Y_real_coupling = YArmRedMode(jY).targetOverlap;

        xq=XArmRedMode(iX).qPropagate(10);
        YArmRedMode(jY).targetq=xq;
        YArmRedMode(jY).targetz=10;
        XY_coupling = YArmRedMode(jY).targetOverlap;
        
        figure
        XArmRedMode(iX).plotBeamWidth(z, 'b');
        YArmRedMode(jY).plotBeamWidth(z, 'r');
        realArmRedMode.plotBeamWidth(zz, 'g');
        title(['Red, ' Xname{iX} '-' Yname{jY}]);
        text(-3700, 0.09, ['X: X-nominal coupling ' num2str(X_real_coupling) ', X-Y coupling ' num2str(XY_coupling)  ', F=' num2str(XArmRedFinesse(iX)), ', T=' num2str(XArmRedTrans(iX))]);
        text(-3700, 0.08, ['Y: Y-nominal coupling ' num2str(Y_real_coupling)  ', F=' num2str(YArmRedFinesse(jY)), ', T=' num2str(YArmRedTrans(jY))]);
        text(-3700, 0.07, 'nominal');
        line([-4000, -3800], [0.09, 0.09], 'Color', 'b');
        line([-4000, -3800], [0.08, 0.08], 'Color', 'r');
        line([-4000, -3800], [0.07, 0.07], 'Color', 'g');
        grid on;
        print('-depsc', ['mode' num2str(figNum) '.eps']);
        figNum=figNum+1;
    end
end

system('cat mode*.eps > modes.eps; ps2pdf modes.eps')