added missing file
This commit is contained in:
parent
8c9273d684
commit
9fd9bf5f9c
287
src/reflections.f90
Normal file
287
src/reflections.f90
Normal file
@ -0,0 +1,287 @@
|
||||
module reflections
|
||||
implicit none
|
||||
private
|
||||
integer, parameter :: r8=selected_real_kind(15,300)
|
||||
real(r8), parameter :: tinyr8=tiny(1._r8)
|
||||
public :: reflect,inters_linewall,inside
|
||||
public :: linecone_coord,interssegm_coord,interssegm
|
||||
contains
|
||||
|
||||
subroutine reflect(ki,nsurf,ko)
|
||||
implicit none
|
||||
real(r8), intent(in), dimension(3) :: ki
|
||||
real(r8), intent(in), dimension(3) :: nsurf
|
||||
real(r8), intent(out), dimension(3) :: ko
|
||||
real(r8) :: twokn,norm2
|
||||
norm2 = dot_product(nsurf,nsurf)
|
||||
if (norm2>0.0_r8) then
|
||||
twokn = 2.0_r8*dot_product(ki,nsurf)/norm2
|
||||
ko=ki-twokn*nsurf
|
||||
else
|
||||
ko=ki
|
||||
end if
|
||||
end subroutine reflect
|
||||
|
||||
subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
|
||||
implicit none
|
||||
real(r8), intent(in), dimension(3) :: xv,kv
|
||||
integer, intent(in) :: nw
|
||||
real(r8), dimension(nw), intent(in) :: rw,zw
|
||||
real(r8), intent(out) :: sint
|
||||
real(r8), dimension(3), intent(out) :: normw
|
||||
integer :: i,j,ni,iint
|
||||
real(r8), dimension(2) :: si,ti
|
||||
real(r8) :: drw,dzw,xint,yint,rint,l,kxy
|
||||
real(r8) :: tol=sqrt(epsilon(1.0_r8))
|
||||
sint=huge(sint)
|
||||
iint=0
|
||||
normw=0.0_r8
|
||||
do i=1,nw-1
|
||||
!search intersections with i-th wall segment
|
||||
call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni)
|
||||
do while (ni>0 .and. si(1)<=tol)
|
||||
!remove solutions with s<=0
|
||||
ni = ni-1
|
||||
si(1) = si(2)
|
||||
ti(1) = ti(2)
|
||||
end do
|
||||
do j=1,ni
|
||||
if ((si(j)<sint .or. iint==0) .and. ti(j)>=0._r8 .and. ti(j)<=1._r8) then
|
||||
!check intersection is in r,z range and keep the closest
|
||||
sint = si(j)
|
||||
iint = i
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
if (iint==0) return
|
||||
!calculate wall normal at intersection point
|
||||
drw = rw(iint+1)-rw(iint)
|
||||
dzw = zw(iint+1)-zw(iint)
|
||||
xint = xv(1) + sint*kv(1)
|
||||
yint = xv(2) + sint*kv(2)
|
||||
rint = sqrt(xint**2+yint**2)
|
||||
l = sqrt(drw**2+dzw**2)
|
||||
kxy = sqrt(kv(1)**2+kv(2)**2)
|
||||
normw(3) = -drw/l
|
||||
if (rint>0.0_r8) then
|
||||
normw(1) = xint/rint*dzw/l
|
||||
normw(2) = yint/rint*dzw/l
|
||||
else
|
||||
normw(1) = kv(1)/kxy*dzw/l
|
||||
normw(2) = kv(2)/kxy*dzw/l
|
||||
end if
|
||||
!reverse normal if k.n>0
|
||||
if (dot_product(normw,kv)>0.0_r8) normw=-normw
|
||||
end subroutine inters_linewall
|
||||
|
||||
subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
|
||||
implicit none
|
||||
real(r8), intent(in), dimension(3) :: xv,kv
|
||||
real(r8), intent(in), dimension(2) :: rs,zs
|
||||
real(r8), dimension(2), intent(out) :: s,t
|
||||
integer, intent(out) :: n
|
||||
real(r8) :: x0,y0,z0,kx,ky,kz
|
||||
real(r8) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin
|
||||
x0=xv(1)
|
||||
y0=xv(2)
|
||||
z0=xv(3)
|
||||
kx=kv(1)
|
||||
ky=kv(2)
|
||||
kz=kv(3)
|
||||
dr = rs(2)-rs(1)
|
||||
dz = zs(2)-zs(1)
|
||||
s = 0
|
||||
t = 0
|
||||
if (abs(dz)<tinyr8) then
|
||||
!surface in horizontal plane
|
||||
if (abs(kz)<tinyr8 .or. abs(dr)<tinyr8) then
|
||||
n = 0
|
||||
else
|
||||
s(1) = (zs(1)-z0)/kz
|
||||
r = sqrt((x0+s(1)*kx)**2+(y0+s(1)*ky)**2)
|
||||
t(1) = (r-rs(1))/dr
|
||||
n = 1
|
||||
end if
|
||||
else
|
||||
a = (kx**2+ky**2) - (dr/dz*kz)**2
|
||||
bhalf = -dr/dz*kz*rs(1) + (kx*x0 + ky*y0) - (dr/dz)**2*kz*(z0-zs(1))
|
||||
c = (x0**2+y0**2) - (rs(1) + dr/dz*(z0-zs(1)))**2
|
||||
if (abs(a)<tinyr8) then
|
||||
!line parallel to cone generator
|
||||
if (abs(dr)<tinyr8) then
|
||||
!cylinder and vertical line
|
||||
n = 0
|
||||
else
|
||||
tvertex = -rs(1)/dr
|
||||
zvertex = zs(1) + tvertex*dz
|
||||
srmin = -(kx*x0 + ky*y0)/(kx**2+ky**2)
|
||||
rmin = sqrt((x0+srmin*kx)**2+(y0+srmin*ky)**2)
|
||||
zrmin = z0 + srmin*kz
|
||||
if (rmin<tinyr8 .and. abs(zrmin-zvertex)<tinyr8) then
|
||||
!line passing by cone vertex
|
||||
!s(1) = srmin
|
||||
!t(1) = tvertex
|
||||
!n = 1
|
||||
n = 0
|
||||
else
|
||||
s(1) = -0.5_r8*c/bhalf
|
||||
t(1) = (kz*s(1)+(z0-zs(1)))/dz
|
||||
n = 1
|
||||
end if
|
||||
end if
|
||||
else
|
||||
delta = bhalf**2 - a*c
|
||||
if (delta<0) then
|
||||
n = 0
|
||||
else
|
||||
s(1) = (-bhalf+sqrt(delta))/a
|
||||
s(2) = (-bhalf-sqrt(delta))/a
|
||||
call bubble(s,2)
|
||||
t(:) = (kz*s(:)+(z0-zs(1)))/dz
|
||||
n = 2
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end subroutine linecone_coord
|
||||
|
||||
subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
|
||||
implicit none
|
||||
real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
|
||||
real(r8), intent(out) :: s,t
|
||||
integer, intent(out) :: ierr
|
||||
real(r8) :: crossprod,dxa,dya,dxb,dyb
|
||||
dxa = xa(2)-xa(1)
|
||||
dya = ya(2)-ya(1)
|
||||
dxb = xb(2)-xb(1)
|
||||
dyb = yb(2)-yb(1)
|
||||
crossprod = dxb*dya - dxa*dyb
|
||||
if (abs(crossprod)<tiny(crossprod)) then
|
||||
s = 0.0_r8
|
||||
t = 0.0_r8
|
||||
ierr = 1
|
||||
else
|
||||
s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod
|
||||
t = (dya*(xa(1)-xb(1)) - dxa*(ya(1)-yb(1)))/crossprod
|
||||
ierr = 0
|
||||
end if
|
||||
end subroutine interssegm_coord
|
||||
|
||||
function interssegm(xa,ya,xb,yb)
|
||||
implicit none
|
||||
real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
|
||||
logical :: interssegm
|
||||
real(r8) :: s,t
|
||||
integer :: ierr
|
||||
interssegm = .false.
|
||||
call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
|
||||
if (ierr==0 .and. s>=0._r8 .and. s<=1._r8 .and. &
|
||||
t>=0._r8 .and. t<=1._r8) interssegm = .true.
|
||||
end function interssegm
|
||||
|
||||
function inside(xc,yc,n,x,y)
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
real(r8), dimension(n), intent(in) :: xc,yc
|
||||
real(r8), intent(in) :: x,y
|
||||
logical :: inside
|
||||
integer, dimension(n) :: jint
|
||||
real(r8), dimension(n) :: xint
|
||||
real(r8), 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)=intlin(yclosed(jint(i)),xclosed(jint(i)), &
|
||||
yclosed(jint(i)+1),xclosed(jint(i)+1),y)
|
||||
end do
|
||||
call bubble(xint,nj)
|
||||
inside=(mod(locate(xint,nj,x),2)==1)
|
||||
end function inside
|
||||
|
||||
function intlin(x1,y1,x2,y2,x) result(y)
|
||||
!linear interpolation
|
||||
!must be x1 != x2
|
||||
implicit none
|
||||
real(r8),intent(in) :: x1,y1,x2,y2,x
|
||||
real(r8) :: y
|
||||
real(r8) :: a
|
||||
a=(x2-x)/(x2-x1)
|
||||
y=a*y1+(1._r8-a)*y2
|
||||
end function intlin
|
||||
|
||||
subroutine locate_unord(a,n,x,j,m,nj)
|
||||
implicit none
|
||||
integer, intent(in) :: n,m
|
||||
integer, intent(out) :: nj
|
||||
real(r8), dimension(n), intent(in) :: a
|
||||
real(r8), intent(in) :: x
|
||||
integer, dimension(m), intent(inout) :: j
|
||||
integer :: i
|
||||
nj=0
|
||||
do i=1,n-1
|
||||
if (x>a(i).neqv.x>a(i+1)) then
|
||||
nj=nj+1
|
||||
if (nj<=m) j(nj)=i
|
||||
end if
|
||||
end do
|
||||
end subroutine locate_unord
|
||||
|
||||
function locate(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)
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
real(r8), dimension(n), intent(in) :: a
|
||||
real(r8), 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 locate
|
||||
|
||||
subroutine order(p,q)
|
||||
!returns p,q in ascending order
|
||||
implicit none
|
||||
real(r8), intent(inout) :: p,q
|
||||
real(r8) :: temp
|
||||
if (p>q) then
|
||||
temp=p
|
||||
p=q
|
||||
q=temp
|
||||
end if
|
||||
end subroutine order
|
||||
|
||||
subroutine bubble(a,n)
|
||||
!bubble sorting of array a
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
real(r8), 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
|
||||
|
||||
end module reflections
|
||||
|
||||
Loading…
Reference in New Issue
Block a user