% Add a beam splitter to the model.  A beam splitter is essentially
% 2 mirrors which share a common position.  The two sets of 4x2 transfer
% matrices, each identical to that of a mirror, produce an 8x4 matrix
% with no cross-terms.  These two sets are simiply called A and B, as
% there is no distinguishing feature to differentiate them with.
%
% [opt, sn] = addBeamSplitter(opt, name, aio, Chr, Thr, Lhr, Rar, Lmd, Nmd)
%
% default parameter values are:
% [opt, sn] = addBeamSplitter(opt, name, 45, 0, 0.5, 1e-4, 1e-4, 0, 1.45);
% (see also addMirror, as it has similar paramters)
%
% A beam splitter has 4 inputs and 8 outputs.
% Inputs:
% 1, frA = front A
% 2, bkA = back A
% 3, frB = front B
% 4, bkB = back B
%
% Outputs:
% 1, frA = front
% 2, bkA = back
% 3, piA = pick-off sampling back input
% 4, poA = pick-off sampling back output
% 5, frB = front
% 6, bkB = back
% 7, piB = pick-off sampling back input
% 8, poB = pick-off sampling back output

function [opt, sn] = addBeamSplitter(opt, name, aoi, Chr, Thr, Lhr, Rar, Lmd, Nmd)

  % assign default arguments
  defs = [45, 0, 0.5, 1e-4, 1e-4, 0, 1.45];
  argNames = {'aoi', 'Chr', 'Thr', 'Lhr', 'Rar', 'Lmd', 'Nmd'};
  for n = 1:length(defs);
    if nargin < n + 2
      eval([argNames{n} ' =  defs(n);']);
    end
  end  

  % ======== Compute Coupling Matrix
  % amplitude reflectivities, transmissivities and phases
  hr = -sqrt(1 - Thr - Lhr);			% HR refl
  ht =  sqrt(Thr);				% HR trans
  ar = -sqrt(Rar);				% AR refl
  at =  sqrt(1 - Rar);				% AR trans
  bt =  sqrt(1 - Lmd);				% bulk trans
  
  % transmission combinations
  hrbt = -hr * bt;
  arbt = -ar * bt;
  htbt =  ht * bt;
  atbt =  at * bt;

  % transfer matrix  
  m = zeros(8, 4);
  
  % from front input node (in1)
  m(1, 1) = hr;
  m(2, 1) = atbt * ht;
  m(3, 1) = 0;
  m(4, 1) = atbt * hrbt * arbt * ht;

  % from back input node (in2)
  m(1, 2) = htbt * at;
  m(2, 2) = atbt * hrbt * at;
  m(3, 2) = ar;
  m(4, 2) = atbt * hrbt * arbt * hrbt * at;

  % dl/dx (non-zero for reflected fields)
  d = zeros(4, 2);
  d(1, 1) = -2 * cos(pi * aoi / 180);
  d(2:4, 2) = 2 * cos(pi * aoi / 180);

  % B-side matrix components
  m(5:8, 3:4) = m(1:4, 1:2);
  d(5:8, 3:4) = d(1:4, 1:2);
  
  % ======== Compute Basis Transform Matrix
  qm = repmat(nullOpHG, 8, 4);
  q0 = planeConvexOpHG(aoi, Chr, Nmd);
  qm(1:4, 1:2) = q0;
  qm(5:8, 3:4) = q0;
  
  % ======== Add Optic
  type = 'BeamSplitter';
  inNames = {{'frA'}, {'bkA'}, {'frB'}, {'bkB'}};
  outNames = {{'frA'}, {'bkA'}, {'piA'}, {'poA'}, ...
              {'frB'}, {'bkB'}, {'piB'}, {'poB'}};
  [opt, sn] = addOptic(opt, type, name, m, d, qm, inNames, outNames);

  % save some info
  opt.optic(sn).p = [aoi, Chr, Nmd];
