function [xIn, Ein, xOut, Eout, alpha, wOut, ROCout, gouy, dotProd] = ...
    focusingCavity( Wfft, Nfft, wIn, ROCin, R1, R2, dz, radius, aIn )
%
% function [xIn, Ein, xOut, Eout, alpha, wOut, ROCout, dotProd] = ...
%    focusingCavity( Wfft, Nfft, wIn, ROCin, R1, R2, dz )
%
%  Field propagation through a strngly focusing cavity.
%  Calculate Eout with TEM00 input which is reflected by mirror 1,
%  propagates to mirror 2 and is reflected by mirrror 2.
%
%  inputs
%    Wfft and Nfft define the FFT grids of the input field
%    wIn and ROCin are input field parameters
%    R1 and R2 are curvatures of mirrors forming the focusing cavity
%    and dz is the spacing of the mirrors.
%    radius is the radius of the first mirror
%    aIn is the ratio forced to be used, 0 to use a optimal value
%
%  outputs
%    (xIn, Ein) : input field coordinates and the complex amplitude
%    (xOut, Eout ) : output field coordinates and the complex amplitude
%    alpha : used scaling factor
%    (wOut,ROCout) : beam size and ROC after propagation
%    gouy : gouy phase through this propagation
%    dotProd : Integral( Eout, TEM00-after-propagation )
%
% by Hiro Yamamoto on January 5, 2010 
% documentation : T0900640

k = 2*pi/1.064e-6;
dx = Wfft / Nfft;
xIn = (-Wfft/2+dx/2):dx:(Wfft/2-dx/2);

% scaling of the coordinate at the end of cavity
if aIn > 0
   alpha = aIn;
else
    alpha = 1 + dz/(-R1/2);
end
xOut = alpha * xIn;

% TEM00 input field
Ein = gaussFld( xIn, wIn, ROCin, radius );
Ein = Ein .* exp( 1i*k/2*((alpha-1)/dz-1/(-R1/2)) * Rsquare(xIn) );

% simple FFT propagation
Eout = fftProp( Ein, 1/Wfft, dz/alpha );

% correction to calculate the field after reflected by mirror 2
Eout = Eout / alpha .*  ...
    exp( -1i*k/2 *((1-1/alpha)/dz-1/(R2/2))*Rsquare(xOut) );

% modal model calculation of the field parameters after mirror 2
[wOut, ROCout, gouy] = MMProp( wIn, ROCin, dz, R1, R2 );

% compare the FFT result and the modal model calculation
EoutMM = gaussFld( xOut, wOut, ROCout, 0 ) * exp( 1i * gouy );
dotProd = sum(sum(Eout .* conj(EoutMM))) * (dx*alpha)^2;

return

function E = gaussFld( x, w, ROC, radius )
% simple TEM00 field generator within radius

k = 2*pi/1.064e-6;
E = sqrt(2/pi) / w * exp( -(1/w^2+1i*k/(2*ROC))*Rsquare(x) );
if radius > 0
    E = E .* (Rsquare(x) < radius^2);
end
return

function Eout = fftProp( Ein, k0, L )
% simple FFT propagator

k = 2*pi/1.064e-6;
Nfft = size( Ein, 1 );
klist = 2*pi*k0*[(0:Nfft/2),(-Nfft/2+1):-1];

Fin = fft2( Ein );
Fout = Fin .* exp( 1i/(2*k)*L*Rsquare(klist) );

Eout = ifft2( Fout );
return

function [wOut, ROCout, gouy] = MMProp( wIn, ROCin, dz, R1, R2 )
% modal model progagation by 
% reflection by a mirror with ROC of R1,
% propagation through dz
% reflection by a mirror with ROC of R2

k = 2*pi/1.064e-6;
qinvIn = 1/ROCin - 1i * 2 / k / wIn^2;

% reflection by mirror 1
CA = 1;
CB = 0;
CC = -2/R1;
CD = 1;

qinvOut = ( CC + CD*qinvIn ) / ( CA + CB*qinvIn );
z = real( 1/qinvOut );
z0 = imag( 1/qinvOut );
gouy = atan2( z+dz, z0 ) - atan2( z, z0 );

% full propagation
CA = 1 - 2*dz/R1;
CB = dz;
CC = -2*(1-2*dz/R2)/R1 - 2/R2;
CD = 1 - 2*dz/R2;

qinvOut = ( CC + CD*qinvIn ) / ( CA + CB*qinvIn );

ROCout = 1/real(qinvOut);
wOut = sqrt( -2/k/imag(qinvOut) );

return

function r2 = Rsquare( x )
if size(x,2) == 1
    x = x';
end
r1 = (x.') * ones( 1, size(x, 2) );
r2 = r1.^2 + (r1').^2;

return

