!! Copyright (C) Stichting Deltares, 2005-2023.
!!
!! This file is part of iMOD.
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see .
!!
!! Contact: imod.support@deltares.nl
!! Stichting Deltares
!! P.O. Box 177
!! 2600 MH Delft, The Netherlands.
!!
MODULE MOD_INTERSECT
USE MOD_POLINT, ONLY : POL1LOCATE
USE MOD_INTERSECT_PAR
USE MOD_QKSORT
CONTAINS
!###======================================================================
SUBROUTINE INTERSECT_EQUI(XMIN1,XMAX1,YMIN1,YMAX1,CSX_IN,CSY_IN,XIN1,XIN2,YIN1,YIN2,N,LHFB,LREFINE)
!###======================================================================
IMPLICIT NONE
REAL(KIND=DP_KIND),INTENT(IN) :: XMIN1,XMAX1,YMIN1,YMAX1,CSX_IN,CSY_IN
REAL(KIND=DP_KIND),INTENT(IN) :: XIN1,XIN2,YIN1,YIN2
INTEGER,INTENT(INOUT) :: N
LOGICAL,INTENT(IN) :: LHFB,LREFINE
REAL(KIND=DP_KIND) :: X,Y,XMN,XMX,YMN,YMX,DX,DY,LENG,TD,X1,X2,Y1,Y2,CSX,CSY,A,B,XMIN,XMAX,YMIN,YMAX
INTEGER :: I,J,ICOL,IROW,ID,N_IN !,NOFFSET
LOGICAL :: LFLIP
XMIN=XMIN1-CSX_IN
XMAX=XMAX1+CSX_IN
YMIN=YMIN1-CSY_IN
YMAX=YMAX1+CSY_IN
X1=XIN1; Y1=YIN1; X2=XIN2; Y2=YIN2
IF(.NOT.ASSOCIATED(XA))ALLOCATE(XA(1000)); IF(.NOT.ASSOCIATED(YA))ALLOCATE(YA(1000))
IF(.NOT.ASSOCIATED(FA))ALLOCATE(FA(1000)); IF(.NOT.ASSOCIATED(LN))ALLOCATE(LN(1000))
IF(.NOT.ASSOCIATED(CA))ALLOCATE(CA(1000)); IF(.NOT.ASSOCIATED(RA))ALLOCATE(RA(1000))
!## initiate all if first call
IF(N.EQ.0)THEN; CA=0; RA=0; ENDIF
CSX=CSX_IN
CSY=CSY_IN
IF(LHFB.OR.LREFINE)THEN
CSX=CSX/2.0D0
CSY=CSY/2.0D0
ENDIF
N_IN=N
IF(.NOT.INTERSECT_EQUATION(X1,X2,Y1,Y2,N,A,B,LFLIP))RETURN
!## find search box - result can be negative, does not matter!
IF(X1.LE.XMIN)THEN
XMN=XMIN
ELSE
! I=0 !1
XMN=XMIN+CSX*(INT((X1-XMIN)/CSX)) !+I)
ENDIF
IF(X2.GE.XMAX)THEN
XMX=XMAX
ELSE
! I=1 !0
! DX=X2-XMIN
! IF(MOD(DX,CSX).EQ.0.0D0)I=-1
XMX=XMIN+CSX*(INT((X2-XMIN)/CSX)) !+I)
ENDIF
Y=MIN(Y1,Y2)
IF(Y.LE.YMIN)THEN
YMN=YMIN
ELSE
I =0
DY=YMAX-Y
IF(MOD(DY,CSY).EQ.0.0D0)I=-1
YMN=YMAX-CSY*(INT((YMAX-Y)/CSY)+I)
ENDIF
Y=MAX(Y1,Y2)
IF(Y.GE.YMAX)THEN
YMX=Y
ELSE
I =1
YMX=YMAX-CSY*(INT((YMAX-Y)/CSY)+I)
ENDIF
!## not for horizontal lines
IF(IHOR.EQ.0)THEN
!## continue seach rest of intersections
!## try intersections with x-axes firstly
! IF(YMN.NE.YMX)THEN
IF(LHFB)THEN
Y=YMN
!## otherwise gaps
IF(MOD(YMN,CSY).EQ.0.0D0)Y=YMN-CSY
ELSE
Y=YMN-CSY
ENDIF
DO
Y=Y+CSY
IF(Y.GT.YMX)EXIT
!## array overwritten
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
IF(IVER.EQ.1)THEN
XA(N)=X1 !## same as xmx
ELSE
XA(N)=(Y-B)/A
ENDIF
YA(N)=Y
ENDDO
! ENDIF
ENDIF
!## not for vertical lines
IF(IVER.EQ.0)THEN
! IF(XMN.NE.XMX)THEN
!## try intersections with y-axes secondly
IF(LHFB)THEN
X=XMN
!## otherwise gaps
IF(MOD(XMN,CSX).EQ.0.0D0)X=XMN-CSX
ELSE
X=XMN-CSX
ENDIF
DO
X=X+CSX
IF(X.GT.XMX)EXIT
!## array overwritten
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
XA(N)=X
IF(IHOR.EQ.1)THEN
YA(N)=Y1 !## same as ymx
ELSE
YA(N)=A*X+B
ENDIF
ENDDO
! ENDIF
ENDIF
!## restore cellsize
CSX=CSX_IN
CSY=CSY_IN
! DX=X1-X2; DY=Y2-Y1
! CALL INTERSECT_SORT(DX,DY,N_IN+1,N,LFLIP)
!## sample each of the point to determine irow/icol, overwrite point with this
!## skip first and last, they represented already by the second and one-last point
! IF(.NOT.LHFB)THEN
XMN=MIN(X1,X2); XMX=MAX(X1,X2)
YMN=MIN(Y1,Y2); YMX=MAX(Y1,Y2)
!## outside current line of the cross-section
J=N_IN
DO I=N_IN+1,N
IF(XA(I).LT.XMN.OR.XA(I).GT.XMX.OR.YA(I).LT.YMN.OR.YA(I).GT.YMX)CYCLE
J =J+1
XA(J)=XA(I)
YA(J)=YA(I)
ENDDO
!## new number of points
N=J
DX=X1-X2; DY=Y2-Y1
CALL INTERSECT_SORT(DX,DY,N_IN+1,N,LFLIP)
! ENDIF
! NOFFSET=2
! DO I=N_IN+NOFFSET,N
DO I=N_IN+2,N
!## mid point
X =(XA(I-1)+XA(I))/2.0D0
Y =(YA(I-1)+YA(I))/2.0D0
!## outside current line of the cross-section
IF(X.LT.XMN.OR.X.GT.XMX.OR.Y.LT.YMN.OR.Y.GT.YMX)THEN
CA(I-1)=0
RA(I-1)=0
LN(I-1)=0.0D0
CYCLE
ENDIF
IF(LHFB)THEN
ICOL=0; IF((X.GT.MIN(X1,X2).AND.X.LT.MAX(X1,X2)).OR.X1.EQ.X2)ICOL=INT((X-XMIN)/CSX)+1
IROW=0; IF((Y.GT.MIN(Y1,Y2).AND.Y.LT.MAX(Y1,Y2)).OR.Y1.EQ.Y2)IROW=INT((YMAX-Y)/CSY)+1
ELSE
ICOL=0
IF(X.GE.XMN.AND.X.LE.XMX)THEN
ICOL=INT((X-XMIN)/CSX)+1
ENDIF
IROW=0
IF(Y.GE.YMN.AND.Y.LE.YMX)THEN
IROW=INT((YMAX-Y)/CSY)+1
ENDIF
ENDIF
!## initiate maximal distance
TD=SQRT(CSX**2.0D0+CSY**2.0D0); ID=0 !## fraction
IF(ICOL.GT.0.AND.IROW.GT.0)THEN
CALL INTERSECT_NCORNER(ID,TD,XMIN+(ICOL-1)*CSX,YMAX-(IROW-1)*CSY,X,Y,1)
CALL INTERSECT_NCORNER(ID,TD,XMIN+ ICOL *CSX,YMAX-(IROW-1)*CSY,X,Y,2)
CALL INTERSECT_NCORNER(ID,TD,XMIN+ ICOL *CSX,YMAX- IROW *CSY,X,Y,3)
CALL INTERSECT_NCORNER(ID,TD,XMIN+(ICOL-1)*CSX,YMAX- IROW *CSY,X,Y,4)
ENDIF
FA(I-1)=DBLE(ID)
DX =XA(I)-XA(I-1)
DY =YA(I)-YA(I-1)
LENG=DX**2.0D0+DY**2.0D0; IF(LENG.GT.0.0D0)LENG=SQRT(LENG)
!## store results in row/column indices - correct as model has been temp. increased
CA(I-1)=ICOL-1
RA(I-1)=IROW-1
LN(I-1)=LENG
END DO
! N=N-2
N=N-1
END SUBROUTINE INTERSECT_EQUI
!###======================================================================
SUBROUTINE INTERSECT_NCORNER(ID,TD,X1,Y1,XC,YC,JD)
!###======================================================================
IMPLICIT NONE
INTEGER,INTENT(INOUT) :: ID
REAL(KIND=DP_KIND),INTENT(INOUT) :: TD
INTEGER,INTENT(IN) :: JD
REAL(KIND=DP_KIND),INTENT(IN) :: X1,Y1,XC,YC
REAL(KIND=DP_KIND) :: D
!## get nearest corner
D=(XC-X1)**2.0D0+(YC-Y1)**2.0D0
IF(D.GT.0.0D0)THEN
D=SQRT(D)
IF(D.LT.TD)THEN
ID=JD
TD=D
ENDIF
ELSE
ID=JD
TD=0.0D0
ENDIF
END SUBROUTINE INTERSECT_NCORNER
!###======================================================================
SUBROUTINE INTERSECT_NONEQUI(DELR,DELC,NROW,NCOL,XIN1,XIN2,YIN1,YIN2,N,LHFB,LREFINE)
!###======================================================================
IMPLICIT NONE
INTEGER,INTENT(IN) :: NROW,NCOL
REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(0:NCOL) :: DELR
REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(0:NROW) :: DELC
REAL(KIND=DP_KIND),INTENT(IN) :: XIN1,XIN2,YIN1,YIN2
INTEGER,INTENT(INOUT) :: N
LOGICAL,INTENT(IN) :: LHFB,LREFINE
REAL(KIND=DP_KIND) :: X,Y,XMN,XMX,YMN,YMX,DX,DY,LENG,XMIN,YMIN,XMAX,YMAX,X1,X2,Y1,Y2,CS,TD,CSX,CSY,A,B
INTEGER :: I,J,ICOL,IROW,IMN,JMN,ID,N_IN
LOGICAL :: LFLIP
X1=XIN1; Y1=YIN1; X2=XIN2; Y2=YIN2
IF(.NOT.ASSOCIATED(XA))ALLOCATE(XA(1000)); IF(.NOT.ASSOCIATED(YA))ALLOCATE(YA(1000))
IF(.NOT.ASSOCIATED(FA))ALLOCATE(FA(1000)); IF(.NOT.ASSOCIATED(LN))ALLOCATE(LN(1000))
IF(.NOT.ASSOCIATED(CA))ALLOCATE(CA(1000)); IF(.NOT.ASSOCIATED(RA))ALLOCATE(RA(1000))
!## initiate all if first call
IF(N.EQ.0)THEN; CA=0; RA=0; ENDIF
XMIN=DELR(0 )
XMAX=DELR(NCOL)
YMIN=DELC(NROW)
YMAX=DELC(0 )
XMIN=XMIN-(DELR(1) -DELR(0) )
XMAX=XMAX+(DELR(NCOL)-DELR(NCOL-1))
YMIN=YMIN-(DELC(NROW)-DELR(NROW-1))
YMAX=YMAX+(DELC(1) -DELC(0) )
N_IN=N
IF(.NOT.INTERSECT_EQUATION(X1,X2,Y1,Y2,N,A,B,LFLIP))RETURN
!## continue search rest of intersections
!## try intersections with y-axes firstly
IMN=0; DO I=1,NROW
IF(DELC(I).GT.MIN(Y1,Y2))THEN; YMN=DELC(I); IMN=I; ENDIF
END DO
JMN=0; DO I=NROW,1,-1
IF(DELC(I).LT.MAX(Y1,Y2))THEN; YMX=DELC(I); JMN=I; ENDIF
END DO
IF(IMN.GT.0.AND.JMN.GT.0)THEN
!## not for horizontal lines
IF(IHOR.EQ.0)THEN
CS=DELC(IMN-1)-DELC(IMN)
IF(LHFB.OR.LREFINE)CS=CS/2.0D0
! IF(LHFB)THEN
Y =YMN !-CS
! ELSE
! Y =YMN-CS
! ENDIF
DO
Y=Y+CS
IF(Y.GT.YMX)EXIT
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
IF(IVER.EQ.1)THEN
XA(N)=X1 !## same as xmx
ELSE
XA(N)=(Y-B)/A
ENDIF
YA(N)=Y
!## double intersections, for better estimate for hfb
IF(LHFB.OR.LREFINE)THEN
Y=Y+CS
!## array overwritten
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
IF(IVER.EQ.1)THEN
XA(N)=X1 !## same as xmx
ELSE
XA(N)=(Y-B)/A
ENDIF
YA(N)=Y
ENDIF
IMN =IMN-1
!## model is not bigger than line-segment
IF(IMN.LE.0)EXIT
CS=DELC(IMN)-DELC(IMN+1)
IF(LHFB.OR.LREFINE)CS=CS/2.0D0
ENDDO
ENDIF
ENDIF
!## try intersections with x-axes secondly
IMN=0; DO I=NCOL,1,-1
IF(DELR(I).GT.X1)THEN; XMN=DELR(I); IMN=I; ENDIF
END DO
JMN=0; DO I=0,NCOL
IF(DELR(I).LT.X2)THEN; XMX=DELR(I); JMN=I; ENDIF
END DO
IF(IMN.GT.0.AND.JMN.GT.0)THEN
!## not for vertical lines
IF(IVER.EQ.0)THEN
CS=DELR(IMN-1)-DELR(IMN)
IF(LHFB.OR.LREFINE)CS=CS/2.0D0
! IF(LHFB)THEN
X =XMN !-CS
! ELSE
! X =XMN-CS
! ENDIF
DO
X=X+CS
IF(X.GT.XMX)EXIT
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
XA(N)=X
IF(IHOR.EQ.1)THEN
YA(N)=Y1 !## same as ymx
ELSE
YA(N)=A*X+B
ENDIF
!## double intersections, for better estimate for hfb
IF(LHFB.OR.LREFINE)THEN
X=X+CS
!## array overwritten
N=N+1; CALL INTERSECT_RESIZEVECTORS(N)
XA(N)=X
IF(IHOR.EQ.1)THEN
YA(N)=Y1 !## same as ymx
ELSE
YA(N)=A*X+B
ENDIF
ENDIF
IMN =IMN+1
!## model is not bigger than line-segment
IF(IMN.GT.NCOL)EXIT
CS=DELR(IMN)-DELR(IMN-1)
IF(LHFB)CS=CS/2.0D0
ENDDO
ENDIF
ENDIF
! !## sort intersections, determined by the one with the largest difference
! DX=X1-X2; DY=Y2-Y1
! CALL INTERSECT_SORT(DX,DY,N_IN+1,N,LFLIP)
!## sample each of the point to determine irow/icol, overwrite point with this
!## skip first and last, they represented already by the second and one-last point
XMN=MIN(X1,X2); XMX=MAX(X1,X2)
YMN=MIN(Y1,Y2); YMX=MAX(Y1,Y2)
!## outside current line of the cross-section
J=N_IN
DO I=N_IN+1,N
IF(XA(I).LT.XMN.OR.XA(I).GT.XMX.OR.YA(I).LT.YMN.OR.YA(I).GT.YMX)CYCLE
J =J+1
XA(J)=XA(I)
YA(J)=YA(I)
ENDDO
!## new number of points
N=J
!## sort intersections, determined by the one with the largest difference
DX=X1-X2; DY=Y2-Y1
CALL INTERSECT_SORT(DX,DY,N_IN+1,N,LFLIP)
DO I=N_IN+2,N
!## mid point
X =(XA(I-1)+XA(I))/2.0D0
Y =(YA(I-1)+YA(I))/2.0D0
!## outside current line of the cross-section
IF(X.LT.XMN.OR.X.GT.XMX.OR.Y.LT.YMN.OR.Y.GT.YMX)THEN
CA(I-1)=0
RA(I-1)=0
LN(I-1)=0.0D0
CYCLE
ENDIF
IF(LHFB)THEN
! ICOL=0; IF((X.GT.XMIN.AND.X.LT.XMAX).OR.X1.EQ.X2)CALL POL1LOCATE(DELR,NCOL+1,DBLE(X),ICOL)
! IROW=0; IF((Y.GT.YMIN.AND.Y.LT.YMAX).OR.Y1.EQ.Y2)CALL POL1LOCATE(DELC,NROW+1,DBLE(Y),IROW)
ICOL=0; IF((X.GT.MIN(X1,X2).AND.X.LT.MAX(X1,X2)).OR.X1.EQ.X2)CALL POL1LOCATE(DELR,NCOL+1,DBLE(X),ICOL)
IROW=0; IF((Y.GT.MIN(Y1,Y2).AND.Y.LT.MAX(Y1,Y2)).OR.Y1.EQ.Y2)CALL POL1LOCATE(DELC,NROW+1,DBLE(Y),IROW)
ELSE
ICOL=0
IF(X.EQ.DELR(0))THEN !XMN)THEN
ICOL=1
ELSEIF(X.EQ.DELR(NCOL))THEN !XMX)THEN
ICOL=NCOL
ELSEIF(X.GT.XMN.AND.X.LT.XMX)THEN
CALL POL1LOCATE(DELR,NCOL+1,X,ICOL)
ENDIF
IROW=0
IF(Y.EQ.DELC(NROW))THEN !YMN)THEN
IROW=NROW
ELSEIF(Y.EQ.DELC(1))THEN !YMX)THEN
IROW=1
ELSEIF(Y.GT.YMN.AND.Y.LT.YMX)THEN
CALL POL1LOCATE(DELC,NROW+1,Y,IROW)
ENDIF
ENDIF
ID=0 !## fraction
IF(ICOL.GE.1.AND.ICOL.LE.NCOL.AND. &
IROW.GE.1.AND.IROW.LE.NROW)THEN
CSX=DELR(ICOL)-DELR(ICOL-1)
CSY=DELC(IROW-1)-DELC(IROW)
!## initial maximal distance
TD=SQRT(CSX**2.0D0+CSY**2.0D0)
CALL INTERSECT_NCORNER(ID,TD,DELR(ICOL-1),DELC(IROW-1),X,Y,1)
CALL INTERSECT_NCORNER(ID,TD,DELR(ICOL) ,DELC(IROW-1),X,Y,2)
CALL INTERSECT_NCORNER(ID,TD,DELR(ICOL) ,DELC(IROW) ,X,Y,3)
CALL INTERSECT_NCORNER(ID,TD,DELR(ICOL-1),DELC(IROW) ,X,Y,4)
ELSE
ICOL=0; IROW=0
ENDIF
FA(I-1)=DBLE(ID)
DX =XA(I)-XA(I-1)
DY =YA(I)-YA(I-1)
LENG=DX**2.0D0+DY**2.0D0; IF(LENG.GT.0.0D0)LENG=SQRT(LENG)
!## store results in row/column indices
CA(I-1)=ICOL
RA(I-1)=IROW
LN(I-1)=LENG
END DO
N=N-2
! N=N-1
END SUBROUTINE INTERSECT_NONEQUI
!###======================================================================
SUBROUTINE INTERSECT_NULLIFY()
!###======================================================================
IMPLICIT NONE
NULLIFY(XA ,YA ,FA ,LN ,CA ,RA )!,ZA)
NULLIFY(XA_DUMMY,YA_DUMMY,FA_DUMMY,LN_DUMMY,CA_DUMMY,RA_DUMMY)!,ZA_DUMMY)
END SUBROUTINE INTERSECT_NULLIFY
!###======================================================================
SUBROUTINE INTERSECT_DEALLOCATE()
!###======================================================================
IMPLICIT NONE
IF(ASSOCIATED(XA))DEALLOCATE(XA)
IF(ASSOCIATED(YA))DEALLOCATE(YA)
! IF(ASSOCIATED(ZA))DEALLOCATE(ZA)
IF(ASSOCIATED(FA))DEALLOCATE(FA)
IF(ASSOCIATED(LN))DEALLOCATE(LN)
IF(ASSOCIATED(CA))DEALLOCATE(CA)
IF(ASSOCIATED(RA))DEALLOCATE(RA)
END SUBROUTINE INTERSECT_DEALLOCATE
!###======================================================================
SUBROUTINE INTERSECT_RESIZEVECTORS(N)
!###======================================================================
IMPLICIT NONE
INTEGER,PARAMETER :: DN=100
INTEGER,INTENT(IN) :: N
INTEGER :: MM,NN,I
IF(N.LE.SIZE(XA))RETURN
NN=SIZE(XA)
MM=MAX(N,NN+DN)
ALLOCATE(XA_DUMMY(MM),YA_DUMMY(MM),FA_DUMMY(MM),LN_DUMMY(MM),CA_DUMMY(MM),RA_DUMMY(MM))
DO I=1,NN
XA_DUMMY(I)=XA(I)
YA_DUMMY(I)=YA(I)
FA_DUMMY(I)=FA(I)
LN_DUMMY(I)=LN(I)
CA_DUMMY(I)=CA(I)
RA_DUMMY(I)=RA(I)
ENDDO
DO I=NN+1,MM
XA_DUMMY(I)=0.0D0
YA_DUMMY(I)=0.0D0
FA_DUMMY(I)=0.0D0
LN_DUMMY(I)=0.0D0
CA_DUMMY(I)=0
RA_DUMMY(I)=0
ENDDO
DEALLOCATE(XA,YA,FA,LN,CA,RA)
XA=>XA_DUMMY
YA=>YA_DUMMY
FA=>FA_DUMMY
LN=>LN_DUMMY
CA=>CA_DUMMY
RA=>RA_DUMMY
END SUBROUTINE INTERSECT_RESIZEVECTORS
!###======================================================================
LOGICAL FUNCTION INTERSECT_EQUATION(X1,X2,Y1,Y2,N,A,B,LFLIP)
!###======================================================================
IMPLICIT NONE
REAL(KIND=DP_KIND),INTENT(INOUT) :: X1,X2,Y1,Y2,A,B
LOGICAL,INTENT(OUT) :: LFLIP
INTEGER,INTENT(INOUT) :: N
REAL(KIND=DP_KIND) :: X,Y,DX,DY
INTERSECT_EQUATION=.FALSE.
! XBEGIN=X1; YBEGIN=Y1
!## arrange x1,x2,y1,y2 such that x1