% Create an Optickle model of the 40m
%
% opt = opt40
 
% lPRC = lPR + (lIX + lIY) / 2;
% lSRC = lSR + (lIX + lIY) / 2;
% lasy = lIX - lIY;

% dTune - the SRC 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 166 MHz SB.  Using
%  -dTune instead resonates the lower 166 MHz SB, which destroys
%  the detectors RSE behaviour.

function opt = opt40

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

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

  lPR = 0.3032;		% distance from PR to BS
  lSR = 0.200;		% distance from SR to BS

  % SRC DeToon in degrees
  dTunePhase = -19.78; % SR micro detuning (see note)
  %dTune = dTunePhase;
  
  lIX = 2.183;		% distance from BS to IX
  lIY = 1.731;		% distance from BS to IY

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

  % input beam parameters
  Pin = 1;			% input power (W)
  f1 = 33.2e6;		% first modulation frequency
  f2 = 5 * f1;			% second modulation frequency
  g1 = 0.1;			% first modulation depth
  g2 = 0.1;			% second modulation depth
  Nmod1 = 2;			% first modulation order
  Nmod2 = 2;			% second modulation order

  % 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, aio, 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, 0.005, 1e-6, 25e-6);
  [opt, nIY] = addMirror(opt, 'IY', 0, 1 / Ri, 0.005, 1e-6, 45e-6);

  [opt, nEX] = addMirror(opt, 'EX', 0, 1 / Re, 1, 1e-6);  % T=1 for DRMI
  [opt, nEY] = addMirror(opt, 'EY', 0, 1 / Re, 1, 1e-6);  %

  % 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
  dTune = dTunePhase * (opt.lambda / 360);
  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, 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 detectors
  opt = addLink(opt, nPR, 'bk', nREFL, 'in', 5);
  opt = addLink(opt, nSR, '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);

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

  % double-demod sum (plus) and difference (minus) frequencies
  fp = f2 + f1;
  fm = f2 - f1;
  
  % REFL signals (reflected or symmetric port)
  REFL33_PHASE = 7;
  REFL133_PHASE = 163.8;
  REFL166_PHASE = 47.5;
  REFL199_PHASE = 178.5;
  %   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 +REFL33_PHASE);	% f1 demod I
  opt = addProbeAt(opt, 'REFL Q1', nREFL, 'in', f1,90+REFL33_PHASE);	% f1 demod Q
  opt = addProbeAt(opt, 'REFL IM', nREFL, 'in', fm,0+REFL133_PHASE);	% fm demod I
  opt = addProbeAt(opt, 'REFL QM', nREFL, 'in', fm,90+REFL133_PHASE);	% fm demod Q
  opt = addProbeAt(opt, 'REFL IP', nREFL, 'in', fp,0+REFL199_PHASE);	% fp demod I
  opt = addProbeAt(opt, 'REFL QP', nREFL, 'in', fp,90+REFL199_PHASE);	% fp demod Q
  opt = addProbeAt(opt, 'REFL I2', nREFL, 'in', f2,0 +REFL166_PHASE);	% f2 demod I
  opt = addProbeAt(opt, 'REFL Q2', nREFL, 'in', f2,90+REFL166_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
  AS33_PHASE = 7;
  AS133_PHASE = 11.65;
  AS166_PHASE = 0;
  AS199_PHASE = 24.24;
  opt = addProbeAt(opt, 'AS DC', nAS, 'in',  0,   0);		% DC
  opt = addProbeAt(opt, 'AS I1', nAS, 'in', f1,AS33_PHASE+0);		% f1 demod I
  opt = addProbeAt(opt, 'AS Q1', nAS, 'in', f1,AS33_PHASE+90);		% f1 demod Q
  opt = addProbeAt(opt, 'AS IM', nAS, 'in', fm,AS133_PHASE+0);		% fm demod I
  opt = addProbeAt(opt, 'AS QM', nAS, 'in', fm,AS133_PHASE+90);		% fm demod Q
  opt = addProbeAt(opt, 'AS IP', nAS, 'in', fp,AS199_PHASE+0);		% fp demod I
  opt = addProbeAt(opt, 'AS QP', nAS, 'in', fp,AS199_PHASE+90);		% fp demod Q
  opt = addProbeAt(opt, 'AS I2', nAS, 'in', f2,AS166_PHASE+0);		% f2 demod Q
  opt = addProbeAt(opt, 'AS Q2', nAS, 'in', f2,AS166_PHASE+90);		% f2 demod Q

  % POX signals (IX pick-off)
  %   POX DD phases set for SRC
  POX33_PHASE = 7;
  POX133_PHASE = -76.152;
  POX199_PHASE = 117.894;
  opt = addProbeAt(opt, 'POX DC', nPOX, 'in',  0,   0);  		% DC
  opt = addProbeAt(opt, 'POX I1', nPOX, 'in', f1,POX33_PHASE+0);	% f1 demod I
  opt = addProbeAt(opt, 'POX Q1', nPOX, 'in', f1,POX33_PHASE+90);	% f1 demod Q
  opt = addProbeAt(opt, 'POX IM', nPOX, 'in', fm,POX133_PHASE+0);	% fm demod I
  opt = addProbeAt(opt, 'POX QM', nPOX, 'in', fm,POX133_PHASE+90);	% fm demod Q
  opt = addProbeAt(opt, 'POX IP', nPOX, 'in', fp,POX199_PHASE+0);	% fp demod I
  opt = addProbeAt(opt, 'POX QP', nPOX, 'in', fp,POX199_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

  % temporary probe
  opt = addProbeAt(opt, 'TMP', nPR, 'bk', 0, 0);
