% Create an Optickle model of Advanced LIGO
%
% opt = optLIGO2
 
% lPRC = lPR + (lIX + lIY) / 2;
% lSRC = lSR + (lIX + lIY) / 2;
% lasy = lIX - lIY;

% dTune - the SRCL detuning distance
%
% The detuning phase used in the Bench is the round trip phase:
%   ifo.Optics.SRM.Tunephase = 2 * k * dTune
% where k = 2 * pi / lambda
%
% Buonanno & Chen gr-qc/0102012 use the one-way phase, given by:
%   Chen.Tunephase = pi / 2 - k * dTune
%  which is similar to the total offset on the SR set with:
%   opt = setPosOffset(opt, nSR, opt.lambda / 4 + dTune);
%  since
%   Chen.Tunephase = k * (lambda / 4 + dTune)
%
%  The lambda / 4 is there to get close to the resonance, and 
%  dTune to set the detuning frequency.  For the LIGO2 detector,
%  the detuning frequency is roughly given by:
%    detuneFreq =?= tan(k * (lambda / 4 + dTune)) * (-13 Hz)
%
%  (lPR - lSR) is set to resonate the upper 180 MHz SB.  Using
%  -dTune instead resonates the lower 180 MHz SB, which destroys
%  the detectors RSE behaviour.

function opt = optLIGO2

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Detector Parameters

  % detector geometry (distances in meters)
  Ri = 2076;		% radius of curvature of input mirrors (IX and IY)
  Re = 2076;		% radius of curvature of end mirrors (EX and EY)
  Rr = 1442;		% radius of curvature of recycling mirrors (PR and SR)

  lPR = 3.4455;		% distance from PR to BS
  lSR = 3.498;		% distance from SR to BS
  dTune = 1e-8;		% SR micro detuning (see note)
  
  lIX = 5.106;		% distance from BS to IX
  lIY = 4.63;		% distance from BS to IY

  lEX = 3995.00;	% length of the X arm
  lEY = 3995.01;        % length of the Y arm
  
  % mechanical parameters
  wI = 2 * pi * 0.45;		% resonance frequency of input mirror (rad/s)
  mI = 39.9;			% mass of input mirror (kg)
  wE = 2 * pi * 0.45;		% resonance frequency of end mirror (rad/s)
  mE = mI;			% mass of end mirror (kg)

  % input beam parameters
  Pin = 125;			% input power (W)
  f1 = 9e6;			% first modulation frequency
  f2 = 180e6;			% second modulation frequency
  g1 = 0.1;			% first modulation depth
  g2 = 0.1;			% second modulation depth
  Nmod1 = 1;			% first modulation order
  Nmod2 = 1;			% second modulation order

  
  % Arm cavity Finesse and imbalance
  Titm = 0.005;
  dTitm = Titm/100;
  
  % the loss on EX and EY is used to set the contrast defect
  % 
  dLoss = 15e-6;


  % construct modulation vectors (Mach Z)
  n1 = (-Nmod1:Nmod1)';
  n2 = (-Nmod2:Nmod2)';
  N1 = length(n1);
  N2 = length(n2);

  vMod1 = zeros(N1, 2);        % f1 sideband vector
  for n = 1:N1
    ff = n1(n) * f1;
    aa = sqrt(Pin/2) * besselj(n1(n),g1);
    vMod1(n,:) = [ff, aa/sqrt(2)];   % sqrt(2) for the MZ output mirror
  end
  
  vMod2 = zeros(N2, 2);        % f2 sideband vector
  for n = 1:N2
    ff = n2(n) * f2;
    aa = sqrt(Pin/2) * besselj(n2(n),g2);
    vMod2(n,:) = [ff, aa/sqrt(2)];   % sqrt(2) for the MZ output mirror
  end
  
  % add them together constructively
  vMod = vMod1;
  for n = 1:N2
    % if the N2 frequency is already in the list, add to the exisiting field
    if sum(ismember(vMod(:,1),vMod2(n,1)))
      gg = find(vMod(:,1) == vMod2(n,1));
      vMod(gg,2) = vMod(gg,2) + vMod2(n,2);  % constructive interference
    % else add the new frequency and field to the list of modulations
    else
      vMod = [vMod;vMod2(n,:)];
    end
  end
  
  % sort them according to frequency
  vMod = sortrows(vMod, 1);    
  
  % create an empty model
  opt = Optickle;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add a Field Source

  % make the modulation vectors and add a source
  [opt, nSrc] = addSource(opt, 'Source', vMod);

  % add an AM modulator
  [opt, nAM] = addModulator(opt, 'AM', 1);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add Mirrors
  % The argument list for addMirror is:
  % [opt, sn] = addMirror(opt, name, AoI, Chr, Thr, Lhr, Rar)
  % type "help addMirror" for more information
  
  % add recycling mirrors
  [opt, nPR] = addMirror(opt, 'PR', 0, 1 / Rr, 0.07, 10e-6);
  [opt, nSR] = addMirror(opt, 'SR', 0, 1 / Rr, 0.07, 10e-6);

  % add beam-splitter
  [opt, nBS] = addBeamSplitter(opt, 'BS', 45, 0, 0.5, 10e-6, 50e-6);

  % add arm mirrors (a.k.a. test-masses)
  [opt, nIX] = addMirror(opt, 'IX', 0, 1 / Ri, Titm + dTitm, 35e-6, 500e-6);
  [opt, nIY] = addMirror(opt, 'IY', 0, 1 / Ri, Titm - dTitm, 35e-6, 500e-6);
  [opt, nEX] = addMirror(opt, 'EX', 0, 1 / Re, 10e-6, 35e-6 + dLoss);
  [opt, nEY] = addMirror(opt, 'EY', 0, 1 / Re, 10e-6, 35e-6 - dLoss);

  % OMC mirrors
  [opt, nOMCa] = addMirror(opt, 'OMCa', 0, 1, pi/500);
  [opt, nOMCb] = addMirror(opt, 'OMCb', 0, 1, pi/500);

  % set mechanical transfer-functions for the test-masses (m/N)
  dampedRes = [-0.1 + 1i, -0.1 - 1i];

  opt = setMechTF(opt, nIX, zpk([], wI * dampedRes, 1 / mI));
  opt = setMechTF(opt, nEX, zpk([], wE * dampedRes, 1 / mE));
  opt = setMechTF(opt, nIY, zpk([], wI * dampedRes, 1 / mI));
  opt = setMechTF(opt, nEY, zpk([], wE * dampedRes, 1 / mE));

  % set position offsets for SR
  opt = setPosOffset(opt, nSR, opt.lambda / 4 + dTune);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add Field Sinks (needed by the probes)

  [opt, nREFL] = addSink(opt, 'REFL');
  [opt, nAS] = addSink(opt, 'AS');
  [opt, nOMC] = addSink(opt, 'OMC');
  [opt, nPOX] = addSink(opt, 'POX');
  [opt, nPOY] = addSink(opt, 'POY');
  [opt, nTRX] = addSink(opt, 'TRX');
  [opt, nTRY] = addSink(opt, 'TRY');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add Links

  % link input to PR
  opt = addLink(opt, nSrc, 'out', nAM, 'in', 1);
  opt = addLink(opt, nAM, 'out', nPR, 'bk', 10);

  % link BS A-side inputs to PR and SR front outputs
  opt = addLink(opt, nPR, 'fr', nBS, 'frA', lPR);
  opt = addLink(opt, nSR, 'fr', nBS, 'bkA', lSR);

  % link BS A-side outputs to and IX and IY back inputs
  opt = addLink(opt, nBS, 'frA', nIY, 'bk', lIY);
  opt = addLink(opt, nBS, 'bkA', nIX, 'bk', lIX);

  % link BS B-side inputs to and IX and IY back outputs
  opt = addLink(opt, nIY, 'bk', nBS, 'frB', lIY);
  opt = addLink(opt, nIX, 'bk', nBS, 'bkB', lIX);

  % link BS B-side outputs to PR and SR front inputs
  opt = addLink(opt, nBS, 'frB', nPR, 'fr', lPR);
  opt = addLink(opt, nBS, 'bkB', nSR, 'fr', lSR);

  % link the arms
  opt = addLink(opt, nIX, 'fr', nEX, 'fr', lEX);
  opt = addLink(opt, nEX, 'fr', nIX, 'fr', lEX);

  opt = addLink(opt, nIY, 'fr', nEY, 'fr', lEY);
  opt = addLink(opt, nEY, 'fr', nIY, 'fr', lEY);

  % link the OMC
  opt = addLink(opt, nSR, 'bk', nOMCa, 'bk', 5);
  opt = addLink(opt, nOMCa, 'fr', nOMCb, 'fr', 0.22);
  opt = addLink(opt, nOMCb, 'fr', nOMCa, 'fr', 0.22);

  % link the detectors
  opt = addLink(opt, nPR,'bk', nREFL, 'in', 5);        % REFL
  opt = addLink(opt, nOMCb,'bk', nAS, 'in', 5);        % After the OMC
  opt = addLink(opt, nOMCa,'bk', nOMC, 'in', 5);       % OMC Reflection 
  opt = addLink(opt, nIX, 'po', nPOX, 'in', 5);
  opt = addLink(opt, nIY, 'po', nPOY, 'in', 5);
  opt = addLink(opt, nEX, 'bk', nTRX, 'in', 5);
  opt = addLink(opt, nEY, 'bk', nTRY, 'in', 5);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add Probes

  % double-demod sum (plus) and difference (minus) frequencies
  fp = f2 + f1;
  fm = f2 - f1;
  
  % REFL signals (reflected or symmetric port)
  REFL9_PHASE = 80;
  REFL171_PHASE = 92.2;
  REFL180_PHASE = 0;
  REFL189_PHASE = -100;
  %   REFL_I1 phase is set to maximize CARM
  %   REFL DD phases set for PRCL
  opt = addProbeAt(opt, 'REFL DC', nREFL, 'in',  0,   0);	% DC
  opt = addProbeAt(opt, 'REFL I1', nREFL, 'in', f1,0 + REFL9_PHASE);	% f1 demod I
  opt = addProbeAt(opt, 'REFL Q1', nREFL, 'in', f1,90 +REFL9_PHASE);	% f1 demod Q
  opt = addProbeAt(opt, 'REFL IM', nREFL, 'in', fm,0 + REFL171_PHASE);	% fm demod I
  opt = addProbeAt(opt, 'REFL QM', nREFL, 'in', fm,90 +REFL171_PHASE);	% fm demod Q
  opt = addProbeAt(opt, 'REFL IP', nREFL, 'in', fp,0 + REFL189_PHASE);	% fp demod I
  opt = addProbeAt(opt, 'REFL QP', nREFL, 'in', fp,90 +REFL189_PHASE);	% fp demod Q
  opt = addProbeAt(opt, 'REFL I2', nREFL, 'in', f2,0 + REFL180_PHASE);	% f2 demod I
  opt = addProbeAt(opt, 'REFL Q2', nREFL, 'in', f2,90 +REFL180_PHASE);	% f2 demod Q

  % AS signals (anti-symmetric port)
  %   AS_Q2 phase is set to maximize DARM above the optical spring
  %   AS_I2 phase is set 90 degrees away from AS_Q2
  %   AS DD phases set for MICH
  AS9_PHASE = 7;
  AS171_PHASE = 11.65;
  AS180_PHASE = 0;
  AS189_PHASE = 24.24;
  opt = addProbeAt(opt, 'AS DC', nAS, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'AS I1', nAS, 'in', f1,AS9_PHASE + 0);		% f1 demod I
  opt = addProbeAt(opt, 'AS Q1', nAS, 'in', f1,AS9_PHASE + 90);		% f1 demod Q
  opt = addProbeAt(opt, 'AS IM', nAS, 'in', fm,AS171_PHASE + 0);	% fm demod I
  opt = addProbeAt(opt, 'AS QM', nAS, 'in', fm,AS171_PHASE + 90);	% fm demod Q
  opt = addProbeAt(opt, 'AS IP', nAS, 'in', fp,AS189_PHASE + 0);	% fp demod I
  opt = addProbeAt(opt, 'AS QP', nAS, 'in', fp,AS189_PHASE + 90);	% fp demod Q
  opt = addProbeAt(opt, 'AS I2', nAS, 'in', f2,AS180_PHASE + 0);	% f2 demod Q
  opt = addProbeAt(opt, 'AS Q2', nAS, 'in', f2,AS180_PHASE + 90);	% f2 demod Q

  % POX signals (IX pick-off)
  %   POX DD phases set for SRC
  POX9_PHASE = 37;
  POX171_PHASE = -32.4;
  POX180_PHASE = 0;
  POX189_PHASE = -45.2;
  opt = addProbeAt(opt, 'POX DC', nPOX, 'in',  0,   0);  		% DC
  opt = addProbeAt(opt, 'POX I1', nPOX, 'in', f1,POX9_PHASE+0);	        % f1 demod I
  opt = addProbeAt(opt, 'POX Q1', nPOX, 'in', f1,POX9_PHASE+90);	% f1 demod Q
  opt = addProbeAt(opt, 'POX IM', nPOX, 'in', fm,POX171_PHASE+0);	% fm demod I
  opt = addProbeAt(opt, 'POX QM', nPOX, 'in', fm,POX171_PHASE+90);	% fm demod Q
  opt = addProbeAt(opt, 'POX IP', nPOX, 'in', fp,POX189_PHASE+0);	% fp demod I
  opt = addProbeAt(opt, 'POX QP', nPOX, 'in', fp,POX189_PHASE+90);	% fp demod Q
  opt = addProbeAt(opt, 'POX IP', nPOX, 'in', f2,POX180_PHASE+0);	% fp demod I
  opt = addProbeAt(opt, 'POX QP', nPOX, 'in', f2,POX180_PHASE+90);	% fp demod Q

  % POY signals (IY pick-off)
  opt = addProbeAt(opt, 'POY DC', nPOY, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'POY I1', nPOY, 'in', f1,   0);		% f1 demod I
  opt = addProbeAt(opt, 'POY Q1', nPOY, 'in', f1,  90);		% f1 demod Q

  % transmitted DC signals
  opt = addProbeAt(opt, 'TRX DC', nTRX, 'in', 0, 0);		% DC
  opt = addProbeAt(opt, 'TRY DC', nTRY, 'in', 0, 0);		% DC

  % OMC reflected signals
  OMR9_PHASE = 37;
  OMR171_PHASE = -166.6;
  OMR180_PHASE = 0;
  OMR189_PHASE = 100;
  opt = addProbeAt(opt, 'OMR DC', nOMC, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'OMR I1', nOMC, 'in', f1,OMR9_PHASE+0);	        % f1 demod I
  opt = addProbeAt(opt, 'OMR Q1', nOMC, 'in', f1,OMR9_PHASE+90);	% f1 demod Q
  opt = addProbeAt(opt, 'OMR IM', nOMC, 'in', fm,OMR171_PHASE+0);	% fm demod I
  opt = addProbeAt(opt, 'OMR QM', nOMC, 'in', fm,OMR171_PHASE+90);	% fm demod Q
  opt = addProbeAt(opt, 'OMR IP', nOMC, 'in', fp,OMR189_PHASE+0);	% fp demod I
  opt = addProbeAt(opt, 'OMR QP', nOMC, 'in', fp,OMR189_PHASE+90);	% fp demod Q
  opt = addProbeAt(opt, 'OMR IP', nOMC, 'in', f2,OMR180_PHASE+0);	% fp demod I
  opt = addProbeAt(opt, 'OMR QP', nOMC, 'in', f2,OMR180_PHASE+90);	% fp demod Q

  % direct AS probe (at OMC input)
  ASp9_PHASE = 37;
  ASp171_PHASE = -166.6;
  ASp180_PHASE = 0;
  ASp189_PHASE = 100;
  opt = addProbeAt(opt, 'ASp DC', nOMCa, 'bk',  0,   0);	% DC
  opt = addProbeAt(opt, 'ASp I1', nOMCa, 'bk', f1,  80);	% f1 demod I
  opt = addProbeAt(opt, 'ASp Q1', nOMCa, 'bk', f1, 170);	% f1 demod Q
  opt = addProbeAt(opt, 'ASp I2', nOMCa, 'bk', f2,  -30);	% f1 demod I
  opt = addProbeAt(opt, 'ASp Q2', nOMCa, 'bk', f2,  60);	% f1 demod Q

  
