% mop = planeConvexOpHG(aoi, crv, ior)
%
% aoi = angle of incidence
% crv = 1 / radius of curvature
% ior = index of refraction

function mop = planeConvexOpHG(aoi, crv, ior)

  % operator matrix  
  mop = repmat(nullOpHG, 4, 2);
  
  % inverse focal lengths
  sinFront = sin(pi * aoi / 180);
  cosFront = cos(pi * aoi / 180);
  cosBack = sqrt(1.0 - (sinFront / ior)^2);

  crx = crv * 2.0 / cosFront;
  cry = crv * 2.0 * cosFront;
  ctx = crv * (ior - cosFront / cosBack);
  cty = crv * cosFront * (ior - 1.0);

  % basis operators
  nop.x = eye(2);				% no op
  fop.x = focusOpHG(crx);			% front refl, x-axis
  bop.x = focusOpHG(-ior * crx);		% back refl, x-axis
  top.x = focusOpHG(-ctx);			% transmission, x-axis

  nop.y = eye(2);				% no op
  fop.y = focusOpHG(cry);			% front refl, y-axis
  bop.y = focusOpHG(-ior * cry);		% back refl, y-axis
  top.y = focusOpHG(-cty);			% transmission, y-axis
  
  % from front input node (in1)
  mop(1, 1) = fop;
  mop(2, 1) = top;

  mop(4, 1).x = bop.x * top.x;
  mop(4, 1).y = bop.y * top.y;

  % from back input node (in2)
  mop(1, 2) = top;
  mop(2, 2) = bop;
  mop(3, 2) = nop;

  mop(4, 2).x = bop.x * bop.x;
  mop(4, 2).y = bop.y * bop.y;
