The program ACORD is taken, with a slight revision, from the article
"ACORD: AUTOMATIC COUNTOURING OF RAW DATA, Computers & Geosciences,
vol. 8, no. 1, p. 97-101, 1982", by D.F. Watson.  It was written in
1979 in an ancient version of FORTRAN to run on the Control Data
Cyber main-frame.  It accepts a set of bivariate data {x, y, f(x,y)}
and finds the Delaunay triangulation.  Then the trace of each contour
level on each triangle is output, that is, piecewise linear triangle-
based interpolation.  This program is notable only for the concise
triangulation algorithm, which is better than O(N**2).  The arrays
will take a maximum of 1200 data plus allowing for three control
vertices.  The data set that follows, included for reference, is
from "Davis, J.C., 1986, Statistics and data analysis in geology,
2ed., Wiley, p.362" and is prefaced by the number of data, the number
of contours, the input format, and the contour levels, in that order.
Most data sets will require perturbation to avoid degenerate
configuations so a small random value is randomly added or subtracted
to or from each x and y coordinate.
Discussion or (gentle) criticsm is invited.

Dave Watson                   Internet: watson@maths.uwa.oz.au
Department of Mathematics
The University of Western Australia       Tel: (61 9) 380 3359
Nedlands, WA 6009  Australia.             FAX: (61 9) 380 1028
--------------------------- snipity-snip -------------------------------
      PROGRAM ACORD
C.....ACORD: Automatic Contouring Of Raw Data - by D.F.Watson
C.....READ DATA POINTS AND FORM ALL 3-TUPLES S.T. NO OTHER POINT
C     LIES WITHIN THAT 3-TUPLE'S CIRCUMCIRCLE
C.....PNT HOLDS THE DATA POINTS AND TETR CARRIES THE CIRCUMCIRCLE
C     CENTRE AND RADIUS SQUARED FOR EACH 3-TUPLE
C.....ITETR HOLDS THE DATA POINT INDICES IN INPUT ORDER OF EACH 3-TUPLE
C.....ISTACK IS A LAST-IN-FIRST-OUT STACK OF INDICES OF VACANT 3-TUPLES
C.....KTETR IS A TEMPORARY LIST OF EDGES OF DELETED 3-TUPLES
C.....ID IS POINTER TO ISTACK, AND JT IS POINTER TO TETR AND ITETR
      REAL PNT(1203,3), TETR(2401,3), XPNT(3,3), CONT(50), DET(2,3)
      INTEGER ITETR(2401,3), ISTACK(2401), KTETR(50,2), ITEMP(3,2)
      CHARACTER*32 FORM
      DATA ITEMP, XPNT/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./
      DATA NTETR,ID,XMIN,YMIN,XMAX,YMAX/1,2,1.E37,1.E37,-1.E37,-1.E37/
C......................................................................|
C.....OPEN FILE FOR DATA INPUT
      OPEN(UNIT=7, FILE="jdavis.dat", ERR=99)
C.....READ NO. OF DATA POINTS, NO. OF CONTOURS, DATA FORMAT
      READ(7,200) NDATA, LC, FORM
  200 FORMAT(I8,I5,A13)
C.....READ THE CONTOUR VALUES
      READ(7,FORM) (CONT(I),I=1,LC)
      RNUM = RAND(1)
C.....READ THE DATA
      DO 5 I = 1,NDATA
        READ(7,FORM) PNT(I,1), PNT(I,2), PNT(I,3)
        RNUM = RAND(RNUM)
        PNT(I,1) = PNT(I,1) + 0.00001 * (RNUM - 0.5)
        RNUM = RAND(RNUM)
        PNT(I,2) = PNT(I,2) + 0.00001 * (RNUM - 0.5)
        IF(PNT(I,1) .GT. XMAX) XMAX = PNT(I,1)
        IF(PNT(I,1) .LT. XMIN) XMIN = PNT(I,1)
        IF(PNT(I,2) .GT. YMAX) YMAX = PNT(I,2)
    5   IF(PNT(I,2) .LT. YMIN) YMIN = PNT(I,2)
      CLOSE(7)
      DATAX = XMAX - XMIN
      Y = YMAX - YMIN
      IF(Y .GT. DATAX) DATAX = Y
C.....NORMALIZE THE DATA
      DO 7 I = 1,NDATA
        PNT(I,1) = (PNT(I,1) - XMIN)/DATAX
    7   PNT(I,2) = (PNT(I,2) - YMIN)/DATAX
      DO 8 I = 1,3
        ITETR(1,I) = I + NDATA
        TETR(1,I) = XPNT(I,3)
        DO 8 J = 1,2
    8     PNT(I + NDATA,J) = XPNT(I,J)
      J = NDATA * 2 + 1
      DO 9 I = 2,J
    9   ISTACK(I) = I
C......................................................................|
      DO 50 NUC = 1,NDATA
        KM = 0
C.......LOOP THRU THE ESTABLISHED 3-TUPLES
        DO 30 JT = 1,NTETR
C.........TEST IF NEW DATA POINT IS WITHIN THE JT CIRCUMCIRCLE
          DSQ = TETR(JT,3) - (PNT(NUC,1) - TETR(JT,1))**2
          IF(DSQ .LT. 0.0) GO TO 30
          DSQ = DSQ - (PNT(NUC,2) - TETR(JT,2))**2
          IF(DSQ .LT. 0.0) GO TO 30
C.........DELETE THIS 3-TUPLE BUT SAVE ITS EDGES
          ID = ID - 1
          ISTACK(ID) = JT
C.........ADD EDGES TO KTETR BUT DELETE IF ALREADY LISTED
          DO 28 I = 1,3
            L1 = ITEMP(I,1)
            L2 = ITEMP(I,2)
            IF(KM .LE. 0) GO TO 26
            KMT = KM
            DO 24 J = 1,KMT
              IF(ITETR(JT,L1) .NE. KTETR(J,1)) GO TO 24
              IF(ITETR(JT,L2) .NE. KTETR(J,2)) GO TO 24
              KM = KM - 1
              IF(J .GT. KM) GO TO 28
              DO 20 K = J,KM
                K1 = K + 1
                DO 20 L = 1,2
   20             KTETR(K,L) = KTETR(K1,L)
                GO TO 28
   24         CONTINUE
   26       KM = KM + 1
            KTETR(KM,1) = ITETR(JT,L1)
            KTETR(KM,2) = ITETR(JT,L2)
   28       CONTINUE
   30     CONTINUE
C
C.......FORM NEW 3-TUPLE3
        DO 48 I = 1,KM
          KT = ISTACK(ID)
          ID = ID + 1
C.........CALCULATE THE CIRCUMCIRCLE CENTER AND RADIUS
C         SQUARED OF POINTS KTETR(I,*) AND PLACE IN TETR(KT,*)
          DO 44 JZ = 1,2
            I2 = KTETR(I,JZ)
            DET(JZ,1) = PNT(I2,1) - PNT(NUC,1)
            DET(JZ,2) = PNT(I2,2) - PNT(NUC,2)
   44       DET(JZ,3) = DET(JZ,1) * (PNT(I2,1) + PNT(NUC,1))/2.0
     +                + DET(JZ,2) * (PNT(I2,2) + PNT(NUC,2))/2.0
          DD = DET(1,1) * DET(2,2) - DET(1,2) * DET(2,1)
          TETR(KT,1) = (DET(1,3) * DET(2,2) - DET(2,3) * DET(1,2))/DD
          TETR(KT,2) = (DET(1,1) * DET(2,3) - DET(2,1) * DET(1,3))/DD
          TETR(KT,3) = (PNT(NUC,1) - TETR(KT,1))**2
     +               + (PNT(NUC,2) - TETR(KT,2))**2
          ITETR(KT,1) = KTETR(I,1)
          ITETR(KT,2) = KTETR(I,2)
   48     ITETR(KT,3) = NUC
   50   NTETR = NTETR + 2
C......................................................................|
C.....SET OUTPUT SIZE FACTOR
      SIZE = 1.0
      OPEN(UNIT=8, FILE="contsegm.dat", ERR=99)
C
      DO 90 JT = 1,NTETR
        IF(ITETR(JT,1) .GT. NDATA .OR. TETR(JT,3) .GT. 1) GO TO 90
C.......FIND CONTOUR INTERSECTIONS
        ABIT = 0.0
        IF(PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,2),3) .OR.
     +     PNT(ITETR(JT,1),3) .EQ. PNT(ITETR(JT,3),3) .OR.
     +     PNT(ITETR(JT,2),3) .EQ. PNT(ITETR(JT,3),3)) ABIT = 1.E-10
        TOP = AMAX1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3),
     +              PNT(ITETR(JT,3),3))
        BOT = AMIN1(PNT(ITETR(JT,1),3),PNT(ITETR(JT,2),3),
     +              PNT(ITETR(JT,3),3))
        DO 87 JC = 1,LC
          IF(CONT(JC) .GT. TOP .OR. CONT(JC) .LT. BOT) GO TO 87
          CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,2),3)
     +                   - PNT(ITETR(JT,1),3) + ABIT)
          IF(CZ .LE. 0.0 .OR. CZ .GE. 1.0) GO TO 83
          X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,2),1)
     +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
          Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,2),2)
     +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
   82     CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3)
     +                   - PNT(ITETR(JT,1),3) + ABIT)
          IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 84
          X2 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1)
     +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
          Y2 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2)
     +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
          GO TO 85
   83     CZ = (CONT(JC) - PNT(ITETR(JT,1),3))/(PNT(ITETR(JT,3),3)
     +                   - PNT(ITETR(JT,1),3) + ABIT)
          IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87
          X1 = (PNT(ITETR(JT,1),1) + (PNT(ITETR(JT,3),1)
     +       -  PNT(ITETR(JT,1),1)) * CZ) * SIZE
          Y1 = (PNT(ITETR(JT,1),2) + (PNT(ITETR(JT,3),2)
     +       -  PNT(ITETR(JT,1),2)) * CZ) * SIZE
   84     CZ = (CONT(JC) - PNT(ITETR(JT,2),3))/(PNT(ITETR(JT,3),3)
     +                   - PNT(ITETR(JT,2),3) + ABIT)
          IF(CZ .LT. 0.0 .OR. CZ .GT. 1.0) GO TO 87
          X2 = (PNT(ITETR(JT,2),1) + (PNT(ITETR(JT,3),1)
     +       -  PNT(ITETR(JT,2),1)) * CZ) * SIZE
          Y2 = (PNT(ITETR(JT,2),2) + (PNT(ITETR(JT,3),2)
     +       -  PNT(ITETR(JT,2),2)) * CZ) * SIZE
   85     CONTINUE
C         WRITE A CONTOUR SEGMENT TO FILE
C         "3" means 'pen-up' and "2" means 'pen-down"
          WRITE(8,400) X1, Y1, "3", X2, Y2, "2"
  400     FORMAT(2F10.4,A5,2F10.4,A5)
   87     CONTINUE
   90   CONTINUE
      CLOSE(8)
   99 STOP
      END
------------------------------- snipity-snip  ------------------------------
      52   11 (3F5.0)
     700  725  750
     775  800  825
     850  875  900
     925  950
     0.3  6.1  870
     1.4  6.2  793
     2.4  6.1  755
     3.6  6.2  690
     5.7  6.2  800
     1.6  5.2  800
     2.9  5.1  730
     3.4  5.3  728
     3.4  5.7  710
     4.8  5.6  780
     5.3  5.0  804
     6.2  5.2  855
     0.2  4.3  830
     0.9  4.2  813
     2.3  4.8  762
     2.5  4.5  765
     3.0  4.5  740
     3.5  4.5  765
     4.1  4.6  760
     4.9  4.2  790
     6.3  4.3  820
     0.9  3.2  855
     1.7  3.8  812
     2.4  3.8  773
     3.7  3.5  812
     4.5  3.2  827
     5.2  3.2  805
     6.3  3.4  840
     0.3  2.4  890
     2.0  2.7  820
     3.8  2.3  873
     6.3  2.2  875
     0.6  1.7  873
     1.5  1.8  865
     2.1  1.8  841
     2.1  1.1  862
     3.1  1.1  908
     4.5  1.8  855
     5.5  1.7  850
     5.7  1.0  882
     6.2  1.0  910
     0.4  0.5  940
     1.4  0.6  915
     1.4  0.1  890
     2.1  0.7  880
     2.3  0.3  870
     3.1  0.0  880
     4.1  0.8  960
     5.4  0.4  890
     6.0  0.1  860
     5.7  3.0  830
     3.6  6.0  705
-------------------------------------------------