module Sph_utils implicit none private public :: InitSphUtils public :: ArcsOvrlap, Pnt_in_Arc public :: Distance, PntsEqual, ArcsEqual, cartsEqual public :: ArcMidpnt, cross_prod, geod2cart, cart2geod public :: MatSlv, GCXsect, CommonGC, getS, getXYZ public :: sTriArea, sQuadArea, Centroid, Orientation, Pnt_in_triangle public :: Pnt_in_SphTri, Pnt_in_SphPoly public :: cartPnt, geodeticPnt, simpleArc, arc, spolygon public :: diagArc, diagPnts public :: ROUNDOFF public :: ZERO, ONEHALF, ONE, TWO, FOUR, FIVE, TEN public :: NINETY, ONE80, THREE60, NOTSET, PI, D2R, R2D real(8), parameter :: ZERO = 0._8 real(8), parameter :: ONEHALF = .5_8 real(8), parameter :: ONE = 1._8 real(8), parameter :: TWO = 2._8 real(8), parameter :: FOUR = 4._8 real(8), parameter :: FIVE = 5._8 real(8), parameter :: TEN = 10._8 real(8), parameter :: NINETY = 90._8 real(8), parameter :: ONE80 = 180._8 real(8), parameter :: THREE60 = 360._8 real(8), parameter :: NOTSET = -1000._8 ! real(8), parameter :: PI = 3.14159265358979323846 ! real(8), parameter :: D2R = PI/ONE80 ! real(8), parameter :: R2D = ONE80/PI real(8) :: ROUNDOFF = 1.e-10_8 real(8) :: SMALL = epsilon(SMALL) real(8) :: PI real(8) :: D2R real(8) :: R2D !--------------------------------------------------- ! define the types !--------------------------------------------------- type cartPnt real(8) :: x, y, z end type cartPnt type geodeticPnt real(8) :: lon real(8) :: lat end type geodeticPnt type simpleArc type(geodeticPnt) :: sPnt type(geodeticPnt) :: ePnt end type simpleArc type arc type(geodeticPnt) :: sPnt type(cartPnt) :: sxyz type(geodeticPnt) :: ePnt type(cartPnt) :: exyz real(8) :: Angle = NOTSET real(8) :: Normal(3) = (/ ZERO, ZERO, ZERO /) logical :: RHS = .true. logical :: IsLatArc = .false. end type arc type spolygon type(arc), allocatable :: arcs(:) end type spolygon !--------------------------------------------------- ! operator functions !--------------------------------------------------- interface operator (==) module procedure ArcsEqual end interface !--------------------------------------------------- ! overloaded functions !--------------------------------------------------- interface geod2cart module procedure ftngeod2cart, ftntypgeod2cart end interface interface cart2geod module procedure ftncart2geod end interface interface cartNorm module procedure ftncartNorm, typcartNorm end interface interface cross_prod module procedure cross_prod_xyz, cross_prod_typ end interface interface getS module procedure getSxyz, getSpnt end interface CONTAINS subroutine InitSphUtils !-------------------------------------------------------- ! initialize spherical utilities module !-------------------------------------------------------- PI = FOUR*atan(ONE) D2R = PI/ONE80 R2D = ONE80/PI end subroutine InitSphUtils function ArcMidpnt( Arc1 ) result(midpnt) !-------------------------------------------------------- ! get arc midpoint (lon,lat) !-------------------------------------------------------- type(arc), intent(inout) :: Arc1 type(geodeticPnt) :: midpnt real(8) :: dist, dist2 real(8) :: wght real(8) :: x, y, z real(8) :: Pnt(2) dist = Distance( Arc1%sPnt,Arc1%ePnt ) dist2 = ONEHALF*dist Arc1%sxyz = geod2cart( Arc1%sPnt ) Arc1%exyz = geod2cart( Arc1%ePnt ) wght = sin(dist2)/sin(dist) x = wght*(Arc1%sxyz%x + Arc1%exyz%x) y = wght*(Arc1%sxyz%y + Arc1%exyz%y) z = wght*(Arc1%sxyz%z + Arc1%exyz%z) Pnt(:) = cart2geod( x, y, z ) midpnt%lon = R2D*Pnt(1) midpnt%lat = R2D*Pnt(2) end function ArcMidpnt function ArcsOvrlap( Arc1, Arc2 ) result(Ovrlap) !-------------------------------------------------------- ! Do arcs overlap? !-------------------------------------------------------- !-------------------------------------------------------- ! dummy arguments !-------------------------------------------------------- type(arc), intent(in) :: Arc1, Arc2 logical :: Ovrlap Ovrlap = Arc1 == Arc2 if( .not. Ovrlap ) then Ovrlap = Pnt_in_Arc( Arc1,Arc2%sPnt ) if( .not. Ovrlap ) then Ovrlap = Pnt_in_Arc( Arc1,Arc2%ePnt ) if( .not. Ovrlap ) then Ovrlap = Pnt_in_Arc( Arc2,Arc1%sPnt ) if( .not. Ovrlap ) then Ovrlap = Pnt_in_Arc( Arc2,Arc1%ePnt ) endif endif endif endif end function ArcsOvrlap function CommonGC( Arc1, Arc2 ) result(sameGC) !-------------------------------------------------------- ! Are arcs in same GC? !-------------------------------------------------------- !-------------------------------------------------------- ! dummy arguments !-------------------------------------------------------- type(arc), intent(in) :: Arc1, Arc2 logical :: sameGC !-------------------------------------------------------- ! local variables !-------------------------------------------------------- integer :: n real(8) :: XYZ(3,2) real(8) :: Norm(3,2) XYZ(:,1) = geod2cart( Arc1%sPnt%lon,Arc1%sPnt%lat ) XYZ(:,2) = geod2cart( Arc1%ePnt%lon,Arc1%ePnt%lat ) Norm(:,1) = cross_prod( XYZ(:,1),XYZ(:,2) ) Norm(:,1) = cartNorm( Norm(:,1) ) XYZ(:,1) = geod2cart( Arc2%sPnt%lon,Arc2%sPnt%lat ) XYZ(:,2) = geod2cart( Arc2%ePnt%lon,Arc2%ePnt%lat ) Norm(:,2) = cross_prod( XYZ(:,1),XYZ(:,2) ) Norm(:,2) = cartNorm( Norm(:,2) ) do n = 1,3 sameGC = PntsEqual( norm(n,1),norm(n,2) ) if( .not. sameGC ) then exit endif enddo end function CommonGC function Pnt_in_Arc( Arc_in, Pnt ) result(InArc) !-------------------------------------------------------- ! Is point in arc? !-------------------------------------------------------- !-------------------------------------------------------- ! dummy arguments !-------------------------------------------------------- type(arc), intent(in) :: Arc_in type(geodeticPnt), intent(in) :: Pnt logical :: InArc !-------------------------------------------------------- ! local variables !-------------------------------------------------------- type(cartPnt ) :: tstcart, chkcart type(arc) :: theArc real(8) :: tstDist real(8) :: arcDist real(8) :: factor1, factor2 real(8) :: Cart(3,3) real(8) :: X(3,2) logical :: onGC !-------------------------------------------------------- ! Does Pnt match either beg or end pnt of arc? !-------------------------------------------------------- theArc = Arc_in InArc = PntsEqual( theArc%sPnt%lon,Pnt%lon ) .and. & PntsEqual( theArc%sPnt%lat,Pnt%lat ) if( .not. InArc ) then InArc = PntsEqual( theArc%ePnt%lon,Pnt%lon ) .and. & PntsEqual( theArc%ePnt%lat,Pnt%lat ) endif fullTest: & if( .not. InArc ) then arcDist = Distance( theArc%sPnt,theArc%ePnt ) tstDist = Distance( theArc%sPnt,Pnt ) InArc = .not. (tstDist > arcDist) if( InArc ) then theArc%sxyz = geod2cart( theArc%sPnt ) theArc%exyz = geod2cart( theArc%ePnt ) tstcart = geod2cart( Pnt ) factor1 = sin(arcDist-tstDist)/sin(arcDist) factor2 = sin(tstDist)/sin(arcDist) chkcart%x = factor1*theArc%sxyz%x + factor2*theArc%exyz%x chkcart%y = factor1*theArc%sxyz%y + factor2*theArc%exyz%y chkcart%z = factor1*theArc%sxyz%z + factor2*theArc%exyz%z InArc = .false. if( PntsEQual( chkcart%x,tstcart%x ) ) then if( PntsEQual( chkcart%y,tstcart%y ) ) then if( PntsEQual( chkcart%z,tstcart%z ) ) then InArc = .true. endif endif endif endif endif fullTest end function Pnt_in_Arc function Distance( pnt1, pnt2 ) result(Angle) !-------------------------------------------------------- ! Get distance between to geodetic pnts on unit sphere !-------------------------------------------------------- implicit none type(geodeticPnt), intent(in) :: pnt1 type(geodeticPnt), intent(in) :: pnt2 real(8) :: Angle real(8) :: dprod real(8) :: X(3,2) X(:,1) = geod2cart( pnt1%lon,pnt1%lat ) X(:,2) = geod2cart( pnt2%lon,pnt2%lat ) dprod = dot_product(X(:,1),X(:,2)) Angle = acos( max( min(ONE,dprod),-ONE ) ) if( abs(Angle) < SMALL ) then write(*,'(''Distance '',1pg15.7,'' is too small'')') Angle Stop 'Distance-ERR' endif end function Distance function cross_prod_xyz( A, B ) result(C) implicit none real(8), intent(in) :: A(3), B(3) real(8) :: C(3) C(1) = A(2)*B(3) - A(3)*B(2) C(2) = A(3)*B(1) - A(1)*B(3) C(3) = A(1)*B(2) - A(2)*B(1) end function cross_prod_xyz function cross_prod_typ( theArc ) result(C) implicit none type(arc), intent(inout) :: theArc real(8) :: C(3) theArc%sxyz = geod2cart( theArc%sPnt ) theArc%exyz = geod2cart( theArc%ePnt ) C(1) = theArc%sxyz%y*theArc%exyz%z - theArc%sxyz%z*theArc%exyz%y C(2) = theArc%sxyz%z*theArc%exyz%x - theArc%sxyz%x*theArc%exyz%z C(3) = theArc%sxyz%x*theArc%exyz%y - theArc%sxyz%y*theArc%exyz%x end function cross_prod_typ function ftngeod2cart( lon, lat ) result( Cart ) !--------------------------------------------------- ! map geodetic (lambda,phi) -> cartesian (x,y,z) !--------------------------------------------------- implicit none real(8), intent(in) :: lon, lat real(8) :: Cart(3) real(8) :: tmp real(8) :: rlon, rlat real(8) :: x, y, z rlon = D2R*lon ; rlat = D2R*lat tmp = cos(rlat) z = sin(rlat) x = tmp*cos(rlon) y = tmp*sin(rlon) Cart(:) = (/ x,y,z /) end function ftngeod2cart function ftntypgeod2cart( geodPnt ) result(xyzPnt) !--------------------------------------------------- ! map geodetic (lambda,phi) -> cartesian (x,y,z) !--------------------------------------------------- implicit none type(geodeticPnt), intent(in) :: geodPnt type(cartPnt) :: xyzPnt real(8) :: tmp real(8) :: rlon, rlat rlon = D2R*geodPnt%lon ; rlat = D2R*geodPnt%lat tmp = cos(rlat) xyzPnt%z = sin(rlat) xyzPnt%x = tmp*cos(rlon) xyzPnt%y = tmp*sin(rlon) end function ftntypgeod2cart function ftncart2geod( x, y, z ) result( Pnt ) implicit none real(8), intent(in) :: x, y, z real(8) :: Pnt(2) real(8) :: tmp real(8) :: lon, lat lat = asin( z ) lon = atan2(y,x) Pnt(:) = (/ lon,lat /) end function ftncart2geod function ArcsEqual( arc1, arc2 ) result( Equal ) !--------------------------------------------------- ! check if arcs are identical !--------------------------------------------------- type(arc), intent(in) :: arc1, arc2 logical :: Equal Equal = PntsEqual( arc1%sPnt%lon,arc2%sPnt%lon ) if( Equal ) then Equal = PntsEqual( arc1%sPnt%lat,arc2%sPnt%lat ) endif end function ArcsEqual function cartsEqual( P1, P2 ) result( Equal ) !--------------------------------------------------- ! check if arcs are identical !--------------------------------------------------- real(8), intent(in) :: P1(3), P2(3) logical :: Equal integer :: n do n = 1,3 Equal = PntsEqual( P1(n),P2(n) ) if( .not. Equal ) then exit endif enddo end function cartsEqual function PntsEqual( pnt1, pnt2 ) result( Equal ) !--------------------------------------------------- ! check if two numbers are identical !--------------------------------------------------- real(8), intent(in) :: pnt1, pnt2 logical :: Equal real(8), parameter :: REL_ROUNDOFF = 1.e-6_8 real(8) :: tmp if( abs(pnt1) <= SMALL .and. abs(pnt2) <= SMALL ) then Equal = .true. elseif( pnt1 /= ZERO .and. pnt2 /= ZERO ) then Equal = abs(pnt2 - pnt1) <= REL_ROUNDOFF*max( abs(pnt1),abs(pnt2) ) else tmp = max(abs(pnt1),abs(pnt2)) Equal = tmp <= REL_ROUNDOFF endif end function PntsEqual function ftncartNorm( x, y, z ) result( Norm ) !--------------------------------------------------- ! Form normed 3d vector !--------------------------------------------------- real(8), intent(in) :: x, y, z real(8) :: Norm(3) real(8) :: Factor Factor = sqrt( x*x + y*y + z*z ) Norm(:) = (/ x/Factor, y/Factor, z/Factor /) end function ftncartNorm function typcartNorm( X ) result( Norm ) !--------------------------------------------------- ! Form normed 3d vector !--------------------------------------------------- real(8), intent(in) :: X(3) real(8) :: Norm(3) real(8) :: Factor Factor = sqrt( sum(X(:)*X(:)) ) Norm(:) = X(:)/Factor end function typcartNorm function MatSlv( A, b ) result(x) !--------------------------------------------------- ! Solve 2x2 linear system !--------------------------------------------------- real(8), intent(in) :: A(2,2) real(8), intent(in) :: b(2) real(8) :: x(2) real(8) :: det real(8) :: cofac(2,2) det = A(1,1)*A(2,2) - A(2,1)*A(1,2) if( Det == ZERO ) then write(*,*) 'MatSlv: zero determinant' write(*,*) 'MatSlv: A' write(*,*) A(1,1),A(1,2) write(*,*) A(2,1),A(2,2) Stop 'Singular matrix' endif cofac(1,1) = A(2,2) ; cofac(2,2) = A(1,1) cofac(1,2) = -A(1,2) ; cofac(2,1) = -A(2,1) x(:) = matmul( cofac,b ) x(:) = x(:)/det end function MatSlv subroutine GCXsect( ArcsIn, geodXsect, cartXsect ) !--------------------------------------------------- ! get great cirle intersection !--------------------------------------------------- !--------------------------------------------------- ! dummy arguments !--------------------------------------------------- type(arc), intent(in) :: ArcsIn(2) type(geodeticPnt), intent(out), optional :: geodXsect(2) type(cartPnt), intent(inout), optional :: cartXsect(2) !--------------------------------------------------- ! local variables !--------------------------------------------------- integer :: n real(8) :: gP(2,2) real(8) :: E1(3,2), E2(3,2) real(8) :: X(3,2) real(8) :: XYZ(3) logical :: Set2Zero type(arc) :: Arcs(2) Arcs(:) = ArcsIn(:) !--------------------------------------------------- ! get x,y,z of arc1 start,end points !--------------------------------------------------- E1(:,1) = geod2cart( Arcs(1)%sPnt%lon,Arcs(1)%sPnt%lat ) E1(:,2) = geod2cart( Arcs(1)%ePnt%lon,Arcs(1)%ePnt%lat ) !--------------------------------------------------- ! get x,y,z of arc2 start,end points !--------------------------------------------------- E2(:,1) = geod2cart( Arcs(2)%sPnt%lon,Arcs(2)%sPnt%lat ) E2(:,2) = geod2cart( Arcs(2)%ePnt%lon,Arcs(2)%ePnt%lat ) !--------------------------------------------------- ! cross for arc1, arc2; then normalize !--------------------------------------------------- X(:,1) = cross_prod( E1(:,1), E1(:,2) ) X(:,2) = cross_prod( E2(:,1), E2(:,2) ) X(:,1) = cartNorm( X(:,1) ) X(:,2) = cartNorm( X(:,2) ) !--------------------------------------------------- ! x,y,z intersect point !--------------------------------------------------- XYZ(:) = cross_prod( X(:,1),X(:,2) ) !--------------------------------------------------- ! check x,z intersect for proximity to zero !--------------------------------------------------- ! if( PntsEqual( XYZ(3),ZERO ) ) then ! XYZ(3) = ZERO ! endif ! if( .not. Set2Zero ) then ! gP(2,1) = R2D*atan2( XYZ(3),sqrt(XYZ(1)*XYZ(1)+XYZ(2)*XYZ(2)) ) ! gP(1,1) = R2D*atan2( XYZ(1),XYZ(2) ) ! gP(1,1) = mod( gp(1,1),THREE60 ) ! gP(2,2) = -gP(2,1) ! gP(1,2) = gP(1,1) + ONE80 ! gP(1,2) = mod( gp(1,2),THREE60 ) ! else ! gP(:,:) = ZERO ! endif ! gP(2,1) = R2D*atan2( XYZ(3),sqrt(XYZ(1)*XYZ(1)+XYZ(2)*XYZ(2)) ) ! gP(1,1) = R2D*atan2( XYZ(2),XYZ(1) ) !--------------------------------------------------- ! gP in radians !--------------------------------------------------- gP(2,1) = atan2( XYZ(3),sqrt(XYZ(1)*XYZ(1)+XYZ(2)*XYZ(2)) ) gP(1,1) = atan2( XYZ(2),XYZ(1) ) do n = 1,2 if( PntsEqual( gP(n,1),ZERO ) ) then gP(n,1) = ZERO endif enddo gP(2,2) = -gP(2,1) if( gP(1,1) > ZERO ) then gP(1,2) = gP(1,1) - PI elseif( gP(1,1) < ZERO ) then gP(1,2) = gP(1,1) + PI else gP(1,2) = gP(1,1) endif do n = 1,2 if( PntsEqual( gP(n,2),ZERO ) ) then gP(n,2) = ZERO endif enddo gP(:,:) = R2D*gP(:,:) if( present(geodXsect) ) then geodXsect(1)%lon = gP(1,1) ; geodXsect(1)%lat = gP(2,1) geodXsect(2)%lon = gP(1,2) ; geodXsect(2)%lat = gP(2,2) endif if( present(cartXsect) ) then cartXsect(1)%x = XYZ(1) ; cartXsect(1)%y = XYZ(2) ; cartXsect(1)%z = XYZ(3) cartXsect(2)%x = -XYZ(1) ; cartXsect(2)%y = -XYZ(2) ; cartXsect(2)%z = -XYZ(3) endif end subroutine GCXsect function getSxyz( AnArc, xyzPnt ) result( S ) !--------------------------------------------------- ! Calculate the distance from 3 pnts on the GC !--------------------------------------------------- type(arc), intent(inout) :: AnArc real(8), intent(in) :: xyzPnt(3) real(8) :: S real(8) :: A(2,2) real(8) :: B(2), X(2) AnArc%sxyz = geod2cart( AnArc%spnt ) AnArc%exyz = geod2cart( AnArc%epnt ) AnArc%Angle = Distance( AnArc%sPnt,AnArc%ePnt ) A(1,1) = AnArc%sxyz%x ; A(1,2) = AnArc%exyz%x A(2,1) = AnArc%sxyz%y ; A(2,2) = AnArc%exyz%y B(:) = xyzPnt(1:2) X(:) = MatSlv( A, B ) S = asin( X(2)*sin(AnArc%Angle) ) end function getSxyz function getSpnt( AnArc, Pnt ) result( S ) !--------------------------------------------------- ! Calculate the distance from 3 pnts on the GC !--------------------------------------------------- type(arc), intent(in) :: AnArc type(geodeticPnt), intent(in) :: Pnt real(8) :: S real(8) :: dist dist = Distance( AnArc%sPnt,AnArc%ePnt ) S = dist*abs(Pnt%lat - AnArc%sPnt%lat)/abs(AnArc%ePnt%lat - AnArc%sPnt%lat) end function getSpnt function getXYZ( AnArc, S ) result(XYZ) !--------------------------------------------------- ! Parametric point on arc GC at "distance" = S !--------------------------------------------------- type(arc), intent(inout) :: AnArc real(8), intent(in) :: S real(8) :: XYZ(3) real(8) :: sind, sins, Factor real(8) :: Sxyz(3), Exyz(3) if( AnArc%Angle == NOTSET ) then AnArc%Angle = Distance( AnArc%sPnt,AnArc%ePnt ) endif sind = sin( AnArc%Angle ) sins = sin( S ) Factor = sin(AnArc%Angle - S) if( S /= ZERO .and. S /= ONE ) then Sxyz(:) = geod2cart( AnArc%sPnt%lon,AnArc%sPnt%lat ) Exyz(:) = geod2cart( AnArc%ePnt%lon,AnArc%ePnt%lat ) XYZ(:) = (Factor*Sxyz(:) + sins*Exyz(:))/sind elseif( S == ZERO ) then XYZ(:) = Sxyz(:) else XYZ(:) = Exyz(:) endif end function getXYZ function sTriArea( TriArcs ) result( Area ) !--------------------------------------------------- ! Compute area of spherical triange !--------------------------------------------------- type(arc), intent(in) :: TriArcs(3) real(8) :: Area integer :: n real(8), pointer :: a, b, c real(8), target :: ang(3) real(8) :: s, E a => ang(1) ; b => ang(2) ; c => ang(3) do n = 1,3 if( TriArcs(n)%Angle == NOTSET ) then ang(n) = Distance( TriArcs(n)%sPnt,TriArcs(n)%ePnt ) else ang(n) = TriArcs(n)%Angle endif enddo s = ONEHALF*sum( ang(:) ) E = tan(ONEHALF*s)*tan(ONEHALF*(s-a))*tan(ONEHALF*(s-b))*tan(ONEHALF*(s-c)) E = sqrt(abs(E)) Area = FOUR*atan(E) end function sTriArea function sQuadArea( QuadArcs ) result( Area ) !--------------------------------------------------- ! Compute area of spherical triange !--------------------------------------------------- type(arc), intent(in) :: QuadArcs(4) real(8) :: Area integer :: n type(arc) :: sTriArcs(3) sTriArcs(1:2) = QuadArcs(1:2) sTriArcs(3)%sPnt = QuadArcs(2)%ePnt sTriArcs(3)%ePnt = QuadArcs(1)%sPnt Area = sTriArea( sTriArcs ) sTriArcs(1:2) = QuadArcs(3:4) sTriArcs(3)%sPnt = QuadArcs(1)%sPnt sTriArcs(3)%ePnt = QuadArcs(2)%ePnt Area = Area + sTriArea( sTriArcs ) end function sQuadArea subroutine diagPnts( Pnts ) !--------------------------------------------------- ! pnts diagnostics !--------------------------------------------------- type(geodeticPnt), intent(in) :: Pnts(:) integer :: n, nPnts nPnts = size(Pnts) write(*,*) '' do n = 1,nPnts write(*,'(''Pnt('',i2,'') = ('',1p,g15.7,'','',g15.7,'')'')') n,Pnts(n)%lon,Pnts(n)%lat enddo write(*,*) '' end subroutine diagPnts subroutine diagArc( theArc ) !--------------------------------------------------- ! arc diagnostics !--------------------------------------------------- type(arc), intent(in) :: theArc write(*,*) '' write(*,*) 'Arc sPnt = ',theArc%sPnt%lon,theArc%sPnt%lat write(*,*) 'Arc ePnt = ',theArc%ePnt%lon,theArc%ePnt%lat write(*,*) 'Arc sxyz = ',theArc%sxyz%x,theArc%sxyz%y,theArc%sxyz%z write(*,*) 'Arc exyz = ',theArc%exyz%x,theArc%exyz%y,theArc%exyz%z write(*,*) 'Arc angle = ',theArc%Angle write(*,*) '' end subroutine diagArc function Centroid( Pnts ) result( center ) !--------------------------------------------------- ! Compute centroid of spherical polynomial !--------------------------------------------------- type(geodeticPnt), intent(in) :: Pnts(:) type(geodeticPnt) :: center integer :: n, nNext, nPnts real(8) :: x0,x1,y0,y1 real(8) :: A, Area, Mult center%lon = ZERO ; center%lat = ZERO Area = ZERO nPnts = size(Pnts) do n = 1,nPnts nNext = mod(n+1,nPnts+1) nNext = max( nNext,1 ) x0 = Pnts(n)%lon ; y0 = Pnts(n)%lat x1 = Pnts(nNext)%lon ; y1 = Pnts(nNext)%lat A = x0 * y1 - x1 * y0 Area = Area + A center%lon = center%lon + A*(x0 + x1) center%lat = center%lat + A*(y0 + y1) enddo Area = ONEHALF*Area Mult= ONE/(6._8*Area) center%lon = center%lon*Mult center%lat = center%lat*Mult end function Centroid function Orientation( Pnts ) result( Hand ) !--------------------------------------------------- ! Compute centroid of spherical polynomial !--------------------------------------------------- type(geodeticPnt), intent(in) :: Pnts(:) real(8) :: Hand integer :: n, nNext, nPnts Hand = ZERO nPnts = size(Pnts) do n = 1,nPnts nNext = mod(n+1,nPnts+1) nNext = max( nNext,1 ) Hand = Hand + (Pnts(nNext)%lon - Pnts(n)%lon) & * (Pnts(nNext)%lat + Pnts(n)%lat) enddo end function Orientation function pnt_in_triangle( x, y ) result(IsInside) !--------------------------------------------------------------- ! determine whether input point is in triangle !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- real(8), intent(in) :: x(:) real(8), intent(in) :: y(:) logical :: IsInside !--------------------------------------------------------------- ! local variables !--------------------------------------------------------------- real(8) :: a, b, c a = (x(1) - x(4))*(y(2) - y(4)) - (x(2) - x(4))*(y(1) - y(4)) b = (x(2) - x(4))*(y(3) - y(4)) - (x(3) - x(4))*(y(2) - y(4)) if( sign( ONE,a) == sign( ONE,b ) ) then c = (x(3) - x(4))*(y(1) - y(4)) - (x(1) - x(4))*(y(3) - y(4)) IsInside = sign( ONE,b ) == sign( ONE,c ) else IsInside = .false. endif end function pnt_in_triangle function pnt_in_SphTri( Tri, pnt ) result(IsInside) !--------------------------------------------------------------- ! determine whether input spherical point is in spherical triangle !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- type(spolygon), intent(in) :: Tri type(geodeticPnt), intent(in) :: Pnt logical :: IsInside integer :: m, mp1, n type(geodeticPnt) :: Xsect(2) type(arc) :: TstArc logical :: found if( size(Tri%arcs) /= 3 ) then write(*,'('' Pnt_in_SphTri; incoming Tri does not have 3 arcs'')') stop 'ArgErr' endif write(*,*) 'pnt_in_sphtri: test point' call diagPnts( (/pnt/) ) !--------------------------------------------------------------- ! check for pnt in arc of triangle !--------------------------------------------------------------- found = .false. do m = 1,3 if( Pnt_in_Arc( Tri%arcs(m),Pnt ) ) then found = .true. exit endif enddo if( .not. found ) then TstArc%ePnt = Pnt arc_loop: & do m = 1,3 mp1 = max( mod(m+1,4),1 ) TstArc%sPnt = Tri%arcs(m)%sPnt call GCXsect( (/ Tri%arcs(mp1),TstArc /),Xsect ) write(*,*) 'pnt_in_sphtri: TstArc' call diagArc( TstArc ) write(*,*) 'pnt_in_sphtri: 2nd arc' call diagArc( Tri%arcs(mp1) ) write(*,*) 'pnt_in_sphtri: xsects' call diagPnts( Xsect ) found = .false. do n = 1,2 if( Pnt_in_Arc( Tri%arcs(mp1),Xsect(n) ) ) then found = .true. exit endif enddo if( .not. found ) then exit arc_loop endif enddo arc_loop endif IsInside = found end function pnt_in_SphTri function pnt_in_SphPoly( Poly, pnt ) result(IsInside) !--------------------------------------------------------------- ! determine whether input spherical point is in spherical polygon !--------------------------------------------------------------- !--------------------------------------------------------------- ! dummy arguments !--------------------------------------------------------------- type(spolygon), intent(in) :: Poly type(geodeticPnt), intent(in) :: Pnt logical :: IsInside integer :: m, n, nArcs type(geodeticPnt) :: PolyVtx(size(Poly%arcs)) type(geodeticPnt) :: Xsect(2) type(arc) :: TstArc logical :: found nArcs = size(Poly%arcs) !--------------------------------------------------------------- ! check for pnt on poly arc !--------------------------------------------------------------- do m = 1,nArcs if( Pnt_in_Arc( Poly%arcs(m),Pnt ) ) then IsInside = .true. return endif enddo !--------------------------------------------------------------- ! full arc test !--------------------------------------------------------------- PolyVtx(1:nArcs) = Poly%arcs(1:nArcs)%sPnt TstArc%sPnt = Centroid( PolyVtx ) ; TstArc%ePnt = Pnt found = .false. arc_loop: & do m = 1,nArcs call GCXsect( (/ Poly%arcs(m),TstArc /),Xsect ) do n = 1,2 found = Pnt_in_Arc( Poly%arcs(m),Xsect(n) ) if( found ) then found = Pnt_in_Arc( TstArc,Xsect(n) ) if( found ) then found = .false. exit arc_loop endif endif enddo enddo arc_loop IsInside = found end function pnt_in_SphPoly end module Sph_utils