% Compute demod signals at each node
%
% [md, mp] = dcDemod(opt, nIn, fMod, gMod, order, dSpec, z)
%
% opt - Optickle model
% nIn - input node
% fMod - modulation frequency
% gMod - modulation depth (gamma)
% order - bessel expansion order
% dSpec - Nx4 detector specification 
%   col 1 - optic index
%   col 2 - port index
%   col 3 - demodulation phase in radians
%   col 4 - demodulation frequency, as an integer multiple of fMod
% z - optic position offsets
%
% md - demod signals at each node
% mp - power signals at each node

function [md, mp] = dcDemod(opt, nIn, fMod, gMod, order, dSpec, z)

  % check arguments
  if size(z, 2) ~= opt.Nopt
    error('size(z, 2) == %d ~= %d == number of optics', ...
           size(z, 2), opt.Nopt)
  end
  Nz = size(z, 1);
  Nd = size(dSpec, 1);
  
  % compute field propagation
  vo = -order:order;
  f = fMod * vo;			% sideband frequencies
  aSB = diag(besselj(vo, gMod));	% sideband amplitudes
  Nf = length(vo);
  
  % compute signals
  md = zeros(Nz, Nd);
  mp = zeros(Nz, Nd);
  
  phiDemod = exp(-i * dSpec(:, 4));
  for nn = 1:Nz
    % compute node-detector transform
    mdet = zeros(Nd, opt.N);
    phi = exp(-2i * opt.k * z(nn, :));
    for nDet = 1:Nd
      nOpt = dSpec(nDet, 1);
      nPort = dSpec(nDet, 2);
      
      % front side input
      nNode = opt.optic(nOpt).in(1);
      if nNode ~= 0
        if( nPort ~= 1 )
          mdet(nDet, nNode) = opt.node(nNode).m(nPort);
        else
          mdet(nDet, nNode) = opt.node(nNode).m(nPort) * phi(nOpt);
        end
      end
      
      % back side input
      nNode = opt.optic(dSpec(nDet, 1)).in(2);
      if nNode ~= 0
        if( nPort == 1 )
          mdet(nDet, nNode) = opt.node(nNode).m(nPort);
        else
          mdet(nDet, nNode) = opt.node(nNode).m(nPort) * conj(phi(nOpt));
        end
      end
    end  

    % compute detector fields
    mf = propFields(opt, f, z(nn, :));		% propagate fields
    mf = squeeze(mf(:, nIn, :)) * aSB;		% apply input amplitudes
    mf = mdet * mf;				% transform to detector basis

    % demod signals
    mf = mf.';
    for nfa = 1:Nf
      for nDet = 1:Nd
        nfb = nfa + dSpec(nDet, 3);
        if nfb <= Nf
          md(nn, nDet) = md(nn, nDet) + phiDemod(nDet) * ...
              mf(nfa, nDet) .* conj(mf(nfb, nDet));
        end
      end
    end

    % power signals
    mp(nn, :) = sum(abs(mf .* conj(mf)));
  end
