experimenting with a python package

not finished

.gitignore

@@ -36,6 +36,9 @@
 /results.txt
 /env
 
+# python
+__pycache__
+
 # directories
 /build
 /doc

radbeltpy/__init__.py

@@ -0,0 +1,75 @@
+#
+# Python interface to the `radbelt_fortran` shared library.
+#
+
+from . import radbelt_fortran
+
+class RadbeltClass:
+    """
+    Class for using the radbelt model.
+    """
+
+    def __init__(self) -> None:
+        """constructor for the fortran class"""
+
+        #`ip` is an integer that represents a c pointer
+        # to a `radbelt_type` in the Fortran library.
+        self.ip = radbelt_fortran.radbelt_c_module.initialize_c()
+
+        #TODO should allow for passing in the data directories here
+
+    def __del__(self) -> None:
+        """destructor for the fortran class"""
+
+        radbelt_fortran.radbelt_c_module.destroy_c(self.ip)
+
+    def set_data_files_paths(self, aep8_dir : str, igrf_dir : str) -> None:
+        """Set the file paths"""
+
+        #TODO split these up so they can be called separately
+
+        radbelt_fortran.radbelt_c_module.set_data_files_paths_c(self.ip, aep8_dir, igrf_dir,
+                                                                len(aep8_dir), len(igrf_dir))
+
+    def get_flux(self, lon : float, lat : float , height : float, year : float,
+                 e : float, particle : str, solar : str) -> float:
+        """
+        Compute the flux.
+
+        Parameters
+        ----------
+        lon : float
+            Longitude (deg)
+        lat : float
+            Latitude (deg)
+        height : float
+            Altitude (deg)
+        year : float
+            Decimal year (needed to account for drift of the Earth's magnetic field).
+        energy : float
+            Minimum energy.
+        particle : {'e', 'p'}
+            Particle species: 'e' for electrons, 'p' for protons.
+        solar : {'min', 'max'}
+            Solar activity: solar minimum or solar maximum.
+
+        Returns
+        -------
+        flux : float
+            The flux of particles above the given energy, (cm^-2 s^-1).
+        """
+
+        model = _model_index(particle, solar)
+        return radbelt_fortran.radbelt_c_module.get_flux_g_c(self.ip,lon,lat,height,year,e,model)
+
+def _model_index(particle : str, solar : str) -> int:
+    """convert the particle and solar options to model number for the low-level routine"""
+
+    p = particle.lower()
+    s = solar.lower()
+    if p not in ('e', 'p'): raise ValueError('particle must be "e" or "p"')
+    if s not in ('min', 'max'): raise ValueError('solar must be "min" or "max"')
+    return { ('e', 'min'): 1,
+             ('e', 'max'): 2,
+             ('p', 'min'): 3,
+             ('p', 'max'): 4 }[(p, s)]

python/radbeltpy.pyf → radbeltpy/radbeltpy.pyf

@@ -1,12 +1,12 @@
 !    -*- f90 -*-
 ! Note: the context of this file is case sensitive.
 
-python module radbelt ! in
-    interface  ! in :radbelt
-        module radbelt_c_module ! in :radbelt:radbelt_c_module.f90
+python module radbelt_fortran ! in
+    interface  ! in :radbelt_fortran
+        module radbelt_c_module ! in :radbelt_fortran:radbelt_c_module.f90
             use iso_c_binding, only: c_double,c_int,c_char,c_intptr_t
 
-            subroutine get_flux_g_c(ipointer,lon,lat,height,year,e,imname,flux)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+            subroutine get_flux_g_c(ipointer,lon,lat,height,year,e,imname,flux)  ! in :radbelt_fortran:radbelt_c_module.f90:radbelt_c_module
                 integer(c_intptr_t),intent(in) :: ipointer
                 real(c_double),intent(in) :: lon
                 real(c_double),intent(in) :: lat
@@ -17,24 +17,24 @@ python module radbelt ! in
                 real(c_double),intent(out) :: flux
             end subroutine get_flux_g_c
 
-            subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+            subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m)  ! in :radbelt_fortran:radbelt_c_module.f90:radbelt_c_module
                 integer(c_intptr_t),intent(in) :: ipointer
                 character(kind=c_char,len=n),intent(in),depend(n) :: aep8_dir
                 character(kind=c_char,len=m),intent(in),depend(m) :: igrf_dir
                 integer(c_int),intent(in) :: n
                 integer(c_int),intent(in) :: m
             end subroutine set_data_files_paths_c
 
-            subroutine initialize_c(ipointer)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+            subroutine initialize_c(ipointer)  ! in :radbelt_fortran:radbelt_c_module.f90:radbelt_c_module
                 integer(c_intptr_t),intent(out) :: ipointer
             end subroutine initialize_c
 
-            subroutine destroy_c(ipointer)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+            subroutine destroy_c(ipointer)  ! in :radbelt_fortran:radbelt_c_module.f90:radbelt_c_module
                 integer(c_intptr_t),intent(in) :: ipointer
             end subroutine destroy_c
 
         end module radbelt_c_module
     end interface
-end python module radbelt
+end python module radbelt_fortran
 
 ! This file was manually created

test/radbeltpy_test.py

@@ -7,38 +7,10 @@
 from pathlib import Path
 
 dir = Path(__file__).resolve().parents[1] # root directory
-sys.path.insert(0,str(dir / 'python'))   # assuming the radbelt lib is in python directory
-import radbelt  # this is the module being tested
+sys.path.insert(0,str(dir))   # assuming the radbelt lib is in python directory
+from radbeltpy import RadbeltClass  # this is the module being tested
 
-class radbelt_class:
-    """
-    Class for using the radbelt model.
-    """
-
-    def __init__(self) -> None:
-        """constructor for the fortran class"""
-
-        #`ip` is an integer that represents a c pointer
-        # to a `radbelt_type` in the Fortran library.
-        self.ip = radbelt.radbelt_c_module.initialize_c()
-
-    def __del__(self) -> None:
-        """destructor for the fortran class"""
-
-        radbelt.radbelt_c_module.destroy_c(self.ip)
-
-    def set_data_files_paths(self, aep8_dir : str, igrf_dir : str) -> None:
-        """Set the file paths"""
-
-        radbelt.radbelt_c_module.set_data_files_paths_c(self.ip, aep8_dir, igrf_dir, len(aep8_dir), len(igrf_dir))
-
-    def get_flux(self, lon : float, lat : float , height : float, year : float, e : float, imname : int) -> float:
-        """Get the flux"""
-
-        return radbelt.radbelt_c_module.get_flux_g_c(self.ip, lon,lat,height,year,e,imname)
-
-
-model = radbelt_class()
+model = RadbeltClass()
 
 # set location of the data files:
 aep8_dir = str(dir / 'data' / 'aep8')
@@ -50,10 +22,11 @@ def get_flux(self, lon : float, lat : float , height : float, year : float, e :
 lat = -30.0
 height = 500.0
 year = 2021.1616438356164  # decimal year
-imname = 4 # 'p', 'max'
+particle = 'p'
+solar = 'max'
 e = 20.0
 
-flux = model.get_flux(lon,lat,height,year,e,imname)
+flux = model.get_flux(lon,lat,height,year,e,particle,solar)
 
 print(f'flux = {flux}')
 
@@ -74,34 +47,6 @@ def get_flux(self, lon : float, lat : float , height : float, year : float, e :
     for lon in range(-180, 181, 45):
         for alt in range(500, 1001, 100):
             n_cases = n_cases + 1
-            f = model.get_flux(lon,lat,height,year, e, 4)
+            f = model.get_flux(lon,lat,height,year, e, particle,solar)
 tend = time.time()
 print(f'Python version runtime: {tend-tstart} sec. {int(n_cases/(tend-tstart))} (cases/sec)')
-
-# print('... VECTORIZE ...')
-# n_cases = 0
-# tstart = time.time()
-
-# lon_vec = []
-# lat_vec = []
-# height_vec = []
-# year_vec = []
-# e_vec = []
-# ifile_vec = []
-
-# for lat in range(-90, 91, 5):
-#     for lon in range(-180, 181, 45):
-#         for alt in range(500, 1001, 100):
-#             n_cases = n_cases + 1
-#             lon_vec.append(lon)
-#             lat_vec.append(lat)
-#             height_vec.append(height)
-#             year_vec.append(year)
-#             e_vec.append(e)
-#             ifile_vec.append(4)
-#             #f = get_flux(lon,lat,height,year, e, 4)
-
-# f_vec = get_flux(lon_vec,lat_vec,height_vec,year_vec,e_vec,ifile_vec)
-# tend = time.time()
-# print(f'Python version runtime: {tend-tstart} sec. {int(n_cases/(tend-tstart))} (cases/sec)')
-