clear
close all
% read HITRAN data, compute broadbening and aborption coefficient
% computes mass absorption coefficient k=S*f(ni-ni0)
% see http://www.forbrf.lth.se/fileadmin/forbrf/Documents/Exjobb/Old/2011_Zhou_Thesis.pdf
% HITRAN parameters https://hitran.org/docs/definitions-and-units/

%filename='5e79029c.par';
filename='5e823169.par';
filename='5e824189.par';
data=importhitran(filename);
Rs=8.13145;   % gas constant
Na=6.022e23;  % Avagadro.s principle
R=287; %Specific Gas Constant
MoleWeight=18/1000;    % H20
MoleWeightAir=28.96/1000;
R=Rs/MoleWeightAir; %Specific Gas Constant for air
Pressure=1000; %hPa
q=10/1000;   % speific humidity  [kg/kg]
PW=1.5;      % Total content of specific gas at standard conditions
epsilon=MoleWeight/(MoleWeightAir);
PartitialPressure=Pressure*q/epsilon;
Temperature=300;  
Temperature_ref=296;   % Proporties in HITRAN
dens=PartitialPressure/(287/epsilon*Temperature);
N_molec=dens/MoleWeight*Na;     % numer of molecules per m3
% Selection of specific molecule (1-H20,2-CO2,3 O3..) and isotop
ind=find(data.isotopologueNumber==1&data.moleculeNumber==1);
ind=find(data.isotopologueNumber==1&data.moleculeNumber==2);

LineStrength_T=data.lineIntensity(ind);
Wavenumber=data.transitionWavenumber(ind);
AirHWHM=data.airBroadenedWidth(ind);
SelfHWHM=data.selfBroadenedWidth(ind);
Delta_air=data.pressureShift(ind); 
n_air=data.temperatureDependence(ind);
RangeLow=Wavenumber(1);
RangeHigh=Wavenumber(end);
Pressure=Pressure/1013.25;    % Normalization to 1 atm
PartitialPressure=PartitialPressure/1013.25;
[Wavenumber_shift,HWHM_Lorentz,HWHM_Doppler]=LineCorrection(Wavenumber,Delta_air,n_air,AirHWHM,SelfHWHM,Pressure,PartitialPressure,Temperature,Temperature_ref,MoleWeight);
[SumStrength_V,wav]=LineShape(LineStrength_T,HWHM_Doppler,HWHM_Lorentz,Wavenumber_shift,RangeLow,RangeHigh);
lambda=10000./wav;
%SumStrength= S*f(v-vo)=ABS[cm^2/molecules]
C=1e-4;
ABS=SumStrength_V*N_molec*C;     % 1/m
Kabs=ABS./dens;                % mass absroptio coeffient m2/kg
% Optical depth for constant (temperature, pressure...)
TAU=Kabs*PW*10;     % 10 from unit change

figure
plot(Wavenumber,LineStrength_T)
hold on
plot(wav,SumStrength_V,'r')

figure
plot(lambda,ABS)
xlabel('Wavelength [\mum]')
ylabel('Absorption coefficient [1/m]')