% compute signals from DC fields as a function of optic positions
%
% [pos, sigDC, fDC] = sweepDC(opt, posStart, posEnd, Npos)
% pos - optic position sweep (Noptic x Npos)
% sigDC - signal vectors (Nprobe x Npos)
% fDC - field matrices (Nlnk x Nrf x Npos)
%
% see also computeFields and computeSignals
%
% Example:
% opt = optFP;
% pos = zeros(opt.Noptic, 1);
% pos(getDriveIndex('EX')) = 1e-10;
% [xpos, sigDC] = sweepDC(opt, -pos, pos, 101);

function [pos, sigDC, fDC] = sweepDC(opt, posStart, posEnd, Npos)

  % === Get Field Source Info
  vFrf = getFieldSource(opt);

  % ==== Sizes of Things
  Nopt = opt.Noptic;				% number of optics
  Nlnk = opt.Nlink;				% number of links
  Nprb = opt.Nprobe;				% number of probes
  Nrf  = length(vFrf);				% number of RF components

  % ==== Compute Fields
  % prepare position sweep
  if nargin == 2
    % posStart is actually pos
    pos = posStart;
    Npos = size(pos, 2);
  else
    % generate position vectors
    x = (0:Npos-1)' / (Npos - 1);
    dpos = posEnd - posStart;
    pos = zeros(Nopt, Npos);			% optic positions
    for n = 1:Npos
      pos(:, n) = dpos * x(n) + posStart;
    end
  end
  
  % loop over positions
  sigDC = zeros(Nprb, Npos);			% signals
  fDC = zeros(Nlnk, Nrf, Npos);			% fields
  for n = 1:Npos
    fDC(:, :, n) = computeFields(opt, pos(:, n), []);
  end

  % ==== Compute DC Signals
  % this is much faster than using computeSignals
  sigDC = zeros(Nprb, Npos);
  fDCp = permute(fDC, [1 3 2]);			% used in the loops (faster)
  
  % loop over probes
  for nPrb = 1:Nprb
    prb = opt.probe(nPrb);
    
    nFld = prb.nField;
    phi = exp(pi * i * prb.phase / 180);
    
    % look for RF frequency matches
    for mRF = 1:Nrf
      for pRF = mRF:Nrf
        if abs(vFrf(pRF) - vFrf(mRF) - prb.freq) < 1e-3
          aDCm = fDCp(nFld, :, mRF) * phi;
          aDCp = conj(fDCp(nFld, :, pRF));
          sigDC(nPrb, :) = sigDC(nPrb, :) + real(aDCp .* aDCm);
        end
      end
    end
  end
  
