Linear energy transfer (LET) model (beta)

Model API name: let

Version: v1
Model developer:
ESTEC
Model provision:
SPENVIS
Category: model/transport
Keywords: stopping-powers
Model flow
External Input: spectra

Input	Value
Quantity	energy_flux_spectrum \| energy_fluence_spectrum
External Input: stopping_power

Input	Value
Quantity	stopping_power
Input group: input_spectra / multiplicity: many

Input	Description	Valid values	Default	Quantity
spectrum	Input spectrum		[]	energy_flux_spectrum
toflux	Agregation period (if fluence values are specified)		spectrum	duration
Name
Quantity
Variables
let_spectrum
LET Spectrum
let_fluence_spectrum
Qualifiers:
Name	Units	Quantity	Variable qualifiers
let	(MeV cm^2 g^-1)	let
integral_fluence axis 1: let	(m^-2 sr^-1)	number_fluence	direction: directional form: integral species: ions
differential_fluence axis 1: let	(m^-2 sr^-1 (MeV cm^2 g^-1)^-1)	number_fluence	direction: directional form: differential species: ions
None

import math
import matplotlib.pyplot as plt

from nom_client.nom_client import NoMClient

plt.rcParams['figure.figsize'] = [15, 10]
plt.rc('font', size=18)


def do_let_calc(input_spectra_set, nom_client_):
    print(f"Doing ... {input_spectra_set}")

    ss_model = nom_client_.get_model('shielded-spectrum')
    ss_model.set_params(shieldThickness=SHIELD_THICKNESS)
    ss_model.set_external_input(external_input_name="spectra", external_input=input_spectra_set)
    ss_results = nom_client_.run_model(ss_model)
    print(ss_results)

    shielded_spectrum = ss_results.get_model_result_by_name(result_name="shielded_spectrum")
    calc_shielded_energy = shielded_spectrum.get_variable_data(variable_name="energy").squeeze()
    calc_shielded_diff_flux = shielded_spectrum.get_variable_data(variable_name="differential_fluence",
                                                                  z=ATOMIC_NUMBER_TO_PLOT).squeeze()
    calc_shielded_int_flux = shielded_spectrum.get_variable_data(variable_name="integral_fluence",
                                                                 z=ATOMIC_NUMBER_TO_PLOT).squeeze()

    # Calculate the stopping powers

    spectrum_energy_min = min(calc_shielded_energy)
    spectrum_energy_max = max(calc_shielded_energy)

    stopping_powers_model = nom_client_.get_model('stopping-powers')
    stopping_powers_model.set_params(energyLow=spectrum_energy_min, energyHigh=spectrum_energy_max,
                                     numberEnergies=66,
                                     material="G4_Si", particle="all", data_source="g4")
    stopping_powers_result = nom_client_.run_model(stopping_powers_model)
    stopping_powers = stopping_powers_result.get_model_result_by_name(result_name="stopping_powers")

    # Calculate the LET spectra

    let_model = nom_client_.get_model('let')
    let_model.set_external_input(external_input_name="spectra", external_input=[shielded_spectrum])
    let_model.set_external_input(external_input_name="stopping_power", external_input=stopping_powers_result)
    let_results = nom_client_.run_model(let_model)

    calc_let_spectrum = let_results.get_model_result_by_name(result_name="let_spectrum")
    calc_let = calc_let_spectrum.get_variable_data(variable_name="let").squeeze()
    calc_diff_flux = calc_let_spectrum.get_variable_data(variable_name="differential_fluence").squeeze()
    calc_int_flux = calc_let_spectrum.get_variable_data(variable_name="integral_fluence").squeeze()

    return calc_let, calc_int_flux, calc_shielded_energy


cm_to_m = 1e4
m_to_cm = 1e-4
to_omni = 4 * math.pi
MIS_DUR_IN_DAYS = 1  # 1820
MIS_DUR_IN_S = MIS_DUR_IN_DAYS * 86400

SHIELD_THICKNESS = 0
SPECIES_NAME_TO_PLOT = "H"
ATOMIC_NUMBER_TO_PLOT = 1


def plot(x_data, y_data, plot_types, line_types, xlabel, ylabel, labels, title):
    if not isinstance(plot_types, list):
        plot_types = [plot_types] * len(x_data)

    if not isinstance(line_types, list):
        line_types = [line_types] * len(x_data)

    if not isinstance(labels, list):
        labels = [labels] * len(x_data)

    for data_index in range(len(x_data)):
        if plot_types[data_index] == 'loglog':
            plt.loglog(x_data[data_index], y_data[data_index], linestyle=line_types[data_index], linewidth=2.0,
                       markersize=12.0,
                       label=labels[data_index])
        elif plot_types[data_index] == 'semilogy':
            plt.semilogy(x_data[data_index], y_data[data_index], linestyle=line_types[data_index], linewidth=2.0,
                         markersize=12.0,
                         label=labels[data_index])
        elif plot_types[data_index] == 'semilogx':
            plt.semilogx(x_data[data_index], y_data[data_index], linestyle=line_types[data_index], linewidth=2.0,
                         markersize=12.0,
                         label=labels[data_index])
        else:
            plt.plot(x_data[data_index], y_data[data_index], linestyle=line_types[data_index], linewidth=2.0,
                     markersize=12.0,
                     label=labels[data_index])

    plt.title(title, fontsize=18)
    plt.grid(True)
    plt.legend(bbox_to_anchor=(0.95, 0.95), loc=1, borderaxespad=0.,
               fancybox=True, fontsize=18, shadow=True)

    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.show()


nom_client = NoMClient(project_name="LET spectrum example", default_server_id="local_server", debug_output=False)

sapre_model_gto = nom_client.get_model('sapre')
sapre_model_gto.set_params(orbitType="NEI", missionDuration=MIS_DUR_IN_DAYS, missionDurationUnits=1)
sapre_gto_result = nom_client.run_model(sapre_model_gto)
print(sapre_gto_result)

LOWEST_ION_SPECIES = 1
HIGHEST_ION_SPECIES = 92

# CREME-96 GCR quiet time

creme96gcr_model = nom_client.get_model('creme-96-gcr')
creme96gcr_model.set_params(useGeomagneticShielding=0,
                            lowestIonSpecies=LOWEST_ION_SPECIES, highestIonSpecies=HIGHEST_ION_SPECIES)
creme96gcr_model.set_external_input(external_input_name="trajectory", external_input=sapre_gto_result)
creme96gcr_results = nom_client.run_model(creme96gcr_model)
print(creme96gcr_results)

creme96gcr_spectrum_m2_s_sr = creme96gcr_results.get_model_result_by_name(result_name="particle_spectrum")
creme96gcr_spectrum_energies = creme96gcr_spectrum_m2_s_sr.get_variable_data(variable_name="energy").squeeze()
creme96gcr_spectrum_diff_flux_m2_s_sr = creme96gcr_spectrum_m2_s_sr.get_variable_data(variable_name="differential_flux",
                                                                                      species=SPECIES_NAME_TO_PLOT).squeeze()
creme96gcr_spectrum_int_flux_m2_s_sr = creme96gcr_spectrum_m2_s_sr.get_variable_data(variable_name="integral_flux",
                                                                                     species=SPECIES_NAME_TO_PLOT).squeeze()

# CREME-96 SEP peak

creme96sep_model = nom_client.get_model('creme-96-sep')
creme96sep_model.set_params(includeGeomagneticShielding=0, eventCase=3,
                            lowestIonSpecies=LOWEST_ION_SPECIES, highestIonSpecies=HIGHEST_ION_SPECIES)
creme96sep_model.set_external_input(external_input_name="trajectory", external_input=sapre_gto_result)
creme96sep_results = nom_client.run_model(creme96sep_model)
print(creme96sep_results)

creme96sep_spectrum_m2_s_sr = creme96sep_results.get_model_result_by_name(result_name="full_ion_spectrum")
creme96sep_spectrum_energies = creme96sep_spectrum_m2_s_sr.get_variable_data(variable_name="energy").squeeze()
creme96sep_peak_spectrum_diff_flux_m2_s_sr = creme96sep_spectrum_m2_s_sr.get_variable_data(
    variable_name="differential_flux",
    species=SPECIES_NAME_TO_PLOT).squeeze()
creme96sep_peak_spectrum_int_flux_m2_s_sr = creme96sep_spectrum_m2_s_sr.get_variable_data(variable_name="integral_flux",
                                                                                          species=SPECIES_NAME_TO_PLOT).squeeze()

# CREME-96 SEP worst week

creme96sepww_model = nom_client.get_model('creme-96-sep')
creme96sepww_model.set_params(includeGeomagneticShielding=0, eventCase=1,
                              lowestIonSpecies=LOWEST_ION_SPECIES, highestIonSpecies=HIGHEST_ION_SPECIES)
creme96sepww_model.set_external_input(external_input_name="trajectory", external_input=sapre_gto_result)
creme96sepww_results = nom_client.run_model(creme96sepww_model)
print(creme96sepww_results)

creme96sepww_spectrum_m2_s_sr = creme96sepww_results.get_model_result_by_name(result_name="full_ion_spectrum")
creme96sepww_spectrum_energies = creme96sepww_spectrum_m2_s_sr.get_variable_data(variable_name="energy").squeeze()
creme96sep_ww_spectrum_diff_flux_m2_s_sr = creme96sepww_spectrum_m2_s_sr.get_variable_data(
    variable_name="differential_flux",
    species=SPECIES_NAME_TO_PLOT).squeeze()
creme96sepww_spectrum_int_flux_m2_s_sr = creme96sepww_spectrum_m2_s_sr.get_variable_data(variable_name="integral_flux",
                                                                                         species=SPECIES_NAME_TO_PLOT).squeeze()

# Your mission

iso_15390_model = nom_client.get_model('iso_15390')
iso_15390_model.set_params(lowestIonSpecies=LOWEST_ION_SPECIES, highestIonSpecies=HIGHEST_ION_SPECIES)
iso_15390_model.set_external_input(external_input_name="trajectory", external_input=sapre_gto_result)
iso_15390_results = nom_client.run_model(iso_15390_model)
print(iso_15390_results)

mission_gcr_spectrum_m2_s_sr = iso_15390_results.get_model_result_by_name(result_name="particle_spectrum")
mission_gcr_spectrum_energies = mission_gcr_spectrum_m2_s_sr.get_variable_data(variable_name="energy").squeeze()
mission_gcr_spectrum_diff_flux_m2_s_sr = mission_gcr_spectrum_m2_s_sr.get_variable_data(
    variable_name="differential_flux",
    species=SPECIES_NAME_TO_PLOT).squeeze()
mission_gcr_spectrum_int_flux_m2_s_sr = mission_gcr_spectrum_m2_s_sr.get_variable_data(variable_name="integral_flux",
                                                                                       species=SPECIES_NAME_TO_PLOT).squeeze()

sapphire_model = nom_client.get_model('sapphire-total-fluence')
sapphire_model.set_params(lowestIonSpecies=LOWEST_ION_SPECIES, highestIonSpecies=HIGHEST_ION_SPECIES)
sapphire_model.set_external_input(external_input_name="trajectory", external_input=sapre_gto_result)
sapphire_results = nom_client.run_model(sapphire_model)
print(sapphire_results)

sep_spectrum_cm2 = sapphire_results.get_model_result_by_name(result_name="full_ion_spectrum")
sep_spectrum_energies = sep_spectrum_cm2.get_variable_data(variable_name="energy").squeeze()
sep_spectrum_diff_fluence_cm2 = sep_spectrum_cm2.get_variable_data(variable_name="differential_fluence",
                                                                   species=SPECIES_NAME_TO_PLOT).squeeze()
sep_spectrum_int_fluence_cm2 = sep_spectrum_cm2.get_variable_data(variable_name="integral_fluence",
                                                                  species=SPECIES_NAME_TO_PLOT).squeeze()

plot(x_data=[mission_gcr_spectrum_energies, sep_spectrum_energies],
     y_data=[m_to_cm * to_omni * mission_gcr_spectrum_int_flux_m2_s_sr, sep_spectrum_int_fluence_cm2 / MIS_DUR_IN_S],
     plot_types=["loglog", "loglog"],
     line_types=["solid", "solid"],
     xlabel="Energy (MeV)",
     ylabel="DIfferential flux (cm^-2 s^-1 MeV^-1)",
     labels=["GCR", "SAPPHIRE total fluence"],
     title="Unshielded GCR / SEP spectra")
#
mission_let = do_let_calc(input_spectra_set=[sep_spectrum_cm2], nom_client_=nom_client)

#mission_let = do_let_calc(input_spectra_set=[sep_spectrum_cm2, mission_gcr_spectrum_m2_s_sr], nom_client_=nom_client)
quiet_let = do_let_calc(input_spectra_set=[creme96gcr_spectrum_m2_s_sr], nom_client_=nom_client)
creme_peak_let = do_let_calc(input_spectra_set=[creme96sep_spectrum_m2_s_sr], nom_client_=nom_client)
creme_ww_let = do_let_calc(input_spectra_set=[creme96sepww_spectrum_m2_s_sr], nom_client_=nom_client)


plot(x_data=[creme_peak_let[0],
             creme_ww_let[0],
             mission_let[0],
             quiet_let[0]],
     y_data=[creme_peak_let[1],
             creme_ww_let[1],
             cm_to_m * mission_let[1] / (MIS_DUR_IN_S * to_omni),
             quiet_let[1]],
     plot_types='loglog',
     line_types="solid",
     xlabel="LET (MeV cm^2 g^-1)",
     ylabel="Differential flux cm^-2 s^-1 (MeV cm^2 g^-1)^-1",
     labels=["Peak 5-min", "Worst week", "mission", "Quiet, GCR only"],
     title="Title")

d = 3
#
# plot(x_data=[mission_let[0]],
#      y_data=[mission_let[1]],
#      plot_types='loglog',
#      line_types="solid",
#      xlabel="LET (MeV cm^2 g^-1)",
#      ylabel="Differential flux cm^-2 s^-1 (MeV cm^2 g^-1)^-1",
#      labels=["Peak 5-min", "Worst week", "mission", "Quiet, GCR only"],
#      title="Title 2")
#
# plot(x_data=[mission_let[2]],
#      y_data=[mission_let[0]],
#      plot_types='loglog',
#      line_types="solid",
#      xlabel="Energy (MeV)",
#      ylabel="LET (MeV cm^2 g^-1)",
#      labels=["peak", "ww", "mission", "quiet"],
#      title="Title 3")
No references
Linear energy transfer (LET) model (beta)

Model API name: let

Model inputs

Model outputs

Documentation

Examples

References

External Input: spectra

External Input: stopping_power

Input group: input_spectra / multiplicity: many
