Initial commit

1 files changed, 432 insertions, 0 deletions

diff --git a/sys/gio/ncarutil/conlib/contng.f b/sys/gio/ncarutil/conlib/contng.f
new file mode 100644
index 00000000..7ebad596
--- /dev/null
+++ b/sys/gio/ncarutil/conlib/contng.f
@@ -0,0 +1,432 @@
+      SUBROUTINE CONTNG (NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
+C                                                                               
+C +-----------------------------------------------------------------+           
+C |                                                                 |           
+C |                Copyright (C) 1986 by UCAR                       |           
+C |        University Corporation for Atmospheric Research          |           
+C |                    All Rights Reserved                          |           
+C |                                                                 |           
+C |                 NCARGRAPHICS  Version 1.00                      |           
+C |                                                                 |           
+C +-----------------------------------------------------------------+           
+C                                                                               
+C                                                                               
+C
+C THIS SUBROUTINE PERFORMS TRIANGULATION.  IT DIVIDES THE X-Y
+C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA
+C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE
+C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS
+C CORRESPONDING TO THE BORDER LINE SEGMENTS.
+C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE
+C ARE LISTED COUNTER-CLOCKWISE.  POINT NUMBERS OF THE END POINTS
+C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE,
+C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.
+C THE INPUT PARAMETERS ARE
+C     NDP = NUMBER OF DATA POINTS,
+C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE
+C           X COORDINATES OF THE DATA POINTS,
+C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE
+C           Y COORDINATES OF THE DATA POINTS.
+C THE OUTPUT PARAMETERS ARE
+C     NT  = NUMBER OF TRIANGLES,
+C     IPT = ARRAY OF DIMENSION 6*NDP-15, WHERE THE POINT
+C           NUMBERS OF THE VERTEXES OF THE (IT)TH TRIANGLE
+C           ARE TO BE STORED AS THE (3*IT-2)ND, (3*IT-1)ST,
+C           AND (3*IT)TH ELEMENTS, IT=1,2,...,NT,
+C     NL  = NUMBER OF BORDER LINE SEGMENTS,
+C     IPL = ARRAY OF DIMENSION 6*NDP, WHERE THE POINT
+C           NUMBERS OF THE END POINTS OF THE (IL)TH BORDER
+C           LINE SEGMENT AND ITS RESPECTIVE TRIANGLE NUMBER
+C           ARE TO BE STORED AS THE (3*IL-2)ND, (3*IL-1)ST,
+C           AND (3*IL)TH ELEMENTS, IL=1,2,..., NL.
+C THE OTHER PARAMETERS ARE
+C     IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
+C           INTERNALLY AS A WORK AREA,
+C     IWP = INTEGER ARRAY OF DIMENSION NDP USED
+C           INTERNALLY AS A WORK AREA,
+C     WK  = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
+C           WORK AREA.
+C DECLARATION STATEMENTS
+C
+      SAVE
+C
+      INTEGER         CONXCH
+      COMMON /CONRA3/ IREC
+      DIMENSION       XD(*)      ,YD(*)      ,IPT(*)     ,IPL(*)     ,
+     1                IWL(*)     ,IWP(*)     ,WK(*)
+      DIMENSION       ITF(2)
+        CHARACTER*4  IP1C, IP2C
+        CHARACTER*64  ITEMP
+      DATA  RATIO/1.0E-6/, NREP/100/
+C
+C STATEMENT FUNCTIONS
+C
+      DSQF(U1,V1,U2,V2) = (U2-U1)**2+(V2-V1)**2
+      SIDE(U1,V1,U2,V2,U3,V3) = (V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
+C
+C PRELIMINARY PROCESSING
+C
+      NDPM1 = NDP-1
+C
+C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
+C
+      DSQMN = DSQF(XD(1),YD(1),XD(2),YD(2))
+      IPMN1 = 1
+      IPMN2 = 2
+      DO  140 IP1=1,NDPM1
+         X1 = XD(IP1)
+         Y1 = YD(IP1)
+         IP1P1 = IP1+1
+         DO  130 IP2=IP1P1,NDP
+            DSQI = DSQF(X1,Y1,XD(IP2),YD(IP2))
+            IF (DSQI .NE. 0.) GO TO  120
+C
+C  ERROR, IDENTICAL INPUT DATA POINTS
+C
+          ITEMP = ' CONTNG-IDENTICAL INPUT DATA POINTS FOUND
+     1 AT      AND      '
+C 
+C + NOAO - FTN internal writes rewritten as calls to encode for IRAF
+C
+C         WRITE(IP1C,'(I4)')IP1
+C         WRITE(IP2C,'(I4)')IP2
+          call encode (4, '(I4)', ip1c, ip1)
+	  call encode (4, '(I4)', ip2c, ip2)
+C - NOAO
+C
+            CALL SETER (ITEMP,1,1)
+            ITEMP(46:49) = IP1C
+            ITEMP(55:58) = IP2C
+            RETURN
+  120       IF (DSQI .GE. DSQMN) GO TO  130
+            DSQMN = DSQI
+            IPMN1 = IP1
+            IPMN2 = IP2
+  130    CONTINUE
+  140 CONTINUE
+      DSQ12 = DSQMN
+      XDMP = (XD(IPMN1)+XD(IPMN2))/2.0
+      YDMP = (YD(IPMN1)+YD(IPMN2))/2.0
+C
+C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
+C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
+C NUMBERS IN THE IWP ARRAY.
+C
+      JP1 = 2
+      DO  150 IP1=1,NDP
+         IF (IP1.EQ.IPMN1 .OR. IP1.EQ.IPMN2) GO TO  150
+         JP1 = JP1+1
+         IWP(JP1) = IP1
+         WK(JP1) = DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
+  150 CONTINUE
+      DO  170 JP1=3,NDPM1
+         DSQMN = WK(JP1)
+         JPMN = JP1
+         DO  160 JP2=JP1,NDP
+            IF (WK(JP2) .GE. DSQMN) GO TO  160
+            DSQMN = WK(JP2)
+            JPMN = JP2
+  160    CONTINUE
+         ITS = IWP(JP1)
+         IWP(JP1) = IWP(JPMN)
+         IWP(JPMN) = ITS
+         WK(JPMN) = WK(JP1)
+  170 CONTINUE
+C
+C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
+C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
+C
+      AR = DSQ12*RATIO
+      X1 = XD(IPMN1)
+      Y1 = YD(IPMN1)
+      DX21 = XD(IPMN2)-X1
+      DY21 = YD(IPMN2)-Y1
+      DO  180 JP=3,NDP
+         IP = IWP(JP)
+         IF (ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21) .GT. AR) GO TO  190
+  180 CONTINUE
+      CALL SETER (' CONTNG - ALL COLLINEAR DATA POINTS',1,1)
+  190 IF (JP .EQ. 3) GO TO  210
+      JPMX = JP
+      JP = JPMX+1
+      DO  200 JPC=4,JPMX
+         JP = JP-1
+         IWP(JP) = IWP(JP-1)
+  200 CONTINUE
+      IWP(3) = IP
+C
+C FORMS THE FIRST TRIANGLE.  STORES POINT NUMBERS OF THE VER-
+C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
+C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
+C THE IPL ARRAY.
+C
+  210 IP1 = IPMN1
+      IP2 = IPMN2
+      IP3 = IWP(3)
+      IF (SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3)) .GE.
+     1    0.0) GO TO  220
+      IP1 = IPMN2
+      IP2 = IPMN1
+  220 NT0 = 1
+      NTT3 = 3
+      IPT(1) = IP1
+      IPT(2) = IP2
+      IPT(3) = IP3
+      NL0 = 3
+      NLT3 = 9
+      IPL(1) = IP1
+      IPL(2) = IP2
+      IPL(3) = 1
+      IPL(4) = IP2
+      IPL(5) = IP3
+      IPL(6) = 1
+      IPL(7) = IP3
+      IPL(8) = IP1
+      IPL(9) = 1
+C
+C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
+C
+      DO  400 JP1=4,NDP
+         IP1 = IWP(JP1)
+         X1 = XD(IP1)
+         Y1 = YD(IP1)
+C
+C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
+C
+         IP2 = IPL(1)
+         JPMN = 1
+         DXMN = XD(IP2)-X1
+         DYMN = YD(IP2)-Y1
+         DSQMN = DXMN**2+DYMN**2
+         ARMN = DSQMN*RATIO
+         JPMX = 1
+         DXMX = DXMN
+         DYMX = DYMN
+         DSQMX = DSQMN
+         ARMX = ARMN
+         DO  240 JP2=2,NL0
+            IP2 = IPL(3*JP2-2)
+            DX = XD(IP2)-X1
+            DY = YD(IP2)-Y1
+            AR = DY*DXMN-DX*DYMN
+            IF (AR .GT. ARMN) GO TO  230
+            DSQI = DX**2+DY**2
+            IF (AR.GE.(-ARMN) .AND. DSQI.GE.DSQMN) GO TO  230
+            JPMN = JP2
+            DXMN = DX
+            DYMN = DY
+            DSQMN = DSQI
+            ARMN = DSQMN*RATIO
+  230       AR = DY*DXMX-DX*DYMX
+            IF (AR .LT. (-ARMX)) GO TO  240
+            DSQI = DX**2+DY**2
+            IF (AR.LE.ARMX .AND. DSQI.GE.DSQMX) GO TO  240
+            JPMX = JP2
+            DXMX = DX
+            DYMX = DY
+            DSQMX = DSQI
+            ARMX = DSQMX*RATIO
+  240    CONTINUE
+         IF (JPMX .LT. JPMN) JPMX = JPMX+NL0
+         NSH = JPMN-1
+         IF (NSH .LE. 0) GO TO  270
+C
+C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
+C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
+C
+         NSHT3 = NSH*3
+         DO  250 JP2T3=3,NSHT3,3
+            JP3T3 = JP2T3+NLT3
+            IPL(JP3T3-2) = IPL(JP2T3-2)
+            IPL(JP3T3-1) = IPL(JP2T3-1)
+            IPL(JP3T3) = IPL(JP2T3)
+  250    CONTINUE
+         DO  260 JP2T3=3,NLT3,3
+            JP3T3 = JP2T3+NSHT3
+            IPL(JP2T3-2) = IPL(JP3T3-2)
+            IPL(JP2T3-1) = IPL(JP3T3-1)
+            IPL(JP2T3) = IPL(JP3T3)
+  260    CONTINUE
+         JPMX = JPMX-NSH
+C
+C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
+C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
+C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
+C
+  270    JWL = 0
+         DO  310 JP2=JPMX,NL0
+            JP2T3 = JP2*3
+            IPL1 = IPL(JP2T3-2)
+            IPL2 = IPL(JP2T3-1)
+            IT = IPL(JP2T3)
+C
+C - - ADDS A TRIANGLE TO THE IPT ARRAY.
+C
+            NT0 = NT0+1
+            NTT3 = NTT3+3
+            IPT(NTT3-2) = IPL2
+            IPT(NTT3-1) = IPL1
+            IPT(NTT3) = IP1
+C
+C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
+C
+            IF (JP2 .NE. JPMX) GO TO  280
+            IPL(JP2T3-1) = IP1
+            IPL(JP2T3) = NT0
+  280       IF (JP2 .NE. NL0) GO TO  290
+            NLN = JPMX+1
+            NLNT3 = NLN*3
+            IPL(NLNT3-2) = IP1
+            IPL(NLNT3-1) = IPL(1)
+            IPL(NLNT3) = NT0
+C
+C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
+C - - LINE SEGMENTS.
+C
+  290       ITT3 = IT*3
+            IPTI = IPT(ITT3-2)
+            IF (IPTI.NE.IPL1 .AND. IPTI.NE.IPL2) GO TO  300
+            IPTI = IPT(ITT3-1)
+            IF (IPTI.NE.IPL1 .AND. IPTI.NE.IPL2) GO TO  300
+            IPTI = IPT(ITT3)
+C
+C - - CHECKS IF THE EXCHANGE IS NECESSARY.
+C
+  300       IF (CONXCH(XD,YD,IP1,IPTI,IPL1,IPL2) .EQ. 0) GO TO  310
+C
+C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
+C
+            IPT(ITT3-2) = IPTI
+            IPT(ITT3-1) = IPL1
+            IPT(ITT3) = IP1
+            IPT(NTT3-1) = IPTI
+            IF (JP2 .EQ. JPMX) IPL(JP2T3) = IT
+            IF (JP2.EQ.NL0 .AND. IPL(3).EQ.IT) IPL(3) = NT0
+C
+C - - SETS FLAGS IN THE IWL ARRAY.
+C
+            JWL = JWL+4
+            IWL(JWL-3) = IPL1
+            IWL(JWL-2) = IPTI
+            IWL(JWL-1) = IPTI
+            IWL(JWL) = IPL2
+  310    CONTINUE
+         NL0 = NLN
+         NLT3 = NLNT3
+         NLF = JWL/2
+         IF (NLF .EQ. 0) GO TO  400
+C
+C - IMPROVES TRIANGULATION.
+C
+         NTT3P3 = NTT3+3
+         DO  390 IREP=1,NREP
+            DO  370 ILF=1,NLF
+               ILFT2 = ILF*2
+               IPL1 = IWL(ILFT2-1)
+               IPL2 = IWL(ILFT2)
+C
+C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
+C - - THE FLAGGED LINE SEGMENT.
+C
+               NTF = 0
+               DO  320 ITT3R=3,NTT3,3
+                  ITT3 = NTT3P3-ITT3R
+                  IPT1 = IPT(ITT3-2)
+                  IPT2 = IPT(ITT3-1)
+                  IPT3 = IPT(ITT3)
+                  IF (IPL1.NE.IPT1 .AND. IPL1.NE.IPT2 .AND.
+     1                IPL1.NE.IPT3) GO TO  320
+                  IF (IPL2.NE.IPT1 .AND. IPL2.NE.IPT2 .AND.
+     1                IPL2.NE.IPT3) GO TO  320
+                  NTF = NTF+1
+                  ITF(NTF) = ITT3/3
+                  IF (NTF .EQ. 2) GO TO  330
+  320          CONTINUE
+               IF (NTF .LT. 2) GO TO  370
+C
+C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
+C - - ON THE LINE SEGMENT.
+C
+  330          IT1T3 = ITF(1)*3
+               IPTI1 = IPT(IT1T3-2)
+               IF (IPTI1.NE.IPL1 .AND. IPTI1.NE.IPL2) GO TO  340
+               IPTI1 = IPT(IT1T3-1)
+               IF (IPTI1.NE.IPL1 .AND. IPTI1.NE.IPL2) GO TO  340
+               IPTI1 = IPT(IT1T3)
+  340          IT2T3 = ITF(2)*3
+               IPTI2 = IPT(IT2T3-2)
+               IF (IPTI2.NE.IPL1 .AND. IPTI2.NE.IPL2) GO TO  350
+               IPTI2 = IPT(IT2T3-1)
+               IF (IPTI2.NE.IPL1 .AND. IPTI2.NE.IPL2) GO TO  350
+               IPTI2 = IPT(IT2T3)
+C
+C - - CHECKS IF THE EXCHANGE IS NECESSARY.
+C
+  350          IF (CONXCH(XD,YD,IPTI1,IPTI2,IPL1,IPL2) .EQ. 0)
+     1             GO TO  370
+C
+C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
+C
+               IPT(IT1T3-2) = IPTI1
+               IPT(IT1T3-1) = IPTI2
+               IPT(IT1T3) = IPL1
+               IPT(IT2T3-2) = IPTI2
+               IPT(IT2T3-1) = IPTI1
+               IPT(IT2T3) = IPL2
+C
+C - - SETS NEW FLAGS.
+C
+               JWL = JWL+8
+               IWL(JWL-7) = IPL1
+               IWL(JWL-6) = IPTI1
+               IWL(JWL-5) = IPTI1
+               IWL(JWL-4) = IPL2
+               IWL(JWL-3) = IPL2
+               IWL(JWL-2) = IPTI2
+               IWL(JWL-1) = IPTI2
+               IWL(JWL) = IPL1
+               DO  360 JLT3=3,NLT3,3
+                  IPLJ1 = IPL(JLT3-2)
+                  IPLJ2 = IPL(JLT3-1)
+                  IF ((IPLJ1.EQ.IPL1 .AND. IPLJ2.EQ.IPTI2) .OR.
+     1                (IPLJ2.EQ.IPL1 .AND. IPLJ1.EQ.IPTI2))
+     2                IPL(JLT3) = ITF(1)
+                  IF ((IPLJ1.EQ.IPL2 .AND. IPLJ2.EQ.IPTI1) .OR.
+     1                (IPLJ2.EQ.IPL2 .AND. IPLJ1.EQ.IPTI1))
+     2                IPL(JLT3) = ITF(2)
+  360          CONTINUE
+  370       CONTINUE
+            NLFC = NLF
+            NLF = JWL/2
+            IF (NLF .EQ. NLFC) GO TO  400
+C
+C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
+C
+            JWL = 0
+            JWL1MN = (NLFC+1)*2
+            NLFT2 = NLF*2
+            DO  380 JWL1=JWL1MN,NLFT2,2
+               JWL = JWL+2
+               IWL(JWL-1) = IWL(JWL1-1)
+               IWL(JWL) = IWL(JWL1)
+  380       CONTINUE
+            NLF = JWL/2
+  390    CONTINUE
+  400 CONTINUE
+C
+C REARRANGE THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
+C ARE LISTED COUNTER-CLOCKWISE.
+C
+      DO  410 ITT3=3,NTT3,3
+         IP1 = IPT(ITT3-2)
+         IP2 = IPT(ITT3-1)
+         IP3 = IPT(ITT3)
+         IF (SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3)) .GE.
+     1       0.0) GO TO  410
+         IPT(ITT3-2) = IP2
+         IPT(ITT3-1) = IP1
+  410 CONTINUE
+      NT = NT0
+      NL = NL0
+      RETURN
+      END