% Interactive routine for computing line absorption coefficient for
% molecules in the atmosphere.  Operates in physical
% variables.  Karp's method (JQSRT 20, 379, 1978) is used.
%
% Syntax: [kr, xnu] = abscof
%
% by D. Holmgren holmgren@phobos.astro.uwo.ca

function [kr, xnu] = abscof
% Note that both coefficient and wavenumber scale are returned.
% input data...
% generate wavenumber scale...
xlo = input('Lower wavenumber limit: ');
xhi = input('Upper wavenumber limit: ');
nnu = input('Number of points over interval: ');
xnu = linspace(xlo,xhi,nnu)';
% get physical parameters...
T = input('Temperature (K): ');
p = input('Pressure (atm): ');
% Lorentz width...
dl = 0.1 * p;
nu0 = input('Line center: ');
% Doppler width...
dnu = 6.e-08 * nu0 * sqrt(T);
disp(' Doppler width: ')
disp( dnu )
disp(' Lorentz width: ')
disp( dl )
sk = input('Line strength: ');
% damping parameter...
a = dl / dnu;
c = 2/pi;
twopi = 2 * pi;
c1 = twopi * dnu;
c2 = twopi * nu0;
% need correct t-scale...
dt = 0.5 / ( xnu(nnu) - xnu(1) );
deltanu = ( xnu(nnu) - xnu(1) ) / nnu;
ndt = 0.5 / deltanu;
t = (0:dt:ndt)';
% generate Voigt profile transform...
f = exp(-(a * c1) .* t - (c1 * c1) .* t .* t ./ 4);
fc = f .* cos(c2 .* t) .* sk;
fs = f .* sin(c2 .* t) .* sk;
% compute transform...
kr = fft(fc,nnu);
ki = fft(fs,nnu);
% note form of knu - due to variety of FFT used.
knu = c .* (real(kr) - imag(ki));
% kr = absorption coefficient.
kr = real(knu);
plot(xnu,kr);
title('Atmospheric Absorption Spectrum');
xlabel('Frequency (wavenumbers)');
ylabel('Absorption Coefficient');
shg