%% IFO LO transfer function

%% Call all spec
run 'Marc2024ALLSQZSPEC.m'
%run CLFTF_highCLF_newVCO_March122024.m

%% load data

LOTF_amp = load('LO_highCLF_March132024_amp.TXT'); %
LOTF_phase = load('LO_highCLF_March132024_phase.TXT');
LOTF_ampx = LOTF_amp(:,1); %freq
LOTF_ampy = LOTF_amp(:,2); %amp
LOTF_phy = LOTF_phase(:,2); %phase

LOTFamp = 10.^(LOTF_ampy/20); %dB to power ratio
LOcmplx = LOTFamp.*exp(i.*deg2rad(LOTF_phy)); 
LOOLmag = abs(LOcmplx);
LOOLphase = angle(LOcmplx);

%% Generic common mode board components
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);
CM_boost_1kHz = zpk(-2.*pi.*[1e3  ], -2*pi.*[20 ], 1); %this doesn't really exist. Used to prject 1kHz LO BW. 2kHz boost has too much gain peaking
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*20); %PZTdriver gain = 20
%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.6e6); %Mephisto PZT calibration Hz/V
PM = zpk(-2*pi*[60],-2*pi*[10], (1/60)*10*.00136); 
%flat = zpk([],[], VCXOV2Hz); %Hz/V
oneonf = zpk([],-2*pi*[0], 2*pi ); %phase detection frequency actuation

VCOV2Hz = 1.1e6; %Hz/V calibration for CLF AOM
flat = zpk([],[], VCOV2Hz); %Hz/V factor of 2 because calibrated at 3MHz but used for 6MHz

%VCO = zpk(-2*pi*[ 40],-2*pi*[1.6], 2.*281479.*(1/40).*1.6); %Hz/V 80MHz VCO calibration at DC. Factor of 2 from freq doubler
%VCO = zpk(-2.*pi.*[40],-2*pi*[1.6], 1.6/40); 
VCO = zpk(-2*pi*[ 40],-2*pi*[1.6], 2.*281479.*(1/40).*1.6); %Hz/V 80MHz VCO calibration at DC. Factor of 2 from freq doubler

%% OMC cavity param alog 47540

OMCFSR = 261e6; 
OMCFin = 390;
OMCFWHM = OMCFSR/OMCFin;
OMCcav_pole = OMCFWHM/2;



%%

OMCcavLP = zpk(-2*pi*[],-2.*pi.*[OMCFWHM./2], 2.*pi.*OMCFWHM/2 );

OMCTF_ampx_lf = logspace(0,7,1e5);

[OMCcavLPamp, OMCcavLPphase] = bode(OMCcavLP, 2.*pi.*LOIMON_x);
%[OMCcavLPamp, OMCcavLPphase] = bode(OMCcavLP, 2.*pi.*OMCTF_ampx_lf);
OMCcavLPamp = squeeze(OMCcavLPamp);
OMCcavLPphase = squeeze(OMCcavLPphase);
OMCcavLPfiltercomplex = OMCcavLPamp.*exp(i.*deg2rad(OMCcavLPphase));
OMCcavLPOLmagzpk = abs(OMCcavLPfiltercomplex);
OMCcavLPOLphasezpk = angle(OMCcavLPfiltercomplex);

%% LO stuff

LO_CMgain = 0.14; %-17dB
LO_fastgain = 1; %0dB

LOcommon_path = LO_CMgain.*CM_boost1.*CM_boost2; 
LOfast_path = LO_fastgain.*VCO; 

%OMC RF3 -14.66dBm
LOcal = 7.48/2; %V/2rad -- Only /2rad for transfer function (slope = derivative of VpkSin(theta) where theta=0)
CM_FCboost = zpk(-2.*pi.*[30 30 30 20], -2*pi.*[0 0 0 0], 1);

fudgefactor = 1/3;
LOcommon_path_no2kHzboost = LO_CMgain.*CM_boost1.*CM_boost_1kHz;
LO_TF = LOcal.*LOcommon_path.* ( (LOfast_path) ).*oneonf.*fudgefactor; %Fudge factor of 1/3
%LO_TF = VCO;
LOcommon_path_2boost1 = CM_FCboost;

%LO_TF_1kHz = fudgefactor.*LOcal.*LOcommon_path_no2kHzboost.* ( (LOfast_path) ).*oneonf;
%FC_LO_TF_100Hz = fudgefactor2.*LOcommon_path_2boost1.*oneonf;

LOCMboard = LOcommon_path;
LOAct = LOfast_path.*oneonf; %Hz/V * rad/Hz

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

% [LO1kHzamp, LO1kHzphase] = bode(LO_TF_1kHz, 2*pi*LOIMON_x); 
% LO1kHzamp = squeeze(LO1kHzamp);
% LO1kHzphase = squeeze(LO1kHzphase);
% LO1kHzfilter = LO1kHzamp.*exp(i.*deg2rad(LO1kHzphase));
% LO1kHzOLmagzpk = abs(LO1kHzfilter);
% LO1kHzOLphasezpk = angle(LO1kHzfilter);
% 
% [LO1kCLGamp, LO1kCLGphase] = bode(LO_TF_1kHz/(1+LO_TF_1kHz), 2*pi*LOIMON_x); 
% LO1kCLGamp = squeeze(LO1kCLGamp);
% LO1kCLGphase = squeeze(LO1kCLGphase);
% LO1kCLGfilter = LO1kCLGamp.*exp(i.*deg2rad(LO1kCLGphase));
% LO1kCLGmagzpk = abs(LO1kCLGfilter);
% LO1kCLGphasezpk = angle(LO1kCLGfilter);



% [LO100CLGamp, LO100CLGphase] = bode(FC_LO_TF_100Hz/(1+FC_LO_TF_100Hz), 2*pi*LOIMON_x); 
% LO100CLGamp = squeeze(LO100CLGamp);
% LO100CLGphase = squeeze(LO100CLGphase);
% LO100CLGfilter = LO100CLGamp.*exp(i.*deg2rad(LO100CLGphase));
% LO100CLGmagzpk = abs(LO100CLGfilter);
% LO100CLGphasezpk = angle(LO100CLGfilter);

% [oneononeplusLO1kHzamp, oneononeplusLO1kHzphase] = bode(1/(1+LO_TF_1kHz), 2*pi*LOIMON_x); 
% oneononeplusLO1kHzamp = squeeze(oneononeplusLO1kHzamp);
% oneononeplusLO1kHzphase = squeeze(oneononeplusLO1kHzphase);
% oneononeplusLO1kHzfilter = oneononeplusLO1kHzamp.*exp(i.*deg2rad(oneononeplusLO1kHzphase));
% oneononeplusLO1kHzOLmagzpk = abs(oneononeplusLO1kHzfilter);
% oneononeplusLO1kHzOLphasezpk = angle(oneononeplusLO1kHzfilter);

[LOoneplusamp, LOoneplusphase] = bode(1+LO_TF, 2*pi*LOIMON_x); 
LOoneplusamp = squeeze(LOoneplusamp);
LOoneplusphase = squeeze(LOoneplusphase);
LOoneplusfilter = LOoneplusamp.*exp(i.*deg2rad(LOoneplusphase));
LOoneplusOLmagzpk = abs(LOoneplusfilter);
LOoneplusOLphasezpk = angle(LOoneplusfilter);


[LOoneononeplusamp, LOoneononeplusphase] = bode(1/(1+LO_TF), 2*pi*LOIMON_x); 
LOoneononeplusamp = squeeze(LOoneononeplusamp);
LOoneononeplusphase = squeeze(LOoneononeplusphase);
LOoneononeplusfilter = LOoneononeplusamp.*exp(i.*deg2rad(LOoneononeplusphase));
LOoneononeplusOLmagzpk = abs(LOoneononeplusfilter);
LOoneononeplusOLphasezpk = angle(LOoneononeplusfilter);


[LOGononeplusamp, LOGononeplusphase] = bode(LO_TF/(1+LO_TF), 2*pi*LOIMON_x); 
LOGononeplusamp = squeeze(LOGononeplusamp);
LOGononeplusphase = squeeze(LOGononeplusphase);
LOGononeplusfilter = LOGononeplusamp.*exp(i.*deg2rad(LOGononeplusphase));
LOGononeplusOLmagzpk = abs(LOGononeplusfilter);
LOGononeplusOLphasezpk = angle(LOGononeplusfilter);

[LOCMamp, LOCMphase] = bode(LOCMboard, 2*pi*LOIMON_x); 
LOCMamp = squeeze(LOCMamp);
LOCMphase = squeeze(LOCMphase);
LOCMfilter = LOCMamp.*exp(i.*deg2rad(LOCMphase));
LOCMOLmagzpk = abs(LOCMfilter);
LOCMOLphasezpk = angle(LOCMfilter);

[LOAamp, LOAphase] = bode(LOAct, 2*pi*LOIMON_x); 
LOAamp = squeeze(LOAamp);
LOAphase = squeeze(LOAphase);
LOAfilter = LOAamp.*exp(i.*deg2rad(LOAphase));
LOAOLmagzpk = abs(LOAfilter);
LOAOLphasezpk = angle(LOAfilter);



OMCLP = OMCcavLPOLmagzpk;

LOCM = LOCMOLmagzpk;
LOA = LOAOLmagzpk; 
LOG = LOOLmagzpk;
oneplusLOG = LOoneplusOLmagzpk;
oneononeplusLOG = LOoneononeplusOLmagzpk;
LO_CLG = oneononeplusLOG;
LO_RLG = LOGononeplusOLmagzpk;
LOS = LOcal; %V/rad
% LO_1kOLG = LO1kHzOLmagzpk;
% LO_1kCLG = oneononeplusLO1kHzOLmagzpk;
% LO_1kRLG_long = LO1kCLGmagzpk;
% LO_RLG_100Hz_long = LO100CLGmagzpk;


%% plot LOTF

figure(234);
subplot(2,1,1);
semilogx(LOTF_ampx, 20*log10(LOOLmag),...
    LOIMON_x, 20*log10(LOG),'--', 'LineWidth', 1.5,'Color', [0.4940 0.1840 0.5560]);
xlabel('Frequency, Hz'), ylabel('Gain, dB');
%ylim([-80,60])
xlim([1e2,1e4])
legend('Measured', 'Modeled', 'LO RLG 100Hz BW', 'LO RLG 4kHz BW')
title('LO TF 6dBm RF6 March 13, 2024');
grid on;
subplot(2,1,2);
semilogx(LOTF_ampx, rad2deg(LOOLphase),...
    LOIMON_x, rad2deg(LOOLphasezpk),'--', 'LineWidth', 1.5, 'Color', [0.4940 0.1840 0.5560]);
xlabel('Frequency, Hz'), ylabel('Phase, deg.');
xlim([1e2,1e4])
grid on;


%% OMC/DCPC sensing noise
OMCDN = 1.45e-6; %V
OMCSN = 1.75e-6; %how did I get this number?

OMCsens = sqrt(OMCSN^2+OMCDN^2)/1.3;
OMCsens = ones(size(LOIMON_x)).*OMCsens;

%% phi LO
%% aLIGO case
% OMCLO makes no different. Cavity pole too far out
LOIMON = LOIMON_y;
NLO = sqrt( (LOIMON./LO_CLG).^2 - OMCsens.^2)./(OMCLP.*LOS);
NLO_inloop = NLO.*LO_CLG.*OMCLP;
LOsens_inloop = OMCsens.*LO_RLG./LOS;
%LOsens_inloop1kHz = OMCsens.*LO_1kRLG./LOS./1.5;

%LOtotal = sqrt(NLO_inloop.^2+LOsens_inloop.^2);

%%
flipdiff1 = flip(diff(LOIMON_x));
flipPSD1 = flip(NLO_inloop.^2);
meansquare1 = cumsum(flipPSD1(1:end-1).*flipdiff1);
LOinlooprms1 = flip(sqrt(meansquare1));

flipdiff2 = flip(diff(LOIMON_x));
flipPSD2 = flip(OMCsens.^2);
meansquare2 = cumsum(flipPSD2(1:end-1).*flipdiff2);
LOsensrms2 = flip(sqrt(meansquare2));


flipdiff8 = flip(diff(LOIMON_x));
flipPSD8 = flip(NLO.^2);
meansquare8 = cumsum(flipPSD8(1:end-1).*flipdiff8);
NLOrms = flip(sqrt(meansquare8));


totalnewLOrms = sqrt(LOinlooprms1.^2+LOsensrms2.^2);

%% plot N_LO in LO loop no FC case (for thesis)

figure(2222)

plotNLO = loglog(LOIMON_x, 1e3.*NLO ,'LineWidth', 2, 'Color', [0.4940 0.1840 0.5560]);
plotNLO.Color(4) = 0.5

xlabel('Frequency, Hz'), ylabel('mrad/\surd Hz');
grid on;
hold on

%loglog(LO_IMON_short_x, 1e3.*LO_IMON_short_y./(2.74/2) ,'LineWidth', 2, 'Color', [0.6350 0.0780 0.1840]);
loglog(LOIMON_x, 1e3.*NLO_inloop ,'LineWidth', 2, 'Color', [0.4940 0.1840 0.5560]);
loglog(LOIMON_x, 1e3.*LOsens_inloop , '-.', 'LineWidth', 2, 'Color', [0.4940 0.1840 0.5560]);
%plotoldLOOOL.Color(4) = 0.5
%plot1ksensing = loglog(LOshort_x, 1e3.*LOsens_inloop1kHz , '-.', 'LineWidth', 2, 'Color', [0 0.4470 0.7410]);
%plot1ksensing.Color(4) = 0.5


loglog(LOIMON_x(1:end-1), 1e3.*totalnewLOrms,  '--' ,'LineWidth', 2, 'Color',[0.4940 0.1840 0.5560]);

%plotNLOrms = loglog(LOshort_x(1:end-1), 1e3.*NLOrms,  '--' ,'LineWidth', 2, 'Color',[0.4940 0.1840 0.5560]);
%plotNLOrms.Color(4) = 0.5




legend({'\deltaN_{LO}','\deltaN_{LO}*LO\_CLG*OMCLP',...
    '\deltaN_{sensLO}*LO\_RLG*(1/S)','\phi_{LO} rms'},'FontSize',14, 'Location', 'Northeast');

set(gca,'FontSize',18);
hTitle = title('LO Loop (6dBm RF6) March 13, 2024');
set(hTitle,'FontSize',16);

xlim([10,1e4])
%ylim([1e-5,1e4])
hold off