% Create an Optickle initial LIGO
%
% opt = iLIGO
 
% lPRC = lPR + (lIX + lIY) / 2;

function opt = iLIGO

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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.9655;		% distance from BS to IX
  lIY = 4.5975;		% distance from BS to IY

  lEX = 3995.14;	% length of the X arm
  lEY = 3995.15;        % length of the Y arm
  
  % mechanical parameters
  wIX = 2 * pi * 0.75;		% resonance frequency of input mirror (rad/s)
  wIY = wIX;
  wEX = wIX;
  wEY = wEX;
  wPR = wIX;
  wBS = wPR;
  
  % Make buh=1 for rad press, or big for none
  buh = 1e9;
  
  mEX = 10.5 * buh;			% mass of ETMX (kg)
  mEY = 10.5 * buh;			% mass of ETMY (kg)
  mIX = 10.5 * buh;			% mass of ITMX (kg)
  mIY = 10.5 * buh;			% mass of ITMY (kg)

  mPR = 10.3 * buh;
  mBS = 6 * buh;   

  % input beam parameters
  Pin = 6;			% 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;
  
  % T_ITMX - T_ITMY = 2*dTitm, 5e-4 gives delta_fc = 2%
  dTitm = -5e-4;

  % 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.027, 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.028+dTitm, 70e-6, 100e-6);
  [opt, nIY] = addMirror(opt, 'IY', 0, 1 / Ri, 0.028-dTitm, 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, name, angle, 1 / ROC, T)
  [opt, nOMCa] = addMirror(opt, 'OMCa', 0, 1, 1);
  [opt, nOMCb] = addMirror(opt, 'OMCb', 0, 1, 1);
  
  % set mechanical transfer-functions for the test-masses (m/N)
  dampedRes = [-0.1 + 1i, -0.1 - 1i];

  opt = setMechTF(opt, nIX, zpk([], wIX * dampedRes, 1 / mIX));
  opt = setMechTF(opt, nEX, zpk([], wEX * dampedRes, 1 / mEX));
  opt = setMechTF(opt, nIY, zpk([], wIY * dampedRes, 1 / mIY));
  opt = setMechTF(opt, nEY, zpk([], wEY * dampedRes, 1 / mEY));

  opt = setMechTF(opt, nPR, zpk([], wPR * dampedRes, 1 / mPR));
  opt = setMechTF(opt, nBS, zpk([], wBS * dampedRes, 1 / mBS));

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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
  REFL1_PHASE = 15.73;
  opt = addProbeAt(opt, 'REFL DC', nREFL, 'in',  0,   0);	% DC
  opt = addProbeAt(opt, 'REFL I1', nREFL, 'in', f1,REFL1_PHASE + 0);	% f1 demod I
  opt = addProbeAt(opt, 'REFL Q1', nREFL, 'in', f1,REFL1_PHASE +90);	% f1 demod Q

  % AS signals (anti-symmetric port, after OMC)
  %   AS_Q1 phase should be set to maximize DARM
  AS_PHASE = -134.5708;
  opt = addProbeAt(opt, 'AS DC', nAS, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'AS I1', nAS, 'in', f1,AS_PHASE);		% f1 demod I
  opt = addProbeAt(opt, 'AS Q1', nAS, 'in', f1,AS_PHASE+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
  
