% compute signals from audio frequency fields
%
% vSig = sigAC(opt, pos, f)
% pos - positions of each optic (Nopt x 1)
% f   - excitation/response frequencies
% vSig - signals for each probe (Nprb x Nopt x Naf)
%
% Example:
% opt = optFP;
% pos = zeros(opt.Noptic, 1);
% f = logspace(0, 4, 300)';
% vSig = sigAC(opt, pos, f);

function vSig = sigAC(opt, pos, f)

  % === Get Field Source Info
  [vSrc, vFrf] = getFieldSource(opt);
  
  % ==== Sizes of Things
  Nopt = opt.Noptic;				% number of optics
  Nprb = opt.Nprobe;				% number of probes
  Nfld = opt.Nfield;				% number of fields
  Nrf  = length(vFrf);				% number of RF components
  Naf  = length(f);				% number of audio frequencies

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ==== Signals
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  rslt = compute(opt, pos, f);
  vDC = rslt.vDC;
  tAF = rslt.tAF;
  vSig = zeros(Nprb, Nopt, Naf);
  
  % 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
          % ==== Field Components
          % for the lower RF component
          aDCm = vDC(nFld, mRF) * phi;
          aAFmm = tAF(2 * Nfld * (mRF - 0.5) + nFld, :, :);
          aAFmp = tAF(2 * Nfld * (mRF - 1) + nFld, :, :);
           
          % for the upper RF component
          aDCp = vDC(nFld, pRF) * conj(phi);
          aAFpm = tAF(2 * Nfld * (pRF - 0.5) + nFld, :, :);
          aAFpp = tAF(2 * Nfld * (pRF - 1) + nFld, :, :);
          
          % ==== Add this RF beat to the signal
          vSig(nPrb, :, :) = vSig(nPrb, :, :) + 0.5 * ...
              (conj(aDCm) * aAFpm + aDCm * aAFpp + ...
               aDCp * aAFmp + conj(aDCp) * aAFmm);

          %[nPrb, mRF, pRF]
          %[squeeze(conj(aDCm) * aAFpm)', ...
          % squeeze(aDCm * aAFpp)', ...
          % squeeze(aDCp * aAFmp)', ...
          % squeeze(conj(aDCp) * aAFmm)']
        end
      end
    end
  end
