Merge pull request #5 from jacobwilliams/develop

Develop

.gitignore

@@ -34,6 +34,7 @@
 # misc
 .DS_Store
 /results.txt
+/env
 
 # directories
 /build

README.md

@@ -29,17 +29,17 @@ The latest API documentation can be found [here](https://jacobwilliams.github.io
 
 ### Original source
 
-* The original sourcecode was hosted at GSFC "Modelweb", which no longer exists, but an archive can be found at the [Internet Archive](https://web.archive.org/web/20210318113325/https://ccmc.gsfc.nasa.gov/pub/modelweb/). It is presumed to be in the public domain. Reference: National Space Science Data Center, Data set PT-11B, Mar 1996. Dieter Bilitza, GSFC/NSSDC code 633, Greenbelt, MD 20771, tel. (301) 286-0190, [email protected]
+* The original sourcecode was hosted at GSFC "Modelweb", an archive of which can be found [here](https://git.smce.nasa.gov/ccmc-share/modelwebarchive). It is presumed to be in the public domain. Reference: National Space Science Data Center, Data set PT-11B, Mar 1996. Dieter Bilitza, GSFC/NSSDC code 633, Greenbelt, MD 20771.
 
 ### See also
 
-* [NASA ModelWebArchive](https://git.smce.nasa.gov/ccmc-share/modelwebarchive)
+* [NASA ModelWebArchive Archive (IGRF)](https://git.smce.nasa.gov/ccmc-share/modelwebarchive/-/tree/main/IGRF)
+* [NASA ModelWebArchive Archive (RADBELT)](https://git.smce.nasa.gov/ccmc-share/modelwebarchive/-/tree/main/RADBELT)
+* [International Geomagnetic Reference Field](https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
 * [An Astropy-friendly wrapper for the AE-8/AP-8 Van Allen belt model](https://github.com/nasa/radbelt)
 * [pyIGRF](https://github.com/rilma/pyIGRF)
 * https://github.com/lanl/RAM-SCB/blob/master/srcExternal/igrf.f
 * https://github.com/space-physics/igrf/blob/main/src/igrf/fortran/igrf13.f
-* https://web.archive.org/web/20210318113325/https://ccmc.gsfc.nasa.gov/pub/modelweb/
-* Model parameters can be computed and plotted online at http://nssdc.gsfc.nasa.gov/space/model/ [broken link]
 
 ### Test case
 
@@ -94,26 +94,13 @@ two lines had been exchanged.
 
 ### References
 
-* G.W. Singley, and J.I. Vette, The AE-4 Model of the Outer Radiation
-  Zone Electron Environment, NSSDC/WDC-A-R&S 72-06, 1972.
-* M.J. Teague, and J.I. Vette, A Model of the Trapped Electron
-  Population for Solar Minimum (AE-5), NSSDC/WDC-A-R&S 74-03, 1974.
-* M.J. Teague, K.W. Chan, and J.I. Vette, AE-6: A Model Environment
-  of Trapped Electrons for Solar Maximum, NSSDC/WDC-A-R&S 76-04, 1976
-* D.W. Sawyer, and J.I. Vette, AP-8 Trapped Proton Environment for
-  Solar Maximum and Minimum, NSSDC/WDC-A-R&S 76-06, 1976.
-* J.I. Vette, K.W. Chan, and M.J. Teague, Problems in Modeling the
-  Earth's Trapped Radiation Environment, AFGL-TR-78-0130, 1978.
-* K.W. Chan, M.J. Teague, N.J. Schofield, and J.I. Vette, Modeling of
-  Electron Time Variation in the Radiation Belts, p. 121-149, in:
-  Quantitative Modeling of Magnetospheric Processes, W.P. Olson
-  (ed.), geophysical monograph 21, American Geophysical Union, 1979.
-* M.T. Teague, N.J. Schofield, K.W. Chan, and J.I. Vette, A Study of
-  Inner Zone Electron Data and their Comparison with Trapped
-  Radiation Models, NSSDC/WDC-A-R&S 79-06, 1979.
-* J.I. Vette, The AE-8 Trapped Electron Model Environment,
-  NSSDC/WDC-A-R&S 91-24, 1991.
-* J.I. Vette, The NASA/National Space Science Data Center Trapped
-  Radiation Environment Model Program (1964-1991), NSSDC/WDC-A-R&S
-  91-29, 1991.
-* D. Heynderickx and A. Beliaev, J. Spacecraft and Rockets 32, 190-192, 1995.
+* G.W. Singley, and J.I. Vette, [The AE-4 Model of the Outer Radiation Zone Electron Environment](https://ntrs.nasa.gov/api/citations/19740012390/downloads/19740012390.pdf), NSSDC/WDC-A-R&S 72-06, 1972.
+* M.J. Teague, and J.I. Vette, [A Model of the Trapped Electron Population for Solar Minimum (AE-5)](https://ntrs.nasa.gov/api/citations/19740018161/downloads/19740018161.pdf), NSSDC/WDC-A-R&S 74-03, 1974.
+* M.J. Teague, K.W. Chan, and J.I. Vette, [AE-6: A Model Environment of Trapped Electrons for Solar Maximum](https://ntrs.nasa.gov/api/citations/19760016051/downloads/19760016051.pdf), NSSDC/WDC-A-R&S 76-04, 1976
+* D.W. Sawyer, and J.I. Vette, AP-8 Trapped Proton Environment for Solar Maximum and Minimum, NSSDC/WDC-A-R&S 76-06, 1976.
+* J.I. Vette, K.W. Chan, and M.J. Teague, [Problems in Modeling the Earth's Trapped Radiation Environment](https://apps.dtic.mil/sti/pdfs/ADA059273.pdf), AFGL-TR-78-0130, 1978.
+* K.W. Chan, M.J. Teague, N.J. Schofield, and J.I. Vette, Modeling of Electron Time Variation in the Radiation Belts, p. 121-149, in: Quantitative Modeling of Magnetospheric Processes, W.P. Olson (ed.), geophysical monograph 21, American Geophysical Union, 1979.
+* M.T. Teague, N.J. Schofield, K.W. Chan, and J.I. Vette, [A Study of Inner Zone Electron Data and their Comparison with Trapped Radiation Models](https://ntrs.nasa.gov/api/citations/19790025500/downloads/19790025500.pdf), NSSDC/WDC-A-R&S 79-06, 1979.
+* J.I. Vette, [The AE-8 Trapped Electron Model Environment](https://ntrs.nasa.gov/api/citations/19920014985/downloads/19920014985.pdf), NSSDC/WDC-A-R&S 91-24, 1991.
+* J.I. Vette, [The NASA/National Space Science Data Center Trapped Radiation Environment Model Program (1964-1991)](https://ntrs.nasa.gov/api/citations/19930001815/downloads/19930001815.pdf), NSSDC/WDC-A-R&S 91-29, 1991.
+* D. Heynderickx and A. Beliaev, [Identification of an error in the distribution of the NASA model AP-8 MIN](https://arc.aiaa.org/doi/10.2514/3.26595), J. Spacecraft and Rockets 32, 190-192, 1995.

ford.md

@@ -18,5 +18,6 @@ graph: true
 externalize: true
 preprocessor: gfortran -E
 extra_mods: iso_fortran_env: https://gcc.gnu.org/onlinedocs/gfortran/ISO_005fFORTRAN_005fENV.html
+            iso_c_binding: https://gcc.gnu.org/onlinedocs/gfortran/ISO_005fC_005fBINDING.html
 
 {!README.md!}

meson.build

@@ -0,0 +1,26 @@
+project('radbelt', ['c', 'fortran'],
+        version: '0.0.1',
+        default_options : [
+          'warning_level=everything',
+          'buildtype=release',]
+        )
+
+add_languages('fortran')
+
+_args = [] # Extra arguments
+_deps = [] # Dependencies
+# _deps += dependency('lapack')
+
+cc = meson.get_compiler('c')
+
+radbeltlib = library('RadBelt',
+              sources: [
+                'src/radbelt_kinds_module.F90' ,
+                'src/trmfun.f90' ,
+                'src/shellig.f90' ,
+                'src/radbelt_module.f90' ,
+                'src/radbelt_c_module.f90' ,
+              ],
+              dependencies: _deps,
+              install: true,
+)

python/.f2py_f2cmap

@@ -0,0 +1,14 @@
+dict(
+    real=dict(
+        real64="double",
+        real32="float",
+        c_double="double",
+        c_float="float"),
+    integer=dict(
+        int64="long long",
+        int32="int",
+        c_int="int",
+        c_long="long",
+        c_intptr_t="long long",
+        c_long_long="long long"),
+)

python/build_python.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+#
+# Build the python extension module using f2py
+#
+
+# this is to eliminate the "ld: warning: could not create compact unwind" on MacOS
+export LDFLAGS=-Wl,-no_compact_unwind
+
+f2py -c --backend meson ./radbeltpy.pyf \
+    ../src/radbelt_kinds_module.F90 \
+    ../src/trmfun.f90 \
+    ../src/shellig.f90 \
+    ../src/radbelt_module.f90 \
+    ../src/radbelt_c_module.f90

python/environment.yml

@@ -0,0 +1,18 @@
+# To use:
+#
+#   $ conda env create --prefix ./env -f ./python/environment.yml
+#   $ conda activate ./env
+
+channels:
+  - nodefaults
+  - conda-forge
+
+dependencies:
+  - fortran-compiler
+  - meson
+  - pkg-config
+  - fpm
+  - ford
+  - graphviz
+  - fprettify
+  - numpy

python/radbeltpy.pyf

@@ -0,0 +1,38 @@
+!    -*- 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
+            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
+                real(c_double),intent(in) :: year
+                real(c_double),intent(in) :: e
+                integer(c_int),intent(in) :: imname
+                real(c_double),intent(out) :: flux
+            end subroutine get_flux_g_c
+
+            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=4096),intent(in) :: aep8_dir
+                character(kind=c_char,len=4096),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
+
+! This file was manually created

src/radbelt_c_module.f90

@@ -0,0 +1,155 @@
+!*****************************************************************************************
+!>
+!  Experimental C interface to the radbelt module.
+
+    module radbelt_c_module
+
+    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
+
+    ! can we eliminate this ?
+    integer,parameter :: STR_LEN = 4096 !! string length for paths
+
+    contains
+!*****************************************************************************************
+
+!*****************************************************************************************
+!>
+!  Convert C string to Fortran
+
+function c2f_str(cstr) result(fstr)
+
+    character(kind=c_char,len=1),dimension(:),intent(in) :: cstr !! string from C
+    character(len=:),allocatable :: fstr !! fortran string
+
+    integer :: i !! counter
+
+    fstr = ''
+    do i = 1, size(cstr)
+        fstr = fstr//cstr(i)
+    end do
+    fstr = trim(fstr)
+
+end function c2f_str
+
+!*****************************************************************************************
+!>
+!  Convert an integer pointer to a [[radbelt_type]] pointer.
+
+subroutine int_pointer_to_f_pointer(ipointer, p)
+
+    integer(c_intptr_t),intent(in) :: ipointer !! integer pointer from C
+    type(radbelt_type),pointer :: p !! fortran pointer
+
+    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
+
+!*****************************************************************************************
+!>
+!  create a [[radbelt_type]] from C
+
+subroutine initialize_c(ipointer) bind(C, name="initialize_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
+
+!*****************************************************************************************
+!>
+!  destroy a [[radbelt_type]] from C
+
+subroutine destroy_c(ipointer) bind(C, name="destroy_c")
+
+    integer(c_intptr_t),intent(in) :: ipointer
+    type(radbelt_type),pointer :: p
+
+    call int_pointer_to_f_pointer(ipointer,p)
+    if (associated(p)) deallocate(p)
+
+end subroutine destroy_c
+
+!*****************************************************************************************
+!>
+!  C interface for setting the data file paths
+
+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
+
+    call int_pointer_to_f_pointer(ipointer, p)
+
+    if (associated(p)) then
+
+        aep8_dir_ = c2f_str(aep8_dir)
+        igrf_dir_ = c2f_str(igrf_dir)
+
+        call p%set_data_files_paths(aep8_dir_, igrf_dir_)
+
+    else
+        error stop 'error in set_data_files_paths_c: class is not associated'
+    end if
+
+ end subroutine set_data_files_paths_c
+!*****************************************************************************************
+
+!*****************************************************************************************
+!>
+!  C interface to [[get_flux_g]].
+
+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
+    real(c_double),intent(in) :: year !! decimal year for which geomagnetic field is to
+                                      !! be calculated (e.g.:1995.5 for day 185 of 1995)
+    real(c_double),intent(in) :: e !! minimum energy
+    integer(c_int),intent(in) :: imname !! which method to use:
+                                        !!
+                                        !! * 1 -- particle species: electrons, solar activity: min
+                                        !! * 2 -- particle species: electrons, solar activity: max
+                                        !! * 3 -- particle species: protons, solar activity: min
+                                        !! * 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.
+
+    type(radbelt_type),pointer :: p
+
+    call int_pointer_to_f_pointer(ipointer, p)
+
+    if (associated(p)) then
+
+        flux = p%get_flux(lon,lat,height,year,e,imname)
+
+    else
+        error stop 'error in get_flux_g_c: class is not associated'
+    end if
+
+end subroutine get_flux_g_c
+
+!*****************************************************************************************
+    end module radbelt_c_module
+!*****************************************************************************************

src/redbelt_module.f90 → src/radbelt_module.f90

@@ -6,7 +6,7 @@
 !   * https://ccmc.gsfc.nasa.gov/pub/modelweb/geomagnetic/igrf/fortran_code/bilcal.for
 !   * https://ccmc.gsfc.nasa.gov/pub/modelweb/radiation_belt/radbelt/fortran_code/radbelt.for
 
-module redbelt_module
+module radbelt_module
 
    use radbelt_kinds_module
    use trmfun_module
@@ -171,4 +171,4 @@ function get_flux_c(v,year,e,imname) result(flux)
 
    end function get_flux_c
 
-end module redbelt_module
+end module radbelt_module