etalon(nu,e0,R,Dfp,A,defect,n0,lambda,etendue,source,tol)

function T = etalon(nu,e0,R,Dfp,A,defect,n0,lambda,etendue,source,tol)

% A realistic model of the etalon frequency response is computed
% numerically based on the model provided in McKay, J. A. (1998). 
% Modeling of direct detection Doppler wind lidar. I. The edge
% technique. Applied optics, 37(27), 6480-6486. equations 7-8

% Inputs
% "nu"       Frequency in units of [GHz]
% "e0"       Spacing in units of [microns]
% "R"        Reflectivity 0-1
% "Dfp"      FPI clear aperture in units of [mm]
% "defect"   Can be a structure containing defect parameters in #lambdas
%            of the design wavelength
%            ex. "struct('sphere',1/100,'parallel',1/90,'rough',1/141)"
%            or a number equal to the defect finesse
% "n0"       Index of refraction of spacer material (air~1.0003)
% "lambda"   Center wavelength in units of [microns]
% "etendue"  System etendue in units of [sr-m^2]
% "source"   1/e half-width of the source spectrum in units of [GHz]
% "tol"      tolerance of model

% Outputs
% "T"       Transmission 0-1

% Written by John A. Smith, CIRES, University of Colorado at Boulder

e0 = e0/1e6;
Dfp = Dfp/1000;
lambda = lambda/1e6;

% fundamentals
c = 299792458; % speed of light

FSR = c/(2*n0*e0)/1e9; % free spectral range in GHz

FR = pi/(2*asin((1-R)/(2*sqrt(R)))); % reflective finesse
FA = (pi*Dfp/2)^2*(lambda/e0)*1/etendue; % aperture finesse
FL = Dfp^2/(2*e0*lambda); % diffraction-limited finesse

% read "defects" and approximate defect finesse if structure or assign value if not
if isstruct(defect)
    ds = defect.sphere; % spherical defect in lambdas
    dp = defect.parallel; % parallelism error in lambdas
    di = defect.rough; % plate roughness in lambdas
    FD = 1/sqrt(4*ds^2+3*dp^2+22*di^2); % defect finesse estimate
else
    FD = defect;
end

FE = 1/sqrt(1/FR^2+1/FD^2+1/FL^2); % effective finesse
eta = cos(pi/FE);
RE = 2-sqrt((eta-3)*(eta-1))-eta; % effective reflectance

k = 1:ceil(log(tol)/log(RE)); % determine # iterations to meet tolerance
[K,NU] = meshgrid(k,nu);

% compute model
T = (1-A/(1-R))^2*(1-R)/(1+R)*(1+FA/pi*sum(RE.^K./K.*exp(-(pi*K*source/FSR).^2) ... 
    .*(sin(2*pi*K.*(NU/FSR+1/2/FA))-sin(2*pi*K.*(NU/FSR-1/2/FA))),2));

end