Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
spec: Spectrum
'''
unit = None
def get_radiance_unit(unit_emisscoeff):
''' get radiance_noslit unit from emisscoeff unit'''
if '/cm3' in unit_emisscoeff:
return unit_emisscoeff.replace('/cm3', '/cm2')
else:
return unit_emisscoeff + '*cm'
# case where we recomputed it already (somehow... ex: no_change signaled)
if 'radiance_noslit' in rescaled:
if __debug__:
printdbg('... rescale: radiance_noslit was scaled already')
assert 'radiance_noslit' in units
return rescaled, units
# Rescale!
if 'emisscoeff' in rescaled and true_path_length and optically_thin:
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2 * L2 ' +
'(optically thin)')
emisscoeff = rescaled['emisscoeff'] # x already scaled
radiance_noslit = emisscoeff*new_path_length # recalculate L
unit = get_radiance_unit(units['emisscoeff'])
elif ('emisscoeff' in rescaled and 'transmittance_noslit' in rescaled
and 'abscoeff' in rescaled and true_path_length and not optically_thin): # not optically thin
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2*(1-T2)/k2')
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2*(1-T2)/k2')
emisscoeff = rescaled['emisscoeff'] # x already scaled
abscoeff = rescaled['abscoeff'] # x already scaled
# mole_fraction, path_length already scaled
transmittance_noslit = rescaled['transmittance_noslit']
b = (transmittance_noslit == 1) # optically thin mask
radiance_noslit = np.empty_like(emisscoeff) # calculate L
radiance_noslit[~b] = emisscoeff[~b] / abscoeff[~b]*(1-transmittance_noslit[~b])
radiance_noslit[b] = emisscoeff[b] * new_path_length # optically thin limit
unit = get_radiance_unit(units['emisscoeff'])
elif ('emisscoeff' in rescaled and 'abscoeff' in rescaled and true_path_length
and not optically_thin): # not optically thin
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2*(1-exp(-k2*L2))/k2')
emisscoeff = rescaled['emisscoeff'] # x already scaled
abscoeff = rescaled['abscoeff'] # x already scaled
b = (abscoeff == 0) # optically thin mask
radiance_noslit = np.empty_like(emisscoeff) # calculate
radiance_noslit[~b] = emisscoeff[~b]/abscoeff[~b] *(1-exp(-abscoeff[~b]*new_path_length))
radiance_noslit[b] = emisscoeff[b] * new_path_length # optically thin limit
unit = get_radiance_unit(units['emisscoeff'])
elif 'radiance_noslit' in initial and optically_thin:
if __debug__:
printdbg('... rescale: radiance_noslit I2 = I1*N2/N1*L2/L1 ' +
'(optically thin)')
_, radiance_noslit = spec.get(
'radiance_noslit', wunit=waveunit, Iunit=units['radiance_noslit'])
radiance_noslit *= new_mole_fraction / old_mole_fraction # rescale
radiance_noslit *= new_path_length / old_path_length # rescale
assert 'radiance_noslit' in units
return rescaled, units
# Rescale!
if 'emisscoeff' in rescaled and true_path_length and optically_thin:
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2 * L2 ' +
'(optically thin)')
emisscoeff = rescaled['emisscoeff'] # x already scaled
radiance_noslit = emisscoeff*new_path_length # recalculate L
unit = get_radiance_unit(units['emisscoeff'])
elif ('emisscoeff' in rescaled and 'transmittance_noslit' in rescaled
and 'abscoeff' in rescaled and true_path_length and not optically_thin): # not optically thin
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2*(1-T2)/k2')
emisscoeff = rescaled['emisscoeff'] # x already scaled
abscoeff = rescaled['abscoeff'] # x already scaled
# mole_fraction, path_length already scaled
transmittance_noslit = rescaled['transmittance_noslit']
b = (transmittance_noslit == 1) # optically thin mask
radiance_noslit = np.empty_like(emisscoeff) # calculate L
radiance_noslit[~b] = emisscoeff[~b] / abscoeff[~b]*(1-transmittance_noslit[~b])
radiance_noslit[b] = emisscoeff[b] * new_path_length # optically thin limit
unit = get_radiance_unit(units['emisscoeff'])
elif ('emisscoeff' in rescaled and 'abscoeff' in rescaled and true_path_length
and not optically_thin): # not optically thin
if __debug__:
printdbg('... rescale: radiance_noslit I2 = j2*(1-exp(-k2*L2))/k2')
emisscoeff = rescaled['emisscoeff'] # x already scaled
abscoeff = rescaled['abscoeff'] # x already scaled
assert 'transmittance_noslit' in units
return rescaled, units
# Calculate rescaled value directly
# ---------------------------------
# Rescale
if 'absorbance' in rescaled:
if __debug__:
printdbg('... rescale: transmittance_noslit T2 = exp(-A2)')
absorbance = rescaled['absorbance'] # x and L already scaled
transmittance_noslit = exp(-absorbance) # recalculate
unit = get_unit()
elif 'abscoeff' in rescaled and true_path_length:
if __debug__:
printdbg('... rescale: transmittance_noslit T2 = exp(-j2*L2)')
abscoeff = rescaled['abscoeff'] # x already scaled
absorbance = abscoeff*new_path_length # calculate
transmittance_noslit = exp(-absorbance) # recalculate
unit = get_unit()
elif 'transmittance_noslit' in initial:
if __debug__:
printdbg('... rescale: transmittance_noslit T2 = ' +
'exp( ln(T1) * N2/N1 * L2/L1)')
# get transmittance from initial transmittance
_, T1 = spec.get('transmittance_noslit', wunit=waveunit,
Iunit=units['transmittance_noslit'])
# We'll have a problem if the spectrum is optically thick
b = (T1 == 0) # optically thick mask
if b.sum() > 0 and (new_mole_fraction < old_mole_fraction or
new_path_length < old_path_length):
Default ``True``
Other Parameters
----------------
assume_equilibrium: boolean
if ``True``, only absorption coefficient ``abscoeff`` is recomputed
and all values are derived from a calculation under equilibrium,
using Kirchoff's Law. Default ``False``
'''
optically_thin = spec.is_optically_thin()
initial = spec.get_vars() # quantities initially in spectrum
if __debug__:
printdbg('... rescale: optically_thin: {0}'.format(optically_thin))
printdbg('... rescale: initial quantities: {0}'.format(initial))
# Check that inputs are valid names
_check_quantity = quantity
if isinstance(quantity, list):
_check_quantity = quantity
else:
_check_quantity = [quantity]
for k in _check_quantity:
assert k in CONVOLUTED_QUANTITIES + \
NON_CONVOLUTED_QUANTITIES + ['all', 'same']
# ... make sure we're not trying to rescale a Spectrum that has non scalable
# ... quantities
if any_in(initial, non_rescalable_keys):
raise NotImplementedError('Trying to rescale a Spectrum that has non scalable ' +
'quantities: {0}'.format([k for k in initial if k in non_rescalable_keys]) +
if greedy:
# ... let's be greedy: recompute all possible quantities. The list of
# all spectral quantities is calculated by parsing a tree in get_reachable
reachable = get_reachable(spec)
extra = [k for k, v in reachable.items() if v]
wanted = set(wanted+extra)
# There are two cases: either we are actually rescaling to another length /
# mole fraction, or we are just updating() without changing length / mole fraction
no_change = (new_mole_fraction ==
old_mole_fraction and new_path_length == old_path_length)
# Quickly stop if no change
if no_change and all_in(wanted,initial):
if __debug__:
printdbg('... rescale: no change')
# Stop here
return
# list of quantities that are needed to recompute what we want
# ... (we're just analysing how to compute them here, the actual calculation
# ... will be done laters)
try:
recompute = get_recompute(
spec, wanted, no_change, true_path_length=true_path_length)
except KeyError as err:
import sys
print(sys.exc_info())
raise KeyError('Error in get_recompute (see above). Quantity `{0}` cannot be recomputed '.format(
err.args[0])+'from available quantities in Spectrum ({0}) with '.format(spec.get_vars()) +
' conditions: optically thin ({0}), true_path_length ({1}), thermal equilibrium ({2})'.format(
optically_thin, true_path_length, assume_equilibrium) +
----------
df: DataFrame
list of transitions
Other Parameters
----------------
calc_Evib_harmonic_anharmonic: boolean
if ``True``, calculate harmonic and anharmonic components of
vibrational energies (for Treanor distributions)
'''
if __debug__:
printdbg(
'called _add_EvibErot_CDSD_pcJN(calc_Evib_harmonic_anharmonic={0})'.format(calc_Evib_harmonic_anharmonic))
if calc_Evib_harmonic_anharmonic:
raise NotImplementedError
molecule = self.input.molecule
state = self.input.state # electronic state
# TODO: for multi-molecule mode: add loops on molecules and states too
assert molecule == 'CO2'
if self.verbose>=2:
printg('Fetching vib / rot energies for all {0} transitions'.format(len(df)))
t0 = time()
# Check energy levels are here
for iso in self._get_isotope_list(molecule):
derives_from('transmittance', ['transmittance_noslit'])
derives_from('emissivity', ['emissivity_noslit'])
if equilibrium:
if __debug__:
printdbg('... build_graph: equilibrium > all keys derive from one')
# Anything can be recomputed from anything
for key in all_keys:
if key in NON_CONVOLUTED_QUANTITIES:
all_but_k = [[k] for k in all_keys if k != key]
derives_from(key, *all_but_k)
# ------------------------------------------------------------
if __debug__:
printdbg('_build_update_graph: dependence/equivalence tree:')
printdbg(derivation)
return derivation
Evib_last = ElecState.Erovib(vmax, 0)
# difference of the 2 last vib level for which Dunham expansion is valid
delta_E_last = Evib_last - ElecState.Erovib(vmax-1, 0)
v_inc = ElecState.get_Morse_inc()
# ... Start loop on all Morse potential levels
Evib = Evib_last
for v in range(vmax+1, vmax_morse+1):
delta_E = delta_E_last-(v+1-vmax)*v_inc
Evib = Evib + delta_E
if Evib > Ediss:
warn('Energy above dissociation threshold: {0}>{1}'.format(
Evib, Ediss))
if __debug__:
printdbg('Calculating Evib for ' +
'v={0}: {1:.2f}cm-1 (Morse Potential)'.format(v, Evib))
# ... Rotational loop
J = 0
E = Evib
# (passes if Jmax is nan)
while 0 <= E < Ediss and not J >= Jmax:
# store rovib level
levels.append([v, J, E, Evib])
# calculate new one:
J += 1
# no Zero-point-energy
Erot = ElecState.Erovib(0, J, remove_ZPE=True)
E = Evib + Erot
Jmaxcalc = max(Jmaxcalc, J)
df = pd.DataFrame(levels, columns=['v', 'j', 'E', 'Evib'])
def _add_Evib123Erot_RADIS_cls5_harmonicanharmonic(self, df):
''' Fetch Evib & Erot in dataframe for HITRAN class 5 (linear triatomic
with Fermi degeneracy... i.e CO2 ) molecules
Parameters
----------
df: DataFrame
list of transitions
'''
if __debug__:
printdbg(
'called _add_Evib123Erot_RADIS_cls5_harmonicanharmonic()')
molecule = self.input.molecule
state = self.input.state # electronic state
# TODO: for multi-molecule mode: add loops on molecules and states too
if self.verbose>=2:
printg('Fetching vib / rot energies for all {0} transitions'.format(len(df)))
t0 = time()
# Check energy levels are here
for iso in self._get_isotope_list(molecule):
if not iso in self.parsum_calc[molecule]:
raise AttributeError('No Partition function calculator defined for isotope {0}'.format(iso)
+ '. You need energies to calculate a non-equilibrium spectrum!'
+ ' Fill the levels parameter in your database definition, '