experimenting with a python class wrapper for the fortran class

python/.f2py_f2cmap

@@ -9,5 +9,6 @@ dict(
         int32="int",
         c_int="int",
         c_long="long",
+        c_intptr_t="long long",
         c_long_long="long long"),
 )

python/radbeltpy.pyf

@@ -4,8 +4,10 @@
 python module radbelt ! in
     interface  ! in :radbelt
         module radbelt_c_module ! in :radbelt:radbelt_c_module.f90
-            use iso_c_binding, only: c_double,c_int,c_char
-            subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+            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
+                integer(c_intptr_t),intent(in) :: ipointer
                 real(c_double),intent(in) :: lon
                 real(c_double),intent(in) :: lat
                 real(c_double),intent(in) :: height
@@ -14,12 +16,23 @@ python module radbelt ! in
                 integer(c_int),intent(in) :: imname
                 real(c_double),intent(out) :: flux
             end subroutine get_flux_g_c
-            subroutine set_data_files_paths_c(aep8_dir, igrf_dir)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+
+            subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir)  ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
+                integer(c_intptr_t),intent(in) :: ipointer
                 character(kind=c_char,len=256),intent(in) :: aep8_dir
                 character(kind=c_char,len=256),intent(in) :: igrf_dir
                 !character(kind=c_char,len=1),dimension(*),intent(in) :: aep8_dir ! doesn't work
                 !character(kind=c_char,len=1),dimension(*),intent(in) :: igrf_dir
             end subroutine set_data_files_paths_c
+
+            subroutine initialize_c(ipointer)  ! in :radbelt: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
+                integer(c_intptr_t),intent(in) :: ipointer
+            end subroutine destroy_c
+
         end module radbelt_c_module
     end interface
 end python module radbelt

src/radbelt_c_module.f90

@@ -7,13 +7,12 @@
 
 module radbelt_c_module
 
-use iso_c_binding, only: c_double, c_int, c_char, c_null_char
+use iso_c_binding, only: c_double, c_int, c_char, c_null_char, &
+                         c_intptr_t, c_ptr, c_loc, c_f_pointer, c_null_ptr, c_associated
 use radbelt_module, only: radbelt_type
 
 implicit none
 
-type(radbelt_type) :: radbelt !! global variable
-
 integer,parameter :: STR_LEN = 256 !! string length for paths
 
 contains
@@ -31,16 +30,57 @@ module radbelt_c_module
 !     end do
 ! end function c2f_str
 
+subroutine initialize_c(ipointer) bind(C, name="initialize_c")
+    !! create a [[radbelt_type]] from C
+    integer(c_intptr_t),intent(out) :: ipointer
+    type(radbelt_type),pointer :: p
+    type(c_ptr) :: cp
+
+    allocate(p)
+    cp = c_loc(p)
+    ipointer = transfer(cp, 0_c_intptr_t)
+
+end subroutine initialize_c
+
+subroutine destroy_c(ipointer) bind(C, name="destroy_c")
+    !! destroy a [[radbelt_type]] from C
+    integer(c_intptr_t),intent(in) :: ipointer
+    type(radbelt_type),pointer :: p
+    type(c_ptr) :: cp
+
+    call int_pointer_to_f_pointer(ipointer,p)
+    if (associated(p)) deallocate(p)
+
+end subroutine destroy_c
+
+subroutine int_pointer_to_f_pointer(ipointer, p)
+
+    integer(c_intptr_t),intent(in) :: ipointer
+    type(radbelt_type),pointer :: p
+
+    type(c_ptr) :: cp
+
+    cp = transfer(ipointer, c_null_ptr)
+    if (c_associated(cp)) then
+        call c_f_pointer(cp, p)
+    else
+        p => null()
+    end if
+
+end subroutine int_pointer_to_f_pointer
+
 !*****************************************************************************************
 !>
 !  C interface for setting the data file paths
 
-subroutine set_data_files_paths_c(aep8_dir, igrf_dir) bind(C, name="set_data_files_paths_c")
+subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir) bind(C, name="set_data_files_paths_c")
 
+    integer(c_intptr_t),intent(in) :: ipointer
     character(kind=c_char,len=1),dimension(STR_LEN),intent(in) :: aep8_dir
     character(kind=c_char,len=1),dimension(STR_LEN),intent(in) :: igrf_dir
 
     character(len=:),allocatable :: aep8_dir_, igrf_dir_
+    type(radbelt_type),pointer :: p
 
     c2f : block
         ! convert to c strings to fortran
@@ -55,13 +95,15 @@ subroutine set_data_files_paths_c(aep8_dir, igrf_dir) bind(C, name="set_data_fil
         igrf_dir_ = trim(igrf_dir_)
     end block c2f
 
+    call int_pointer_to_f_pointer(ipointer, p)
+
     ! !... doesn't work this way... ??
     ! aep8_dir_ = c2f_str(aep8_dir)
     ! igrf_dir_ = c2f_str(igrf_dir)
     ! write(*,*) trim(aep8_dir_)
     ! write(*,*) trim(igrf_dir_)
 
-    call radbelt%set_data_files_paths(aep8_dir_, igrf_dir_)
+    call p%set_data_files_paths(aep8_dir_, igrf_dir_)
 
  end subroutine set_data_files_paths_c
 !*****************************************************************************************
@@ -70,8 +112,9 @@ end subroutine set_data_files_paths_c
 !>
 !  C interface to [[get_flux_g]].
 
-subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux) bind(C, name="get_flux_g_c")
+subroutine get_flux_g_c(ipointer,lon,lat,height,year,e,imname,flux) bind(C, name="get_flux_g_c")
 
+    integer(c_intptr_t),intent(in) :: ipointer
     real(c_double),intent(in) :: lon !! geodetic longitude in degrees (east)
     real(c_double),intent(in) :: lat !! geodetic latitude in degrees (north)
     real(c_double),intent(in) :: height !! altitude in km above sea level
@@ -86,7 +129,11 @@ subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux) bind(C, name="get_flu
                                  !! * 4 -- particle species: protons, solar activity: max
     real(c_double),intent(out) :: flux !! The flux of particles above the given energy, in units of cm^-2 s^-1.
 
-    flux = radbelt%get_flux(lon,lat,height,year,e,imname)
+    type(radbelt_type),pointer :: p
+
+    call int_pointer_to_f_pointer(ipointer, p)
+
+    flux = p%get_flux(lon,lat,height,year,e,imname)
 
 end subroutine get_flux_g_c
 

test/radbeltpy_test.py

@@ -10,6 +10,33 @@
 sys.path.insert(0,str(dir / 'python'))   # assuming the radbelt lib is in python directory
 import radbelt  # 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)
+
+    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)
+
 #########################################################################################################
 def set_data_files_paths(aep8_dir : str, igrf_dir : str) -> None:
     """Python function to set the file paths"""
@@ -21,11 +48,12 @@ def get_flux(lon : float, lat : float , height : float, year : float, e : float,
     """Python function to get the flux"""
     return radbelt.radbelt_c_module.get_flux_g_c(lon,lat,height,year,e,imname)
 
+model = radbelt_class()
 
 # set location of the data files:
 aep8_dir = str(dir / 'data' / 'aep8')
 igrf_dir = str(dir / 'data' / 'igrf')
-set_data_files_paths(aep8_dir, igrf_dir)
+model.set_data_files_paths(aep8_dir, igrf_dir)
 
 EPS = sys.float_info.epsilon # machine precision for error checking
 lon = -45.0
@@ -35,7 +63,7 @@ def get_flux(lon : float, lat : float , height : float, year : float, e : float,
 imname = 4 # 'p', 'max'
 e = 20.0
 
-flux = get_flux(lon,lat,height,year,e,imname)
+flux = model.get_flux(lon,lat,height,year,e,imname)
 
 print(f'flux = {flux}')
 
@@ -56,7 +84,7 @@ def get_flux(lon : float, lat : float , height : float, year : float, e : float,
     for lon in range(-180, 181, 45):
         for alt in range(500, 1001, 100):
             n_cases = n_cases + 1
-            f = get_flux(lon,lat,height,year, e, 4)
+            f = model.get_flux(lon,lat,height,year, e, 4)
 tend = time.time()
 print(f'Python version runtime: {tend-tstart} sec. {int(n_cases/(tend-tstart))} (cases/sec)')