Integrate ESTAR density effect corrections into EGSnrc

[rtownson]

The ESTAR fortran code from NIST has been translated into C++ and integrated into EGSnrc, so that density effect corrections are now calculated on-the-fly during the initialization of a simulation. Materials can be defined directly in the input file, so .pegs4dat and .density (density correction files) are no longer needed.

Example:

:start media definition:
  :start water:
    elements = H, O
    number of atoms = 2, 1
    bulk density = 1.0 # g/cm^3 used for transport
    local density = 1.0 # g/cm^3 used for density corrections
    ivalue = 0.000075 # MeV, optional
    output density file = /mypath/water.density # full path to output density corrections to file, optional
  :stop water:
:stop media definition:

For mixtures of compounds:

    :start saltWater:
        mixture compounds = H2O, NaCl
        mass fractions = 0.95, 0.05
        bulk density = 1.0
    :stop saltWater:

For comparison, there's an option to output a density file when the simulation is run.

This adds the options for bulk and local densities to be specified separately. The bulk density is the density used for transport, and the local density is used for the density effect corrections (when a density file is not provided). Density correction files, if provided, will always override these inputs. If the local density is not provided, bulk density (or rho) is used for everything. The old rho input is still accepted, and if provided is used both as the bulk and local density.

Thanks to our co-op student @Sehmimul for this significant contribution!

[ftessier]

This will be merged in develop after 2022 release.

[ftessier]

Is there an option in the integrated estar code to compute density corrections for a density that is different from the bulk density, e.g., for the grain density in graphite? Perhaps we can have inputs for bulk density (which is simply passed through, and local density which is used for the density correction calculation?

[rtownson]

Is there an option in the integrated estar code to compute density corrections for a density that is different from the bulk density, e.g., for the grain density in graphite? Perhaps we can have inputs for bulk density (which is simply passed through, and local density which is used for the density correction calculation?

I like your suggestion, yes we could add those options. Currently density is just used for both, so you can't have a different bulk density.

[rtownson]

I'm also planning to change the I-value to be specified in MeV, like all other energy parameters..

[ftessier]

@rtownson What is the status of this draft? Should I start working on this branch towards merging, or should I hold off?

[rtownson]

I think hold off, I'm planning to add a few more features.

[ftessier]

Sounds good. I will start reviewing and "cleaning" it on the side, please add stuff on top of the existing estar-cpp-integration f21697d, so that I can integrate your new commits cleanly at the end. Don't forget to update to f21697d, I synced it with develop last week.

[ftessier]

Note that I started cleaning up this branch on the side on branch cleanup-estar-cpp-integration. CLeanup strategy is to keep a null diff between the two branches:

https://github.com/nrc-cnrc/EGSnrc/compare/estar-cpp-integration..cleanup-estar-cpp-integration

[rtownson]

Currently the gas pressure is applied only when a density correction file is not used. It has only one line of action, in the subroutine MIX which "calculates MS parameters", it modifies the density:

IF (GASPP.NE.0.0) [RHOP=GASPP*RHOP;]

Is there a reason that the density itself is not just immediately modified after arguments are read in? For the ESTAR implementation, this would mean that the density change from gas pressure is not considered in the density corrections (just like it is ignored when you use a density correction file). Since a user can just modify the density directly in-line, why have the gas pressure separated?

[rtownson]

This PR is ready for review / testing :)

[rtownson]

There seems to be an outstanding bug where elements input with 2 letters (e.g. Al) results in a crash. Mixture compounds with 2 letter elements work fine.