%% OPO transfer function

%% constants
L = 2*0.175; %OPO cavity length 
FSR = c/(L); %Hz c/round trip
Fin = pi*(.98*.999*.999*.999)^(1/5)/( 1-sqrt(.98*.999*.999*.999) ); %green Fin

OPOFWHM = FSR / Fin;

%% load OPO TF
OPOTF_amp = load('OPO_TFamp_Oct312019.TXT');
OPOTF_phase = load('OPO_TFph_Oct312019.TXT');
OPOTF_ampx = OPOTF_amp(:,1); %freq
OPOTF_ampy = OPOTF_amp(:,2); %amp
OPOTF_phy = OPOTF_phase(:,2); %phase

%% TF model

OPOTFamp = 10.^(OPOTF_ampy/20); %dB to power ratio
OPOcmplx = OPOTFamp.*exp(i.*deg2rad(OPOTF_phy)); 
OPOOLmag = abs(OPOcmplx);
OPOOLphase = angle(OPOcmplx);

%% 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);
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);
PZTcal = zpk([],[], FSR/40);  %Hz/V   
oneonf = zpk(-2*pi*[],-2*pi*[0], 2*pi ); %phase detection frequency actuation
       
% optional daughter board
pztnotch = zpk([-694.29+1j*41651.73, -694.29-1j*41651.73 -2*pi*265],...
            [-20828.759+1j*36076.469, -20828.759-1j*36076.469 -2*pi*66],1);
        

%% OPO stuff
        
OPOcommongain = 0.22; %-13dB
%what's the OPO input mW during this measurement? 3.57mW
OPOcal = 3.57*164e-3/0.4/OPOFWHM; %V/Hz 164mVpp measure at 0.4mW to OPO, 1mW launch. 
%OPOcalm = l*L/(2*c); %V/m        
OPOcavLP = zpk(-2*pi*[],-2.*pi.*[OPOFWHM./2], 2.*pi.*OPOFWHM/2 ); 
OPOcavHP = zpk(-2*pi*[0],-2.*pi.*[OPOFWHM./2], 1); 


OPOAfilter = PZTdriver.*PZTcal;       
OPOfilter = OPOcal.*OPOcommongain.*CM_boost1.*SLOW_boost.*pztnotch.*SLOW_filters.*OPOAfilter;
OPOCM = OPOcommongain.*CM_boost1.*SLOW_boost.*pztnotch.*SLOW_filters;
%OPOfilter = pztnotch;       
%Daneil's notch: complex zeroes at 6.63kHz, poles at 6.63 and 6.56 kHz,
%boost: pole 66Hz, zero 265Hz

%OPOTF_ampx_lf = logspace(0,7,1e5);

[OPOamp, OPOphase] = bode(OPOfilter, 2.*pi.*OPOTF_ampx);
OPOamp = squeeze(OPOamp);
OPOphase = squeeze(OPOphase);
OPOfiltercomplex = OPOamp.*exp(i.*deg2rad(OPOphase));
OPOOLmagzpk = abs(OPOfiltercomplex);
OPOOLphasezpk = angle(OPOfiltercomplex);

%% Plot OPO TF
figure(229);
subplot(2,1,1);
semilogx(OPOTF_ampx, 20*log10(OPOOLmag),...
    OPOTF_ampx, 20*log10(OPOOLmagzpk), '--', 'LineWidth', 1.5, 'Color', [0 0.4470 0.7410]);
xlabel('Frequency, Hz'), ylabel('Gain, dB');
ylim([-80,60])
xlim([1e2,1e5])
legend('Measured', 'Calculated')
title('OPO TF Oct 31, 2019');
grid on;
subplot(2,1,2);
semilogx(OPOTF_ampx, rad2deg(OPOOLphase),...
    OPOTF_ampx, rad2deg(OPOOLphasezpk),'--', 'LineWidth', 1.5, 'Color', [0 0.4470 0.7410]);
xlabel('Frequency, Hz'), ylabel('Phase, deg.');
xlim([1e2,1e5])
grid on;


