% Display the some DARM noises.
%
% respDARM(f, offset)
% f - frequency vector for plots
% offDC - DARM offset (m) for DC readout
% offRF - DARM offset (m) for RF readout
%
% Example:
% f = logspace(0, 4, 200)';
% respDARM(f, 1e-11, 1e-13)

function lasernoisecoupling(f, offDC, offRF)

  
  % Make the frequency vector into a column vector
  f = f(:);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % RF detection scheme
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  opt = mLIGO;

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

  nRF = strmatch('ASp Q1', prbNames);	% which signal to use for RF
  nDC = strmatch('AS DC', prbNames);	% which signal to use for RF
  nFreq = getDriveIndex(opt, 'Source'); % phase modulation input
  nAM = getDriveIndex(opt, 'AM');
  nEX = getDriveIndex(opt, 'EX');
  nEY = getDriveIndex(opt, 'EY');
  
  % run the small-offset simulation for RF locking
  pos = zeros(opt.Noptic, 1);
  pos(nEX) = offRF / 2;
  pos(nEY) = -offRF / 2;

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

  % ==== Optical transfer functions
  % hRF and hDC are from mirror position to sensor signal (W/m)
  hAX = squeeze(sigAC(nRF, nEX, :));
  hAY = squeeze(sigAC(nRF, nEY, :));
  hRF = (hAX - hAY);

  % intensity noise
  hAMRF = squeeze(sigAC(nRF, nAM, :));
  aAMRF = abs(hAMRF ./ hRF);

  % frequency noise
  hFreqRF = squeeze(sigAC(nRF, nFreq, :)) ./ (i * f);
  aFreqRF = abs(hFreqRF ./ hRF);

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

  % run the offset simulation for RF locking
  opt = optLIGOp;
  [sigDC, sigAC, mechTF, fDC] = compute(opt, pos, f);

  % hDC is from mirror position to sensor signal (W/m)
  hAX = squeeze(sigAC(nDC, nEX, :));
  hAY = squeeze(sigAC(nDC, nEY, :));
  hDC = (hAX - hAY);

  % intensity noise
  hAMDC = squeeze(sigAC(nDC, nAM, :));
  aAMDC = abs(hAMDC ./ hDC);

  % frequency noise
  hFreqDC = squeeze(sigAC(nDC, nFreq, :)) ./ (i * f);
  aFreqDC = abs(hFreqDC ./ hDC);

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Plots
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot('Position',[0.09,0.52,0.89,0.43]);
  loglog(f, aAMDC, 'b', f, aAMRF, 'r');
  %title('Amplitude Noises Coupling')
  legend('AM DC', 'AM RF', 'Location', 'SouthWest')
  %xlabel('Hz')
  set(gca,'XTickLabel',[])
  ylabel('Amplitude Noise Coupling [m/RIN]')
  grid on
  axis([9 10000 1e-14 1e-9])

  subplot('Position',[0.09,0.06,0.89,0.43]);
  loglog(f, aFreqDC, 'b', f, aFreqRF, 'r');
  %title('Frequency Noises Coupling')
  legend('Freq DC', 'Freq RF', 'Location', 'SouthWest')
  xlabel('Hz')
  ylabel('Frequency Noise Coupling [m/Hz]')
  grid on
  axis([9 10000 1e-15 1e-12])
