% 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 [varargout] = senseMat(varargin)

  if nargin < 1
    f = logspace(-1,4,200)';
  elseif nargin == 1
    f = varargin{1};
  elseif nargin > 1
    error('Poo')
  end
  
  f = f(:);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % RF detection scheme
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  opt = mLIGO;

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

  %nASQ = strmatch('AS Q1', prbNames);	% which signal to use for RF
  nASDC = strmatch('AS DC', prbNames);	% which signal to use for DC
  
  nPOXI = strmatch('POX I1', prbNames);
  nPOXQ = strmatch('POX Q1', prbNames);
  
  nREFLI = strmatch('REFL I1', prbNames);
  
  nFreq = getDriveIndex(opt, 'Source'); % phase modulation input
  nAM = getDriveIndex(opt, 'AM');
  nEX = getDriveIndex(opt, 'EX');
  nEY = getDriveIndex(opt, 'EY');
  nPR = getDriveIndex(opt, 'PR');
  nBS = getDriveIndex(opt, 'BS');
  
  % LSC fringe offsets in meters
  dLm = 17e-12;
  dLp = 1e-16;
  dlm = 1e-19;
  dlp = 1e-19;
  
  offDC = dLm;
  
  % run the small-offset simulation for RF locking
  pos = zeros(opt.Noptic, 1);
  pos(nEX) = offDC / 2;
  pos(nEY) = -offDC / 2;

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

  % ==== Optical transfer functions
  % hFred and hBarney are from mirror position to sensor signal (W/m)
  hAX = squeeze(sigAC(nASDC, nEX, :));
  hAY = squeeze(sigAC(nASDC, nEY, :));
  hDARM = (hAX - hAY);
  
  hAX = squeeze(sigAC(nREFLI, nEX, :));
  hAY = squeeze(sigAC(nREFLI, nEY, :));
  hCARM = (hAX + hAY);
  
  hPR1 = squeeze(sigAC(nPOXI, nPR, :));
  hPRC = hPR1;

  hBS2 = squeeze(sigAC(nPOXI, nBS, :));
  hBS = squeeze(sigAC(nPOXQ, nBS, :));
  
 % -- - -- Other Couplings - --- - - -- --- ----- ---- -- - - - 
 
  hPR2AS = squeeze(sigAC(nASDC, nPR, :));
 
  % Set the MICH output matrix at 100 Hz
  nnn = min(find(f > 100));
  MICH2PR = abs(hBS2(nnn)/hPR1(nnn));
  hnull = hBS2 - MICH2PR * hPR1;
  hMICH = hBS - MICH2PR * hPR1;
  
  hMICH2AS = squeeze(sigAC(nASDC, nBS, :)) - ...
      MICH2PR * squeeze(sigAC(nASDC, nPR, :));
  
  % intensity noise
  hAMRF = squeeze(sigAC(nASDC, nAM, :));
  aAMRF = abs(hAMRF ./ hDARM);

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

  pendeTF = squeeze(mechTF(nEX,nEY,:));
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Plots
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot('Position',[0.08,0.52,0.89,0.43]);
  loglog(f, abs(hDARM), 'b',...
         f, abs(hCARM), 'r',...
         f, abs(hPRC), 'm',...
         f, abs(hMICH), 'k',...
         f, abs(hMICH2AS), 'k--',...
         f, abs(hPR2AS), 'm--');
  title('Sensing Matrix')
  legend('DARM', 'CARM', 'PRC', 'MICH', 'MICH->AS','PRC->AS')
  %xlabel('Hz')
  set(gca,'XTickLabel',[])          
  ylabel('Mag [W/m]',...
         'FontWeight','bold')
  grid on
  grid minor
  axis([min(f) max(f) 1e5 1e11])

  subplot('Position',[0.08,0.05,0.89,0.43]);
  semilogx(f, 180/pi*angle(hDARM), 'b',...
         f, 180/pi*angle(hCARM), 'r',...
         f, 180/pi*angle(hPRC), 'm',...
         f, 180/pi*angle(hMICH), 'k',...
         f, 180/pi*angle(hMICH2AS), 'k--',...
         f, 180/pi*angle(hPR2AS), 'm--');
  axis([min(f) max(f) -181 181])
  set(gca,'YTick',[-270:45:270])
  xlabel('Hz')
  ylabel('Phase [deg]',...
         'FontWeight','bold')
  grid on
  grid minor
  orient tall
  print -dpdf mLIGO-SensingMatrix.pdf
  
  
  if nargout > 0
    varargout{1} = [f hDARM hCARM hPRC hMICH];
  end
  
