Atmosphere not updating after coordinates are set

[wjwoodley]

I am trying to set custom coordinates using the set_location_coord() function in the cNRLMSISE00 class, and it looks like the atmosphere is not being updated, even when _calculate_integration_path() is called with force = True. The code below produces the same results for both sets of coordinates and months in the loop.

# Import packages

import numpy as np

from MCEq.core import MCEqRun
import mceq_config
import crflux.models as pm
from MCEq.geometry.nrlmsise00_mceq import cNRLMSISE00

# Configure MCEq

mceq_config.debug_level = 2
mceq_config.enable_default_tracking = False

msis_obj = cNRLMSISE00()
mceq_run = MCEqRun(interaction_model = "SIBYLL-2.3c", primary_model = (pm.GlobalSplineFitBeta, None), theta_deg = 0.0, **mceq_config.config)

# Loop over coordinates and months

coordinates = ((-100, 50), (110, 30))
months = ("January", "July")

for (coords, month) in zip(coordinates, months):
    
    longitude, latitude = coords
    
    # Change the atmosphere
    
    msis_obj.set_location_coord(longitude, latitude)
    msis_obj.set_season(month)
    mceq_run.set_theta_deg(0)
    mceq_run._calculate_integration_path(int_grid = None, grid_var = "X", force = True)
    
    # Run MCEq
    
    mceq_run.solve()
    
    print(mceq_run.get_solution("total_mu+", mag = 0) + mceq_run.get_solution("total_mu-", mag = 0))

[afedynitch]

Hi,

while were working on solving this in a more general fashion, there is one problem in your code that prevents these settings to propagate. You're creating an independent object of the atmosphere class, which is not the one used by your MCEq instance msis_obj = cNRLMSISE00(). Instead of creating the object, you need to use the MSIS object from MCEq:

msis_obj = mceq_run.density_profile._msis

and then proceed with the calculation.

[jncots]

Hi,

with the update in cfb08ab your example should work if the following changes are applied:

• Change a density model because coordinates can be changed only in MSIS00:
  mceq_config.density_model = ("MSIS00", ("SouthPole", "January"))

• get the reference to the object of the density model used by the mceq_run object
  density_model_MSIS00 = mceq_run.density_model

• Use new wrapper methods (around methods of _msis object):
  density_model_MSIS00.update_parameters(location_coord = (longitude, latitude), season = month)
  or alternatively:

density_model_MSIS00.set_location_coord(longitude, latitude)
density_model_MSIS00.set_season(month)

If the methods of _msis object are used instead of new wrapper function, put the line:

density_model_MSIS00.theta_deg = None

to force recalculation of the density profile.

The modification of your example that should work:

from MCEq.core import MCEqRun
import mceq_config
import crflux.models as pm

# Configure MCEq

mceq_config.debug_level = 0
mceq_config.enable_default_tracking = False


# coordinates can be changed only in MSIS00:  
mceq_config.density_model = ("MSIS00", ("SouthPole", "January"))
mceq_run = MCEqRun(interaction_model = "SIBYLL-2.3c", 
                   primary_model = (pm.GlobalSplineFitBeta, None), 
                   theta_deg = 0.0, **mceq_config.config)

# get the density model used by the mceq_run object
density_model_MSIS00 = mceq_run.density_model


# Loop over coordinates and months

coordinates = ((-100, 50), (110, 30))
months = ("January", "July")

results = []

for (coords, month) in zip(coordinates, months):
    
    longitude, latitude = coords
    
    # Change the atmosphere
    
    # The parameters can be updated this way:
    density_model_MSIS00.update_parameters(location_coord = (longitude, latitude), 
                                           season = month,
                                           doy = 219)
    
    # Or this way:
    # density_model_MSIS00.set_location_coord(longitude, latitude)
    # density_model_MSIS00.set_season(month)
    
    mceq_run.set_theta_deg(0)
    
    print("longitude = {0}, latitude = {1}, month = {2}".format(longitude, latitude, month))
    print("Day of year = {0}, latitude = {1}, longitude = {2}".format(
        mceq_run.density_model._msis.inp.doy, 
        mceq_run.density_model._msis.inp.g_lat, 
        mceq_run.density_model._msis.inp.g_long))
    
    # Run MCEq
    mceq_run.solve()
    
    results.append(mceq_run.get_solution("total_mu+", mag = 3) + mceq_run.get_solution("total_mu-", mag = 3))


import matplotlib.pyplot as plt
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "sans-serif",
    "font.sans-serif": ["Helvetica"]})

plt.figure(figsize=(4.2, 3))
plt.semilogx(mceq_run.e_grid, results[0]/results[1], label=r'$\mathrm{Flux}_{\mu}$(0)/$\mathrm{Flux}_{\mu}$(1)')

plt.xlim(1., 1e9)
plt.xlabel('Kinetic energy (GeV)')
plt.legend()
plt.tight_layout()

[jncots]

Fixed by #47