I think I have found a bug with how HITRAN line intensities are weighted for different isotopologues. HITRAN line intensities at the reference temperature (298 K) are already weighted by the terrestrial isotopologue abundance. However, HELIOS-K will weight by abundance again in host.h:
|
Abundance *= m.ISO[j].m / meanMass; |
So, in other words, I think HELIOS-K is effectively weighting by abundance^2, which will significantly under predict the relative line intensity of non-main isotopologues (because a minor abundace of 1e-2 will be squared resulting in a weight of 1e-4)
I first realized there was an issue because I computed some CO2 opacities weighted by terrestrial isotopologue abundances. However, the opacities seemed to under predict the contribution of non-main isotopologues.
I could not figure out the issue, until I ran HELIOS-K through the Jules AI tool, and it discovered the problem. I will paste some of our conversation below. Note that Jules also found a related bug for when HELIOS-K is computing opacities for a single HITRAN isotopologue.
Conversation
Nick:
I have done some calculations recently with HELIOS-K to produce opacities. When I did these, I did not pass the -ISO flag, so apparently there should not be an issue. However, when compare these opacities to other opacities computed with different programs, the HELIOS-K opacities seems to under-estimate the contribution of the non-main isotopologue. The main isotopologue appears to have the correct opacity. But by eye, the other more minor isotopologues have an underesimtated opacity. I am suspicious that this has to do with some isotope weighting problem in the code. Can you find no issue?
Jules:
Thank you for that crucial information. You are right to be suspicious. Based on what you've described, I re-examined the code for the case where the -ISO flag is not used, and I have found a second, more subtle bug that explains exactly what you're seeing.
My initial analysis was incomplete, and I apologize for that. Here is the new finding:
The Issue with Minor Isotopologues (no -ISO flag):
The problem is in the readFile function in host.h. When calculating the final line strength, the code incorrectly uses the isotopic abundance (Abundance) in the scaling. The line strength S read from the data file is already proportional to the natural abundance, and the code effectively multiplies it by a factor that includes the abundance again. This results in the final opacity being proportional to the square of the isotopic abundance (abundance^2).
- For the main isotopologue, the abundance is close to 1, so abundance^2 is also close to 1. The error is small, which is why it appears correct by eye.
- For a minor isotopologue, the abundance is a small number (e.g., 0.01). Squaring it makes it much smaller (e.g., 0.0001), causing a significant underestimation of its opacity. This perfectly matches your observation.
Two Separate Bugs
This means there are actually two separate bugs in the code:
- The one we just discussed in host.h that affects calculations with default natural abundances.
- The original issue I found in hitran.cpp that only affects calculations when you specify a single isotopologue with the -ISO flag.
I can fix both of these issues. The plan would be:
- Correct the calculation in host.h to remove the incorrect scaling factor that causes the abundance^2 dependency.
- Apply the original fix I proposed to hitran.cpp to handle the -ISO flag case correctly.
This should resolve both the underestimation of minor isotopologues and the issue with single-isotopologue calculations. Would you like me to proceed with fixing both bugs?
I think I have found a bug with how HITRAN line intensities are weighted for different isotopologues. HITRAN line intensities at the reference temperature (298 K) are already weighted by the terrestrial isotopologue abundance. However, HELIOS-K will weight by abundance again in host.h:
HELIOS-K/host.h
Line 629 in 1245778
So, in other words, I think HELIOS-K is effectively weighting by abundance^2, which will significantly under predict the relative line intensity of non-main isotopologues (because a minor abundace of 1e-2 will be squared resulting in a weight of 1e-4)
I first realized there was an issue because I computed some CO2 opacities weighted by terrestrial isotopologue abundances. However, the opacities seemed to under predict the contribution of non-main isotopologues.
I could not figure out the issue, until I ran HELIOS-K through the Jules AI tool, and it discovered the problem. I will paste some of our conversation below. Note that Jules also found a related bug for when HELIOS-K is computing opacities for a single HITRAN isotopologue.
Conversation
Nick:
I have done some calculations recently with HELIOS-K to produce opacities. When I did these, I did not pass the -ISO flag, so apparently there should not be an issue. However, when compare these opacities to other opacities computed with different programs, the HELIOS-K opacities seems to under-estimate the contribution of the non-main isotopologue. The main isotopologue appears to have the correct opacity. But by eye, the other more minor isotopologues have an underesimtated opacity. I am suspicious that this has to do with some isotope weighting problem in the code. Can you find no issue?
Jules:
Thank you for that crucial information. You are right to be suspicious. Based on what you've described, I re-examined the code for the case where the -ISO flag is not used, and I have found a second, more subtle bug that explains exactly what you're seeing.
My initial analysis was incomplete, and I apologize for that. Here is the new finding:
The Issue with Minor Isotopologues (no -ISO flag):
The problem is in the readFile function in host.h. When calculating the final line strength, the code incorrectly uses the isotopic abundance (Abundance) in the scaling. The line strength S read from the data file is already proportional to the natural abundance, and the code effectively multiplies it by a factor that includes the abundance again. This results in the final opacity being proportional to the square of the isotopic abundance (abundance^2).
Two Separate Bugs
This means there are actually two separate bugs in the code:
I can fix both of these issues. The plan would be:
This should resolve both the underestimation of minor isotopologues and the issue with single-isotopologue calculations. Would you like me to proceed with fixing both bugs?