% Display the DARM loop and a few noises.
%
% [sigDC, sigAC, mechTF] = respDARM(f, sigAC, mechTF)
% f - frequency vector for plots
% sigAC and mechTF - from compute
%
% This is an EXAMPLE loop for DARM and some noise projections.  The loop
% is certainly not good, and the noise projections are almost certainly
% wrong, but this should show how Optickle can be used.
%
% Example:
% f = logspace(0, 4, 1000)';
% opt = optLIGO2;
% [sigDC, sigAC, mechTF, fDC] = compute(opt, [], f);
% respDARM(opt, f, sigAC, mechTF, fDC);

function shotDARM(f,offDC)

  opt = mLIGO;

  nEX = getDriveIndex(opt, 'EX');
  nEY = getDriveIndex(opt, 'EY');

  pos = zeros(opt.Noptic, 1);
  pos(nEX) = offDC / 2;
  pos(nEY) = -offDC / 2;

  f = f(:);
  
  [sigDC, sigAC, mechTF, fDC] = compute(opt, pos, f);

  % probe names and numbers
  oprb = opt.probe;
  prbNames = {oprb.name};

  % some parameters
  wE = 2 * pi * 0.75;			% resonance of end mirror (rad/s)
  mE = 10;				% mass of end mirror (kg)
  gAct = 1e-3;				% actuator DC gain (N/V)
  gSensor = 1e10;				% sensor DC gain (W/m)
  fUGF = 200;				% unity gain of servo (Hz)
  nSig = strmatch('AS DC', prbNames);	% which signal to use
  nEX = getDriveIndex(opt, 'EX');
  nEY = getDriveIndex(opt, 'EY');
  nAS = getFieldAt(opt, 'AS', 'in');
  
  dampedRes = [-0.1 + 1i, -0.1 - 1i];
  
  % ==== Optical transfer functions
  % hSens is from mirror position to sensor signal (W/m)
  hAX = squeeze(sigAC(nSig, nEX, :));
  hAY = squeeze(sigAC(nSig, nEY, :));
  hSig = (hAX-hAY);

  % ==== Technical transfer functions
  w = 2 * pi * f;

  % tranfer function from actuator voltage to free mirror position (m/V)
  %tfAct = zpk([], wE * dampedRes, gAct / mE);
  %hAct = squeeze(freqresp(tfAct, w));
  tfAct = zpk([], wE * dampedRes, gAct / mE);
  hAct = squeeze(freqresp(tfAct,w));

  % tranfer function from optical RF to sensor voltage (V/W)
  tfSensor = zpk([], 1e6 * dampedRes, gSensor);
  hSensor = squeeze(freqresp(tfSensor, w));

  % servo transfer function (V/V)
  tfServo = zpk(-2 * pi * [15+30;15-30i;90;10], -2 * pi * [0;0;4e3;4e3], 1);
  hServo = squeeze(freqresp(tfServo,w));
  
  % set servo gain to make the desired UGF.
  hLoop = hSensor .* hSig .* hAct .* hServo;
  mugf = max(find(f < fUGF));
  gugf = abs(hLoop(mugf));
  hLoop = hLoop / gugf;
  
  % hLoop is the open loop TF
  figure(101)
  subplot(2, 1, 1)
  mybodeplot(f,hLoop)
  %title('DARM open loop gain')
  %legend('loop','sensor','sig','act','servo')
  %xlabel('Hz')
  %axis([min(f) max(f) 0.01 1e4])
  %grid on
  
  % === Noise propagation to the error point

  % a flat "shot" noise
  % see also Buonanno & Chen gr-qc/0102012
  %   http://www.iop.org/EJ/article/0264-9381/18/15/102/q115l2.html
  % or Bench
  %   http://www.its.caltech.edu/~rana/aLIGO/Bench/bench33.m
  Phv = 2 * opt.h * opt.c / opt.lambda;
  Pas = sum(abs(fDC(nAS, :))).^2;
  aSensor = sqrt(Phv * Pas) ./ abs(hSig);

  % a flat actuation noise (force) on then mirror
  aAct = 5e-12 * abs(hAct);

  % incoherent sum
  aSum = sqrt(aSensor.^2 + aAct.^2);
  
  figure(103)
  loglog(f, [aSum, aSensor, aAct]);
  title('A few noises')
  legend('Sum', 'Sensor', 'Actuator')
  xlabel('Hz')
  ylabel('m/rtHz')
  a0 = axis;
  a0(3:4) = [1e-20 1e-17];
  axis(a0);
  grid on
