From 652f45fa0b60870f5ac41263357261fd30b5b911 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 11 Feb 2024 21:31:45 -0600 Subject: [PATCH 01/10] some refactoring --- src/shellig.f90 | 75 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 25 deletions(-) diff --git a/src/shellig.f90 b/src/shellig.f90 index 1fb501c..20846ba 100644 --- a/src/shellig.f90 +++ b/src/shellig.f90 @@ -359,10 +359,9 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v) !! instead of glat,glon,alt. See [[shellc]]. real(wp) :: arg1 , arg2 , bequ , bq1 , bq2 , bq3 , c0 , c1 , c2 , c3 , & - ct , d , d0 , d1 , d2, dimob0 , e0 , e1 , e2 , ff , fi , gg , & + d0 , d1 , d2, dimob0 , e0 , e1 , e2 , ff , fi , gg , & hli , oradik , oterm , p(8,100) , r , r1 , r2 , r3 , r3h , radik , & - rlat , rlon , rq , st , step12 , step2 , & - stp , t , term , xx , z , zq , zz + rq , step12 , step2 , stp , t , term , xx , z , zq , zz integer :: i , iequ , n real(wp),parameter :: rmin = 0.05_wp !! boundaries for identification of `icode=2 and 3` @@ -376,15 +375,7 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v) me%xi(2) = v(2) me%xi(3) = v(3) else - rlat = glat*umr - rlon = glon*umr - ct = sin(rlat) - st = cos(rlat) - d = sqrt(aquad-(aquad-bquad)*ct*ct) - me%xi(1) = (alt+aquad/d)*st/era - me%xi(3) = (alt+bquad/d)*ct/era - me%xi(2) = me%xi(1)*sin(rlon) - me%xi(1) = me%xi(1)*cos(rlon) + me%xi = geo_to_cart(glat,glon,alt) end if !*****convert to dipol-oriented co-ordinates @@ -642,19 +633,21 @@ subroutine feldg(me,glat,glon,alt,bnorth,beast,bdown,babs) x , xxx , y , yyy , z , zzz integer :: i , ih , ihmax , il , imax , k , last , m - rlat=glat*umr - ct=sin(rlat) - st=cos(rlat) - d=sqrt(aquad-(aquad-bquad)*ct*ct) - rlon=glon*umr - cp=cos(rlon) - sp=sin(rlon) - zzz=(alt+bquad/d)*ct/era - rho=(alt+aquad/d)*st/era - xxx=rho*cp - yyy=rho*sp - - rq=1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz) + ! same calculation as geo_to_cart, but not used here + ! because the intermediate variables are also used below. + rlat = glat*umr + ct = sin(rlat) + st = cos(rlat) + d = sqrt(aquad-(aquad-bquad)*ct*ct) + rlon = glon*umr + cp = cos(rlon) + sp = sin(rlon) + zzz = (alt+bquad/d)*ct/era + rho = (alt+aquad/d)*st/era + xxx = rho*cp + yyy = rho*sp + + rq = 1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz) me%xi = [xxx,yyy,zzz] * rq ihmax=me%nmax*me%nmax+1 @@ -1130,4 +1123,36 @@ subroutine extrashc(me,date,dte1,nmax1,gh1,nmax2,gh2,nmax,gh) end subroutine extrashc +!***************************************************************************************** +!> +! geodetic to scaled cartesian coordinates + +pure function geo_to_cart(glat,glon,alt) result(x) + + real(wp),intent(in) :: glat !! geodetic latitude in degrees (north) + real(wp),intent(in) :: glon !! geodetic longitude in degrees (east) + real(wp),intent(in) :: alt !! altitude in km above sea level + real(wp),dimension(3) :: x !! cartesian coordinates in earth radii (6371.2 km) + !! + !! * x-axis pointing to equator at 0 longitude + !! * y-axis pointing to equator at 90 long. + !! * z-axis pointing to north pole + + real(wp) :: rlat !! latitude in radians + real(wp) :: rlon !! longitude in radians + real(wp) :: d, rho + + ! deg to radians: + rlat = glat*umr + rlon = glon*umr + + ! JW : it's weird that ct is sin, and st is cos...it was like that in the original code + associate (ct => sin(rlat), st => cos(rlat), cp => cos(rlon), sp => sin(rlon)) + d = sqrt(aquad-(aquad-bquad)*ct*ct) + rho = (alt+aquad/d)*st/era + x = [rho*cp, rho*sp, (alt+bquad/d)*ct/era] + end associate + +end function geo_to_cart + end module SHELLIG_module \ No newline at end of file From 4d73bd1d28fe2c7dabf0332bb4ff3bf15a74d53b Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 13 Feb 2024 09:09:49 -0600 Subject: [PATCH 02/10] fix module name --- src/redbelt_module.f90 | 4 ++-- test/radbelt_test.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/redbelt_module.f90 b/src/redbelt_module.f90 index 930848f..8a9ae1f 100644 --- a/src/redbelt_module.f90 +++ b/src/redbelt_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 diff --git a/test/radbelt_test.f90 b/test/radbelt_test.f90 index 7151572..1fdcac6 100644 --- a/test/radbelt_test.f90 +++ b/test/radbelt_test.f90 @@ -2,7 +2,7 @@ program radbelt_test !! comparison to the python radbelt example - use redbelt_module + use radbelt_module use radbelt_kinds_module implicit none From 41a1417b30a98675adca22489942cb1ecf6c7907 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 13 Feb 2024 09:10:01 -0600 Subject: [PATCH 03/10] gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c2607a5..918a3a2 100644 --- a/.gitignore +++ b/.gitignore @@ -34,6 +34,7 @@ # misc .DS_Store /results.txt +/env # directories /build From 86f5c1771e7597a30a50ac8ccf33cb8517f45101 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 13 Feb 2024 10:24:36 -0600 Subject: [PATCH 04/10] renamed file --- src/{redbelt_module.f90 => radbelt_module.f90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{redbelt_module.f90 => radbelt_module.f90} (100%) diff --git a/src/redbelt_module.f90 b/src/radbelt_module.f90 similarity index 100% rename from src/redbelt_module.f90 rename to src/radbelt_module.f90 From bf0d1d91dc91f0914064dfdf2c2284c83578febd Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 13 Feb 2024 16:08:56 -0600 Subject: [PATCH 05/10] experimenting with a python interface work in progress --- build_python.sh | 15 +++++++ meson.build | 26 +++++++++++ python/.f2py_f2cmap | 13 ++++++ python/environment.yml | 20 +++++++++ python/radbeltpy.pyf | 27 ++++++++++++ src/radbelt_c_module.f90 | 93 +++++++++++++++++++++++++++++++++++++++ src/shellig.f90 | 1 + src/trmfun.f90 | 5 +++ test/radbeltpy_test.py | 94 ++++++++++++++++++++++++++++++++++++++++ 9 files changed, 294 insertions(+) create mode 100755 build_python.sh create mode 100644 meson.build create mode 100644 python/.f2py_f2cmap create mode 100644 python/environment.yml create mode 100644 python/radbeltpy.pyf create mode 100644 src/radbelt_c_module.f90 create mode 100644 test/radbeltpy_test.py diff --git a/build_python.sh b/build_python.sh new file mode 100755 index 0000000..336a614 --- /dev/null +++ b/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 ./python/radbeltpy.pyf \ + ./src/radbelt_kinds_module.F90 \ + ./src/trmfun.f90 \ + ./src/shellig.f90 \ + ./src/radbelt_module.f90 \ + ./src/radbelt_c_module.f90 \ No newline at end of file diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..faeab68 --- /dev/null +++ b/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, +) diff --git a/python/.f2py_f2cmap b/python/.f2py_f2cmap new file mode 100644 index 0000000..4f5336c --- /dev/null +++ b/python/.f2py_f2cmap @@ -0,0 +1,13 @@ +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_long_long="long long"), +) diff --git a/python/environment.yml b/python/environment.yml new file mode 100644 index 0000000..e7990ee --- /dev/null +++ b/python/environment.yml @@ -0,0 +1,20 @@ +# To use: +# +# $ conda env create --prefix ./env -f ./python/environment.yml +# $ conda activate ./env + +channels: + - nodefaults + - conda-forge + +dependencies: + #- compilers + #- gfortran + - fortran-compiler + - meson + - pkg-config + - fpm + - ford + - graphviz + - fprettify + - numpy diff --git a/python/radbeltpy.pyf b/python/radbeltpy.pyf new file mode 100644 index 0000000..6967b34 --- /dev/null +++ b/python/radbeltpy.pyf @@ -0,0 +1,27 @@ +! -*- 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 + subroutine get_flux_g_c(lon,lat,height,year,e,imname,flux) ! in :radbelt:radbelt_c_module.f90:radbelt_c_module + 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(aep8_dir, igrf_dir) ! in :radbelt:radbelt_c_module.f90:radbelt_c_module + 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 + end module radbelt_c_module + end interface +end python module radbelt + +! This file was manually created diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 new file mode 100644 index 0000000..5d7d754 --- /dev/null +++ b/src/radbelt_c_module.f90 @@ -0,0 +1,93 @@ +!***************************************************************************************** +!> +! Experimental C interface to the radbelt module. +! +!@note Currently, this makes use of a global variable to store the state, +! so should not be considered threadsafe. + +module radbelt_c_module + +use iso_c_binding, only: c_double, c_int, c_char, c_null_char +use radbelt_module, only: radbelt_type + +implicit none + +type(radbelt_type) :: radbelt !! global variable + +integer,parameter :: STR_LEN = 256 !! string length for paths + +contains + +! function c2f_str(cstr) result(fstr) +! character(kind=c_char,len=1),dimension(*),intent(in) :: cstr +! character(len=:),allocatable :: fstr +! integer :: i +! fstr = '' +! i = 0 +! do +! i = i + 1 +! if (cstr(i)==c_null_char) exit +! fstr = fstr // cstr(i) +! end do +! end function c2f_str + +!***************************************************************************************** +!> +! 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") + + 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_ + + c2f : block + ! convert to c strings to fortran + integer :: i + aep8_dir_ = '' + igrf_dir_ = '' + do i = 1, STR_LEN + aep8_dir_ = aep8_dir_//aep8_dir(i) + igrf_dir_ = igrf_dir_//igrf_dir(i) + end do + aep8_dir_ = trim(aep8_dir_) + igrf_dir_ = trim(igrf_dir_) + end block c2f + + ! !... 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_) + + 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") + + 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. + + flux = radbelt%get_flux(lon,lat,height,year,e,imname) + +end subroutine get_flux_g_c + +end module radbelt_c_module diff --git a/src/shellig.f90 b/src/shellig.f90 index 20846ba..68c65ce 100644 --- a/src/shellig.f90 +++ b/src/shellig.f90 @@ -214,6 +214,7 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0) step=stps irun=0 + rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings main : do irun=irun+1 diff --git a/src/trmfun.f90 b/src/trmfun.f90 index 3c784e9..5565c31 100644 --- a/src/trmfun.f90 +++ b/src/trmfun.f90 @@ -152,6 +152,11 @@ subroutine trara1(me,descr,map,fl,bb0,e,f,n) integer :: i0 , i1 , i2 , i3 , ie , l3 , nb , nl logical :: s0 , s1 , s2 + e0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings + f0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings + i0 = 0 ! to avoid -Wmaybe-uninitialized warnings + s0 = .false. ! to avoid -Wmaybe-uninitialized warnings -- but not sure what default value here should be ! -JW + bb0_ = bb0 me%fistep = descr(7)/descr(2) escale = descr(4) diff --git a/test/radbeltpy_test.py b/test/radbeltpy_test.py new file mode 100644 index 0000000..bfc940b --- /dev/null +++ b/test/radbeltpy_test.py @@ -0,0 +1,94 @@ +# +# Test of the f2py-created Python interface. +# + +from pathlib import Path + +dir = Path(__file__).resolve().parents[1] # root directory + +import sys +sys.path.insert(0,'.') # assuming the radbelt lib is in the main directory + +import radbelt +import time +import os +# import numpy as np + +######################################################################################################### +def set_data_files_paths(aep8_dir : str, igrf_dir : str) -> None: + """Python function to set the file paths""" + radbelt.radbelt_c_module.set_data_files_paths_c(aep8_dir, igrf_dir) + +######################################################################################################### +#@np.vectorize <-- makes it run slower for scalar inputs?? +def get_flux(lon : float, lat : float , height : float, year : float, e : float, imname :int) -> float: + """Python function to get the flux""" + return radbelt.radbelt_c_module.get_flux_g_c(lon,lat,height,year,e,imname) + + +# set location of the data files: +# dir = Path(__file__).resolve().parents[1] # root directory +aep8_dir = str(dir / 'data' / 'aep8') +igrf_dir = str(dir / 'data' / 'igrf') +set_data_files_paths(aep8_dir, igrf_dir) + +EPS = sys.float_info.epsilon # machine precision for error checking +lon = -45.0 +lat = -30.0 +height = 500.0 +year = 2021.1616438356164 # decimal year +imname = 4 # 'p', 'max' +e = 20.0 + +flux = get_flux(lon,lat,height,year,e,imname) + +print(f'flux = {flux}') + +error = flux - 2642.50370051985726336559603128948869 #difference from real128 version +relerror = abs(error/flux) + +print(f'Flux = {flux}') +print(f'Error = {error}') +print(f'Rel Error = {relerror}') +if relerror>10*EPS: + raise Exception('error') + +print('') + +n_cases = 0 +tstart = time.time() +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 + f = 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)') + +# 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)') + From 933fa5bdede4187bf26286e4be3125d5d95ccd6d Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Tue, 13 Feb 2024 17:10:51 -0600 Subject: [PATCH 06/10] python updates --- build_python.sh | 15 --------------- python/build_python.sh | 15 +++++++++++++++ python/environment.yml | 2 -- test.sh | 15 +++++++++++++++ test/radbeltpy_test.py | 13 ++++--------- 5 files changed, 34 insertions(+), 26 deletions(-) delete mode 100755 build_python.sh create mode 100755 python/build_python.sh create mode 100755 test.sh diff --git a/build_python.sh b/build_python.sh deleted file mode 100755 index 336a614..0000000 --- a/build_python.sh +++ /dev/null @@ -1,15 +0,0 @@ -#!/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 ./python/radbeltpy.pyf \ - ./src/radbelt_kinds_module.F90 \ - ./src/trmfun.f90 \ - ./src/shellig.f90 \ - ./src/radbelt_module.f90 \ - ./src/radbelt_c_module.f90 \ No newline at end of file diff --git a/python/build_python.sh b/python/build_python.sh new file mode 100755 index 0000000..c81cfae --- /dev/null +++ b/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 \ No newline at end of file diff --git a/python/environment.yml b/python/environment.yml index e7990ee..5a2c322 100644 --- a/python/environment.yml +++ b/python/environment.yml @@ -8,8 +8,6 @@ channels: - conda-forge dependencies: - #- compilers - #- gfortran - fortran-compiler - meson - pkg-config diff --git a/test.sh b/test.sh new file mode 100755 index 0000000..1c272f8 --- /dev/null +++ b/test.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +# +# basic script to build a conda environment, build the python extension, and test it +# + +rm -rf ./env +conda env create --prefix ./env -f ./python/environment.yml +conda activate ./env + +cd python +chmod +x ./build_python.sh +./build_python.sh +cd .. +python ./test/radbeltpy_test.py diff --git a/test/radbeltpy_test.py b/test/radbeltpy_test.py index bfc940b..0ea1076 100644 --- a/test/radbeltpy_test.py +++ b/test/radbeltpy_test.py @@ -2,17 +2,13 @@ # Test of the f2py-created Python interface. # +import sys +import time from pathlib import Path dir = Path(__file__).resolve().parents[1] # root directory - -import sys -sys.path.insert(0,'.') # assuming the radbelt lib is in the main directory - -import radbelt -import time -import os -# import numpy as np +sys.path.insert(0,str(dir / 'python')) # assuming the radbelt lib is in python directory +import radbelt # this is the module being tested ######################################################################################################### def set_data_files_paths(aep8_dir : str, igrf_dir : str) -> None: @@ -27,7 +23,6 @@ def get_flux(lon : float, lat : float , height : float, year : float, e : float, # set location of the data files: -# dir = Path(__file__).resolve().parents[1] # root directory aep8_dir = str(dir / 'data' / 'aep8') igrf_dir = str(dir / 'data' / 'igrf') set_data_files_paths(aep8_dir, igrf_dir) From 7b29933c6cb380848e66e70f7e5a4b3c0652f76d Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Wed, 14 Feb 2024 00:15:17 -0600 Subject: [PATCH 07/10] experimenting with a python class wrapper for the fortran class --- python/.f2py_f2cmap | 1 + python/radbeltpy.pyf | 19 +++++++++++-- src/radbelt_c_module.f90 | 61 +++++++++++++++++++++++++++++++++++----- test/radbeltpy_test.py | 34 ++++++++++++++++++++-- 4 files changed, 102 insertions(+), 13 deletions(-) diff --git a/python/.f2py_f2cmap b/python/.f2py_f2cmap index 4f5336c..a96534b 100644 --- a/python/.f2py_f2cmap +++ b/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"), ) diff --git a/python/radbeltpy.pyf b/python/radbeltpy.pyf index 6967b34..73f1cdb 100644 --- a/python/radbeltpy.pyf +++ b/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 diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 index 5d7d754..6981638 100644 --- a/src/radbelt_c_module.f90 +++ b/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 diff --git a/test/radbeltpy_test.py b/test/radbeltpy_test.py index 0ea1076..0dc9e5e 100644 --- a/test/radbeltpy_test.py +++ b/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)') From 27f4593141914724d3925381a5aeda17021318b1 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Wed, 14 Feb 2024 09:34:00 -0600 Subject: [PATCH 08/10] commenting --- src/radbelt_c_module.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 index 6981638..aa7127f 100644 --- a/src/radbelt_c_module.f90 +++ b/src/radbelt_c_module.f90 @@ -1,9 +1,6 @@ !***************************************************************************************** !> ! Experimental C interface to the radbelt module. -! -!@note Currently, this makes use of a global variable to store the state, -! so should not be considered threadsafe. module radbelt_c_module @@ -13,6 +10,7 @@ module radbelt_c_module implicit none +!... how can we eliminate this ? integer,parameter :: STR_LEN = 256 !! string length for paths contains From 87e9a990bfbf02188e6968d759a94c7789a44691 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Wed, 14 Feb 2024 17:13:22 -0600 Subject: [PATCH 09/10] some reafactoring elimination of some warnings --- ford.md | 1 + python/radbeltpy.pyf | 6 +- src/radbelt_c_module.f90 | 147 ++++++++------- src/shellig.f90 | 19 +- src/trmfun.f90 | 376 ++++++++++++++++++++++----------------- test/radbelt_test.f90 | 2 +- 6 files changed, 308 insertions(+), 243 deletions(-) diff --git a/ford.md b/ford.md index 5a194d3..2fc5007 100644 --- a/ford.md +++ b/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!} diff --git a/python/radbeltpy.pyf b/python/radbeltpy.pyf index 73f1cdb..4fc1198 100644 --- a/python/radbeltpy.pyf +++ b/python/radbeltpy.pyf @@ -19,10 +19,8 @@ python module radbelt ! in 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 + 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 diff --git a/src/radbelt_c_module.f90 b/src/radbelt_c_module.f90 index aa7127f..31c038d 100644 --- a/src/radbelt_c_module.f90 +++ b/src/radbelt_c_module.f90 @@ -2,59 +2,48 @@ !> ! Experimental C interface to the radbelt module. -module radbelt_c_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 + 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 + implicit none -!... how can we eliminate this ? -integer,parameter :: STR_LEN = 256 !! string length for paths + ! can we eliminate this ? + integer,parameter :: STR_LEN = 4096 !! string length for paths -contains + contains +!***************************************************************************************** -! function c2f_str(cstr) result(fstr) -! character(kind=c_char,len=1),dimension(*),intent(in) :: cstr -! character(len=:),allocatable :: fstr -! integer :: i -! fstr = '' -! i = 0 -! do -! i = i + 1 -! if (cstr(i)==c_null_char) exit -! fstr = fstr // cstr(i) -! end do -! end function c2f_str +!***************************************************************************************** +!> +! Convert C string to Fortran -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 +function c2f_str(cstr) result(fstr) - allocate(p) - cp = c_loc(p) - ipointer = transfer(cp, 0_c_intptr_t) + character(kind=c_char,len=1),dimension(:),intent(in) :: cstr !! string from C + character(len=:),allocatable :: fstr !! fortran string -end subroutine initialize_c + integer :: i !! counter -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 + fstr = '' + do i = 1, size(cstr) + fstr = fstr//cstr(i) + end do + fstr = trim(fstr) - call int_pointer_to_f_pointer(ipointer,p) - if (associated(p)) deallocate(p) +end function c2f_str -end subroutine destroy_c +!***************************************************************************************** +!> +! Convert an integer pointer to a [[radbelt_type]] pointer. subroutine int_pointer_to_f_pointer(ipointer, p) - integer(c_intptr_t),intent(in) :: ipointer - type(radbelt_type),pointer :: p + integer(c_intptr_t),intent(in) :: ipointer !! integer pointer from C + type(radbelt_type),pointer :: p !! fortran pointer type(c_ptr) :: cp @@ -67,6 +56,36 @@ subroutine int_pointer_to_f_pointer(ipointer, p) 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 @@ -80,28 +99,18 @@ subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir) bind(C, name="se character(len=:),allocatable :: aep8_dir_, igrf_dir_ type(radbelt_type),pointer :: p - c2f : block - ! convert to c strings to fortran - integer :: i - aep8_dir_ = '' - igrf_dir_ = '' - do i = 1, STR_LEN - aep8_dir_ = aep8_dir_//aep8_dir(i) - igrf_dir_ = igrf_dir_//igrf_dir(i) - end do - aep8_dir_ = trim(aep8_dir_) - 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_) + if (associated(p)) then - call p%set_data_files_paths(aep8_dir_, igrf_dir_) + 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 !***************************************************************************************** @@ -120,19 +129,27 @@ subroutine get_flux_g_c(ipointer,lon,lat,height,year,e,imname,flux) bind(C, name !! 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 + !! + !! * 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) - flux = p%get_flux(lon,lat,height,year,e,imname) + 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 +!***************************************************************************************** + end module radbelt_c_module +!***************************************************************************************** \ No newline at end of file diff --git a/src/shellig.f90 b/src/shellig.f90 index 68c65ce..35eed16 100644 --- a/src/shellig.f90 +++ b/src/shellig.f90 @@ -75,7 +75,7 @@ module shellig_module procedure, public :: feldg, feldc procedure, public :: shellg, shellc procedure, public :: findb0 - procedure :: stoer, getshc, intershc, extrashc, feldi + procedure :: stoer, feldi procedure,public :: set_data_file_dir, get_data_file_dir end type shellig_type @@ -873,19 +873,19 @@ subroutine feldcof(me,year,dimo) if (read_file) then ! get igrf coefficients for the boundary years ! [if they have not ready been loaded] - call me%getshc(me%name,me%nmax1,erad,me%g,ier) + call getshc(me%name,me%nmax1,erad,me%g,ier) if ( ier/=0 ) error stop 'error reading file: '//trim(me%name) me%g_cache = me%g ! because it is modified below, we have to cache the original values from the file - call me%getshc(fil2,me%nmax2,erad,me%gh2,ier) + call getshc(fil2,me%nmax2,erad,me%gh2,ier) if ( ier/=0 ) error stop 'error reading file: '//trim(fil2) else me%g = me%g_cache end if !-- determine igrf coefficients for year if ( l<=numye-1 ) then - call me%intershc(year,dte1,me%nmax1,me%g,dte2,me%nmax2,me%gh2,me%nmax,gha) + call intershc(year,dte1,me%nmax1,me%g,dte2,me%nmax2,me%gh2,me%nmax,gha) else - call me%extrashc(year,dte1,me%nmax1,me%g,me%nmax2,me%gh2,me%nmax,gha) + call extrashc(year,dte1,me%nmax1,me%g,me%nmax2,me%gh2,me%nmax,gha) endif !-- determine magnetic dipol moment and coeffiecients g f0 = 0.0_wp @@ -929,9 +929,8 @@ end subroutine feldcof ! * Version 1.01, A. Zunde, USGS, MS 964, ! Box 25046 Federal Center, Denver, CO 80225 -subroutine getshc(me,Fspec,Nmax,Erad,Gh,Ier) +subroutine getshc(Fspec,Nmax,Erad,Gh,Ier) - class(shellig_type),intent(inout) :: me character(len=*),intent(in) :: Fspec !! File specification integer,intent(out) :: Nmax !! Maximum degree and order of model real(wp),intent(out) :: Erad !! Earth's radius associated with the spherical @@ -1023,9 +1022,8 @@ END subroutine getshc ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -subroutine intershc(me,date,dte1,nmax1,gh1,dte2,nmax2,gh2,nmax,gh) +subroutine intershc(date,dte1,nmax1,gh1,dte2,nmax2,gh2,nmax,gh) - class(shellig_type),intent(inout) :: me real(wp),intent(in) :: date !! Date of resulting model (in decimal year) real(wp),intent(in) :: dte1 !! Date of earlier model integer,intent(in) :: nmax1 !! Maximum degree and order of earlier model @@ -1082,9 +1080,8 @@ end subroutine intershc ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -subroutine extrashc(me,date,dte1,nmax1,gh1,nmax2,gh2,nmax,gh) +subroutine extrashc(date,dte1,nmax1,gh1,nmax2,gh2,nmax,gh) - class(shellig_type),intent(inout) :: me real(wp),intent(in) :: date !! Date of resulting model (in decimal year) real(wp),intent(in) :: dte1 !! Date of base model integer,intent(in) :: nmax1 !! Maximum degree and order of base model diff --git a/src/trmfun.f90 b/src/trmfun.f90 index 5565c31..feacd6d 100644 --- a/src/trmfun.f90 +++ b/src/trmfun.f90 @@ -162,9 +162,9 @@ subroutine trara1(me,descr,map,fl,bb0,e,f,n) escale = descr(4) fscale = descr(7) xnl = min(15.6_wp,abs(fl)) - nl = xnl*descr(5) + nl = int(xnl*descr(5)) if ( bb0_<1.0_wp ) bb0_ = 1.0_wp - nb = (bb0_-1.0_wp)*descr(6) + nb = int((bb0_-1.0_wp)*descr(6)) ! i2 is the number of elements in the flux map for the first energy. ! i3 is the index of the last element of the second energy map. @@ -261,192 +261,244 @@ function trara2(me,map,il,ib) fkbm , fll1 , fll2 , flog , flog1 , flog2 , flogm , & fnb , fnl , sl1 , sl2 integer :: i1 , i2 , itime , j1 , j2 , kt , l1 , l2 - integer :: itask + logical :: dummy fistep = me%fistep - itask = 1 - main: do - select case (itask) - case (1) - fnl = il - fnb = ib - itime = 0 - i2 = 0 - do - - ! find consecutive sub-sub-maps for scaled l-values ls1,ls2, - ! with il less or equal ls2. l1,l2 are lengths of sub-sub-maps. - ! i1,i2 are indeces of first elements minus 1. - - l2 = map(i2+1) - if ( map(i2+2)<=il ) then - i1 = i2 - l1 = l2 - i2 = i2 + l2 - - ! if sub-sub-maps are empty, i. e. length less 4, than trara2=0 - - elseif ( (l1<4) .and. (l2<4) ) then - trara2 = 0.0_wp - return - else + !******** + ! to avoid -Wmaybe-uninitialized warning + dfl = 0.0_wp + fincr1 = 0.0_wp + fincr2 = 0.0_wp + fkb = 0.0_wp + fkb1 = 0.0_wp + fkb2 = 0.0_wp + fkbm = 0.0_wp + flog = 0.0_wp + flog1 = 0.0_wp + flog2 = 0.0_wp + flogm = 0.0_wp + fnb = 0.0_wp + fnl = 0.0_wp + sl2 = 0.0_wp + i1 = 0 + i2 = 0 + itime = 0 + j2 = 0 + l1 = 0 + l2 = 0 + !******** + + ! these are recursive functions that + ! replace the gotos in the original code + call task1(dummy) - ! if flog2 less flog1, than ls2 first map and ls1 second map + contains - if ( map(i2+3)<=map(i1+3) ) exit - itask = 3 - cycle main - endif + recursive subroutine task1(done) + logical,intent(out) :: done + done = .false. + fnl = il + fnb = ib + itime = 0 + i2 = 0 + do + ! find consecutive sub-sub-maps for scaled l-values ls1,ls2, + ! with il less or equal ls2. l1,l2 are lengths of sub-sub-maps. + ! i1,i2 are indeces of first elements minus 1. + l2 = map(i2+1) + if ( map(i2+2)<=il ) then + i1 = i2 + l1 = l2 + i2 = i2 + l2 + ! if sub-sub-maps are empty, i. e. length less 4, than trara2=0 + elseif ( (l1<4) .and. (l2<4) ) then + trara2 = 0.0_wp + done = .true. + return + else + ! if flog2 less flog1, than ls2 first map and ls1 second map + if ( map(i2+3)<=map(i1+3) ) exit + call task3(done) + return + endif + enddo + call task2(done) + end subroutine task1 + recursive subroutine task2(done) + logical,intent(out) :: done + done = .false. + kt = i1 + i1 = i2 + i2 = kt + kt = l1 + l1 = l2 + l2 = kt + call task3(done) + end subroutine task2 + recursive subroutine task3(done) + logical,intent(out) :: done + logical :: check + done = .false. + ! determine interpolate in scaled l-value + fll1 = map(i1+2) + fll2 = map(i2+2) + dfl = (fnl-fll1)/(fll2-fll1) + flog1 = map(i1+3) + flog2 = map(i2+3) + fkb1 = 0.0_wp + fkb2 = 0.0_wp + if ( l1>=4 ) then + ! b/b0 loop + check = .true. + do j2 = 4 , l2 + fincr2 = map(i2+j2) + if ( fkb2+fincr2>fnb ) then + check = .false. + exit + end if + fkb2 = fkb2 + fincr2 + flog2 = flog2 - fistep enddo - itask = 2 - case (2) - kt = i1 - i1 = i2 - i2 = kt - kt = l1 - l1 = l2 - l2 = kt - itask = 3 - case (3) - - ! determine interpolate in scaled l-value - - fll1 = map(i1+2) - fll2 = map(i2+2) - dfl = (fnl-fll1)/(fll2-fll1) - flog1 = map(i1+3) - flog2 = map(i2+3) - fkb1 = 0.0_wp - fkb2 = 0.0_wp - if ( l1>=4 ) then - - ! b/b0 loop - - do j2 = 4 , l2 - fincr2 = map(i2+j2) - if ( fkb2+fincr2>fnb ) goto 10 - fkb2 = fkb2 + fincr2 - flog2 = flog2 - fistep - enddo + if (check) then itime = itime + 1 if ( itime==1 ) then - itask = 2 - cycle main + call task2(done) + return endif trara2 = 0.0_wp + done = .true. return - 10 if ( itime/=1 ) then - if ( j2==4 ) then - itask = 4 - cycle main - endif - sl2 = flog2/fkb2 - do j1 = 4 , l1 - fincr1 = map(i1+j1) - fkb1 = fkb1 + fincr1 - flog1 = flog1 - fistep - fkbj1 = ((flog1/fistep)*fincr1+fkb1)/((fincr1/fistep)*sl2+1.0_wp) - if ( fkbj1<=fkb1 ) goto 15 - enddo + end if + if ( itime/=1 ) then + if ( j2==4 ) then + call task4(done) + return + endif + sl2 = flog2/fkb2 + check = .true. + do j1 = 4 , l1 + fincr1 = map(i1+j1) + fkb1 = fkb1 + fincr1 + flog1 = flog1 - fistep + fkbj1 = ((flog1/fistep)*fincr1+fkb1)/((fincr1/fistep)*sl2+1.0_wp) + if ( fkbj1<=fkb1 ) then + check = .false. + exit + end if + enddo + if (check) then if ( fkbj1<=fkb2 ) then trara2 = 0.0_wp + done = .true. return endif - 15 if ( fkbj1<=fkb2 ) then - fkbm = fkbj1 + (fkb2-fkbj1)*dfl - flogm = fkbm*sl2 - flog2 = flog2 - fistep - fkb2 = fkb2 + fincr2 - sl1 = flog1/fkb1 - sl2 = flog2/fkb2 - itask = 5 - cycle main - else - fkb1 = 0.0_wp - endif + end if + if ( fkbj1<=fkb2 ) then + fkbm = fkbj1 + (fkb2-fkbj1)*dfl + flogm = fkbm*sl2 + flog2 = flog2 - fistep + fkb2 = fkb2 + fincr2 + sl1 = flog1/fkb1 + sl2 = flog2/fkb2 + call task5(done) + return + else + fkb1 = 0.0_wp endif - fkb2 = 0.0_wp endif - j2 = 4 - fincr2 = map(i2+j2) - flog2 = map(i2+3) - flog1 = map(i1+3) - itask = 4 - case (4) - flogm = flog1 + (flog2-flog1)*dfl - fkbm = 0.0_wp - fkb2 = fkb2 + fincr2 - flog2 = flog2 - fistep - sl2 = flog2/fkb2 - if ( l1<4 ) then - fincr1 = 0.0_wp - sl1 = -900000.0_wp - itask = 6 - cycle main + fkb2 = 0.0_wp + endif + j2 = 4 + fincr2 = map(i2+j2) + flog2 = map(i2+3) + flog1 = map(i1+3) + call task4(done) + end subroutine task3 + recursive subroutine task4(done) + logical,intent(out) :: done + done = .false. + flogm = flog1 + (flog2-flog1)*dfl + fkbm = 0.0_wp + fkb2 = fkb2 + fincr2 + flog2 = flog2 - fistep + sl2 = flog2/fkb2 + if ( l1<4 ) then + fincr1 = 0.0_wp + sl1 = -900000.0_wp + call task6(done) + return + else + j1 = 4 + fincr1 = map(i1+j1) + fkb1 = fkb1 + fincr1 + flog1 = flog1 - fistep + sl1 = flog1/fkb1 + endif + call task5(done) + end subroutine task4 + recursive subroutine task5(done) + logical,intent(out) :: done + done = .false. + do while ( sl1>=sl2 ) + fkbj2 = ((flog2/fistep)*fincr2+fkb2)/((fincr2/fistep)*sl1+1.0_wp) + fkb = fkb1 + (fkbj2-fkb1)*dfl + flog = fkb*sl1 + if ( fkb>=fnb ) then + call task7(done) + return + endif + fkbm = fkb + flogm = flog + if ( j1>=l1 ) then + trara2 = 0.0_wp + done = .true. + return else - j1 = 4 + j1 = j1 + 1 fincr1 = map(i1+j1) - fkb1 = fkb1 + fincr1 flog1 = flog1 - fistep + fkb1 = fkb1 + fincr1 sl1 = flog1/fkb1 endif - itask = 5 - case (5) - do while ( sl1>=sl2 ) - fkbj2 = ((flog2/fistep)*fincr2+fkb2)/((fincr2/fistep)*sl1+1.0_wp) - fkb = fkb1 + (fkbj2-fkb1)*dfl - flog = fkb*sl1 - if ( fkb>=fnb ) then - itask = 7 - cycle main - endif - fkbm = fkb - flogm = flog - if ( j1>=l1 ) then - trara2 = 0.0_wp - return - else - j1 = j1 + 1 - fincr1 = map(i1+j1) - flog1 = flog1 - fistep - fkb1 = fkb1 + fincr1 - sl1 = flog1/fkb1 - endif - enddo - itask = 6 - case (6) - fkbj1 = ((flog1/fistep)*fincr1+fkb1)/((fincr1/fistep)*sl2+1.0_wp) - fkb = fkbj1 + (fkb2-fkbj1)*dfl - flog = fkb*sl2 - if ( fkb=l2 ) then - trara2 = 0.0_wp - return - else - j2 = j2 + 1 - fincr2 = map(i2+j2) - flog2 = flog2 - fistep - fkb2 = fkb2 + fincr2 - sl2 = flog2/fkb2 - itask = 5 - cycle main - endif - endif - itask = 7 - case (7) - if ( fkb=l2 ) then trara2 = 0.0_wp + done = .true. + return else - trara2 = flogm + (flog-flogm)*((fnb-fkbm)/(fkb-fkbm)) - trara2 = max(trara2,0.0_wp) + j2 = j2 + 1 + fincr2 = map(i2+j2) + flog2 = flog2 - fistep + fkb2 = fkb2 + fincr2 + sl2 = flog2/fkb2 + call task5(done) return endif - exit main - end select - enddo main + endif + call task7(done) + end subroutine task6 + recursive subroutine task7(done) + logical,intent(out) :: done + if ( fkb Date: Fri, 16 Feb 2024 11:46:20 -0600 Subject: [PATCH 10/10] readme update --- README.md | 41 ++++++++++++++--------------------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index 2d517d8..011b6e8 100755 --- a/README.md +++ b/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, dbilitza@pop600.gsfc.nasa.gov +* 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.