はじめに
*2025/04/10更新 (H,I,D)の説明の更新を行いました
以前の記事でHITRANの使い方を説明したが、当時はproxyを通さなきゃいけなかったのでいちいちraw dataをダウンロードして解析を行っていた。あちらのやり方のほうが実際のデータが見れてわかりやすいですが、hapiを使ったほうがよりほしいデータに近づくので改めて記事を書き直すことにしました。
本編
最初の使い方
hapiをDLします。以下ipynbで実行する場合
!pip install hitran-apicommand promptで実行する場合
pip install hitran-apiここで基本のmoduleは以下の様にinstallを改めてしておく
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from hapi import *本編のコードは以下の様に書ける。
db_begin()
start_wavenumber = 5000 #λ=2000nm
end_wavenumber = 25000 #λ=400nm
#data load
fetch("H2O",1,1,start_wavenumber,end_wavenumber)
#空気中に圧力比1.88%(=a)存在するとき
#1.88%は25degC,湿度60%
a = 1902 / 101325
H2O_components = [(1,1,1)]
axis_hitran, value_hitran = absorptionCoefficient_Voigt(H2O_components,"H2O",OmegaStep=0.01,HITRAN_units=False,GammaL='gamma_self',Environment={'p':1,'T':296.},Diluent={"self":a,"air":1-a})
#l's unit is cm
x,y = transmittanceSpectrum(axis_hitran, value_hitran, Environment={'l': 500.0})
def nu2lambda(nu):
lam = 10**7/nu
return lam
plt.plot(nu2lambda(x),y)
plt.xlabel("wavelength(nm)")
plt.ylabel("Transmittance (/5m)@25degC,60%")
plt.show()
このコードは25degC,湿度60%,距離5mにおける透過率を表している。
初めからコードを解説すると
・db.begin()…データを保存するためのdirectoryを指定している。
ここで何も指定していないので今回は以下の様に実行ファイルと同directoryにファイルが保存される。
・fecth()…データをDLする。このとき波長ではなく波数で取得する必要がある。
ちなみに一回fecth()した後はlocalに保存されるので追加で行う必要はない。
・a…空気中の圧力比でガスがどの程度あるかを計算。今回は以下のサイトを参考にした
化学（理想気体と実在気体）｜技術情報館「SEKIGIN」｜気体の状態方程式などの前提となっている理想化された気体と実際の気体との違いを紹介
・components…[(M,I,D)]を指定。 Mは分子番号、Iは同位体番号、Dは分子の存在比。分子番号等はHITRANに設定されているのでHITRANのこちらとにらめっこしながら記入してください。
・absorptionCoefficient_Voigt…吸収係数を求める。ここでVoigtはvoigt関数広がりでほかにも
> absorptionCoefficient_HT -> absorptionCoefficient_Voigt -> absorptionCoefficient_Lorentz -> absorptionCoefficient_Doppler -> absorptionCoefficient_SDVoigtがサポートされている。
・x,y = transmittanceSpectrum(axis_hitran, value_hitran, Environment={‘l’: 500.0})…最後に得た吸収係数から実際の透過率に計算しなおす。
実測データとの比較
db_begin()
start_wavenumber = 5000 #λ=2000nm
end_wavenumber = 25000 #λ=400nm
#data load
fetch("H2O",1,1,start_wavenumber,end_wavenumber)
#空気中に圧力比1.88%(=a)存在するとき
#1.88%は25degC,湿度60%
a = 1902 / 101325 *11.8/13.8
H2O_components = [(1,1,a)]
axis_hitran, value_hitran = absorptionCoefficient_Voigt(H2O_components,"H2O",OmegaStep=0.01,HITRAN_units=False,GammaL='gamma_self',Environment={'p':1,'T':296.})
#l's unit is cm
x,y = absorptionSpectrum(axis_hitran, value_hitran, Environment={'l': 16.8})
def nu2lambda(nu):
lam = 10**7/nu
return lam
plt.plot(nu2lambda(x),y)
plt.xlabel("wavelength(nm)")
plt.ylabel("Transmittance (/16.8mm)@25degC,60%")
plt.xlim(1300,2000)
plt.show()
上記のコードを実行すると16.8mmの透過率がわかる。実際の論文にある値と比較しても近い値になることがわかる。(Near-infrared imaging of water vapour in air – IOPscience)
追記 水蒸気の湿度をどこに代入すべきか
上記の説明で以下のような説明をしました。
・components…[(M,I,D)]を指定。 Mは分子番号、Iは同位体番号、Dは分子の存在比。分子番号等はHITRANに設定されているのでHITRANのこちらとにらめっこしながら記入してください。
これについて補足となります。こちらの説明はhapiのマニュアルをそのまま翻訳説明したものとなっています。”abondance”というのは通常同位体があるなかでの存在確率を意味しており、このパラメータに気体の質量を代入するのは一見して正しくないように見えます。
Diluent: (optional parameter)
Dictionary defining the gas mixture composition. For instance, Diluent={‘self’:A, ‘air’:B, ‘CO2’:C, …}
defines the mixture where A,B,C are volume mixing ratios of the corresponding gases (A+B+C+…=1).
一方でDilulantという項があり、こちらに体積比を代入するのが正しそうに見えます。
結論から言うと自分としてはDに体積比を代入するほうが確からしいと考えています。
この判断材料として以下の３点があります。
・Dilulantに体積比を代入しても上記に示しているように論文と一致している形状にならず、絶対値も大きく違う値になった
・Githubの大気の透過率を測定したプログラムでもDに代入していることが確認できる。https://github.com/M-Howorucha/hitran_sims/blob/main/Calculate_Spectrum.py
・ソースファイルを読む限りDに入れたほうがよさそう
hapiソースコードからの考察
下記コードがHAPI version: 1.2.2.3の吸収スペクトルを計算するときに呼び出されるコードになります。
ここで
def absorptionCoefficient_Generic(Components=None,SourceTables=None,partitionFunction=PYTIPS,
Environment=None,OmegaRange=None,OmegaStep=None,OmegaWing=None,
IntensityThreshold=DefaultIntensityThreshold,
OmegaWingHW=DefaultOmegaWingHW,
GammaL='gamma_air', HITRAN_units=True, LineShift=True,
File=None, Format=None, OmegaGrid=None,
WavenumberRange=None,WavenumberStep=None,WavenumberWing=None,
WavenumberWingHW=None,WavenumberGrid=None,
Diluent={},LineMixingRosen=False,
profile=None,calcpars=None,exclude=[],
DEBUG=None):
if profile is None: raise Exception('user must provide the line profile function')
if calcpars is None: raise Exception('user must provide the function for calculating profile parameters')
if DEBUG is not None:
VARIABLES['abscoef_debug'] = True
else:
VARIABLES['abscoef_debug'] = False
exclude = set(exclude)
if not LineMixingRosen: exclude.add('YRosen')
if not LineShift: exclude.update({'Delta0','Delta2'})
if WavenumberRange is not None: OmegaRange=WavenumberRange
if WavenumberStep is not None: OmegaStep=WavenumberStep
if WavenumberWing is not None: OmegaWing=WavenumberWing
if WavenumberWingHW is not None: OmegaWingHW=WavenumberWingHW
if WavenumberGrid is not None: OmegaGrid=WavenumberGrid
Components = listOfTuples(Components)
SourceTables = listOfTuples(SourceTables)
Components,SourceTables,Environment,OmegaRange,OmegaStep,OmegaWing,\
IntensityThreshold,Format = \
getDefaultValuesForXsect(Components,SourceTables,Environment,OmegaRange,
OmegaStep,OmegaWing,IntensityThreshold,Format)
if OmegaStep>0.005 and profile is PROFILE_DOPPLER:
warn('Big wavenumber step: possible accuracy decline')
elif OmegaStep>0.1:
warn('Big wavenumber step: possible accuracy decline')
if OmegaGrid is not None:
Omegas = npsort(OmegaGrid)
else:
Omegas = arange_(OmegaRange[0],OmegaRange[1],OmegaStep)
number_of_points = len(Omegas)
Xsect = zeros(number_of_points)
T_ref_default = __FloatType__(296.)
p_ref_default = __FloatType__(1.)
T = Environment['T']
p = Environment['p']
ABUNDANCES = {}
NATURAL_ABUNDANCES = {}
for Component in Components:
M = Component[0]
I = Component[1]
if len(Component) >= 3:
ni = Component[2]
else:
try:
ni = ISO[(M,I)][ISO_INDEX['abundance']]
except KeyError:
raise Exception('cannot find component M,I = %d,%d.' % (M,I))
ABUNDANCES[(M,I)] = ni
NATURAL_ABUNDANCES[(M,I)] = ISO[(M,I)][ISO_INDEX['abundance']]
if HITRAN_units:
factor = __FloatType__(1.0)
else:
factor = volumeConcentration(p,T)
GammaL = GammaL.lower()
if not Diluent:
if GammaL == 'gamma_air':
Diluent = {'air':1.}
elif GammaL == 'gamma_self':
Diluent = {'self':1.}
else:
raise Exception('Unknown GammaL value: %s' % GammaL)
print(Diluent)
for key in Diluent:
val = Diluent[key]
if val < 0 or val > 1:
raise Exception('Diluent fraction must be in [0,1]')
t = time()
CALC_INFO_TOTAL = []
for TableName in SourceTables:
DATA_DICT = LOCAL_TABLE_CACHE[TableName]['data']
parnames_exclude = ['a','global_upper_quanta','global_lower_quanta',
'local_upper_quanta','local_lower_quanta','ierr','iref','line_mixing_flag']
parnames = set(DATA_DICT)-set(parnames_exclude)
nlines = len(DATA_DICT['nu'])
for RowID in range(nlines):
TRANS = CaselessDict({parname:DATA_DICT[parname][RowID] for parname in parnames})
TRANS['T'] = T
TRANS['p'] = p
TRANS['T_ref'] = T_ref_default
TRANS['p_ref'] = p_ref_default
TRANS['Diluent'] = Diluent
TRANS['Abundances'] = ABUNDANCES
if (TRANS['molec_id'],TRANS['local_iso_id']) not in ABUNDANCES: continue
TRANS['SigmaT'] = partitionFunction(TRANS['molec_id'],TRANS['local_iso_id'],TRANS['T'])
TRANS['SigmaT_ref'] = partitionFunction(TRANS['molec_id'],TRANS['local_iso_id'],TRANS['T_ref'])
LineIntensity = calculate_parameter_Sw(None,TRANS)
if LineIntensity < IntensityThreshold: continue
CALC_INFO = {}
PARAMETERS = calcpars(TRANS=TRANS,CALC_INFO=CALC_INFO,exclude=exclude)
try:
GammaD = PARAMETERS['GammaD']
except KeyError:
GammaD = 0
try:
Gamma0 = PARAMETERS['Gamma0']
except KeyError:
Gamma0 = 0
GammaMax = max(Gamma0,GammaD)
if GammaMax==0 and OmegaWingHW==0:
OmegaWing = 10.0
warn('Gamma0 and GammaD are missing; setting OmegaWing to %f cm-1'%OmegaWing)
OmegaWingF = max(OmegaWing,OmegaWingHW*GammaMax)
BoundIndexLower = bisect(Omegas,TRANS['nu']-OmegaWingF)
BoundIndexUpper = bisect(Omegas,TRANS['nu']+OmegaWingF)
PARAMETERS['WnGrid'] = Omegas[BoundIndexLower:BoundIndexUpper]
lineshape_vals = profile(**PARAMETERS)
Xsect[BoundIndexLower:BoundIndexUpper] += factor * lineshape_vals
if VARIABLES['abscoef_debug']: DEBUG.append(CALC_INFO)
print('%f seconds elapsed for abscoef; nlines = %d'%(time()-t,nlines))
if File: save_to_file(File,Format,Omegas,Xsect)
return Omegas,Xsect
上記コードを読み解いていきます。
最終的に出力されるデータのうち、Xsectが吸収スペクトルの強度となります。
その最終更新部分は
Xsect[BoundIndexLower:BoundIndexUpper] += factor * lineshape_valsとなっております。このうちfactor項は以下のように正確な1.0(__FloatType__(1.0))または圧力温度(volumeConcentration(p,T))によって一意的に決まる
# actual temperature and pressure
T = Environment['T'] # K
p = Environment['p'] # atm
# pre-calculation of volume concentration
if HITRAN_units:
factor = __FloatType__(1.0)
else:
factor = volumeConcentration(p,T)よって議論すべきはlineshape_valsとなりますが、これを深堀すると
TRANS['SigmaT'] = partitionFunction(TRANS['molec_id'],TRANS['local_iso_id'],TRANS['T'])
TRANS['SigmaT_ref'] = partitionFunction(TRANS['molec_id'],TRANS['local_iso_id'],TRANS['T_ref'])
LineIntensity = calculate_parameter_Sw(None,TRANS)
if LineIntensity < IntensityThreshold: continue
# calculate profile parameters
CALC_INFO = {}
PARAMETERS = calcpars(TRANS=TRANS,CALC_INFO=CALC_INFO,exclude=exclude)
PARAMETERS['WnGrid'] = Omegas[BoundIndexLower:BoundIndexUpper]
lineshape_vals = profile(**PARAMETERS)以上の計算まとまりのうち、calcpars関数内で"LineIntensity"が非明示的に使用されています。
calcparsには各広がりに応じた関数が代入されVoigt広がりの場合は以下のようになります。
(calculate_parameter_Sw~LinIntesityがつかわれていることがわかります。)
def calculateProfileParametersVoigt(TRANS,CALC_INFO=None,exclude=None):
"""
Get values for abstract profile parameters for Voigt profile.
"""
envdep_presets = [('Voigt','default')]
parameters = [
('Nu', calculate_parameter_Nu),
('Sw', calculate_parameter_Sw),
('GammaD', calculate_parameter_GammaD),
('Gamma0', calculate_parameter_Gamma0),
('Delta0', calculate_parameter_Delta0),
('YRosen', calculate_parameter_YRosen),
]
return calculateProfileParameters(envdep_presets,parameters,CALC_INFO=CALC_INFO,TRANS=TRANS,exclude=exclude)そのLineIntensityをつかさどるcalculate_parameter_Sw関数をのぞいてみると以下のようになっています。
def calculate_parameter_Sw(dummy,TRANS,CALC_INFO=None):
molec_id = TRANS['molec_id']
local_iso_id = TRANS['local_iso_id']
nu = TRANS['nu']
sw = TRANS['sw']
T = TRANS['T']
Tref = TRANS['T_ref']
SigmaT = TRANS['SigmaT']
SigmaTref = TRANS['SigmaT_ref']
elower = TRANS['elower']
Sw_calc = EnvironmentDependency_Intensity(sw,T,Tref,SigmaT,SigmaTref,elower,nu)
if 'Abundances' in TRANS:
Sw_calc *= TRANS['Abundances'][(molec_id,local_iso_id)]/abundance(molec_id,local_iso_id)
if type(CALC_INFO) is dict:
CALC_INFO['Sw'] = {
'value':Sw_calc,
'mixture':{
'generic': {
'args': {
'SigmaT':{'value':SigmaT,'source':''},
'SigmaT_ref':{'value':SigmaTref,'source':''},
'elower':{'value':elower,'source':'elower'},
'nu':{'value':nu,'source':'nu'},
}
}
}
}
return Sw_calc これを要素分解すると
Sw_calc = EnvironmentDependency_Intensity(sw,T,Tref,SigmaT,SigmaTref,elower,nu)
if 'Abundances' in TRANS:
Sw_calc *= TRANS['Abundances'][(molec_id,local_iso_id)]/abundance(molec_id,local_iso_id)
以上2点が最終的な吸収強度に寄与していることがわかります。
これを読んでみると確かにAboundanceが"aboundanceに代入された数値/通常時の存在比"として掛け算されていることがわかり、
吸収係数は吸収断面積×分子数の線形な関係があることからaboundanceに代入するというのは間違っていない計算結果を導くことがわかります。
コメント