function tau_R = rayleigh_OT(wavelengths_nm, surf_p_hPa, map_file)

%% function to calculate the Rayleigh scattering using ggg files for inputs
% This function is based on the work of Bodhaine et al. 1999.
% (http://dx.doi.org/10.1175/1520-0426(1999)016<1854:ORODC>2.0.CO;2)
% If finding a .map file each time is too difficult - may consider
% hard-coding values derived from them as constants...

% example inputs
% wavelengths_nm  = 250:5:1000;
% map_file = 'H:\Research\Data\TCCON\priors\Caltech\all\ci20150101.map'; 
% % surf_p_hPa = [1000:-10:950]; %list of surface pressures
% surf_p_hPa = 1013.25;

%% check inputs
if exist('surf_p_hPa', 'var')~=1;
  surf_p_hPa = 1000; %default surface P of 1000
end
if exist(map_file, 'file')~=2;
  error(['Missing .map file: ', map_file])
end

if size(wavelengths_nm,1)==1; %take transpose
  wavelengths_nm = wavelengths_nm';
end
if size(surf_p_hPa,2)==1; %take transpose
  surf_p_hPa = surf_p_hPa';
end

%% read in map file
map_raw = read_map_2012(map_file);

% calculate approximate pressure weighting function (Connor et al 2008).  Approximate because non-exact surface p used.
p_surf_wgt = map_raw.Pressure_hPa(1); %surface pressure for weighting
h_vals_aprx(1) = abs((-1*map_raw.Pressure_hPa(1) + (map_raw.Pressure_hPa(2)-...
  map_raw.Pressure_hPa(1))/log(map_raw.Pressure_hPa(2)/map_raw.Pressure_hPa(1))))/p_surf_wgt; %first value
h_vals_aprx(length(map_raw.Pressure_hPa)) = abs((map_raw.Pressure_hPa(end) - (map_raw.Pressure_hPa(end)-...
  map_raw.Pressure_hPa(end-1))/log(map_raw.Pressure_hPa(end)/map_raw.Pressure_hPa(end-1))))/p_surf_wgt; %end value

plsts = [map_raw.Pressure_hPa(1:end-2), map_raw.Pressure_hPa(2:end-1), map_raw.Pressure_hPa(3:end)]; %abbreviated variable name of pressure level listings
h_vals_aprx(2:(end-1)) = abs((-1*plsts(:,2)+(plsts(:,3)-plsts(:,2))./log(plsts(:,3)./plsts(:,2))) + ...
  (plsts(:,2)-(plsts(:,2)-plsts(:,1))./log(plsts(:,2)./plsts(:,1))))./p_surf_wgt;

aprx_xgas.co2_ppm = sum(map_raw.co2_ppm.*h_vals_aprx');  %approximate amount of CO2, from prior
g_wgt = sum(map_raw.gravity_m_s2.*h_vals_aprx')*100; %pressure weighted g in cm/s^2

%% estimate the refractive index of dry air
lambda_um = wavelengths_nm/1000;
bodhaine_eq18 = 8060.51 + 2480990./(132.274- lambda_um.^(-2)) + 17455.7./(39.32957 - lambda_um.^(-2));
bodhaine_eq19 = 1 + 0.54*(aprx_xgas.co2_ppm*1E-6 - 0.0003);

nair = bodhaine_eq18*bodhaine_eq19*1E-8+1; %refractive index of DRY air, including current CO2 value. FOR: 288.15 K, 1013.25 hPa

%% calculate depolarization factor
F_N2 = 1.034 + 3.17*1E-4*lambda_um.^(-2); %depolarization of N2
F_O2 = 1.096 + 1.385*1E-3*lambda_um.^(-2) + 1.448*1E-4*lambda_um.^(-4); %depolarization of O2
F_air_co2 = (78.084*F_N2 + 20.946*F_O2 + 0.934*1.00 + aprx_xgas.co2_ppm*1E-6*1.15)./ ...
  (78.084 + 20.946 + 0.934 + aprx_xgas.co2_ppm*1E-6); 
% rho_air = (6*F_air_co2-6)./(7*F_air_co2+3); %depolarization ratio

%% scattering cross section of air
lambda_cm = wavelengths_nm/1E7; %wavelengths in cm
N_s = 2.546899E19; %molecules/cm^-3 at 288.15 K and 1013.25 hPa
xsect = 24*pi^3*(nair.^2-1).^2./(lambda_cm.^4.*N_s.^2.*(nair.^2+2).^2).*F_air_co2; %scattering cross section of air
xsect_rep = repmat(xsect, 1, length(surf_p_hPa));

%% Rayleigh optical depths
surf_p_dyncm2 = surf_p_hPa*1000; %pressure in dyn/cm^2
surf_p_dyncm2_rep = repmat(surf_p_dyncm2, length(xsect), 1);
m_a = 15.0556*aprx_xgas.co2_ppm*1E-6 + 28.9595; %molecular weight dry air

tau_R = xsect_rep.*surf_p_dyncm2_rep*map_raw.Avogadro./(m_a.*g_wgt); 

