Skip to content

Commit

Permalink
eliminated max str len var for c/python interface
Browse files Browse the repository at this point in the history
  • Loading branch information
@jacobwilliams
jacobwilliams committed Feb 16, 2024
1 parent f822a67 commit 2e10bff
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 47 deletions.
8 changes: 5 additions & 3 deletions python/radbeltpy.pyf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,12 @@ 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) ! in :radbelt:radbelt_c_module.f90:radbelt_c_module
subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m) ! 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
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
Expand Down
11 changes: 5 additions & 6 deletions src/radbelt_c_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,6 @@ module radbelt_c_module

implicit none

! can we eliminate this ?
integer,parameter :: STR_LEN = 4096 !! string length for paths

contains
!*****************************************************************************************

Expand Down Expand Up @@ -90,11 +87,13 @@ 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")
subroutine set_data_files_paths_c(ipointer, aep8_dir, igrf_dir, n, m) 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
integer(c_int),intent(in) :: n !! size of `aep8_dir`
character(kind=c_char,len=1),dimension(n),intent(in) :: aep8_dir
integer(c_int),intent(in) :: m !! size of `igrf_dir`
character(kind=c_char,len=1),dimension(m),intent(in) :: igrf_dir

character(len=:),allocatable :: aep8_dir_, igrf_dir_
type(radbelt_type),pointer :: p
Expand Down
73 changes: 36 additions & 37 deletions src/shellig.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0)
value=.false.
exit main
endif
!*********************first three points
! first three points
p(1,2)=me%sp(1)
p(2,2)=me%sp(2)
p(3,2)=me%sp(3)
Expand All @@ -240,7 +240,7 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0)
p(2,3)=p(2,2)+step*(20.0_wp*p(5,3)-3.*p(5,2)+p(5,1))/18.0_wp
p(3,3)=p(3,2)+step
call me%stoer(p(1,3),bq3,r3)
!******************invert sense if required
! invert sense if required
if (bq3>bq1) then
step=-step
r3=r1
Expand All @@ -251,18 +251,18 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0)
p(i,3)=zz
end do
end if
!******************initialization
! initialization
step12=step/12.0_wp
value=.true.
bmin=1.0e4_wp
bold=1.0e4_wp
!******************corrector (field line tracing)
! corrector (field line tracing)
n=0
corrector : do
p(1,3)=p(1,2)+step12*(5.0_wp*p(4,3)+8.0_wp*p(4,2)-p(4,1))
n=n+1
p(2,3)=p(2,2)+step12*(5.0_wp*p(5,3)+8.0_wp*p(5,2)-p(5,1))
!******************predictor (field line tracing)
! predictor (field line tracing)
p(1,4)=p(1,3)+step12*(23.0_wp*p(4,3)-16.0_wp*p(4,2)+5.0_wp*p(4,1))
p(2,4)=p(2,3)+step12*(23.0_wp*p(5,3)-16.0_wp*p(5,2)+5.0_wp*p(5,1))
p(3,4)=p(3,3)+step
Expand Down Expand Up @@ -379,13 +379,13 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
me%xi = geo_to_cart(glat,glon,alt)
end if

!*****convert to dipol-oriented co-ordinates
! convert to dipol-oriented co-ordinates
rq = 1.0_wp/(me%xi(1)*me%xi(1)+me%xi(2)*me%xi(2)+me%xi(3)*me%xi(3))
r3h = sqrt(rq*sqrt(rq))
p(1,2) = (me%xi(1)*u(1,1)+me%xi(2)*u(2,1)+me%xi(3)*u(3,1))*r3h
p(2,2) = (me%xi(1)*u(1,2)+me%xi(2)*u(2,2))*r3h
p(3,2) = (me%xi(1)*u(1,3)+me%xi(2)*u(2,3)+me%xi(3)*u(3,3))*rq
! *****first three points of field line
! first three points of field line
me%step = -sign(me%step,p(3,2))
call me%stoer(p(1,2),bq2,r2)
b0 = sqrt(bq2)
Expand All @@ -401,7 +401,7 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
p(2,3) = p(2,2) + me%step*(20.0_wp*p(5,3)-3.*p(5,2)+p(5,1))/18.0_wp
p(3,3) = p(3,2) + me%step
call me%stoer(p(1,3),bq3,r3)
!*****invert sense if required
! invert sense if required
if ( bq3>bq1 ) then
me%step = -me%step
r3 = r1
Expand All @@ -412,7 +412,7 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
p(i,3) = zz
enddo
endif
!*****search for lowest magnetic field strength
! search for lowest magnetic field strength
if ( bq1<bequ ) then
bequ = bq1
iequ = 1
Expand All @@ -425,7 +425,7 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
bequ = bq3
iequ = 3
endif
!*****initialization of integration loops
! initialization of integration loops
step12 = me%step/12.0_wp
step2 = me%step + me%step
me%steq = sign(me%steq,me%step)
Expand All @@ -438,13 +438,13 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
stp = stp/0.75_wp
p(8,1) = step2*(p(1,1)*p(4,1)+p(2,1)*p(5,1))
p(8,2) = step2*(p(1,2)*p(4,2)+p(2,2)*p(5,2))
!*****main loop (field line tracing)
! main loop (field line tracing)
main: do n = 3 , max_loop_index
!*****corrector (field line tracing)
! corrector (field line tracing)
p(1,n) = p(1,n-1) + step12*(5.0_wp*p(4,n)+8.0_wp*p(4,n-1)-p(4,n-2))
p(2,n) = p(2,n-1) + step12*(5.0_wp*p(5,n)+8.0_wp*p(5,n-1)-p(5,n-2))
!*****prepare expansion coefficients for interpolation
!*****of slowly varying quantities
! prepare expansion coefficients for interpolation
! of slowly varying quantities
p(8,n) = step2*(p(1,n)*p(4,n)+p(2,n)*p(5,n))
c0 = p(1,n-1)**2 + p(2,n-1)**2
c1 = p(8,n-1)
Expand All @@ -457,15 +457,15 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
e1 = (p(7,n)-p(7,n-2))*0.5_wp
e2 = (p(7,n)+p(7,n-2)-e0-e0)*0.5_wp
inner: do
!*****inner loop (for quadrature)
! inner loop (for quadrature)
t = (z-p(3,n-1))/me%step
if ( t>1.0_wp ) then
!*****predictor (field line tracing)
! predictor (field line tracing)
p(1,n+1) = p(1,n) + step12*(23.0_wp*p(4,n)-16.0_wp*p(4,n-1)+5.0_wp*p(4,n-2))
p(2,n+1) = p(2,n) + step12*(23.0_wp*p(5,n)-16.0_wp*p(5,n-1)+5.0_wp*p(5,n-2))
p(3,n+1) = p(3,n) + me%step
call me%stoer(p(1,n+1),bq3,r3)
!*****search for lowest magnetic field strength
! search for lowest magnetic field strength
if ( bq3<bequ ) then
iequ = n + 1
bequ = bq3
Expand All @@ -476,7 +476,7 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
zq = z*z
r = hli + sqrt(hli*hli+zq)
if ( r<=rmin ) then
!*****approximation for high values of l.
! approximation for high values of l.
icode = 3
t = -p(3,n-1)/me%step
fl = 1.0_wp/(abs(((c3*t+c2)*t+c1)*t+c0)+1.0e-15_wp)
Expand Down Expand Up @@ -504,21 +504,19 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v)
me%sp(2) = p(2,iequ-1)
me%sp(3) = p(3,iequ-1)
if ( oradik>=1.0e-15_wp ) fi = fi + stp/0.75_wp*oterm*oradik/(oradik-radik)
!
!-- the minimal allowable value of fi was changed from 1e-15 to 1e-12,
!-- because 1e-38 is the minimal allowable arg. for alog in our envir.
!-- d. bilitza, nov 87.
!

! the minimal allowable value of fi was changed from 1e-15 to 1e-12,
! because 1e-38 is the minimal allowable arg. for alog in our envir.
! d. bilitza, nov 87.
fi = 0.5_wp*abs(fi)/sqrt(b0) + 1.0e-12_wp
!*****compute l from b and i. same as carmel in invar.
!
!-- correct dipole moment is used here. d. bilitza, nov 87.
!

! compute l from b and i. same as carmel in invar.
! correct dipole moment is used here. d. bilitza, nov 87.
dimob0 = dimo/b0
arg1 = log(fi)
arg2 = log(dimob0)
! arg = fi*fi*fi/dimob0
! if(abs(arg)>88.0_wp) arg=88.0_wp
! arg = fi*fi*fi/dimob0
! if(abs(arg)>88.0_wp) arg=88.0_wp
xx = 3*arg1 - arg2
if ( xx>23.0_wp ) then
gg = xx - 3.0460681_wp
Expand Down Expand Up @@ -566,39 +564,40 @@ subroutine stoer(me,p,bq,r)
real(wp) :: dr , dsq , dx , dxm , dy , dym , dz , &
dzm , fli , q , rq , wr , xm , ym , zm

!*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES
! xm,ym,zm are geomagnetic cartesian inverse co-ordinates
zm = P(3)
fli = P(1)*P(1) + P(2)*P(2) + 1.0e-15_wp
R = 0.5_wp*(fli+sqrt(fli*fli+(zm+zm)**2))
rq = R*R
wr = sqrt(R)
xm = P(1)*wr
ym = P(2)*wr
!*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM
! transform to geographic co-ordinate system
me%Xi(1) = xm*u(1,1) + ym*u(1,2) + zm*u(1,3)
me%Xi(2) = xm*u(2,1) + ym*u(2,2) + zm*u(2,3)
me%Xi(3) = xm*u(3,1) + zm*u(3,3)
!*****COMPUTE DERIVATIVES
! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results
! are the same; dkb Feb 1998.
! JW : feb 2024 : xi, h now class variables.
! compute derivatives
! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results
! are the same; dkb Feb 1998.
! JW : feb 2024 : xi, h now class variables.
call me%feldi()
q = me%H(1)/rq
dx = me%H(3) + me%H(3) + q*me%Xi(1)
dy = me%H(4) + me%H(4) + q*me%Xi(2)
dz = me%H(2) + me%H(2) + q*me%Xi(3)
!*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM
! transform back to geomagnetic co-ordinate system
dxm = u(1,1)*dx + u(2,1)*dy + u(3,1)*dz
dym = u(1,2)*dx + u(2,2)*dy
dzm = u(1,3)*dx + u(2,3)*dy + u(3,3)*dz
dr = (xm*dxm+ym*dym+zm*dzm)/R
!*****FORM SLOWLY VARYING EXPRESSIONS
! form slowly varying expressions
P(4) = (wr*dxm-0.5_wp*P(1)*dr)/(R*dzm)
P(5) = (wr*dym-0.5_wp*P(2)*dr)/(R*dzm)
dsq = rq*(dxm*dxm+dym*dym+dzm*dzm)
Bq = dsq*rq*rq
P(6) = sqrt(dsq/(rq+3.0_wp*zm*zm))
P(7) = P(6)*(rq+zm*zm)/(rq*dzm)

end subroutine stoer

!*****************************************************************************************
Expand Down
2 changes: 1 addition & 1 deletion test/radbeltpy_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ def __del__(self) -> None:
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)
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"""
Expand Down

0 comments on commit 2e10bff

Please sign in to comment.