% Compute field propagation matrix
%
% pos - a Nx1 matrix of optic positions where N = opt.Nopt
%  Positions of non-mirror optics are ignored!
%
% ny, nx - HG mode number in y and x (TEM01 => ny = 0, nx = 1)
%
% fldCC = getFieldMatrix(opt, pos, ny, nx)

function fldCC = getFieldMatrix(opt, pos, ny, nx)

  % ==== DC and Distance Matrices
  % mirror position phases
  phi = ones(opt.Nopt, 1);
  if nargin > 1
    for n = 1:opt.Nopt
      if opt.optic(n).tn == opt.types.Mirror
        aoi = opt.optic(n).p(1);		% angle of incidence
        phi(n) = exp(-2i * opt.k * pos(n) * cos(aoi));
      end
    end
  end

  % mirror reflections flip the odd X mode amplitudes
  flipMirror = ones(4, 2);
  flipMirror(1, 1) = (-1)^nx;
  flipMirror(2:4, 2) = (-1)^nx;
  
  % build matrices
  N = opt.Nlink;
  fldCC = sparse(N, N);
  
  for n = 1:N
    % set short-cuts
    lnk = opt.link(n);
    src = opt.optic(lnk.snSource);
    snk = opt.optic(lnk.snSink);

    % find linked inputs
    m = find(src.in ~= 0);
    mm = src.in(m);

    % set static coefficient
    fldCC(n, mm) = src.m(lnk.portSource, m);

    % add mirror reflection phase
    if snk.tn == opt.types.Mirror
      phiMirror = ones(4, 2);
      phiMirror(1, 1) = phi(lnk.snSource);
      phiMirror(2:4, 2) = conj(phi(lnk.snSource));

      fldCC(n, mm) = fldCC(n, mm) .* phiMirror(lnk.portSource, m);
    end
  
    % add guoy phases
    qx = src.qx(lnk.portSource);
    gx = angle(qx + lnk.len) - angle(qx);

    qy = src.qy(lnk.portSource);
    gy = angle(qy + lnk.len) - angle(qy);

    fldCC(n, mm) = fldCC(n, mm) * exp(i * (gx * nx + gy * ny));
  
    % add mirror TEM10 phase flip
    if snk.tn == opt.types.Mirror
      fldCC(n, mm) = fldCC(n, mm) .* flipMirror(lnk.portSource, m);
    end  
  end
