% Offset DARM generates plots related to selecting the DARM
% offset for DC readout.  (This function has no arguments and
% no return values.)

function varargout = offsetDARM(varargin)

  % prepare some variables
  sweepSetup

  nDC = nAS(1);

  % Whether or not to fit the AS_DC plot
  fitstr = 0; % 0 for NO, 1 for YES
  
  if nargin > 0
    npts = varargin{1};
    if nargin > 1
      fitstr = varargin{2};
      if nargin > 2
        error('Too Many input arguments')
      end
     end
  elseif nargin == 0
    npts = 200;
    % Need 2000 pts to get a 0.5% accuracy in the power
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % DARM sweep for AS power and signal as a function of DARM offset
  pos = zeros(opt.Noptic, 1);
  pos(nEX) = 5e-13;
  pos(nEY) = -5e-13;
  pos = pos * logspace(-1, 2, npts);
  
  [xPos, sigDC, fDC] = sweepDC(opt, pos);
  x = 2 * xPos(nEX, :);
  fDC = permute(fDC, [1 3 2]);
  
  %%%%%%%%%%%%%%%% Plot DARM signal

  % DARM signal, just a difference of DC values
  sd = [diff(sigDC(nDC, :)) ./ diff(x)];
  nux = x(2:end);
  
  % shot noise, more or less, with some normalization
  ad = 5.92e10 * sqrt(sigDC(nDC, 2:end));

  % SNR
  snare = sd ./ ad;
  
  % plot signal / noise ratio
  subplot('Position',[0.08,0.53,0.89,0.43]);
  semilogx(nux * 1e12, snare)
  %xlabel('DARM offset [pm]')
  set(gca,'XTickLabel',[])
  title('DARM Shot Noise SNR (arbitrary units)')
  grid on
  axis([1e12*min(x) 1e12*max(x) 0 1])
  
  %%%%%%%%%%%%%%%% Plot AS port power
  Pas = sigDC(nDC,:);
  x12 = x * 1e12;
  
  subplot('Position',[0.08,0.06,0.89,0.43]);
  loglog(x12, Pas,'b',...
         x12, sigDC(nPOX(1), :),'r')
  %loglog(x * 1e12, sigDC(nDC, :))
  title('AS DC (after OMC) and POX DC (100 ppm pick-off)')
  xlabel('DARM offset [pm]')
  ylabel('P_{AS} [W]')
  legend('AS DC', 'POX DC', 'Location', 'NorthWest');
  grid on
  grid minor
  axis([min(x12) max(x12) min(Pas) max(Pas)])

  [yy,ii] = max(snare);
  optoff = nux(ii);
  
  disp(['Optimum offset = ' num2str(optoff*1e12,'%9.1f') ' pm'])
  disp(['AS Power = ' num2str(Pas(ii+1)*1000,'%9.1f') ' mW'])

  
 if fitstr > 0 
 % Fit AS_DC to a polynomial to get the nonlinear terms =============== 
  nn = find(x > optoff*0.8 & x < optoff*1.2);
 
  xnew = x(nn) - optoff;
  Perr = Pas(nn) - Pas(ii+1); 
  
  pfitcoef = polyfit(xnew,Perr,4)'
  pfitmodel = polyval(pfitcoef,xnew);
  %pfitresidual = sqrt(sum((Pas-pfitmodel).^2))/length(Pas);
  %if pfitresidual > 1e-1
  %  error('Bad fit for AS_DC')
  %end
 % =================================================================== 

 clf
 plot(xnew,pfitmodel)

 end
 
 
 
  if nargout > 1
    varargout{1} = x;          % dLm offset
    varargout{2} = Pas;        % AS_DC as a function of x
    if nargout > 2
      varargout{3} = pfitcoef'; % polyfit coefficients to AS_DC(x12)
    end
  end

  
  
return
