% Nutsinee Kijbunchoo Jan 2, 2018
% Intensity noise to freq noise projection


%% constants
c= 299792458; %speed of light
h = 6.626e-34; %J.s
l = 532e-9; %laser wavelength
%l = 1064e-9; %laser wavelength
fl = c/l; %laser frequency

n_ktp = 1.8868; % crystal index of refraction for green. 1.8296 for 1064. 
d_crys = 10e-3; %crystal length
L = 2*0.175; %OPO cavity length 
FSR = c/(L); %Hz c/round trip
%Fin = pi*(.875*.999*.999*.994)^(1/5)/( 1-sqrt(.875*.999*.999*.994) ); %Chua's thesis eq3.53
Fin = pi*(.98*.999*.999*.994)^(1/5)/( 1-sqrt(.98*.999*.999*.994) ); %green Fin
%Fin = 44
FWHM = FSR / Fin;
V2Hz = 40/FSR;
cav_pole = (FSR/Fin)/2;


%% load EOM TF data 
IMON_SHGlaunch_TF_amp = load('LOIMON_SHGlaunch_TF_Amp_Dec21.TXT');
IMON_SHGlaunch_TF_phase = load('LOIMON_SHGlaunch_TF_Phase_Dec21.TXT');

IMON_SHGlaunch_TF_ampx = IMON_SHGlaunch_TF_amp(:,1); %freq
IMON_SHGlaunch_TF_ampy = IMON_SHGlaunch_TF_amp(:,2); %amp
IMON_SHGlaunch_TF_phy = IMON_SHGlaunch_TF_phase(:,2); %phase

%% load SHG launch spec data 

GRPD_spec = load('SHGlaunch_Jan2.TXT');
GRPD_spec_x = GRPD_spec(:,1);
GRPD_spec_y = GRPD_spec(:,2);


%% make it plottable

IMON_SHGlaunch_TFamp = 10.^(IMON_SHGlaunch_TF_ampy/20); %dB to power ratio
IMON_SHGlaunch_cmplx = IMON_SHGlaunch_TFamp.*exp(i.*deg2rad(IMON_SHGlaunch_TF_phy)); 
IMON_SHGlaunch_OLmag = abs(IMON_SHGlaunch_cmplx);
IMON_SHGlaunch_OLphase = angle(IMON_SHGlaunch_cmplx);

%% LO TF model
LO_CMgain = 3.98; %14dB
LO_fastgain = 1.99; %6dB, factor of 2 going from red to green
IQdemod_cal = 3./2; %V/rad

CM_comp = zpk(-2.*pi.*[4e3  ], -2*pi.*[40 ], 1);
CM_boost1 = zpk(-2.*pi.*[10  ], -2*pi.*[0 ], 1);
CM_boost2 = zpk(-2.*pi.*[2e3  ], -2*pi.*[20 ], 1);
SLOW_boost = zpk(-2.*pi.*[400  ], -2*pi.*[4 ], 1);
SLOW_filters = zpk(-2.*pi.*[  ], -2*pi.*[100e3 ], 100e3*2*pi);
PZTdriver = zpk(-2.*pi.*[  10 ],...
            -2*pi.*[ 1 400], (1/10)*400*2*pi*20); %PZTdriver gain = 20 -- not included here, already in the calibration
%PFD = zpk(-2.*pi*[],-2.*pi*[0], 1.6*2.*pi); %1.6V/rad, 32deg/V, has to be measured in the lab
%PZT = zpk([],[], 10.7e6); %PZT calibration Hz/V for red, 21MHz/V for green
PZT = zpk([],[], 21e6); %PZT calibration Hz/V for red, 21MHz/V for green

PM = zpk(-2*pi*[60],-2*pi*[10], (1/60)*10*.0084); 
%flat = zpk([],[], 3.54*72e3/10); %flat from fast path. 11dB fast gain, 1/10 from AO
flat = zpk([],[], 5.9e3); %Hz/V calibration fast path
oneonf = zpk([],-2*pi*[0], 2*pi ); %phase detection frequency actuation
laserPZTpole = zpk([],-2*pi*[0], 2*pi );

common_path = LO_CMgain.*CM_comp.*CM_boost1; %no boost2 today
slow_path = SLOW_boost.*SLOW_filters.*PM.*PZTdriver.*PZT;
fast_path = LO_fastgain.*flat; 

%gain = 1./(0.233.*FSR); %LO rad per Hz


LO_TF = (1/IQdemod_cal).*common_path.* ( (fast_path) + (slow_path) ).*oneonf;


%LOTF_ampx_lf = logspace(0,6,1e5);

[LOamp, LOphase] = bode(LO_TF, 2*pi*IMON_SHGlaunch_TF_ampx); 
LOamp = squeeze(LOamp);
LOphase = squeeze(LOphase);
LOfilter = LOamp.*exp(i.*deg2rad(LOphase));
LOOLmagzpk = abs(LOfilter);
LOOLphasezpk = angle(LOfilter);

LOCLgain = abs(1./(1+LOOLmagzpk));

PDamp = zpk(-2*pi*[ ],-2*pi*[8e3 15e3 ], 1*8e3*15e3*4*pi^2); %have to multiple EOMcalib rad/V 
%PDampboard = zpk(-2*pi*[ ],-2*pi*[8e3 ], 1*8e3*2*pi); %have to multiple EOMcalib rad/V 


[PDamp, PDphase] = bode(PDamp, 2*pi*IMON_SHGlaunch_TF_ampx);
PDamp = squeeze(PDamp);
PDphase = squeeze(PDphase);
PDfilter = PDamp.*exp(i.*deg2rad(PDphase));
PDOLmagzpk = abs(PDfilter);
PDOLphasezpk = angle(PDfilter);

%%
figure(18);
subplot(2,1,1);
semilogx(IMON_SHGlaunch_TF_ampx, 20*log10(PDOLmagzpk),'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Gain, dB');
%xlim([1e2,1e6])
%legend('PD amp board (PDOLmagzpk)')
title('PD amp board (D1200543)');
grid on;
subplot(2,1,2);
semilogx(IMON_SHGlaunch_TF_ampx, rad2deg(PDOLphasezpk), 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Phase, deg.');
grid on;

%% plot both (COMMON)

figure(19);
subplot(2,1,1);
semilogx(IMON_SHGlaunch_TF_ampx, 20*log10(LOOLmagzpk),...
    IMON_SHGlaunch_TF_ampx, 20*log10(LOCLgain),'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Gain, dB');
%xlim([1e2,1e6])
legend('LO OLG model', '1/(1+OLG)')
title('Bode Plot LO OLG');
grid on;
subplot(2,1,2);
semilogx(IMON_SHGlaunch_TF_ampx, rad2deg(LOOLphasezpk),...
    IMON_SHGlaunch_TF_ampx, rad2deg(LOCLgain), 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Phase, deg.');
grid on;


%% plot suppressed, measured
% LO IMON/SHG launch

figure(20);
subplot(2,1,1);
semilogx(IMON_SHGlaunch_TF_ampx, 20*log10(IMON_SHGlaunch_OLmag), 'r','LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Gain, dB');
title('Bode Plot LO IMON/GR launch PD, raw data, Nov 1, 2018');
grid on;
subplot(2,1,2);
semilogx(IMON_SHGlaunch_TF_ampx, rad2deg(IMON_SHGlaunch_OLphase), 'r', 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Phase, deg.');
grid on;

%% plot unsuppressed, calibrated

%TTFSS_calib = 3.4/FWHM; %Vpk/FWHM
LO_calib = 3/2; %V/rad

GR_PD_calib = 0.192*2e3*10*(5.5/100)*1e-3; %V/W gain of 10 from PD amp toggle

figure(21);
loglog(IMON_SHGlaunch_TF_ampx, PDOLmagzpk.*IMON_SHGlaunch_OLmag.*GR_PD_calib.*(1./LO_calib)./(LOCLgain), 'r','LineWidth', 1.5);
%semilogx(EOMTF_ampx, 20*log10(EOMOLmag)./(20*log10(OPOCLgain)), 'r','LineWidth', 1.5);
hold on
interp_radpermW = interp1(IMON_SHGlaunch_TF_ampx, PDOLmagzpk.*IMON_SHGlaunch_TFamp.*GR_PD_calib.*(1./LO_calib)./(LOCLgain),GRPD_spec_x, 'pchip');

semilogx(GRPD_spec_x, interp_radpermW, 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('rad/mW');
xlim([1e1,1e5])
legend({'LO I-MON/GR launch PD', 'LO I-MON/GR launch PD interpolated'},'FontSize',18)
title('TTFSS/GR launch PD, Calibrated, Corrected for Close Loop gain, Nov 1, 2018');
grid on;
hold off


%% GR Launch PD response
% PD to SR785 CH2
% EXC to EOM and T-off to CH1

cal = 0.19*2e3*10*(5.5/100); %(V/A)*(A/W)*PDgain20dB*beam splitter ratio

%% plot GR launch in rad/rtHz

data = (GRPD_spec_y.*1e3./cal).*interp_radpermW;
% RMS
flipdiff = flip(diff(GRPD_spec_x));
flipPSD = flip(data.^2);
meansquare = cumsum(flipPSD(1:end-1).*flipdiff);
rms = flip(sqrt(meansquare));


figure(24);
loglog(GRPD_spec_x, (GRPD_spec_y.*1e3./cal).*interp_radpermW, 'r','LineWidth', 1.5);
%semilogx(EOMTF_ampx, 20*log10(EOMOLmag)./(20*log10(OPOCLgain)), 'r','LineWidth', 1.5);
hold on
loglog(GRPD_spec_x(1:end-1), rms, 'LineWidth' ,2 );
%semilogx(EOMTF_ampx, 20*log10(EOMOLmag)./(20.*log10(OPOCLgain)), 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('rad/\surdHz');
xlim([1e1,1e5])
legend({'Estimated frequency noise contribution to LO from GR launch PD', 'rms'},'FontSize',16)
title('GR launch PD Calibrated, Jan 2, 2018');
grid on;
hold off

%%
figure(25);
loglog(GRPD_spec_x, (GRPD_spec_y.*1e3./cal).*interp_radpermW.*GRPD_spec_x, 'r','LineWidth', 1.5);

%semilogx(EOMTF_ampx, 20*log10(EOMOLmag)./(20.*log10(OPOCLgain)), 'LineWidth', 1.5);
xlabel('Frequency, Hz'), ylabel('Hz/\surdHz');
xlim([1e1,1e5])
legend({'Estimated frequency noise contribution to LO from GR launch PD', 'rms'},'FontSize',16)
title('GR launch PD Calibrated, Jan 2, 2018');
grid on;
hold off
