''' Using ice experiment
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pyana.ice11.ene_spec_file as icespec
def hydrogen_1keV_as_ganymede_cusp():
''' 1keV hydrogen beam to simulate 2 eV incoming flux. Hopeless!
'''
dpath = os.path.join(icespec.get_data_path(), 'sv3/hydrogen/Hp_1keV_1.0nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 2.5e7
enespec[1, :] *= flux_ganymede
energy_factor = 2. / 1000
enespec *= energy_factor
plt.plot(enespec[0, :], enespec[1, :], label='Cusp:H+ (1keV/2eV) -> H')
def hydrogen_1keV_as_ganymede_corotation():
''' 1keV hydrogen beam to simulate 120 eV incoming flux. Hopeless!
'''
dpath = os.path.join(icespec.get_data_path(), 'sv3/hydrogen/Hp_1keV_1.0nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 1.5e7
enespec[1, :] *= flux_ganymede
energy_factor = 120. / 1000
enespec *= energy_factor
plt.plot(enespec[0, :], enespec[1, :], label='Corotate:H+ (1keV/120eV) -> H')
def oxygen_3keV_as_ganymede_corotation():
dpath = os.path.join(icespec.get_data_path(), 'sv3/oxygen/Op_3keV_2nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 4.5e7
enespec[1, :] *= flux_ganymede
energy_factor = 3000. / 3000
enespec *= energy_factor
plt.plot(enespec[0, :], enespec[1, :], label='Corotate:O+ (3keV/3keV) -> H')
def oxygen_3keV_as_ganymede_cusp():
dpath = os.path.join(icespec.get_data_path(), 'sv3/oxygen/Op_3keV_2nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 7.5e7
enespec[1, :] *= flux_ganymede
energy_factor = 50. / 3000
enespec *= energy_factor
plt.plot(enespec[0, :], enespec[1, :], label='Cusp:O+ (3keV/50eV) -> H')
def hydrogen_23keV_as_ganymede_energetic():
dpath = os.path.join(icespec.get_data_path(), 'sv2/hydrogen/Hp_23keV_1.2nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 5e6
enespec[1, :] *= flux_ganymede
energy_factor = 100. / 23
enespec *= energy_factor
plt.plot(enespec[0, :], enespec[1, :], label='Energetic:H+ (23keV/100keV) -> H')
def oxygen_43keV_as_ganymede_energetic():
dpath = os.path.join(icespec.get_data_path(), 'sv2/oxygen/Op_43keV_2nA.dat')
specobj = icespec.EspecFile(dpath)
enespec = specobj.get_energy_spectra_cps(species='H')
# The unit is normalized by ``nA``, meaning # / nA / s
# The beam width is roughly 5mm diameter.
# The corresponding flux is 1e9 / e / D = 3.185e10 #/cm2 s
unit_factor = 3.185e10
enespec[1, :] /= unit_factor # Now the unit is for incoming flux of 1 #/cm2 s.
flux_ganymede = 2e6
enespec[1, :] *= flux_ganymede # This is a count rate.
energy_factor = 100. / 43
enespec *= energy_factor # Both energy and count rate.
plt.plot(enespec[0, :], enespec[1, :], label='Energetic:O+ (43keV/100keV) -> H')
def main():
'''Main script'''
### The lowest energy hydrogen beam
plt.figure(figsize=(10, 10))
hydrogen_1keV_as_ganymede_cusp()
oxygen_3keV_as_ganymede_cusp()
hydrogen_1keV_as_ganymede_corotation()
oxygen_3keV_as_ganymede_corotation()
hydrogen_23keV_as_ganymede_energetic()
oxygen_43keV_as_ganymede_energetic()
plt.axhline(np.sqrt(1.97/(15./8.)), color='orange')
plt.text(5000, np.sqrt(1.97/(15./8.)), 'BG (1 scan=15s)', horizontalalignment='right', color='orange')
plt.axhline(np.sqrt(1.97/(60/8.)), linestyle='--', color='orange') # 60 sec
plt.text(5000, np.sqrt(1.97/(60./8.)), 'BG (4 scan=60s)', horizontalalignment='right', color='orange')
plt.axhline(np.sqrt(1.97/(300/8.)), linestyle='-.', color='orange') # 300 sec
plt.text(5000, np.sqrt(1.97/(300./8.)), 'BG (20 scan=5m)', horizontalalignment='right', color='orange')
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='upper left')
plt.xlim(5, 5000)
plt.ylim(1e-5, 1e2)
plt.savefig('ice_scaling.eps')
if __name__ == "__main__":
main()