% Create an Optickle LIGO+ (a.k.a. Middle LIGO)
%
% opt = optLIGOp
 
% lPRC = lPR + (lIX + lIY) / 2;

function opt = optLIGOp

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

  % detector geometry (distances in meters)
  Rr = 14400;		% radius of curvature of recycling mirror (PR)
  Ri = 14300;		% radius of curvature of input mirrors (IX and IY)
  Re = 8000;		% radius of curvature of end mirrors (EX and EY)

  lPR = 4.397;		% distance from PR to BS
  
  lIX = 4.9765;		% distance from BS to IX
  lIY = 4.5975;		% distance from BS to IY

  lEX = 3995;		% length of the X arm
  lEY = 3995.15;        % length of the Y arm
  
  % mechanical parameters
  wI = 2 * pi * 0.75;		% resonance frequency of input mirror (rad/s)
  wE = 2 * pi * 0.76;		% resonance frequency of end mirror (rad/s)
  mEX = 10.5;			% mass of ETMX (kg)
  mEY = 10.5;			% mass of ETMY (kg)
  mIX = 10.5;			% mass of ITMX (kg)
  mIY = 10.5;			% mass of ITMY (kg)

  % input beam parameters
  Pin = 30;			% input power (W)
  f1 = 24485446;		% modulation frequency (24.5 MHz)
  f2 = 5/2*f1;			% second modulation frequency
  g1 = 0.3;			% first modulation depth
  g2 = 0.1;			% second modulation depth
  Nmod1 = 1;			% first modulation order
  Nmod2 = 0;			% second modulation order

  % the loss on EX and EY is used to set the contrast defect
  % 30.5 ppm gives CD = 5e-6
  dLoss = 30.5e-6;

  % construct modulation vectors
  n1 = (-Nmod1:Nmod1)';
  n2 = (-Nmod2:Nmod2)';
  N1 = length(n1);
  N2 = length(n2);
  vMod = zeros(N1 * N2, 2);
  for n = 1:N1
    nn = N2 * (n - 1) + (1:N2);
    ff = n1(n) * f1 + n2 * f2;
    aa = sqrt(Pin) * bessel(n1(n), g1) * bessel(n2, g2);
    vMod(nn, :) = [ff, aa];
  end
  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, aio, Chr, Thr, Lhr, Rar)
  % type "help addMirror" for more information
  
  % add recycling mirror
  [opt, nPR] = addMirror(opt, 'PR', 0, 1 / Rr, 0.03, 100e-6);

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

  % add arm mirrors (a.k.a. test-masses)
  [opt, nIX] = addMirror(opt, 'IX', 0, 1 / Ri, 0.0277, 70e-6, 100e-6);
  [opt, nIY] = addMirror(opt, 'IY', 0, 1 / Ri, 0.0283, 70e-6, 100e-6);
  [opt, nEX] = addMirror(opt, 'EX', 0, 1 / Re, 10e-6, 70e-6 - dLoss);
  [opt, nEY] = addMirror(opt, 'EY', 0, 1 / Re, 10e-6, 70e-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 / mIX));
  opt = setMechTF(opt, nEX, zpk([], wE * dampedRes, 1 / mEX));
  opt = setMechTF(opt, nIY, zpk([], wI * dampedRes, 1 / mIY));
  opt = setMechTF(opt, nEY, zpk([], wE * dampedRes, 1 / mEY));

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Add Field Sinks (needed by the probes)

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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 input to PR front output
  opt = addLink(opt, nPR, 'fr', nBS, 'frA', lPR);

  % 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 front input
  opt = addLink(opt, nBS, 'frB', nPR, 'fr', lPR);

  % 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, nBS, 'bkB', nOMCa, 'bk', 5);
  opt = addLink(opt, nOMCa, 'fr', nOMCb, 'fr', 0.22);
  opt = addLink(opt, nOMCb, 'fr', nOMCa, 'fr', 0.22);

  % link the sinks
  opt = addLink(opt, nPR, 'bk', nREFL, 'in', 5);
  opt = addLink(opt, nOMCb, 'bk', nAS, 'in', 5);
  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);
  opt = addLink(opt, nOMCa, 'bk', nOMC, 'in', 5);

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

  % REFL signals (reflected or symmetric port)
  %   REFL_I1 phase should be set to maximize CARM
  opt = addProbeAt(opt, 'REFL DC', nREFL, 'in',  0,   0);	% DC
  opt = addProbeAt(opt, 'REFL I1', nREFL, 'in', f1,   0);	% f1 demod I
  opt = addProbeAt(opt, 'REFL Q1', nREFL, 'in', f1,  90);	% f1 demod Q

  % AS signals (anti-symmetric port, after OMC)
  %   AS_Q1 phase should be set to maximize DARM
  opt = addProbeAt(opt, 'AS DC', nAS, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'AS I1', nAS, 'in', f1,   0);		% f1 demod I
  opt = addProbeAt(opt, 'AS Q1', nAS, 'in', f1,  90);		% f1 demod Q

  % POX signals (IX pick-off)
  opt = addProbeAt(opt, 'POX DC', nPOX, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'POX I1', nPOX, 'in', f1,   0);		% f1 demod I
  opt = addProbeAt(opt, 'POX Q1', nPOX, 'in', f1,  90);		% f1 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
  opt = addProbeAt(opt, 'OMC DC', nOMC, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'OMC I1', nOMC, 'in', f1,   0);		% f1 demod I
  opt = addProbeAt(opt, 'OMC Q1', nOMC, 'in', f1,  90);		% f1 demod Q

  % direct AS probe (at OMC input)
  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
  
