diff --git a/src/beams.f90 b/src/beams.f90 index 2bfcd57..2d6a4ce 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -256,8 +256,7 @@ contains ! ellipses in the transverse plane at the launch point (deg) use gray_params, only : antenna_parameters - use utils, only : get_free_unit, intlin, locate - use reflections, only : inside + use utils, only : get_free_unit, intlin, locate, inside use dierckx, only : curfit, splev, surfit, bispev use logger, only : log_error @@ -557,7 +556,7 @@ contains ! ! check if (xcoord0, ycoord0) is out of table range ! incheck = 1 inside / 0 outside - if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) incheck = 1 + if(inside(xpolyg, ypolyg, xcoord0, ycoord0)) incheck = 1 end if deallocate(wrk,iwrk) !####################################################################################### @@ -666,17 +665,17 @@ contains ! v1 A v2 ! (8) | (1) | (2) ! - if(inside(xoutA,youtA,nxcoord+3,xcoord0,ycoord0)) then + if(inside(xoutA, youtA, xcoord0, ycoord0)) then in = 1 if(xcoord0.GT.xvert(2)) then in = 2 end if - else if(inside(xoutB,youtB,nycoord+3,xcoord0,ycoord0)) then + else if(inside(xoutB, youtB, xcoord0, ycoord0)) then in = 3 if(ycoord0.GT.yvert(3)) then in = 4 end if - else if(inside(xoutC,youtC,nxcoord+3,xcoord0,ycoord0)) then + else if(inside(xoutC, youtC, xcoord0, ycoord0)) then in = 5 if(xcoord0.LT.xvert(4)) then in = 6 diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index a97d574..c623829 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -408,8 +408,7 @@ contains use const_and_precisions, only : zero, one use gray_params, only : equilibrium_parameters, equilibrium_data use gray_params, only : iequil - use reflections, only : inside - use utils, only : vmaxmin, vmaxmini + use utils, only : vmaxmin, vmaxmini, inside use logger, only : log_info implicit none @@ -463,7 +462,7 @@ contains do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle + if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle else if(data%psin(i,j).le.0.0d0) cycle end if @@ -477,7 +476,7 @@ contains do j=1,nz do i=1,nr if (nbnd.gt.0) then - if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle + if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle else if(data%psin(i,j).le.0.0d0) cycle end if @@ -910,7 +909,7 @@ contains ! ! Note: all output arguments are optional. use gray_params, only : iequil - use reflections, only : inside + use utils, only : inside use const_and_precisions, only : pi implicit none @@ -1037,7 +1036,7 @@ contains if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz else ! Numerical data - if (inside(psi_domain%R, psi_domain%z, psi_domain%npoints, R, z)) then + if (inside(psi_domain%R, psi_domain%z, R, z)) then ! Within the interpolation range if (present(psi_n)) psi_n = psi_spline%eval(R, z) if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index c1edae8..310a73d 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -20,8 +20,7 @@ contains use beamdata, only : init_btr, dealloc_beam use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab - use utils, only : vmaxmin - use reflections, only : inside + use utils, only : vmaxmin, inside use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & initmultipass, turnoffray, plasma_in, plasma_out, & wall_out @@ -232,7 +231,7 @@ contains rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ if(inside(data%equilibrium%rlim, data%equilibrium%zlim, & - nlim, rrm, zzm)) then ! * start propagation in/outside vessel? + rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside else iow(jk) = 0 ! + outside @@ -295,7 +294,7 @@ contains ins_pl = (psinv>=zero .and. psinvset_globals use const_and_precisions, only : cm @@ -73,7 +73,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Convert to a list of R,z call range2rect(xmin=params%misc%rwall, xmax=R_max, ymin=z(1), ymax=z(mz), & - xv=data%equilibrium%rlim, yv=data%equilibrium%zlim) + x=data%equilibrium%rlim, y=data%equilibrium%zlim) ! Set the global variables of `limiter` call limiter_set_globals(data%equilibrium) diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index a841d20..b78c339 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -206,6 +206,7 @@ contains rhopjp=height psinjp=height*height rhotjp=frhotor(rhopjp) + if (rhotjp /= rhotjp) print *, 'Nan!!!!!!!!!!' psicon(jp)=height call contours_psi(psinjp,rctemp,zctemp,rup,zup,rlw,zlw) diff --git a/src/main.f90 b/src/main.f90 index 95c3de8..8ffb442 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -385,7 +385,7 @@ contains subroutine init_misc(params, data) ! Performs miscellanous initial tasks, before the gray_main subroutine. - use reflections, only : range2rect + use utils, only : range2rect use limiter, only : limiter_set_globals=>set_globals use const_and_precisions, only : cm @@ -440,7 +440,7 @@ contains ! call range2rect(xmin=params%misc%rwall, xmax=R_max, & ymin=z_launch - delta_z, ymax=z_launch + delta_z, & - xv=data%equilibrium%rlim, yv=data%equilibrium%zlim) + x=data%equilibrium%rlim, y=data%equilibrium%zlim) end block end if diff --git a/src/reflections.f90 b/src/reflections.f90 index 7b40f39..8f9c9f4 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -3,9 +3,9 @@ module reflections implicit none private - public :: reflect,inters_linewall,inside - public :: linecone_coord,interssegm_coord,interssegm - public :: wall_refl,range2rect + public :: reflect, inters_linewall + public :: linecone_coord, interssegm_coord, interssegm + public :: wall_refl contains @@ -191,47 +191,9 @@ function interssegm(xa,ya,xb,yb) end function interssegm - - subroutine range2rect(xmin,xmax,ymin,ymax,xv,yv) - implicit none - real(wp_), intent(in) :: xmin,xmax,ymin,ymax - real(wp_), intent(out), dimension(5) :: xv,yv - xv=(/xmin,xmax,xmax,xmin,xmin/) - yv=(/ymin,ymin,ymax,ymax,ymin/) - end subroutine range2rect - - - -function inside(xc,yc,n,x,y) - use utils, only : locatef, locate_unord, intlinf, bubble - implicit none - integer, intent(in) :: n - real(wp_), dimension(n), intent(in) :: xc,yc - real(wp_), intent(in) :: x,y - logical :: inside - integer, dimension(n) :: jint - real(wp_), dimension(n) :: xint - real(wp_), dimension(n+1) :: xclosed,yclosed - integer :: i,nj - xclosed(1:n)=xc(1:n) - yclosed(1:n)=yc(1:n) - xclosed(n+1)=xc(1) - yclosed(n+1)=yc(1) - call locate_unord(yclosed,n+1,y,jint,n,nj) - inside=.false. - if (nj==0) return - do i=1,nj - xint(i)=intlinf(yclosed(jint(i)),xclosed(jint(i)), & - yclosed(jint(i)+1),xclosed(jint(i)+1),y) - end do - call bubble(xint,nj) - inside=(mod(locatef(xint,nj,x),2)==1) -end function inside - - - subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) use limiter, only : rlim,zlim,nlim + use utils, only : inside implicit none ! arguments integer :: irfl @@ -253,7 +215,7 @@ subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) smax=smax*1.0e2_wp_ xvrfl=xv+smax*anv0 irfl=1 - if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + if (.not.inside(rlim,zlim,rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln)a(i).neqv.x>a(i+1)) then - nj=nj+1 - if (nj<=m) j(nj)=i + logical :: larger_than_last + + nlocs = 0 + if (size(array) < 2) return + + larger_than_last = value > array(1) + do i = 2, size(array) + ! Note: the condition is equivalent to + ! (array(i-1) < value < array(i)) + ! .or. (array(i) < value < array(i-1)) + if (array(i) < value .neqv. larger_than_last) then + larger_than_last = value > array(i) + nlocs = nlocs + 1 + if (nlocs <= n) locs(nlocs) = i - 1 end if end do end subroutine locate_unord + function intlinf(x1,y1,x2,y2,x) result(y) !linear interpolation !must be x1 != x2 @@ -247,6 +266,72 @@ contains end subroutine bubble + 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) + ! + + implicit none + + ! 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`. + + implicit none + + ! 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. @@ -299,6 +384,7 @@ contains end if end function dirname + function isrelative(filepath) ! Check if `filepath` is a relative or an absolute path