% function [Fcav, Ein, rList] = HnklFP( N, a, T1, R1, R2, invF, Lcav, phi )
% Field in a cavity formed by RM and ITM. 
% ITM HR is facing away from RM, so cavity length should be lambda/4 off
% 
% T1 : power transmittance of front mirror
% R1/R2 : ROC of front and back mirrro
% invF : inverse of the f number in the second mirror
% Lcav : cavity length
% phi : extra cavity phase. Should resonate at -eta00+pi/2

function [Fcav, Ein, rList] = ...
    HnklFP( N, a, T1, R1, R2, invF, Lcav, phi )

global kLam__ nSubs__ lambda__ oneN__
lambda__ = 1.064e-6;
kLam__ = 2*pi/lambda__;
nSubs__ = 1.44963;
oneN__ = diag(ones(N,1),0);

zetaList = zeros(N,1);
for in = 1 : N
    zetaList(in) = getZhiK( in-1, 1e-6 );
end
[Hplus, J0Norm] = calcHPlus( a, zetaList );

rList = a / zetaList(end) * zetaList;

E00 = HG00( R1, -R2/nSubs__, Lcav, rList );
Ein = sqrt(T1) * E00;

MRef1 = MirrorFrontRef( sqrt(1-T1), 1/R1, 0, rList );
MRef2 = MirrorBackRef( 1, 1/R2, invF, rList );
FProp = Hplus \ diag(Propagator( Lcav, a, zetaList ),0) * Hplus;

Cround = diag(MRef1,0) * FProp * diag(MRef2,0) * FProp;

Nphi = length( phi );
Fcav(1:Nphi) = struct( 'Ecav', zeros(N,1), 'PowerTotal', 0, 'Power00', 0 );

for iphi = 1 : Nphi
    Fcav(iphi).Ecav = ...
        (oneN__ - exp(2i*phi(iphi))*Cround) \ Ein;
    Fcav(iphi).PowerTotal = ...
        scalarProd( Fcav(iphi).Ecav, Fcav(iphi).Ecav, J0Norm );
    Fcav(iphi).Power00 = ...
        abs(scalarProd( Fcav(iphi).Ecav(:,1), E00, J0Norm )).^2;
end

return

% mirror front reflection
% Eref = r * exp( i*k*(2/2ROC + 2*del) * r^2 )
% del is surface deformation
function RefM = MirrorFrontRef( r, invROC, del, rList )
global kLam__
RefM = -r * exp( 1i*kLam__*( invROC + 2*del ).* rList.^2 );
return

% mirror back reflection
% Eref = r * exp( i*k*(-2*n/2ROC + 2/f) * r^2 )
% invF is 1/f, lens effect in the substrate
function RefM = MirrorBackRef( r, invROC, invF, rList )
global kLam__ nSubs__
RefM = r * exp( 1i*kLam__*( -nSubs__*invROC + 2*invF ).*rList.^2 );
return

% mirror transmission
% Eout = t * exp( -i*k*((n-1)/2ROC - 2/f) )
function TraM = MirrorTrans( t, invROC, invF, rList )
global kLam__ nSubs__
TraM = t * exp( -1i*kLam__*( (nSubs__ - 1)*invROC/2 - invF )*rList.^2 );
return

% propagator
function G = Propagator( dist, a, zetaList )
global kLam__
G = exp(1i*dist/(2*kLam__)/a^2*zetaList.^2);
return

function [H, qa] = calcHPlus( a, zetaList )

N = length( zetaList );
H = zeros( N, N );
qa = (zetaList(end) * besselj(0, zetaList)/a).^2 / 2;
for ix = 1 : N
    H(ix,:) = besselj(0, zetaList(ix)/zetaList(end)*zetaList) ./ qa;
end
return

function zhi = getZhiK( k, accuracy )

if k <= 0
    zhi = 0;
else
    zhi = ( k + 0.25 ) * pi;
	delzval = 0.01;
	while delzval > accuracy
        rat10 = besselj(1,zhi) / besselj(0, zhi);
        zhiNew = zhi * (1 - rat10/(zhi - rat10) );
        delzval = abs( zhiNew - zhi );
        zhi = zhiNew;
    end
end
return

function sp = scalarProd( f, g, qa )
sp = 2*pi*sum( conj(f) .* g ./ qa );
return

function E00 = HG00( R1, R2, Lcav, rList )
global kLam__
z = -Lcav*(R2-Lcav)/(R1+R2-2*Lcav);
z0 = sqrt( Lcav*(R1-Lcav)*(R2-Lcav)*(R1+R2-Lcav) / (R1+R2-2*Lcav)^2 );
w0 = sqrt(2*z0/kLam__);
w = w0*sqrt( 1 + (z/z0)^2 );
E00 = sqrt( 2/pi ) / w * exp( -(1/w^2 + 1i*kLam__/(2*(-R1))) * rList.^2 );
return
