% Plot scaled IMC input beam jitter 
% P. Fulda 2013/10/18

clear all

% Mean counts for two transmission sensors in aligned state
alignedIM4T=8987;
alignedMC2T=1180;

% IMC parameters
lambda=1064e-9;
w0=2.1255e-3;
thetaIMC=1/(pi*w0/lambda);
Rc=27.275;
d=0.465/2;
L=sqrt(16.2406^2-d^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate jitter -> RIN couplings for all 4 D.O.Fs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Yaw tilt (MC2 yaw)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mean misaligned transmission sensor counts
meanIM4T_MC2Y=7300;
meanMC2T_MC2Y=900;

% Waist tilt as function of MC2 yaw
MC2Ytilt2waist_tilt=Rc/(Rc-L-d);  % From Kawazoe et al eq. 16
% Calculate average power drop from both sensors
MC2Y_pdrop=1-mean([meanIM4T_MC2Y/alignedIM4T,meanMC2T_MC2Y/alignedMC2T]);
% Calculate MC2 Yaw offset from power drop
MC2Y_offset=thetaIMC*sqrt(MC2Y_pdrop)/MC2Ytilt2waist_tilt;
% Coupling factor from input jitter in radians to RIN for calculated MC2 yaw offset
MC2Y_inputJ2RIN=2*MC2Y_offset*MC2Ytilt2waist_tilt/thetaIMC^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pitch tilt (MC1/3 diff pitch)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mean misaligned transmission sensor counts
meanIM4T_MC13diffP=6700;
meanMC2T_MC13diffP=820;

% Waist tilt as function of MC13 differential pitch
MC13diffP2waist_tilt=1/sqrt(2); % From Kawazoe et al eq. 65
% Calculate average power drop from both sensors
MC13diffP_pdrop=1-mean([meanIM4T_MC13diffP/alignedIM4T,meanMC2T_MC13diffP/alignedMC2T]);
% Calculate MC1/3 diff pitch offset from power drop
MC13diffP_offset=thetaIMC*sqrt(MC13diffP_pdrop)/MC13diffP2waist_tilt;
% Coupling factor from input jitter in radians to RIN for calculated MC1/3 diff pitch offset
MC13diffP_inputJ2RIN=2*MC13diffP_offset*MC13diffP2waist_tilt/thetaIMC^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pitch shift (MC2 pitch)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mean misaligned transmission sensor counts
meanIM4T_MC2P=7500;
meanMC2T_MC2P=1000;

% Waist shift as function of MC2 pitch
MC2P2waist_shift=Rc; % From Kawazoe et al eq. 60
% Calculate average power drop from both sensors
MC2P_pdrop=1-mean([meanIM4T_MC2P/alignedIM4T,meanMC2T_MC2P/alignedMC2T]);
% Calculate MC2 pitch offset from power drop
MC2P_offset=w0*sqrt(MC2P_pdrop)/MC2P2waist_shift;
% Coupling factor from input jitter waist translation in m to RIN for calculated MC2 pitch offset
MC2P_inputJ2RIN=2*MC2P_offset*MC2P2waist_shift/w0^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Yaw shift (MC1/3 diff yaw)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mean misaligned transmission sensor counts
meanIM4T_MC13diffY=5150;
meanMC2T_MC13diffY=700;

% Waist shift as function of MC1/3 diff yaw
MC13diffY2waist_shift=sqrt(L^2+d^2); % From Kawazoe et al eq. 4
% Calculate average power drop from both sensors
MC13diffY_pdrop=1-mean([meanIM4T_MC13diffY/alignedIM4T,meanMC2T_MC13diffY/alignedMC2T]);
% Calculate MC1/3 diff yaw offset from power drop
MC13diffY_offset=w0*sqrt(MC13diffY_pdrop)/MC13diffY2waist_shift;
% Coupling factor from input jitter waist translation in m to RIN for calculated MC13 diff yaw offset
MC13diffY_inputJ2RIN=2*MC13diffY_offset*MC13diffY2waist_shift/w0^2;

% loading data...
data=load('Commiss_MC2Y.txt');
MC2Yfreq=data(:,1);
% Normalize by mean to get RIN
MC2Y_RIN_IM4T=data(:,2)/meanIM4T_MC2Y;
MC2Y_RIN_MC2T=data(:,3)/meanMC2T_MC2Y;

% Convert to rad/rt(Hz)
MC2Y_JIT_IM4T=MC2Y_RIN_IM4T/MC2Y_inputJ2RIN;
MC2Y_JIT_MC2T=MC2Y_RIN_MC2T/MC2Y_inputJ2RIN;
% Convert to HG10 mode amplitude/rt(Hz)
MC2Y_HG10amp_IM4T=MC2Y_RIN_IM4T/MC2Y_inputJ2RIN/thetaIMC;
MC2Y_HG10amp_MC2T=MC2Y_RIN_MC2T/MC2Y_inputJ2RIN/thetaIMC;

data=load('Commiss_MC13diffP.txt');
MC13diffPfreq=data(:,1);
% Normalize by mean to get RIN
MC13diffP_RIN_IM4T=data(:,2)/meanIM4T_MC13diffP;
MC13diffP_RIN_MC2T=data(:,3)/meanMC2T_MC13diffP;
% Convert to rad/rt(Hz)
MC13diffP_JIT_IM4T=MC13diffP_RIN_IM4T/MC13diffP_inputJ2RIN;
MC13diffP_JIT_MC2T=MC13diffP_RIN_MC2T/MC13diffP_inputJ2RIN;
% Convert to HG10 mode amplitude/rt(Hz)
MC13diffP_HG10amp_IM4T=MC13diffP_RIN_IM4T/MC13diffP_inputJ2RIN/thetaIMC;
MC13diffP_HG10amp_MC2T=MC13diffP_RIN_MC2T/MC13diffP_inputJ2RIN/thetaIMC;

data=load('Commiss_MC2P.txt');
MC2Pfreq=data(:,1);
MC2P_RIN_IM4T=data(:,2)/meanIM4T_MC2P;
MC2P_RIN_MC2T=data(:,3)/meanMC2T_MC2P;
% Convert to m/rt(Hz)
MC2P_JIT_IM4T=MC2P_RIN_IM4T/MC2P_inputJ2RIN;
MC2P_JIT_MC2T=MC2P_RIN_MC2T/MC2P_inputJ2RIN;
% Convert to HG10 mode amplitude/rt(Hz)
MC2P_HG10amp_IM4T=MC2P_RIN_IM4T/MC2P_inputJ2RIN/w0;
MC2P_HG10amp_MC2T=MC2P_RIN_MC2T/MC2P_inputJ2RIN/w0;

data=load('Commiss_MC13diffY.txt');
MC13diffYfreq=data(:,1);
MC13diffY_RIN_IM4T=data(:,2)/meanIM4T_MC13diffY;
MC13diffY_RIN_MC2T=data(:,3)/meanMC2T_MC13diffY;
% Convert to m/rt(Hz)
MC13diffY_JIT_IM4T=MC13diffY_RIN_IM4T/MC13diffY_inputJ2RIN;
MC13diffY_JIT_MC2T=MC13diffY_RIN_MC2T/MC13diffY_inputJ2RIN;
% Convert to HG10 mode amplitude/rt(Hz)
MC13diffY_HG10amp_IM4T=MC13diffY_RIN_IM4T/MC13diffY_inputJ2RIN/w0;
MC13diffY_HG10amp_MC2T=MC13diffY_RIN_MC2T/MC13diffY_inputJ2RIN/w0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate RMS vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MC2YRMSjitIM4T=zeros(1,length(MC2Yfreq));
MC2YRMSjitMC2T=zeros(1,length(MC2Yfreq));
MC13diffPRMSjitIM4T=zeros(1,length(MC13diffPfreq));
MC13diffPRMSjitMC2T=zeros(1,length(MC13diffPfreq));
MC2PRMSjitIM4T=zeros(1,length(MC2Pfreq));
MC2PRMSjitMC2T=zeros(1,length(MC2Pfreq));
MC13diffYRMSjitIM4T=zeros(1,length(MC13diffYfreq));
MC13diffYRMSjitMC2T=zeros(1,length(MC13diffYfreq));

fstepMC2Y=MC2Yfreq(2)-MC2Yfreq(1);
fstepMC13diffP=MC13diffPfreq(2)-MC13diffPfreq(1);
fstepMC2P=MC2Pfreq(2)-MC2Pfreq(1);
fstepMC13diffY=MC13diffYfreq(2)-MC13diffYfreq(1);


for j=1:length(MC2YRMSjitIM4T)
    MC2YRMSjitIM4T(j)=sqrt(sum((MC2Y_JIT_IM4T(j:length(MC2Yfreq))).^2)*fstepMC2Y);
    MC2YRMSjitMC2T(j)=sqrt(sum((MC2Y_JIT_MC2T(j:length(MC2Yfreq))).^2)*fstepMC2Y);
    MC13diffPRMSjitIM4T(j)=sqrt(sum((MC13diffP_JIT_IM4T(j:length(MC13diffPfreq))).^2)*fstepMC13diffP);
    MC13diffPRMSjitMC2T(j)=sqrt(sum((MC13diffP_JIT_MC2T(j:length(MC13diffPfreq))).^2)*fstepMC13diffP);
    MC2PRMSjitIM4T(j)=sqrt(sum((MC2P_JIT_IM4T(j:length(MC2Pfreq))).^2)*fstepMC2P);
    MC2PRMSjitMC2T(j)=sqrt(sum((MC2P_JIT_MC2T(j:length(MC2Pfreq))).^2)*fstepMC2P);
    MC13diffYRMSjitIM4T(j)=sqrt(sum((MC13diffY_JIT_IM4T(j:length(MC13diffYfreq))).^2)*fstepMC13diffY);
    MC13diffYRMSjitMC2T(j)=sqrt(sum((MC13diffY_JIT_MC2T(j:length(MC13diffYfreq))).^2)*fstepMC13diffY);
end

%%%%%%%%%%%%%%%%%%%%%%%%
% plot results
%%%%%%%%%%%%%%%%%%%%%%%%

% Define requirements for IMC input beam jitter from T0900142_V2
reqabove10Hz=2.5e-6*sqrt(1+(100./MC2Yfreq).^4);

% Plot results from IM4 trans sensor
figure2=figure()
loglog(MC2Yfreq(5:length(MC2Yfreq)),reqabove10Hz(5:length(MC2Yfreq)),'k--',...
    MC2Yfreq,MC2Y_HG10amp_IM4T,'y',MC2Yfreq,MC2YRMSjitIM4T/thetaIMC,'y--',...
    MC13diffPfreq,MC13diffP_HG10amp_IM4T,'g',MC13diffPfreq,MC13diffPRMSjitIM4T/thetaIMC,'g--',...
    MC13diffYfreq,MC13diffY_HG10amp_IM4T,'r',MC13diffYfreq,MC13diffYRMSjitIM4T/w0,'r--',...
    MC2Pfreq,MC2P_HG10amp_IM4T,'b',MC2Pfreq,MC2PRMSjitIM4T/w0,'b--');
legend('requirement','yaw tilt','yaw tilt RMS','pitch tilt','pitch tilt RMS',...
    'yaw shift','yaw shift RMS','pitch shift','pitch shift RMS','Location','SouthWest')
xlabel('Frequency [Hz]')
ylabel('Relative HG10 amplitude [1 / \surd Hz]')
title('IMC input beam jitter (measured from RIN at IM4 trans)')
axis tight
grid on
%plotalltitle('IMCJitterHG10mode_amp_IM4T')

% Plot results from MC2 trans sensor
figure3=figure()
loglog(MC2Yfreq(5:length(MC2Yfreq)),reqabove10Hz(5:length(MC2Yfreq)),'k--',...
    MC2Yfreq,MC2Y_HG10amp_MC2T,'y',MC2Yfreq,MC2YRMSjitMC2T/thetaIMC,'y--',...
    MC13diffPfreq,MC13diffP_HG10amp_MC2T,'g',MC13diffPfreq,MC13diffPRMSjitMC2T/thetaIMC,'g--',...
    MC13diffYfreq,MC13diffY_HG10amp_MC2T,'r',MC13diffYfreq,MC13diffYRMSjitMC2T/w0,'r--',...
    MC2Pfreq,MC2P_HG10amp_MC2T,'b',MC2Pfreq,MC2PRMSjitMC2T/w0,'b--');
legend('requirement','yaw tilt','yaw tilt RMS','pitch tilt','pitch tilt RMS',...
    'yaw shift','yaw shift RMS','pitch shift','pitch shift RMS','Location','SouthWest')
xlabel('Frequency [Hz]')
ylabel('Relative HG10 amplitude [1 / \surd Hz]')
%title('IMC input beam jitter (measured from RIN at MC2 trans)')
axis tight
grid on
plotalltitle('IMCJitterHG10mode_amp_MC2T')




