istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/utils.f90: clean up

Browse Source 
- Replace the `get_free_unit` subroutine with the built-in
  `newutin` option of the `open` statement.

- Replace `locatex` with just `locate` + an index offset.

- Replace `inside` with `contour%contains`.

- Merge `vmaxmin` and `vmaxmini` into a single subroutine
  with optional arguments.

- Remove unused `range2rect`, `bubble`.
This commit is contained in:
Michele Guerini Rocco 2024-09-11 11:52:42 +02:00 committed by rnhmjoj
parent 751cca3bfc
commit d5c81268de
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
10 changed files with 190 additions and 408 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
src/beams.f90
View File

@ -5,7 +5,7 @@ module beams
contains
subroutine read_beam0(params, err, unit)
subroutine read_beam0(params, err)
! Reads the wave launcher parameters for the simple case
! where w(z) and 1/R(z) are fixed.
!
@ -30,21 +30,17 @@ contains
use const_and_precisions, only : pi, vc=>ccgs_
use gray_params, only : antenna_parameters
use utils, only : get_free_unit
use logger, only : log_error
! subroutine arguments
type(antenna_parameters), intent(inout) :: params
integer, intent(out) :: err
integer, intent(in), optional :: unit
type(antenna_parameters), intent(inout) :: params
integer, intent(out) :: err
! local variables
integer :: u
real(wp_) :: k0, w0(2), d0(2), z_R(2), phi
u = get_free_unit(unit)
open(unit=u, file=params%filenm, status='old', action='read', iostat=err)
open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
if (err /= 0) then
call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
mod='beams', proc="read_beam0")
@ -79,7 +75,7 @@ contains
end subroutine read_beam0
subroutine read_beam1(params, err, unit)
subroutine read_beam1(params, err)
! Reads the wave launcher parameters for the case
! where w(z, α) and 1/R(z, α) depend on the launcher angle α.
!
@ -100,13 +96,12 @@ contains
use gray_params, only : antenna_parameters
use splines, only : spline_simple
use utils, only : get_free_unit,locate
use utils, only : locate
use logger, only : log_error
! subroutine arguments
type(antenna_parameters), intent(inout) :: params
integer, intent(out) :: err
integer, optional, intent(in) :: unit
! local variables
integer :: u, nisteer, i, k, ii
@ -118,9 +113,7 @@ contains
rci1, rci2, phi1, phi2, &
x0, y0, z0
u = get_free_unit(unit)
open(unit=u, file=params%filenm, status='old', action='read', iostat=err)
open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
if (err /= 0) then
call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
mod='beams', proc="read_beam1")
@ -164,7 +157,7 @@ contains
call z0%init(alphastv, z00v, nisteer)
if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
call locate(alphastv, nisteer, params%alpha , k)
call locate(alphastv, params%alpha, k)
dal = params%alpha - alphastv(k)
params%beta = beta%raw_eval(k, dal)
params%pos(1) = x0%raw_eval(k, dal)
@ -211,7 +204,7 @@ contains
end subroutine read_beam1
subroutine read_beam2(params, beamid, err, unit)
subroutine read_beam2(params, beamid, err)
! Reads the wave launcher parameters for the general case
! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β.
!
@ -232,7 +225,8 @@ contains
! 5. Each line stores one record with the same fields as in the 1D format.
use gray_params, only : antenna_parameters
use utils, only : get_free_unit, intlin, locate, inside
use utils, only : linear_interp, locate
use types, only : contour
use dierckx, only : curfit, splev, surfit, bispev
use logger, only : log_error
@ -240,7 +234,6 @@ contains
type(antenna_parameters), intent(inout) :: params
integer, intent(in) :: beamid
integer, intent(out), optional :: err
integer, intent(in), optional :: unit
! local variables
character(len=20) :: beamname
@ -259,17 +252,14 @@ contains
ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2, &
txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, &
txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, &
cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, &
xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, &
ypolygC, xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC
cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2
type(contour) :: polyg, polygA, polygB, polygC, polygD, outA, outB, outC
real(wp_), DIMENSION(4) :: xvert, yvert
real(wp_), dimension(1) :: fi
integer, parameter :: kspl=1
real(wp_), parameter :: sspl=0.01_wp_
u = get_free_unit(unit)
open(unit=u, file=params%filenm, status='old', action='read', iostat=err)
open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
if (err /= 0) then
call log_error('opening beams file ('//trim(params%filenm)//') failed!', &
mod='beams', proc="read_beam1")
@ -470,7 +460,7 @@ contains
cy0(nxest*nyest), txy0(nmax), tyy0(nmax), &
cz0(nxest*nyest), txz0(nmax), tyz0(nmax), &
wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), &
xpolyg(npolyg), ypolyg(npolyg), w(nisteer))
polyg%R(npolyg), polyg%z(npolyg), w(nisteer))
!
w = 1.d0
!
@ -511,26 +501,26 @@ contains
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
! data range polygon
xpolyg(1:nxcoord) = xcoord(1:nxcoord)
ypolyg(1:nxcoord) = ycoord(1:nxcoord)
polyg%R(1:nxcoord) = xcoord(1:nxcoord)
polyg%z(1:nxcoord) = ycoord(1:nxcoord)
!
! c====================================================================================
do i=1,nycoord-2
xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord)
xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1)
ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord)
ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1)
polyg%R(nxcoord+i) = xcoord((i+1)*nxcoord)
polyg%R(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1)
polyg%z(nxcoord+i) = ycoord((i+1)*nxcoord)
polyg%z(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1)
end do
! c====================================================================================
do i=1,nxcoord
xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1)
ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1)
polyg%R(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1)
polyg%z(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1)
end do
! c====================================================================================
!
! check if (xcoord0, ycoord0) is out of table range
! incheck = 1 inside / 0 outside
if(inside(xpolyg, ypolyg, xcoord0, ycoord0)) incheck = 1
if (polyg%contains(xcoord0, ycoord0)) incheck = 1
end if
deallocate(wrk,iwrk)
!#######################################################################################
@ -582,50 +572,52 @@ contains
else
! c====================================================================================
if(incheck.eq.0) then
allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), &
ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), &
xpolygD(nycoord), ypolygD(nycoord), &
xoutA(nxcoord+3), youtA(nxcoord+3), xoutB(nycoord+3), &
youtB(nycoord+3), xoutC(nxcoord+3), youtC(nxcoord+3))
allocate(polygA%R(nxcoord), polygA%z(nxcoord))
allocate(polygC%R(nxcoord), polygC%z(nxcoord))
allocate(polygB%R(nycoord), polygB%z(nycoord))
allocate(polygD%R(nycoord), polygD%z(nycoord))
allocate(outA%R(nxcoord+3), outA%z(nxcoord+3))
allocate(outB%R(nycoord+3), outB%z(nycoord+3))
allocate(outC%R(nxcoord+3), outC%z(nxcoord+3))
! coordinates of vertices v1,v2,v3,v4
xvert(1) = xpolyg(1)
xvert(2) = xpolyg(nxcoord)
xvert(3) = xpolyg(nxcoord+nycoord-1)
xvert(4) = xpolyg(2*nxcoord+nycoord-2)
yvert(1) = ypolyg(1)
yvert(2) = ypolyg(nxcoord)
yvert(3) = ypolyg(nxcoord+nycoord-1)
yvert(4) = ypolyg(2*nxcoord+nycoord-2)
xvert(1) = polyg%R(1)
xvert(2) = polyg%R(nxcoord)
xvert(3) = polyg%R(nxcoord+nycoord-1)
xvert(4) = polyg%R(2*nxcoord+nycoord-2)
yvert(1) = polyg%z(1)
yvert(2) = polyg%z(nxcoord)
yvert(3) = polyg%z(nxcoord+nycoord-1)
yvert(4) = polyg%z(2*nxcoord+nycoord-2)
! coordinates of side A,B,C,D
xpolygA = xpolyg(1:nxcoord)
ypolygA = ypolyg(1:nxcoord)
xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1)
ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1)
xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg)
xpolygD(nycoord) = xpolyg(1)
ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg)
ypolygD(nycoord) = ypolyg(1)
polygA%R = polyg%R(1:nxcoord)
polygA%z = polyg%z(1:nxcoord)
polygB%R = polyg%R(nxcoord:nxcoord+nycoord-1)
polygB%z = polyg%z(nxcoord:nxcoord+nycoord-1)
polygC%R = polyg%R(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
polygC%z = polyg%z(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
polygD%R(1:nycoord-1) = polyg%R(2*nxcoord+nycoord-2:npolyg)
polygD%R(nycoord) = polyg%R(1)
polygD%z(1:nycoord-1) = polyg%z(2*nxcoord+nycoord-2:npolyg)
polygD%z(nycoord) = polyg%z(1)
! contour of outside regions A (1,2), B(3,4), C(5,6)
xoutA = huge(one)
xoutA(1:nxcoord) = xpolygA
xoutA(nxcoord+3) = xvert(1)
youtA = -huge(one)
youtA(1:nxcoord) = ypolygA
youtA(nxcoord+1) = yvert(2)
xoutB = huge(one)
xoutB(1:nycoord) = xpolygB
xoutB(nycoord+1) = xvert(3)
youtB = huge(one)
youtB(1:nycoord) = ypolygB
youtB(nycoord+3) = yvert(2)
xoutC = -huge(one)
xoutC(1:nxcoord) = xpolygC
xoutC(nxcoord+3) = xvert(3)
youtC = huge(one)
youtC(1:nxcoord) = ypolygC
youtC(nxcoord+1) = yvert(4)
outA%R = huge(one)
outA%R(1:nxcoord) = polygA%R
outA%R(nxcoord+3) = xvert(1)
outA%z = -huge(one)
outA%z(1:nxcoord) = polygA%z
outA%z(nxcoord+1) = yvert(2)
outB%R = huge(one)
outB%R(1:nycoord) = polygB%R
outB%R(nycoord+1) = xvert(3)
outB%z = huge(one)
outB%z(1:nycoord) = polygB%z
outB%z(nycoord+3) = yvert(2)
outC%R = -huge(one)
outC%R(1:nxcoord) = polygC%R
outC%R(nxcoord+3) = xvert(3)
outC%z = huge(one)
outC%z(1:nxcoord) = polygC%z
outC%z(nxcoord+1) = yvert(4)
! c----------------------------------------------------------------------------------
! search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid
!
@ -639,17 +631,17 @@ contains
! v1 A v2
! (8) | (1) | (2)
!
if(inside(xoutA, youtA, xcoord0, ycoord0)) then
if (outA%contains(xcoord0, ycoord0)) then
in = 1
if(xcoord0.GT.xvert(2)) then
in = 2
end if
else if(inside(xoutB, youtB, xcoord0, ycoord0)) then
else if (outB%contains(xcoord0, ycoord0)) then
in = 3
if(ycoord0.GT.yvert(3)) then
in = 4
end if
else if(inside(xoutC, youtC, xcoord0, ycoord0)) then
else if (outC%contains(xcoord0, ycoord0)) then
in = 5
if(xcoord0.LT.xvert(4)) then
in = 6
@ -678,9 +670,9 @@ contains
CASE (1)
! params%beta outside table range
! locate position of xcoord0 with respect to x coordinates of side A
call locate(xpolygA,nxcoord,xcoord0,ii)
call locate(polygA%R, xcoord0, ii)
! find corresponding y value on side A for xcoord position
call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0)
ycoord0 = linear_interp(polygA%R(ii), polygA%z(ii), polygA%R(ii+1), polygA%z(ii+1), xcoord0)
incheck = 1
CASE (2)
! params%alpha and params%beta outside table range
@ -690,8 +682,8 @@ contains
ii = nxcoord !indice per assegnare valori waist, rci, phi
CASE (3)
! params%alpha outside table range
call locate(ypolygB,nycoord,ycoord0,ii)
call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0)
call locate(polygB%z, ycoord0, ii)
xcoord0 = linear_interp(polygB%z(ii), polygB%R(ii), polygB%z(ii+1), polygB%R(ii+1), ycoord0)
incheck = 1
CASE (4)
! params%alpha and params%beta outside table range
@ -700,8 +692,8 @@ contains
ii = nxcoord+nycoord-1
CASE (5)
! params%beta outside table range
call locate(xpolygC,nxcoord,xcoord0,ii)
call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0)
call locate(polygC%R, xcoord0, ii)
ycoord0 = linear_interp(polygC%R(ii+1), polygC%z(ii+1), polygC%R(ii), polygC%z(ii), xcoord0)
incheck = 1
CASE (6)
! params%alpha and params%beta outside table range
@ -710,8 +702,8 @@ contains
ii = 2*nxcoord+nycoord-2
CASE (7)
! params%alpha outside table range
call locate(ypolygD,nycoord,ycoord0,ii)
call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0)
call locate(polygD%z, ycoord0, ii)
xcoord0 = linear_interp(polygD%z(ii), polygD%R(ii), polygD%z(ii+1), polygD%R(ii+1), ycoord0)
incheck = 1
CASE (8)
! params%alpha and params%beta outside table range
@ -721,8 +713,8 @@ contains
END SELECT
! c----------------------------------------------------------------------------------
!
deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, &
xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC)
deallocate(polygA%R, polygA%z, polygC%R, polygC%z, polygB%R, polygB%z, &
polygD%R, polygD%z, outA%R, outA%z, outB%R, outB%z, outC%R, outC%z)
end if
! c====================================================================================
!
@ -785,7 +777,7 @@ contains
crci1, txrci1, tyrci1, crci2, txrci2, tyrci2, &
cphi1, txphi1, typhi1, cphi2, txphi2, typhi2, &
cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0, &
wrk2, xpolyg, ypolyg, w)
wrk2, polyg%R, polyg%z, w)
end if
!
!#######################################################################################

6
src/gray_core.f90
View File

@ -23,7 +23,7 @@ contains
use magsurf_data, only : compute_flux_averages, dealloc_surfvec
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use utils, only : vmaxmin, inside
use utils, only : vmaxmin
use multipass, only : initbeam, initmultipass, turnoffray, &
plasma_in, plasma_out, wall_out
use logger, only : log_info, log_debug
@ -508,7 +508,7 @@ contains
end if
! test whether further trajectory integration is unnecessary
call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling
call vmaxmin(tau1+tau0+lgcpl1, taumn, taumx) ! test on tau + coupling
if(is_critical(error)) then ! stop propagation for current beam
istop_pass = istop_pass +1 ! * +1 non propagating beam
@ -531,7 +531,7 @@ contains
if(i>params%raytracing%nstep) i=params%raytracing%nstep
pabs_beam = sum(ppabs(:,i))
icd_beam = sum(ccci(:,i))
call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print
call vmaxmin(tau0, taumn, taumx) ! taumn,taumx for print
if (params%equilibrium%iequil /= EQ_VACUUM) then
! compute power and current density profiles for all rays

5
src/gray_equil.f90
View File

@ -683,7 +683,7 @@ contains
use gray_params, only : equilibrium_parameters
use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
use utils, only : vmaxmin, vmaxmini
use utils, only : vmaxmin
use logger, only : log_debug, log_info
! subroutine arguments
@ -844,7 +844,7 @@ contains
! for the search of O/X points
nbnd = size(data%boundary%R)
if (nbnd > 0) then
call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup)
call vmaxmin(data%boundary%z, zbinf, zbsup, ibinf, ibsup)
rbinf = data%boundary%R(ibinf)
rbsup = data%boundary%R(ibsup)
else
@ -1387,7 +1387,6 @@ contains
! For a description of the G-EQDSK format, see the GRAY user manual.
use const_and_precisions, only : one
use gray_params, only : equilibrium_parameters
use utils, only : get_free_unit
use logger, only : log_error
! subroutine arguments

5
src/gray_params.f90
View File

@ -447,7 +447,6 @@ contains
subroutine read_gray_params(filename, params, err)
! Reads the GRAY parameters from the legacy gray_params.data file
use utils, only : get_free_unit
use logger, only : log_error
! subrouting arguments
@ -458,9 +457,7 @@ contains
! local variables
integer :: u, temp(3)
u = get_free_unit()
open(u, file=filename, status='old', action='read', iostat=err)
open(newunit=u, file=filename, status='old', action='read', iostat=err)
if (err /= 0) then
call log_error('opening gray_params file ('//filename//') failed!', &
mod='gray_params', proc='read_gray_params')

4
src/ini_parser.f90
View File

@ -60,7 +60,6 @@ contains
! ERR_VALUE on invalid values for this property;
! ERR_UNKNOWN on unknown property.
!
use utils, only : get_free_unit
! function argument
character(*), intent(in) :: filepath
@ -74,8 +73,7 @@ contains
character(len=:), allocatable :: section, name, value
! open the INI file
ini = get_free_unit()
open(unit=ini, file=filepath, action='read', iostat=error)
open(newunit=ini, file=filepath, action='read', iostat=error)
if (error /= 0) then
write (msg, '("failed to open INI file: ", a)') filepath
call log_error(msg, proc='parse_ini', mod='ini_parser')

28
src/pec.f90
View File

@ -154,7 +154,7 @@ contains
subroutine pec_tab(xxi, ypt, yamp, ii, xtab1, wdpdv, wajphiv)
! Power and current projected on ψ grid - mid points
use utils, only : locatex, intlin
use utils, only : locate, linear_interp
! subroutine arguments
real(wp_), dimension(:), intent(in) :: xxi, ypt, yamp
@ -222,14 +222,15 @@ contains
indi = ind
if (idecr == -1) indi = ind - 1
rt1 = xtab1(indi)
call locatex(xxi,iise,iise0,iise,rt1,itb1)
if (itb1 >= iise0 .and. itb1 <= iise) then
call locate(xxi(iise0:iise), rt1, itb1)
if (itb1 >= 1 .and. itb1 <= iise) then
itb1 = itb1 + (iise0 - 1)
if(itb1 == iise) then
ppa2=ypt(itb1)
cci2=yamp(itb1)
else
call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2)
call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2)
ppa2 = linear_interp(xxi(itb1), ypt(itb1), xxi(itb1+1), ypt(itb1+1), rt1)
cci2 = linear_interp(xxi(itb1), yamp(itb1), xxi(itb1+1), yamp(itb1+1), rt1)
end if
dppa = ppa2 - ppa1
didst = cci2 - cci1
@ -321,7 +322,7 @@ contains
subroutine profwidth(xx,yy,xpk,ypk,dxxe)
use const_and_precisions, only : emn1
use utils, only : locatex, locate, intlin, vmaxmini
use utils, only : locate, linear_interp, vmaxmin
! subroutine arguments
real(wp_), intent(in) :: xx(:), yy(:)
@ -332,7 +333,7 @@ contains
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
real(wp_) :: ypkp,ypkm
call vmaxmini(yy,size(yy),ymn,ymx,imn,imx)
call vmaxmin(yy, ymn, ymx, imn, imx)
ypk = 0
xmx = xx(imx)
xmn = xx(imn)
@ -355,15 +356,16 @@ contains
xpk = xpkp
ypk = ypkp
yye = ypk*emn1
call locatex(yy,size(yy),1,ipk,yye,ie)
call locate(yy(1:ipk), yye, ie)
if(ie > 0 .and. ie < size(yy)) then
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1)
rte1 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
else
rte1 = 0
end if
call locatex(yy,size(yy),ipk,size(yy),yye,ie)
call locate(yy(ipk:size(yy)), yye, ie)
if(ie > 0 .and. ie < size(yy)) then
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
ie = ie + (ipk - 1)
rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
else
rte2 = 0
end if
@ -373,9 +375,9 @@ contains
ypk=yy(2)
rte1=0.0_wp_
yye=ypk*emn1
call locate(yy,size(yy),yye,ie)
call locate(yy, yye, ie)
if(ie > 0 .and. ie < size(yy)) then
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye)
else
rte2 = 0
end if

4
src/reflections.f90
View File

@ -86,7 +86,7 @@ end subroutine inters_linewall
subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
use utils, only : bubble
use utils, only : sort_pair
real(wp_), intent(in), dimension(3) :: xv,kv
real(wp_), intent(in), dimension(2) :: rs,zs
real(wp_), dimension(2), intent(out) :: s,t
@ -147,7 +147,7 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
else
s(1) = (-bhalf+sqrt(delta))/a
s(2) = (-bhalf-sqrt(delta))/a
call bubble(s,2)
call sort_pair(s)
t(:) = (kz*s(:)+(z0-zs(1)))/dz
n = 2
end if

6
src/splines.f90
View File

@ -121,7 +121,7 @@ contains
integer :: i
real(wp_) :: dx
call locate(self%data, self%ndata, x, i)
call locate(self%data, x, i)
i = min(max(1, i), self%ndata - 1)
dx = x - self%data(i)
y = self%raw_eval(i, dx)
@ -156,7 +156,7 @@ contains
integer :: i
real(wp_) :: dx
call locate(self%data, self%ndata, x, i)
call locate(self%data, x, i)
i = min(max(1, i), self%ndata - 1)
dx = x - self%data(i)
y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4))
@ -691,7 +691,7 @@ contains
! local variables
integer :: i
call locate(self%xdata, self%ndata, x, i)
call locate(self%xdata, x, i)
i = min(max(1, i), self%ndata - 1)
y = self%raw_eval(i, x)
end function linear_1d_eval

12
src/types.f90
View File

@ -190,7 +190,7 @@ contains
! Save the table to a file
!
! Note: this operation consumes the table
use utils, only : get_free_unit, digits
use utils, only : digits
! suborutine arguments
class(table), intent(inout) :: self
@ -206,8 +206,6 @@ contains
if (self%empty()) return
file = get_free_unit()
! set default file name
if (present(filepath)) then
filepath_ = filepath
@ -216,7 +214,7 @@ contains
end if
! open output file
open(unit=file, file=filepath_, action='write', iostat=error)
open(newunit=file, file=filepath_, action='write', iostat=error)
if (error /= 0) return
! write a custom file header
@ -299,7 +297,7 @@ contains
pure function contour_contains(self, R0, z0) result(inside)
! Tests whether the point (`R`, `z`) lies inside the 2D contour
use utils, only : intlinf, locate_unord
use utils, only : linear_interp, locate_unord
! subroutine arguments
class(contour), intent(in) :: self
@ -324,8 +322,8 @@ contains
! means that the point is outside the polygon.
do i = 1, nsegs
! coordinate of the intersection between segment and z=z0
R = intlinf(self%z(seg(i)), self%R(seg(i)), &
self%z(seg(i)+1), self%R(seg(i)+1), z0)
R = linear_interp(self%z(seg(i)), self%R(seg(i)), &
self%z(seg(i)+1), self%R(seg(i)+1), z0)
if (R < R0) inside = .not. inside
end do
end function contour_contains

358
src/utils.f90
View File

@ -6,83 +6,38 @@ module utils
contains
pure function locatef(a,n,x) result(j)
! Given an array a(n), and a value x, with a(n) monotonic, either
! increasing or decreasing, returns a value j such that
! a(j) < x <= a(j+1) for a increasing, and such that
! a(j+1) < x <= a(j) for a decreasing.
! j=0 or j=n indicate that x is out of range (Numerical Recipes)
integer, intent(in) :: n
real(wp_), dimension(n), intent(in) :: a
real(wp_), intent(in) :: x
integer :: j
integer :: jl,ju,jm
logical :: incr
jl=0
ju=n+1
incr=a(n)>a(1)
do while ((ju-jl)>1)
jm=(ju+jl)/2
if(incr.eqv.(x>a(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end function locatef
pure subroutine locate(array, value, index)
! Given an `array`, and a `value` returns the `index`
! such that `array(index) < value < array(index+1)`
!
! Notes:
! 1. `array` must be monotonic, either increasing or decreasing.
! 2. an index equal to 0 or `size(array)` indicate the value was not found
! Source: Numerical Recipes
pure subroutine locate(xx,n,x,j)
integer, intent(in) :: n
real(wp_), intent(in) :: xx(n), x
integer, intent(out) :: j
integer :: jl,ju,jm
! subroutine arguments
real(wp_), intent(in) :: array(:), value
integer, intent(out) :: index
! local variables
integer :: jl, ju, jm
logical :: incr
!
! Given an array xx(n), and a value x
! returns a value j such that xx(j) < x < xx(j+1)
! xx(n) must be monotonic, either increasing or decreasing.
! j=0 or j=n indicate that x is out of range (Numerical Recipes)
!
jl=0
ju=n+1
incr=xx(n)>xx(1)
do while ((ju-jl)>1)
jm=(ju+jl)/2
if(incr .eqv. (x>xx(jm))) then
jl=jm
jl = 0
ju = size(array) + 1
incr = array(size(array)) > array(1)
do while ((ju - jl) > 1)
jm = (ju + jl)/2
if (incr .eqv. (value > array(jm))) then
jl = jm
else
ju=jm
ju = jm
endif
end do
j=jl
index = jl
end subroutine locate
pure subroutine locatex(xx,n,n1,n2,x,j)
integer, intent(in) :: n,n1,n2
real(wp_), intent(in) :: xx(n), x
integer, intent(out) :: j
integer :: jl,ju,jm
!
! Given an array xx(n), and a value x
! returns a value j such that xx(j) < x < xx(j+1)
! xx(n) must be monotonic, either increasing or decreasing.
! j=n1-1or j=n2+1 indicate that x is out of range
! modified from subr. locate (Numerical Recipes)
!
jl=n1-1
ju=n2+1
do while ((ju-jl)>1)
jm=(ju+jl)/2
if((xx(n2)>xx(n1)) .eqv. (x>xx(jm))) then
jl=jm
else
ju=jm
endif
end do
j=jl
end subroutine locatex
pure subroutine locate_unord(array, value, locs, n, nlocs)
! Given an `array` of size `n` and a `value`, finds at most
@ -117,233 +72,74 @@ contains
end subroutine locate_unord
pure function intlinf(x1,y1,x2,y2,x) result(y)
!linear interpolation
!must be x1 != x2
use const_and_precisions, only : one
real(wp_),intent(in) :: x1,y1,x2,y2,x
pure function linear_interp(x1, y1, x2, y2, x) result(y)
! Computes the linear interpolation between
! the points `(x1, y1)` and `(x2, y2)` at `x`.
!
! Note: x1 != x2 is assumed
! function arguments
real(wp_), intent(in) :: x1, y1, x2, y2, x
real(wp_) :: y
! local variables
real(wp_) :: a
a=(x2-x)/(x2-x1)
y=a*y1+(one-a)*y2
end function intlinf
pure subroutine intlin(x1,y1,x2,y2,x,y)
real(wp_), intent(in) :: x1,y1,x2,y2,x
real(wp_), intent(out) :: y
real(wp_) :: dx,aa,bb
!
! linear interpolation
! (x1,y1) < (x,y) < (x2,y2)
!
dx=x2-x1
aa=(x2-x)/dx
bb=1.0_wp_-aa
y=aa*y1+bb*y2
end subroutine intlin
a = (x2 - x)/(x2 - x1)
y = a*y1 + (1 - a)*y2
end function linear_interp
pure subroutine vmax(x,n,xmax,imx)
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmax
integer, intent(out) :: imx
pure subroutine vmaxmin(x, xmin, xmax, imin, imax)
! Computes the maximum and minimum of the array `x`
! and optionally their indices
! subroutine arguments
real(wp_), intent(in) :: x(:)
real(wp_), intent(out) :: xmin, xmax
integer, intent(out), optional :: imin, imax
! local variables
integer :: i
if (n<1) then
imx=0
if (size(x) < 1) then
if (present(imin)) imin = 0
if (present(imax)) imax = 0
return
end if
imx=1
xmax=x(1)
do i=2,n
if(x(i)>xmax) then
xmax=x(i)
imx=i
end if
end do
end subroutine vmax
pure subroutine vmin(x,n,xmin,imn)
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin
integer, intent(out) :: imn
integer :: i
if (present(imin)) imin = 1
if (present(imax)) imax = 1
xmin = x(1)
xmax = x(1)
if (n<1) then
imn=0
return
end if
imn=1
xmin=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
imn=i
end if
end do
end subroutine vmin
pure subroutine vmaxmini(x,n,xmin,xmax,imn,imx)
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin, xmax
integer, intent(out) :: imn, imx
integer :: i
if (n<1) then
imn=0
imx=0
return
end if
imn=1
imx=1
xmin=x(1)
xmax=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
imn=i
else if(x(i)>xmax) then
xmax=x(i)
imx=i
end if
end do
end subroutine vmaxmini
pure subroutine vmaxmin(x,n,xmin,xmax)
integer, intent(in) :: n
real(wp_), intent(in) :: x(n)
real(wp_), intent(out) :: xmin, xmax
integer :: i
if (n<1) then
return
end if
xmin=x(1)
xmax=x(1)
do i=2,n
if(x(i)<xmin) then
xmin=x(i)
else if(x(i)>xmax) then
xmax=x(i)
do i = 2, size(x)
if(x(i) < xmin) then
xmin = x(i)
if (present(imin)) imin = i
else if(x(i) > xmax) then
xmax = x(i)
if (present(imax)) imax = i
end if
end do
end subroutine vmaxmin
pure subroutine order(p,q)
! returns p,q in ascending order
real(wp_), intent(inout) :: p,q
pure subroutine sort_pair(p)
! Sorts the pair `p` in ascending order
! subroutine inputs
real(wp_), intent(inout) :: p(2)
! local variables
real(wp_) :: temp
if (p>q) then
temp=p
p=q
q=temp
if (p(1) > p(2)) then
temp = p(1)
p(1) = p(2)
p(2) = temp
end if
end subroutine order
pure subroutine bubble(a,n)
! bubble sorting of array a
integer, intent(in) :: n
real(wp_), dimension(n), intent(inout) :: a
integer :: i, j
do i=1,n
do j=n,i+1,-1
call order(a(j-1), a(j))
end do
end do
end subroutine bubble
pure subroutine range2rect(xmin, xmax, ymin, ymax, x, y)
! Given two ranges [xmin, xmax], [ymin, ymax] builds
! the x and y vertices of the following rectangle:
!
! (xmin, ymax)╔═════╗(xmax, ymax)
! ║4 3║
! ║ ║
! ║1 2║
! (xmin, ymin)╚═════╝(xmin, ymax)
!
! subroutine arguments
real(wp_), intent(in) :: xmin, xmax, ymin, ymax
real(wp_), intent(out), dimension(5) :: x, y
x = [xmin, xmax, xmax, xmin, xmin]
y = [ymin, ymin, ymax, ymax, ymin]
end subroutine range2rect
function inside(vertx, verty, x0, y0)
! Tests whether the point (`x0`, `y0`) lies inside the
! simple polygon of vertices `vertx`, `verty`.
! subroutine arguments
real(wp_), dimension(:), intent(in) :: vertx, verty
real(wp_), intent(in) :: x0, y0
logical :: inside
! local variables
integer :: seg(size(vertx))
real(wp_) :: x, vertx_(size(vertx)+1), verty_(size(vertx)+1)
integer :: i, nsegs, n
! Ensure the first and last point are the same
n = size(vertx)
vertx_(1:n) = vertx(1:n)
verty_(1:n) = verty(1:n)
vertx_(n+1) = vertx(1)
verty_(n+1) = verty(1)
inside = .false.
! Find the `nsegs` segments that intersect the horizontal
! line y=y0, i.e. `verty(seg(i)) < y < verty(seg(i)+1)`
call locate_unord(verty_, y0, seg, n, nsegs)
! No intersections, it must be outside (above or below)
if (nsegs == 0) return
! Count the number of intersections that lie to the left
! (equivalently, to the right) of the point. An even number
! means that the point is outside the polygon.
do i = 1, nsegs
! coordinate of the intersection between segment and y=y0
x = intlinf(verty_(seg(i)), vertx_(seg(i)), &
verty_(seg(i)+1), vertx_(seg(i)+1), y0)
if (x < x0) inside = .not. inside
end do
end function inside
function get_free_unit(unit) result(i)
! Returns `unit` back or the first free unit
! number `i` if `unit` is absent.
! When no unit is available, returns -1.
! function arguments
integer :: i
integer, intent(in), optional :: unit
! local variables
integer, parameter :: max_allowed = 999
integer :: error
logical :: ex, op
if (present(unit)) then
i = unit
return
end if
do i=0,max_allowed
inquire(unit=i, exist=ex, opened=op, iostat=error)
! if unit i exists and is free
if (error == 0 .and. ex .and. .not. op) return
end do
i = -1
end function get_free_unit
end subroutine sort_pair
function dirname(filepath) result(directory)