% o_O v8.3
% 
% Add Pcal line for calibration
%
% Quadrature + sideband picture with mode-mismatch
% Direct calculation of quantum noise
% 
% Wenxuan Jia
%
%
% Version history
% 
% v1 Direct calculation of quantum noise, using 4x4 matrix for DRFPMI
% 
% v2 Reformulate everything from KLMTV to quadrature + sideband picture
% 
% v3 Fixed dr/dphi(omega) for QRPN
% 
% v4 Package into a function. Add various quantum noise budget curves
% 
% v5 Add diagnostic homodyne curves
% 
% v6 Add mode-mismatch to homodyne
% 
% v7 Add mode-mismatch to everything
% 
% v7.2 Re-organize quantum modes to 2x1 vector
% 
% v7.3 Removed function handle for best performance
% 
% v7.4 Parallel computing enabled
% 
% v8 Add Pcal line for calibration, fixed bug on rl_CC. 
%
% v8.2 Add SQL lines, change _DC (double cavity) to _CC (coupled cavity)
%
% v8.3 Add contrast defect
%
% v9 Separate classical calculation from the main loop. PRC locking enabled



%#ok<*UNRCH> 
%#ok<*SEPEX> 
%#ok<*NASGU> 
%#ok<*AGROW>
%#ok<*PFOUS>
%#ok<*MINV> 


function results = o_O(varargin)
%% Read inputs
% 'param', {'theta'; 0; ...}
% 'omegas', 2*pi*logspace(1, 4, 100)
% 'traces', {'S_h'; 'S_rel'; ...}
% 'printDefaultParam', 1
% 'echo', 1

    echo = false;
    echo2 = false;
    SenseMon = false;
    parallel = false;
    printDefaultParam = false;
    
    for i = 1:2:length(varargin)
        temp = varargin{i+1};
        eval([varargin{i} '= temp;'])
    end
    
    clear i temp;
    


%% Calc classical + quantum fields in LIGO

    %% Default Parameters
    
    % theta is pre-defined angle
    % phi is propagation angle by classical parts, carrier. Affects detuning
    % varphi is propagation angle by quantum-mechanical parts, sideband. No
    % effect on detuning
    
    % constants
    c = 2.99792458e8; hbar = 1.0545718e-34; G = 6.67e-11;
    A2 = (1/sqrt(2))*[1, 1; -1i, 1i]; % we don't need A2 at all in this picture
    
    R2 = @(theta) [cos(theta), -sin(theta); sin(theta), cos(theta)]; % 2x2 rotation matrix
    I2 = eye(2,2);
    Z2 = zeros(2,2);
    
    
    
    %%%% HOM order
    HOMorder = 2; % only affects Gouy phase
    % phi = omega0*L/c - (n+m+1)*psi from Sigmann's Laser. HOM only adds
    % -n*psi to the single-trip phase phi_FC
    
    % single-trip rotation matrix
    R = @(theta, psi) [R2(theta), Z2; Z2, R2(theta-HOMorder*psi)]; % 4x4 rotation matrix for 00+HOM
    I = eye(4,4);
    Z = zeros(4,4);
    
    % Mode-mismatch scatter matrix
    U_mm = @(modemismatch, mismatchPsi) [I2*sqrt(1-modemismatch),              R2(-mismatchPsi)*sqrt(modemismatch); ...
                                         -R2(mismatchPsi)*sqrt(modemismatch),  I2*sqrt(1-modemismatch)];
    
                                     
    %%%%%%%%%%%%%%%%%%%%%%%%% input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Advanced LIGO paper: https://arxiv.org/pdf/1411.4547.pdf
    % O3 data: https://journals.aps.org/prd/pdf/10.1103/PhysRevD.102.062003
    % O4 SQZ design: https://dcc.ligo.org/DocDB/0156/T1800447/007/FC_DRD_T1800447-v7.pdf
    
    
    lambda0 = 1064e-9;
    
    
    %%%% Input mode cleaner
    P_IM4 = 39.8; % power [W]
    theta_LO = 90; % 0 for SR, 90 deg for RSE
    mismatch_IMC = 0; % between IMC to PRC
    mismatchPsi_IMC = 0; % extra mode-mismatch phase

    eta_IM4 = 0.97*0.977; % calibration
    
    
    %%%% ADF positive sidebands injected
%     A_ADF = 0; % amplitude
%     theta_ADF = pi/4; 
%     f_ADF = 1000; % [Hz] 
    
    
    %%%% Squeezer
    r_OPO1_r = sqrt(0.845);
%     r1_OPO_g = sqrt(0.7);
    loss_OPO = 0.0012; % gives 98% OPO efficiency
    T_OPO2_r = loss_OPO; 
    z_eff_dB = 0; % squeezing
    theta_SQZ = 0;
    
    loss_INJ = 0.046; % https://alog.ligo-la.caltech.edu/aLOG/index.php?callRep=62872
    
    
    %%%% Filter cavity
    L_FC = 297.85;
    L_FC_rms = 0.5e-12; % length rms fluctuation ~ pm
    deltaOmega_FC = 2*pi*25; % detuning frequency
    loss_FC = 60e-6;
    T_FC1 = 1000e-6; % 
%     T_FC1 = 1.2e-3; % gwinc A+
    T_FC1 = sqrt(deltaOmega_FC^2 + (loss_FC*c/4/L_FC)^2)*4*L_FC/c; % optimal detuning

    psi_FC = 49.7; % single-trip Gouy phase shift
    mismatch_FC = 0; % mode-matching to FC
    mismatchPsi_FC = 0; % extra mode-mismatch phase

    mismatch_FC2SR = 0; % mode-matching from FC to SRC
    mismatch_FC2OMC = 0; % mode-matching from FC to OMC
    mismatchPsi_FC2SR = 0; % extra mode-mismatch phase
    mismatchPsi_FC2OMC = 0; % extra mode-mismatch phase
    
    
    %%%% from FC to SRM
%     L_FC2SRM = 0; % path length from OPO -> FC -> SRM
%     phi_OPO2SRM = 0; % phase accumulated from propagation to OPO -> FC -> SRM
    
    
    %%%% Arm cavity
    m = 39.9; % [kg], test mass
    T_ITMX = 0.0148; T_ETMX = 4e-6; 
    T_ITMY = 0.0148; T_ETMY = 3.9e-6;
%     loss_arm = 81e-6; % set to have 200 kW at 35 W input
%     loss_arm = 74.5e-6; % set to have PRG = 47 
%     PRG = 43.4;
    loss_xarm = 77.4e-6;
    loss_yarm = 77.4e-6;
    deltaLoss_YMinusX = 0; 
    L_arm = 3994.5;
    DARM_offset = 10e-12;
%     deltaphi = 0; % Schnupp differential phase P = (1-cos(deltaphi))/2
    psi_xarm = 155.7; % single-trip Gouy phase
    psi_yarm = 155.7; 
    
    mismatch_XY = 0; % mode matching between X arm and Y arm (affect contrast)
    mismatchPsi_X = 0; % extra mismatch phase

    %%%% Michelson 
    L_MICH = 0; %(length absorbed to recycling cavity)
    
    H_BS = (1/sqrt(2))*[I I; I -I];
    
    % Schnupp asymmetry ignored
    phi_y = 0; % phi_y is the single-trip phase from BS to ITMY
    phi_x = 0; % degenerate with PRC and SRC
    
    psi_y = 0; % single-trip Gouy phase along y Michelson part (absorbed to recycling cavity)
    psi_x = 0; 
    
    
    %%%% Power-recycling
    loss_PR = 500e-6;
    T_PR = 0.031;
%     T_PR = 1-loss_PR;
    L_PR = 57.7;
    % single-trip extra phase from PRM to ITM
    % the phase from PRM to BS is assumed to be 0
    phi_PR = 0; 
    psi_PR = 25; % single-trip Gouy phase 
    mismatch_PR = 0; % mode-matching from PRC to X-arm (XY mismatch is handled at Y-arm)
    mismatchPsi_PR = 0; % extra mode-mismatch phase
    
    %%%% Resonant-sideband extraction cavity
    loss_SR = 500e-6;
    T_SR = 0.325; 
%     T_SR = 1-loss_SR;
    L_SR = 56.0;
    % single-trip extra phase from SRM to ITM
    % phase from SRM to BS is assumed to be 0
    % (deltaOmega>0 -> smaller wavelength -> shorter cavity to maintain the
    % same number of nodes in cavity -> less phase to travel)
    phi_SR = 90; % 0 for SR, 90 deg for RSE
    SRCL_detuning = 0; % in degree
    psi_SR = 19; % single-trip Gouy phase
    
    mismatch_SR = 0; % mode-matching from SRC to X-arm
    mismatchPsi_SR = 0; % extra mode-mismatch phase
    
    %%%% Output mode cleaner
    T_OMC1 = 8221e-6;
    T_OMC2 = 8221e-6;
    loss_OMCtrans = 0.025; % alog.ligo-la.caltech.edu/aLOG/index.php?callRep=61552
    L_OMC = 1.146/2; % single-trip cavity length
    
    phi_OMC = 0; % on resonance
    HOMspacing_OMC = 2*pi*57.6e6; % OMC002 https://dcc.ligo.org/DocDB/0117/T1500060/002/T1500060_OMC_Optical_Testing_Procedure.pdf
    psi_OMC = 180/pi* (HOMspacing_OMC*2*L_OMC/c/2);
    
    mismatch_OMC = 0;
    mismatchPsi_OMC = 0; % extra mode-mismatch phase
    
    %%%% Readout loss and PD quantum efficiency
    
    loss_RO = 0.01 + 800e-6 + 0.015 + 0.007; % OFI + OM1 + OM3 + OMC QPD pickoff
    loss_SQZUnknown = 0; % Unknown SQZ loss before OFI
    loss_QE = 0.02; % nominal, need reference?
    
    
    %%%% Balanced Homodyne
    theta_HD_LO = 90; % homodyne LO readout phase
    P_HD_LO = 1e-3; % homodyne LO power [W]
    r_HD_BS = (sqrt(1/2));
    t_HD_BS = sqrt(1-r_HD_BS^2);
    nu_HD = 1; % homodyne visibility (sqrt of mode-matching), not a BS-like loss
    mismatch_LO = 1-nu_HD^2;
    mismatchPsi_LO = 0; % extra mode-mismatch phase
    
    
    %%%% Photon calibrator Pcal
    omega_Pcal = 2*pi*434.9;
    
    
    if echo 
        fprintf('\n');
    end
    
    % print default arguments
    
    if printDefaultParam

        fprintf('Default Parameters: \n');
        defaultParam = who;
        for i = 1:length(defaultParam)
            temp = eval([defaultParam{i} ';']);
            try
                if isscalar(temp)
                    fprintf([defaultParam{i} ' = ']);
                    fprintf('%g; \n', temp);
                else
                    fprintf([defaultParam{i} ' = \n']);
                    disp(temp);
                end
            catch
                fprintf('Not displayable\n')
            end
        end
        if echo fprintf('\n\n'); end


    end



    

    
    
    
    
    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Read arguments and calculate other parameters.
    
    if exist('param', 'var')
        if echo fprintf('Input Parameters: \n'); end
        for i = 1:2:length(param) %#ok<USENS>
            temp = param{i+1};
            eval([param{i} '= temp;'])
            if echo fprintf([param{i} ' = ' num2str(param{i+1}) ';\n']); end
        end
        if echo fprintf('\n\n'); end
    end
    
    
    if echo fprintf('Derived Parameters: \n'); end
    
    omega0 = c*2*pi/lambda0;
    
    theta_LO = theta_LO*pi/180;
    phi_SR = (phi_SR + SRCL_detuning)*pi/180;
    theta_SQZ = theta_SQZ*pi/180; 
    theta_HD_LO = theta_HD_LO*pi/180;
    
    A_IM4 = sqrt(P_IM4*eta_IM4)*R(theta_LO, 0)*[1 0 0 0]'; % all in 00 mode of IMC basis

    mismatchPsi_IMC = pi/2 + sqrt(mismatch_IMC) + mismatchPsi_IMC*pi/180; 
    U_IMC = U_mm(mismatch_IMC, mismatchPsi_IMC); % convert IMC basis to PRC
%     omega_ADF = 2*pi*f_ADF;
    A_IM4 = U_IMC*A_IM4;
    

    
    
    %%%% OPO
    T_OPO1_r = 1-r_OPO1_r.^2; t1_OPO_r = sqrt(T_OPO1_r);
    t_OPO2_r = sqrt(T_OPO2_r);
    r_OPO2_r = sqrt(1-T_OPO2_r);
    finesse_OPO_r = pi*sqrt(r_OPO1_r.*r_OPO2_r)./(1-r_OPO1_r.*r_OPO2_r);
    
    z = abs(z_eff_dB)/(20*log10(exp(1)));
    expz1 = ((exp(z)+r_OPO1_r)./(r_OPO1_r.*exp(z)+1)); % single-pass antisqz
    
    eta_p = 1 - (2*r_OPO1_r.*(1-r_OPO2_r).*expz1 - r_OPO1_r.^2.*T_OPO2_r.*expz1.^2)./(1-r_OPO1_r.*expz1).^2;
    eta_m = 1 - (2*r_OPO1_r.*(1-r_OPO2_r)./expz1 - r_OPO1_r.^2.*T_OPO2_r./expz1.^2)./(1-r_OPO1_r./expz1).^2;
    if echo fprintf('OPO efficiency (asqz) = %g\n', eta_p); end
    if echo fprintf('OPO efficiency  (sqz) = %g\n', eta_m); end
    
    loss_INJ = loss_INJ + loss_SQZUnknown;
    
    % no frequency dependence if the detuning is ignored (not a function of omega)
    antisqz = (r_OPO1_r - (r_OPO2_r*t1_OPO_r^2.*expz1)./(1-r_OPO1_r*r_OPO2_r.*expz1)).^2 + ((t1_OPO_r*t_OPO2_r.*expz1)./(1-r_OPO1_r*r_OPO2_r.*expz1)).^2;
    sqz =  (r_OPO1_r - (r_OPO2_r*t1_OPO_r^2./expz1)./(1-r_OPO1_r*r_OPO2_r./expz1)).^2 + ((t1_OPO_r*t_OPO2_r./expz1)./(1-r_OPO1_r*r_OPO2_r./expz1)).^2;
    
    h_OPO = R2(theta_SQZ)*[sqrt(antisqz) 0; 0 sqrt(sqz)]*R2(theta_SQZ)'; % lossy
    
    % a_SQZ is 4x1, and <aa†> is [RS^2R†, 0; 0, 1] 
    % Treat a_SQZ = [RSR†, 0; 0, 1] = m_OPO_SQZ
    % and (m'*A)'*(m'*A) = A'*mm'*A = A'<aa†>A, equivalently
    m_OPO_SQZ = sqrt(1-loss_INJ)*[h_OPO, Z2; Z2, I2]; 
    
    m_qINJ = sqrt(loss_INJ)*I; % injection loss vacuum
    


    
    %%%% Filter cavity
%     T_FC1 = sqrt(deltaOmega_FC^2 + (loss_FC*c/4/L_FC)^2)*4*L_FC/c;
    r_FC1 = sqrt(1-T_FC1);
    T_FC2 = loss_FC;
    
    t_FC1 = sqrt(T_FC1);
    t_FC2 = sqrt(T_FC2);
    r_FC2 = sqrt(1-T_FC2);
    if T_FC1 == 0
        finesse_FC = 0;
    else
        finesse_FC = pi*sqrt(r_FC1.*r_FC2)./(1-r_FC1.*r_FC2);
    end
    
    phi_FC = -deltaOmega_FC*L_FC/c; % single-trip phase delay due to detuning
    psi_FC = psi_FC*pi/180;
    R1_FC = R(phi_FC, psi_FC); % single-trip phase in FC
    
    
    if echo fprintf('FC finesse = %g\n', finesse_FC); end
    
    mismatchPsi_FC = pi/2 + sqrt(mismatch_FC) + mismatchPsi_FC*pi/180; 
    U_FC = U_mm(mismatch_FC, mismatchPsi_FC); % convert OPO basis to FC 
    
    
    %%%% Arm cavity
    
%     loss_data = [65    66    67    68    69    70    71    72    73    74    75    76    77    78    79    80    81    82    83    84    85]*1e-6;
%     PRG_data =  [51.9692   51.4167   50.8731   50.3382   49.8116   49.2934   48.7833   48.2812   47.7868   47.3000   46.8207   46.3487   45.8839   45.4260   44.9751   44.5309   44.0933   43.6622   43.2374   42.8189   42.4065];
%     loss_arm = interp1(PRG_data, loss_data, PRG, 'linear', 'extrap');
     
    if deltaLoss_YMinusX ~= 0
        loss_yarm = loss_xarm + deltaLoss_YMinusX;
    end
    r_xarmloss = sqrt(1-loss_xarm); t_xarmloss = sqrt(loss_xarm);
    r_yarmloss = sqrt(1-loss_yarm); t_yarmloss = sqrt(loss_yarm);
    t_ITMX = sqrt(T_ITMX); r_ITMX = sqrt(1-T_ITMX);
    t_ETMX = sqrt(T_ETMX); r_ETMX = sqrt(1-T_ETMX);
    t_ITMY = sqrt(T_ITMY); r_ITMY = sqrt(1-T_ITMY);
    t_ETMY = sqrt(T_ETMY); r_ETMY = sqrt(1-T_ETMY);
    
    deltaL_xarm = DARM_offset/2; % DARM offset/2 (detuning fringe offset)
    deltaL_yarm = -deltaL_xarm;
    
    PRG = [];
    P_xarm = [];
    P_yarm = [];
    
    phi_xarm = omega0*deltaL_xarm/c;
    phi_yarm = omega0*deltaL_yarm/c;
    
    psi_xarm = psi_xarm*pi/180;
    psi_yarm = psi_yarm*pi/180;
    
    R1_xarm = R(phi_xarm, psi_xarm);
    R1_yarm = R(phi_yarm, psi_yarm);
    
    mismatch_X  = mismatch_XY/2;
    mismatchPsi_X = pi/2 + sqrt(mismatch_X) + mismatchPsi_X*pi/180; 
    U_X = U_mm(mismatch_X, mismatchPsi_X); % Convert to X arm basis

    mismatch_Y = mismatch_X;
    mismatchPsi_Y = pi - mismatchPsi_X;
    U_Y = U_mm(mismatch_Y, mismatchPsi_Y); % Convert to Y arm basis
    
    % R(pi/2) converts amplitude quadrature to phase quadrature
    RpiOver2 = [R2(pi/2), Z2;
                Z2, R2(pi/2)];

    % Classical part (quantum part below)
    g_xarm0 = t_ITMX*r_xarmloss*R1_xarm*exp(-1i*0*L_arm/c) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*0*L_arm/c));
    g_yarm0 = t_ITMY*r_yarmloss*R1_yarm*exp(-1i*0*L_arm/c) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*0*L_arm/c));
    r_xarm0 = (r_ITMX*I - r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-2i*0*L_arm/c)) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*0*L_arm/c));  
    r_yarm0 = (r_ITMY*I - r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-2i*0*L_arm/c)) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*0*L_arm/c)); 

    


    %%%% Michelson
    psi_y = psi_y*pi/180; % 0, because it's absorbed to recycling cavity
    psi_x = psi_x*pi/180;
    
    R1_MICHx = R(phi_x, psi_x);
    R1_MICHy = R(phi_y, psi_y);

    % Michelson for classical field, no optomechanical coupling
    % Classical part (quantum part below)
    h_MICH_NoOpto = H_BS*[R1_MICHy*U_Y'*r_yarm0*U_Y*R1_MICHy, Z;...
                    Z, R1_MICHx*U_X'*r_xarm0*U_X*R1_MICHx]*H_BS; % X arm and Y arm mode-mismatch handled here

    
    
    
    %%%% Double cavity
    t_PR = sqrt(T_PR);
    r_PR = sqrt(1-T_PR);
    t_SR = sqrt(T_SR);
    r_SR = sqrt(1-T_SR);
    t_PRloss = sqrt(loss_PR);
    t_SRloss = sqrt(loss_SR);
    
    psi_PR = psi_PR*pi/180;
    psi_SR = psi_SR*pi/180;
    
    mismatchPsi_PR = pi/2 + sqrt(mismatch_PR) + mismatchPsi_PR*pi/180;
    mismatchPsi_SR = pi/2 + sqrt(mismatch_SR) + mismatchPsi_SR*pi/180;
    U_PR = U_mm(mismatch_PR, mismatchPsi_PR); % convert PRC basis to X, Y arm on average
    U_SR = U_mm(mismatch_SR, mismatchPsi_SR); % convert SRC basis to X, Y arm on average

    r_CC = [r_PR*I     Z;...
            Z   r_SR*I];

    t_CC = [t_PR*I       Z;...
            Z    t_SR*I];

    tl_CC = [t_PRloss*I       Z;...
            Z      t_SRloss*I];

    rl_CC = [sqrt(1-t_PRloss^2)*I       Z;...
            Z                  sqrt(1-t_SRloss^2)*I];

    % mode-matching from IMC to PRC
    U_CC = [U_PR, Z;
            Z, U_SR];

    % Classical part (quantum part below)
    R1_PR = R(phi_PR, psi_PR); % single-trip phase rotation
    R1_SR = R(phi_SR, psi_SR);

    e_CC0 = [R1_PR*exp(-1i*0*L_PR/c), Z;...
            Z, R1_SR*exp(-1i*0*L_SR/c)];

    % Field amplitude entering Y arm and X arm
    F_y = [I, Z]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)]; % replace last zeros(4,1) with field injected from AS and mode-matched to SRC
    F_x = [Z, I]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)];

    PRG = (F_y'*F_y + F_x'*F_x)/(A_IM4'*A_IM4);

    % Locking PRC if needed
%     if 0
    if phi_PR == 0 && mismatch_PR ~= 0 % if its default value, optimize phi_PR to lock PRC
%         fprintf('Locking PRC. Consider inputting a known phi_PR \n');
        err = Inf;
        perturbPhi_PR = 0.001; % gradient descent step size
        j = 0;

        PRG0 = 4*T_PR/(T_PR+loss_PR+mismatch_PR+(1-r_yarm0(1)^2))^2;

%         err_ = []; slope_ = []; delta_ = [];
        while abs(err) > 0.001
            oldPRG = PRG;
            
            phi_PR = phi_PR + perturbPhi_PR;
            R1_PR = R(phi_PR, psi_PR); % single-trip phase rotation
            R1_SR = R(phi_SR, psi_SR);
            e_CC0 = [R1_PR*exp(-1i*0*L_PR/c), Z;...
                    Z, R1_SR*exp(-1i*0*L_SR/c)];
            % Field amplitude entering Y arm and X arm
            F_y = [I, Z]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)]; % replace last zeros(4,1) with field injected from AS and mode-matched to SRC
            F_x = [Z, I]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)];
            perturbPRG = (F_y'*F_y + F_x'*F_x)/(A_IM4'*A_IM4); 
            slope = (perturbPRG - PRG)/perturbPhi_PR;

            deltaPhi_PR = slope/(2*PRG0*16*PRG0/(4*T_PR));
            phi_PR = phi_PR + deltaPhi_PR;

            R1_PR = R(phi_PR, psi_PR); % single-trip phase rotation
            R1_SR = R(phi_SR, psi_SR);
            e_CC0 = [R1_PR*exp(-1i*0*L_PR/c), Z;...
                    Z, R1_SR*exp(-1i*0*L_SR/c)];
            % Field amplitude entering Y arm and X arm
            F_y = [I, Z]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)]; % replace last zeros(4,1) with field injected from AS and mode-matched to SRC
            F_x = [Z, I]*H_BS*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC*[A_IM4; zeros(4,1)];
            PRG = (F_y'*F_y + F_x'*F_x)/(A_IM4'*A_IM4); 
            
            err = (PRG - oldPRG)/PRG;
%             err_ = [err_; err];
%             slope_ = [slope_; slope];
%             delta_ = [delta_; deltaPhi_PR];
            j = j+1;
        end

    end

    
    


    P_xarm = (g_xarm0*F_x)'*g_xarm0*F_x;
    P_yarm = (g_yarm0*F_y)'*g_yarm0*F_y;

    h_CC0 = r_CC + t_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0 / (eye(8,8)+r_CC*e_CC0*U_CC'*h_MICH_NoOpto*U_CC*rl_CC*e_CC0)*t_CC;
    

    % Classical field PRC REFL
    A_REFL = [I Z]*h_CC0*[A_IM4; zeros(4,1)]; % replace last zeros(4,1) with field injected from AS and mode-matched to SRC
    P_REFL = A_REFL'*A_REFL;
    % Classical field leaving SRC
    B = [Z I]*h_CC0*[A_IM4; zeros(4,1)]; % replace last zeros(4,1) with field injected from AS and mode-matched to SRC

    
    %%%% Output mode cleaner
    r_OMC1 = sqrt(1-T_OMC1);
    t_OMC1 = sqrt(T_OMC1);
    t_OMC2 = sqrt(T_OMC2);
    r_OMC2 = sqrt(1-T_OMC2);
    throughput_OMC = 1 - loss_OMCtrans;
    loss_OMCRTL = 1 - throughput_OMC/(t_OMC1*t_OMC2+sqrt(throughput_OMC)*r_OMC1*r_OMC2)^2;
    r_OMCloss = sqrt(1-loss_OMCRTL); t_OMCloss = sqrt(loss_OMCRTL);
    finesse_OMC = pi*sqrt(r_OMC1*r_OMC2*r_OMCloss)./(1-r_OMC1*r_OMC2*r_OMCloss);
    
    psi_OMC = psi_OMC*pi/180;
    R1_OMC = R(phi_OMC, psi_OMC);
    
    mismatchPsi_OMC = pi/2 + sqrt(mismatch_OMC) + mismatchPsi_OMC*pi/180; 
    U_OMC = U_mm(mismatch_OMC, mismatchPsi_OMC); % convert SRC basis to OMC
    
    P_DCPD = []; % DC readout classical field
    
    % FC to IFO modematching, if we know FC to OMC or not
    mismatchPsi_FC2SR = pi/2 + sqrt(mismatch_FC2SR) + mismatchPsi_FC2SR*pi/180; 
    mismatchPsi_FC2OMC = pi/2 + sqrt(mismatch_FC2OMC) + mismatchPsi_FC2OMC*pi/180; 
    if mismatch_FC2SR ~= 0
        U_FC2SR = U_mm(mismatch_FC2SR, mismatchPsi_FC2SR); % convert OPO basis to FC
    else
        U_FC2OMC = U_mm(mismatch_FC2OMC, mismatchPsi_FC2OMC); % convert OPO basis to OMC, single bounce
        U_FC2SR = U_FC2OMC*(U_OMC)'; 
    end

    
    t_OMC0 = t_OMC1*r_OMCloss*t_OMC2*R1_OMC*exp(-1i*0*L_OMC/c) / (I-r_OMC1*r_OMC2*r_OMCloss*R1_OMC*R1_OMC*exp(-1i*2*0*L_OMC/c));

    %%%% Readout PD quantum efficiency
    r_ROloss = sqrt(1-loss_RO);
    r_QEloss = sqrt(1-loss_QE);
%     r_SQZUnknownloss = sqrt(1-loss_SQZUnknown);

    B = r_QEloss*t_OMC0*U_OMC*r_ROloss*B; % OMC transmission

    P_DCPD = B'*B;
    
    
    
    
    %%%% Balanced homodyne detector
    A_HD_LO = sqrt(P_HD_LO)*R(theta_HD_LO, 0)*[1 0 0 0]'; % in LO basis
    
    h_HD_BS = [r_HD_BS*I t_HD_BS*I; ...
               t_HD_BS*I -r_HD_BS*I];
    R_HDLOvac = R(theta_HD_LO, 0);
    
    mismatchPsi_LO = pi/2 + sqrt(mismatch_LO) + mismatchPsi_LO*pi/180; 
    % Mode-mismatch to homodyne
    U_HD = U_mm(mismatch_LO, mismatchPsi_LO); % convert seed/SQZ to LO basis


    
    
    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                               % main loop %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
        

    if ~exist('traces', 'var')
        flag_traces = false;
        omegas = 1;
        results = [];
    else
        flag_traces = true;
        if ~iscell(traces) && strcmp(traces, 'SenseMon') %#ok<USENS>
            SenseMon = true;
            f_range = logspace(1, 4, 100);
            omegas = 2*pi*f_range(1:end-1);
        else
            if ~exist('omegas', 'var')
%                 omegas = logspace(log10(20), log10(6000), 100); % default frequency
                omegas = [];
            end
            
            results = cell(length(traces), 1);
        end
    end
    
    % initialize all traces
    v = zeros(1, length(omegas));
%     S = struct('PD', v, 'rel', v, 'rel_AS', v, 'rel_qAS', v, 'rel_qFC', v, 'rel_qSY', v, 'rel_qPR', v, 'rel_qSR', v, ...
%         'rel_qETMY', v, 'rel_qETMX', v, 'rel_qy', v, 'rel_qx', v, 'rel_qRO', v, 'rel_qOMC', v, 'rel_qQEta', v, ...
%         'h', v, 'h_AS', v, 'h_qAS', v, 'h_qFC', v, 'h_qSY', v, 'h_qPR', v, 'h_qSR', v, 'h_qETMY', v, 'h_qETMX', v, ...
%         'h_qy', v, 'h_qx', v, 'h_qRO', v, 'h_qOMC', v, 'h_qQEta', v, 'h_shot', v, ...
%         'HD_diff', v, 'rel_HD', v, 'rel_HD_SQZ', v, 'rel_HD_qINJ', v, 'rel_HD_qFC', v, 'rel_HD_LOsubtract', v );

    S_PD = v; S_rel = v; S_rel_AS = v; S_rel_qAS = v; S_rel_qFC = v; S_rel_qSY = v; S_rel_qPR = v; S_rel_qSR = v; 
    S_rel_qETMY = v; S_rel_qETMX = v; S_rel_qy = v; S_rel_qx = v; S_rel_qRO = v; S_rel_qOMC = v; S_rel_qQEta = v; 
    S_h = v; S_h_AS = v; S_h_qAS = v; S_h_qFC = v; S_h_qSY = v; S_h_qPR = v; S_h_qSR = v; S_h_qETMY = v; S_h_qETMX = v;
    S_h_qy = v; S_h_qx = v; S_h_qRO = v; S_h_qOMC = v; S_h_qQEta = v; S_h_shot = v; S_SQL = v;
    S_HD_diff = v; S_rel_HD = v; S_rel_HD_SQZ = v; S_rel_HD_qINJ = v; S_rel_HD_qFC = v; S_rel_HD_LOsubtract = v;
    S_rel_SQL = v; 
    opticalGain = v; g_Pcal = v; Pcal_m = v; f_CC = v; spare = v;
    
    cores = 0;
    if parallel
        cores = feature('numcores');
    end


parfor (id = 1:length(omegas), cores)
% for id = 1:length(omegas)
    omega = omegas(id);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%% Passing OPO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % this is the two-photon already, acting on the pump vacuum [a1, a2]
%     H_opo_cdf = R(theta_SQZ)*(t1_OPO_r*t2_OPO_r*S1/(I-r1_OPO_r*r2_OPO_r*S1))*R(theta_SQZ)';


%     antisqz = (r1_OPO_r^2 - 2*r1_OPO_r*r2_OPO_r*expz1 + r2_OPO_r^2*expz1^2 + t1_OPO_r^2*t2_OPO_r^2*expz1^2)/(1-r1_OPO_r*r2_OPO_r*expz1)^2;
%     sqz = (r1_OPO_r^2 - 2*r1_OPO_r*r2_OPO_r/expz1 + r2_OPO_r^2/expz1^2 + t1_OPO_r^2*t2_OPO_r^2/expz1^2)/(1-r1_OPO_r*r2_OPO_r/expz1)^2;





    % = H_opo_cdf*\vec{a}
    % OPO transfer function, ignoring detuning
%     
%     S1 = [expz1, 0; 0, 1/expz1];
%     
%     % no frequency-dependence if the detuning is ignored
%     H_ADF = A2'*R(theta_SQZ)*(t1_OPO_r*t2_OPO_r*S1/(I-r1_OPO_r*r2_OPO_r*S1))*R(theta_SQZ)'*A2;
%     
%     A_OPO_trans = H_ADF*[A_ADF*exp(1i*theta_ADF); 0];
%     A_ADFp = A_OPO_trans(1); % complex amplitude of negative ADF sideband generated by OPO
%     A_ADFn = A_OPO_trans(2)'; % Conjugate because the sideband is [a+, a_']
%   


    %%%%%%%%%%%%%%%%%% Passing through filter cavity %%%%%%%%%%%%%%%%%%%%%%



    % r_FC2 is too close to 1 and this gives about 2e-13 floating error.  
    
    r_FC = (r_FC1*I - r_FC2*R1_FC*R1_FC*exp(-1i*2*omega*L_FC/c)) / (I-r_FC1*r_FC2*R1_FC*R1_FC*exp(-1i*2*omega*L_FC/c));       
    t_FC = t_FC1*t_FC2*R1_FC*exp(-1i*omega*L_FC/c) / (I-r_FC1*r_FC2*R1_FC*R1_FC*exp(-1i*2*omega*L_FC/c));

    m_FC_SQZ = r_FC*U_FC*m_OPO_SQZ;

    m_FC_qINJ = r_FC*U_FC*m_qINJ;

    m_FC_qFC = t_FC*I;



%         
%     R_opo2srm = @(omega) R(phi_OPO2SRM + omega*L_FC2SRM/c);
%     R_FC2srm = @(omega) R(phi_OPO2SRM + omega*L_FC2SRM/c);

    % = R_FC2srm(omega)*H_FC(omega)*H_OPO_ADF*\vec{a}


%     A_ADFp = r_FC(omega_ADF)*A_ADFp;
%     A_ADFn = r_FC(-omega_ADF)*A_ADFn;
%     
%     A_ADFp = exp(-1i*phi_OPO2SRM)*A_ADFp;
%     A_ADFn = exp(-1i*phi_OPO2SRM)*A_ADFn;



    %%%%%%%%%%%%%%%%%%% Passing IFO, DC-readout scheme %%%%%%%%%%%%%%%%%%%%

    
    g_xarm = t_ITMX*r_xarmloss*R1_xarm*exp(-1i*omega*L_arm/c) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*omega*L_arm/c));
    
    g_yarm = t_ITMY*r_yarmloss*R1_yarm*exp(-1i*omega*L_arm/c) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*omega*L_arm/c));
    

    r_xarm = (r_ITMX*I - r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-2i*omega*L_arm/c)) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*omega*L_arm/c));  
    
    r_yarm = (r_ITMY*I - r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-2i*omega*L_arm/c)) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*omega*L_arm/c));    
    
    t_xarm = t_ITMX*t_ETMX*R1_xarm*exp(-1i*omega*L_arm/c) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*omega*L_arm/c));
    t_yarm = t_ITMY*t_ETMY*R1_yarm*exp(-1i*omega*L_arm/c) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*omega*L_arm/c));


    l_xarm = -t_ITMX*r_ETMX*t_xarmloss*R1_xarm*exp(-1i*omega*L_arm/c) / (I-r_ITMX*r_ETMX*r_xarmloss*R1_xarm*R1_xarm*exp(-1i*2*omega*L_arm/c));
    l_yarm = -t_ITMY*r_ETMY*t_yarmloss*R1_yarm*exp(-1i*omega*L_arm/c) / (I-r_ITMY*r_ETMY*r_yarmloss*R1_yarm*R1_yarm*exp(-1i*2*omega*L_arm/c));




    drxarm_dphi = -r_ETMX*r_xarmloss*(1+exp(1i*omega*L_arm/c))*RpiOver2*g_xarm0*g_xarm;
    dryarm_dphi = -r_ETMY*r_yarmloss*(1+exp(1i*omega*L_arm/c))*RpiOver2*g_yarm0*g_yarm;

%     drxarm_dphi = @(omega) -r_ETMX*r_xarmloss*(1+exp(1i*varphi(omega)))*( g_xarm(0)*g_xarm(omega)*R(pi/2) - 1i*(omega/omega0)*exp(-1i*varphi(omega))*R(2*phi_xarm)*t_ITMX^2*(I*exp(1i*2*omega*L_arm/c)-r_ITMX*r_ETMX*r_xarmloss*R(2*phi_xarm)*exp(-1i*2*omega*L_arm/c))^(-1) );
%     dryarm_dphi = @(omega) -r_ETMY*r_yarmloss*(1+exp(1i*varphi(omega)))*( g_yarm(0)*g_yarm(omega)*R(pi/2) - 1i*(omega/omega0)*exp(-1i*varphi(omega))*R(2*phi_yarm)*t_ITMY^2*(I*exp(1i*2*omega*L_arm/c)-r_ITMY*r_ETMY*r_yarmloss*R(2*phi_yarm)*exp(-1i*2*omega*L_arm/c))^(-1) );


    kappa_rxarm = -4*omega0/(m*c^2*omega^2)*(1+r_ETMX^2*exp(-1i*omega*L_arm/c))*drxarm_dphi*(F_x*F_x')*g_xarm0'*g_xarm;
    kappa_txarm = 4*omega0/(m*c^2*omega^2)*(t_ETMX/t_ITMX)*(r_ITMX*R1_xarm*exp(-1i*omega*L_arm/c)+(r_ETMX/r_xarmloss)*inv(R1_xarm))*drxarm_dphi*(F_x*F_x')*g_xarm0'*g_xarm; 
    kappa_lxarm = 4*omega0/(m*c^2*omega^2)*(exp(1i*omega*L_arm/c)+r_ETMX^2)*(t_xarmloss/(t_ITMX*r_xarmloss))*inv(R1_xarm)*drxarm_dphi*(F_x*F_x')*g_xarm0'*g_xarm;


    kappa_ryarm = -4*omega0/(m*c^2*omega^2)*(1+r_ETMY^2*exp(-1i*omega*L_arm/c))*dryarm_dphi*(F_y*F_y')*g_yarm0'*g_yarm;
    kappa_tyarm = 4*omega0/(m*c^2*omega^2)*(t_ETMY/t_ITMY)*(r_ITMY*R1_yarm*exp(-1i*omega*L_arm/c)+(r_ETMY/r_yarmloss)*inv(R1_yarm))*dryarm_dphi*(F_y*F_y')*g_yarm0'*g_yarm;
    kappa_lyarm = 4*omega0/(m*c^2*omega^2)*(exp(1i*omega*L_arm/c)+r_ETMY^2)*(t_yarmloss/(t_ITMY*r_yarmloss))*inv(R1_yarm)*dryarm_dphi*(F_y*F_y')*g_yarm0'*g_yarm;


    R_optox = r_xarm + kappa_rxarm;
    R_optoy = r_yarm + kappa_ryarm;


    T_optox = t_xarm + kappa_txarm;
    T_optoy = t_yarm + kappa_tyarm;


    L_optox = l_xarm + kappa_lxarm;
    L_optoy = l_yarm + kappa_lyarm;

    %%%% Michelson

    % Michelson for quantum field, including optomechanical coupling
    h_MICH = H_BS*[R1_MICHy*U_Y'*R_optoy*U_Y*R1_MICHy*exp(-1i*2*omega*L_MICH/c) Z;...
        Z R1_MICHx*U_X'*R_optox*U_X*R1_MICHx*exp(-1i*2*omega*L_MICH/c)]*H_BS; % X arm and Y arm mode-mismatch handled here

    l_MICHq = H_BS*[R1_MICHy*T_optoy*exp(-1i*omega*L_MICH/c) Z;...
        Z R1_MICHx*T_optox*exp(-1i*omega*L_MICH/c)];

    l_MICHp = H_BS*[R1_MICHy*L_optoy*exp(-1i*omega*L_MICH/c) Z;...
        Z R1_MICHx*L_optox*exp(-1i*omega*L_MICH/c)];

    l_MICHh = H_BS*[R1_MICHy*exp(-1i*omega*L_MICH/c) Z;...
        Z R1_MICHx*exp(-1i*omega*L_MICH/c)]; % GW strain


    %%%% double cavity, aka recycling cavity

    e_CC = [R1_PR*exp(-1i*omega*L_PR/c), Z;...
            Z, R1_SR*exp(-1i*omega*L_SR/c)];


    h_CC  = r_CC + t_CC*e_CC*U_CC'*h_MICH*U_CC*rl_CC*e_CC / (eye(8,8)+r_CC*e_CC*U_CC'*h_MICH*U_CC*rl_CC*e_CC)*t_CC;
    

    l_CCq = t_CC*e_CC*( eye(8,8)-U_CC'*h_MICH*U_CC*rl_CC*e_CC / (eye(8,8)+r_CC*e_CC*U_CC'*h_MICH*U_CC*rl_CC*e_CC)*r_CC*e_CC )*U_CC'*h_MICH*U_CC*tl_CC;

    l_CCm = t_CC*e_CC*( eye(8,8)-U_CC'*h_MICH*U_CC*rl_CC*e_CC / (eye(8,8)+r_CC*e_CC*U_CC'*h_MICH*U_CC*rl_CC*e_CC)*r_CC*e_CC );


    %%%% OMC
    
    t_OMC = t_OMC1*r_OMCloss*t_OMC2*R1_OMC*exp(-1i*omega*L_OMC/c) / (I-r_OMC1*r_OMC2*r_OMCloss*R1_OMC*R1_OMC*exp(-1i*2*omega*L_OMC/c));
    
    
    % OMC loss entered from lossy mirror
    tl_OMC = t_OMCloss*t_OMC2*R1_OMC*exp(-1i*omega*L_OMC/c) / (I-r_OMC1*r_OMC2*r_OMCloss*R1_OMC*R1_OMC*exp(-1i*2*omega*L_OMC/c)); 
    % OMC loss entered from input mirror
    r_OMC = (r_OMC1*I - r_OMC2*r_OMCloss*R1_OMC*R1_OMC*exp(-2i*omega*L_OMC/c)) / (I-r_OMC1*r_OMC2*r_OMCloss*R1_OMC*R1_OMC*exp(-1i*2*omega*L_OMC/c)); 

    % OMC effective loss (floating error ~ 1e-14)
    l_OMC = sqrt((tl_OMC'*tl_OMC)+(r_OMC'*r_OMC));

    


    % calculate the output transfer matrices
    % It's painful to write each quantum operator one-by-one, because I
    % don't know how to zero all cross-spectral density values when
    % taking inner product of the sum of all of them
    % Making a giant matrix with diagonal elements for each m isn't great

    % SQZ vacuum injected from ASymmetric port vacuum at readout PD
    m_AS = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*h_CC*[Z; I]*U_FC2SR*m_FC_SQZ;

    % Unsqz vacuum injected from AS port
    m_qAS = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*h_CC*[Z; I]*U_FC2SR*m_FC_qINJ;

    % FC loss from AS port
    m_qFC = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*h_CC*[Z; I]*U_FC2SR*m_FC_qFC;

    % vacuum injected from SYmmetric port at readout PD
    m_qSY = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*h_CC*[I; Z];

    % vaccum leaked from power-recycling
    m_qPR = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCq*[I; Z];

    % vacuum leaked from signal-recycling
    m_qSR = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCq*[Z; I];

    % vaccum leaked from ETMY
    m_qETMY = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHq*[I; Z];

    % vaccum leaked from ETMX
    m_qETMX = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHq*[Z; I];

    % vacuum leaked from Y-arm round-trip loss
    m_qy = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHp*[I; Z];

    % vacuum leaked from X-arm round-trip loss
    m_qx = r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHp*[Z; I];

    % Readout loss
    m_qRO = r_QEloss*t_OMC*sqrt(loss_RO);

    % OMC loss
    m_qOMC = r_QEloss*l_OMC;

    % Quantum efficiency loss
    m_qQEta = sqrt(loss_QE);




    % Optical gain
    g_GWy = sqrt(2/(hbar*omega0))*omega0*L_arm/(2*c)*dryarm_dphi*F_y; % times h to give sideband
    g_GWx = -sqrt(2/(hbar*omega0))*omega0*L_arm/(2*c)*drxarm_dphi*F_x;

    g_opticalG = sqrt(2*hbar*omega0)*B'*r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHh*[g_GWy; g_GWx]; % [W/strain] h*g_cal = PD power
    opticalGain(id) = g_opticalG;

    % Pcal
    g_Pcaly = sqrt(2/(hbar*omega0))*(omega0/c)*dryarm_dphi*F_y; % times [m] to give sideband
    g_Pcal(id) = (sqrt(2*hbar*omega0)*B'*r_QEloss*t_OMC*U_OMC*r_ROloss*[Z I]*l_CCm*l_MICHh*[g_Pcaly; zeros(4,1)]); % [W/m] L*g_Pcal = PD power
    
    Pcal_m(id) = (B'*B)/(g_Pcal(id));

    f_CC(id) = -omega/2/pi*real(g_Pcal(id))/imag(g_Pcal(id));
    
    

    % PSD of the DC readout power P [W^2/Hz]
    S_PD(id) = 2*hbar*omega0*(...
        ((m_AS'*B)'*(m_AS'*B)) + ...
        ((m_qAS'*B)'*(m_qAS'*B)) + ...
        ((m_qFC'*B)'*(m_qFC'*B)) + ...
        ((m_qSY'*B)'*(m_qSY'*B)) + ...
        ((m_qPR'*B)'*(m_qPR'*B)) + ...
        ((m_qSR'*B)'*(m_qSR'*B)) + ...
        ((m_qETMY'*B)'*(m_qETMY'*B)) + ...
        ((m_qETMX'*B)'*(m_qETMX'*B)) + ...
        ((m_qy'*B)'*(m_qy'*B)) + ...
        ((m_qx'*B)'*(m_qx'*B)) +...
        ((m_qRO'*B)'*(m_qRO'*B)) +...
        ((m_qOMC'*B)'*(m_qOMC'*B)) +...
        ((m_qQEta'*B)'*(m_qQEta'*B)) ...
        ); 



    % PSD relative quantum noise (normalized by shot noise) [W^2/W^2]
    % Current PSD isn't accurate at low f due to DARM loop
    S_rel(id) = S_PD(id)/(2*hbar*omega0*(B'*B));

    % Breakdown of relative quantum noise

    S_rel_AS(id) = 2*hbar*omega0*((m_AS'*B)'*(m_AS'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qAS(id) = 2*hbar*omega0*((m_qAS'*B)'*(m_qAS'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qFC(id) = 2*hbar*omega0*((m_qFC'*B)'*(m_qFC'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qSY(id) = 2*hbar*omega0*((m_qSY'*B)'*(m_qSY'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qPR(id) = 2*hbar*omega0*((m_qPR'*B)'*(m_qPR'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qSR(id) = 2*hbar*omega0*((m_qSR'*B)'*(m_qSR'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qETMY(id) = 2*hbar*omega0*((m_qETMY'*B)'*(m_qETMY'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qETMX(id) = 2*hbar*omega0*((m_qETMX'*B)'*(m_qETMX'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qy(id) = 2*hbar*omega0*((m_qy'*B)'*(m_qy'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qx(id) = 2*hbar*omega0*((m_qx'*B)'*(m_qx'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qRO(id) = 2*hbar*omega0*((m_qRO'*B)'*(m_qRO'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qOMC(id) = 2*hbar*omega0*((m_qOMC'*B)'*(m_qOMC'*B))/(2*hbar*omega0*(B'*B));
    S_rel_qQEta(id) = 2*hbar*omega0*((m_qQEta'*B)'*(m_qQEta'*B))/(2*hbar*omega0*(B'*B));


    % Unitarity test
%     omega = 10000*2*pi;
%     temp = [func_S_rel_AS;
%         func_S_rel_qAS;
%         func_S_rel_qFC;
%         func_S_rel_qSY;
%         func_S_rel_qPR;
%         func_S_rel_qSR;
%         func_S_rel_qETMY;
%         func_S_rel_qETMX;
%         func_S_rel_qy;
%         func_S_rel_qx;
%         func_S_rel_qRO;
%         func_S_rel_qOMC;
%         func_S_rel_qQEta;
%         ]
% 
%     1-sum(temp)





    % PSD of GW signal in strain [1/Hz]
    S_h(id) = S_PD(id)./abs(g_opticalG)^2;
    %   added an extra parenthesis to make the result perfectly Hermitian 

    % Breakdown of strain noise
    S_h_AS(id) = 2*hbar*omega0*((m_AS'*B)'*(m_AS'*B))./abs(g_opticalG)^2;
    S_h_qAS(id) = 2*hbar*omega0*((m_qAS'*B)'*(m_qAS'*B))./abs(g_opticalG)^2;
    S_h_qFC(id) = 2*hbar*omega0*((m_qFC'*B)'*(m_qFC'*B))./abs(g_opticalG)^2;
    S_h_qSY(id) = 2*hbar*omega0*((m_qSY'*B)'*(m_qSY'*B))./abs(g_opticalG)^2;
    S_h_qPR(id) = 2*hbar*omega0*((m_qPR'*B)'*(m_qPR'*B))./abs(g_opticalG)^2;
    S_h_qSR(id) = 2*hbar*omega0*((m_qSR'*B)'*(m_qSR'*B))./abs(g_opticalG)^2;
    S_h_qETMY(id) = 2*hbar*omega0*((m_qETMY'*B)'*(m_qETMY'*B))./abs(g_opticalG)^2;
    S_h_qETMX(id) = 2*hbar*omega0*((m_qETMX'*B)'*(m_qETMX'*B))./abs(g_opticalG)^2;
    S_h_qy(id) = 2*hbar*omega0*((m_qy'*B)'*(m_qy'*B))./abs(g_opticalG)^2;
    S_h_qx(id) = 2*hbar*omega0*((m_qx'*B)'*(m_qx'*B))./abs(g_opticalG)^2;
    S_h_qRO(id) = 2*hbar*omega0*((m_qRO'*B)'*(m_qRO'*B))./abs(g_opticalG)^2;
    S_h_qOMC(id) = 2*hbar*omega0*((m_qOMC'*B)'*(m_qOMC'*B))./abs(g_opticalG)^2;     
    S_h_qQEta(id) = 2*hbar*omega0*((m_qQEta'*B)'*(m_qQEta'*B))./abs(g_opticalG)^2;

    % Shot noise in DARM
    S_h_shot(id) = 2*hbar*omega0*(B'*B)./abs(g_opticalG)^2;

    % SQL
    S_SQL(id) = 8*hbar./(m*omega.^2*L_arm^2);
    S_rel_SQL(id) = S_SQL(id).*abs(g_opticalG)^2/(2*hbar*omega0*(B'*B));

%     func_S_h_ADF = (func_S_ADF + func_S_ADF_BA)./(2*hbar*omega0*abs(A_0)^2*((B'*c_h)'*(B'*c_h)));
%     
%     func_S_h_BA = func_S_ADF_BA./(2*hbar*omega0*abs(A_0)^2*((B'*c_h)'*(B'*c_h)));




    

    A_seed_00 = 0*[1 0]'; % No DC light from OPO for balanced-homodyne readout
    A_seed_HOM = 0*[1 0]'; % it's always zero because why do you choose a seed basis where the most mode isn't zero? Well maybe other modes without mode-cleaner.
    A_HD_seed = U_HD*[A_seed_00; A_seed_HOM]; % HD mode is defined by LO beam mode

    % Classical field at homodyne PD 1 and PD2
    % [ PD1 ] = [r  t] [ SQZ ]
    % [ PD2 ] = [t -r] [ LO  ]
    % A_HD_PD1 - A_HD_PD2 gives about 2e-15 floating error even if the BS is
    % perfectly balanced 
    A_HD_PD1 = [I Z]*h_HD_BS*[A_HD_seed; A_HD_LO]; % other modes that's mismatch to LO are ignored, replace the zeros(2,1) with the classical field from squeezer (seed)
    A_HD_PD2 = [Z I]*h_HD_BS*[A_HD_seed; A_HD_LO];

    
    m_HD_PD1_SQZ = [I Z]*h_HD_BS*[I; Z]*U_HD*m_FC_SQZ; % the SQZ mode that's matched to LO
    m_HD_PD1_qINJ = [I Z]*h_HD_BS*[I; Z]*U_HD*m_FC_qINJ; % injection loss
    m_HD_PD1_qFC = [I Z]*h_HD_BS*[I; Z]*U_HD*m_FC_qFC; % FC loss
    m_HD_PD1_LO = [I Z]*h_HD_BS*[Z; I]*R_HDLOvac; % LO vacuum (aka symmetric vacuum for IFO)

    m_HD_PD2_SQZ = [Z I]*h_HD_BS*[I; Z]*U_HD*m_FC_SQZ; % the SQZ mode that's matched to LO
    m_HD_PD2_qINJ = [Z I]*h_HD_BS*[I; Z]*U_HD*m_FC_qINJ; % injection loss
    m_HD_PD2_qFC = [Z I]*h_HD_BS*[I; Z]*U_HD*m_FC_qFC; % FC loss
    m_HD_PD2_LO = [Z I]*h_HD_BS*[Z; I]*R_HDLOvac; % LO vacuum (aka symmetric vacuum for IFO)


    % PSD of homodyne detector [W^2/Hz]

    S_HD_diff(id) = 2*hbar*omega0*(...
        ((m_HD_PD1_SQZ'*A_HD_PD1 - m_HD_PD2_SQZ'*A_HD_PD2)'*(m_HD_PD1_SQZ'*A_HD_PD1 - m_HD_PD2_SQZ'*A_HD_PD2)) + ...
        ((m_HD_PD1_qINJ'*A_HD_PD1 - m_HD_PD2_qINJ'*A_HD_PD2)'*(m_HD_PD1_qINJ'*A_HD_PD1 - m_HD_PD2_qINJ'*A_HD_PD2)) + ...
        ((m_HD_PD1_qFC'*A_HD_PD1 - m_HD_PD2_qFC'*A_HD_PD2)'*(m_HD_PD1_qFC'*A_HD_PD1 - m_HD_PD2_qFC'*A_HD_PD2)) + ...
        ((m_HD_PD1_LO'*A_HD_PD1 - m_HD_PD2_LO'*A_HD_PD2)'*(m_HD_PD1_LO'*A_HD_PD1 - m_HD_PD2_LO'*A_HD_PD2)) ...
        );

    S_rel_HD(id) = S_HD_diff(id)/(2*hbar*omega0*(A_HD_PD1'*A_HD_PD1 + A_HD_PD2'*A_HD_PD2));

    S_rel_HD_SQZ(id) = ((m_HD_PD1_SQZ'*A_HD_PD1 - m_HD_PD2_SQZ'*A_HD_PD2)'*(m_HD_PD1_SQZ'*A_HD_PD1 - m_HD_PD2_SQZ'*A_HD_PD2))/...
        ((A_HD_PD1'*A_HD_PD1 + A_HD_PD2'*A_HD_PD2));
    S_rel_HD_qINJ(id) = ((m_HD_PD1_qINJ'*A_HD_PD1 - m_HD_PD2_qINJ'*A_HD_PD2)'*(m_HD_PD1_qINJ'*A_HD_PD1 - m_HD_PD2_qINJ'*A_HD_PD2))/...
        ((A_HD_PD1'*A_HD_PD1 + A_HD_PD2'*A_HD_PD2));
    S_rel_HD_qFC(id) = ((m_HD_PD1_qFC'*A_HD_PD1 - m_HD_PD2_qFC'*A_HD_PD2)'*(m_HD_PD1_qFC'*A_HD_PD1 - m_HD_PD2_qFC'*A_HD_PD2))/...
        ((A_HD_PD1'*A_HD_PD1 + A_HD_PD2'*A_HD_PD2));
    S_rel_HD_LOsubtract(id) = ((m_HD_PD1_LO'*A_HD_PD1 - m_HD_PD2_LO'*A_HD_PD2)'*(m_HD_PD1_LO'*A_HD_PD1 - m_HD_PD2_LO'*A_HD_PD2))/...
        ((A_HD_PD1'*A_HD_PD1 + A_HD_PD2'*A_HD_PD2));

%     % Unitarity test
%     omega = 1;
%     temp = [func_S_rel_HD_SQZ;
%     func_S_rel_HD_qINJ;
%     func_S_rel_HD_qFC;
%     func_S_rel_HD_LOsubtract]
% 
%     1-sum(temp)

    
    % collect results

    
end % main loop
        
    if echo2 fprintf('PRG = %g\n', mean(PRG)); end
    if echo2 fprintf('X arm power = %g kW \n', mean(P_xarm)/1e3); end
    if echo fprintf('Y arm power = %g kW \n', mean(P_yarm)/1e3); end 
    if echo fprintf('Optical gain = %g W/strain\n', mean(abs(opticalGain))); end % L1 specific
    if echo2 fprintf('IFO REFL = %g W \n', mean(P_REFL)); end
    P_DCPD_mA =  mean(P_DCPD)*0.86*1000;
    if echo2 fprintf('OMC DCPD = %g mA\n', P_DCPD_mA); end % L1 specific
    
    if ~flag_traces
        return;
    end
    
    if SenseMon
        % SenseMon range https://dcc.ligo.org/DocDB/0121/T1500491/002/rangeNote.pdf
        rho_0 = 8; % SNR ratio
        Msun = 2e30;
        Mpc = 1e6*3.0857e16;
        results = sqrt( 5*(G*1.2*Msun)^(5/3)/(96*c^3*pi^(4/3)*rho_0^2)* sum(diff(f_range).*f_range(1:end-1).^(-7/3)./S_h) )/Mpc;
        return;
    else
        if iscell(traces)
            for j = 1:length(traces)
                eval(['results{j} = ' traces{j} ';']); %#ok<EVLEQ> 
            end
        else
            eval(['results = ' traces ';']); %#ok<EVLEQ> 
        end
    end
    



end
