% Add a mirror to the model.
%
% [opt, sn] = addMirror(opt, name, aio, Chr, Thr, Lhr, Rar, Lmd, Nmd)
%
% default parameter values are:
% [opt, sn] = addMirror(opt, name, 0, 0, 0, 0, 1e-4, 0, 1.45);
%
% Paramters:
% aio - angle of incidence (in degrees)
% Chr - curvature of HR surface (Chr = 1 / radius of curvature)
% Thr - power transmission of HR suface
% Lhr - power loss on reflection from HR surface
% Rar - power reflection of AR surface
% Nmd - refractive index of medium (1.45 for fused silica, SiO2)
% Lmd - power loss in medium (one pass)
%
% Chr is the mirror "curvature" which is the inverse of the radius
% of curvature (c = 1/roc = 1/(2 * f), where f is the focal length).
% Positive curvatures are used for concave, zero for a flat mirror,
% and negative values for convex mirrors.
% 
% A mirror has 2 inputs and 4 outputs.
% Inputs:
% 1, fr, HR = front (HR side)
% 2, bk, AR = back (AR side)
%
% Outputs:
% 1, fr = front
% 2, bk = back
% 3, pi = pick-off sampling back input
% 4, po = pick-off sampling back output

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

  % assign default arguments
  defs = [0, 0, 0, 0, 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(4, 2);
  
  % 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);

  % ======== Compute Basis Transform Matrix
  qm = planeConvexOpHG(aoi, Chr, Nmd);
  
  % ======== Add Optic
  type = 'Mirror';
  inNames = {{'fr' 'HR'}, {'bk' 'AR'}};
  outNames = {{'fr' 'HR'}, {'bk' 'AR'}, {'pi'}, {'po'}};
  [opt, sn] = addOptic(opt, type, name, m, d, qm, inNames, outNames);

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