C zcont.f
C
C VIS-5D version 3.0
C
C
C vis5d program for visualizing five dimensional gridded data sets
C Copyright (C) 1990, 1991, 1992  Bill Hibbard, Brian Paul, Dave Santek,
C and Andre Battaiola.
C 
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 1, or (at your option)
C any later version.
C 
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C 
C You should have received a copy of the GNU General Public License
C along with this program; if not, write to the Free Software
C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C
C
C Compute contour lines over a 2-D grid.
C
C

      SUBROUTINE ZCONT WXYZ (G, NR,NC, CINT, MARK,
     *  VX1,VY1, MAXV1, NUMV1,
     *  VX2,VY2, MAXV2, NUMV2,
     *  VX3,VY3, MAXV3, NUMV3 )
C
C Compute contour lines in a 2-D grid at a given interval.
C
C G = (input) grid to be contoured
C NR, NC = (input) dimensions of G
C CINT = (input) contour interval
C MARK = (input) work array for marking grid squares for labels
C        an array of size nr x nc
C VX1,VY1 = (output) array of contour line segment vertices
C MAXV1 = (input) max size of vx1 and vy1 arrays
C NUMV1 = (output) number of vertices computed and stored in VX1 and VY1
C VX2,VY2 = (output) array of 'hidden' line segment vertices
C MAXV2 = (input) max size of VX2 and VY2 arrays
C NUMV2 = (output) number of vertices computed and stored in VX2 and VY2
c VX3,VY3 = (output) array of label line segment vertices
C MAXV3 = (input) max size of VX3 and VY3 arrays
C NUMV3 = (output) number of vertices computed and stord in VX3 and VY3
C
C
      REAL*4 G(NR,NC)
      INTEGER*4 MARK(NR,NC)
      REAL*4 VX1(MAXV1),VY1(MAXV1), VX2(MAXV2),VY2(MAXV2),
     *       VX3(MAXV3),VY3(MAXV3)
C
      PARAMETER (NLINES=20000)
      INTEGER*4 IPNT(NLINES)
      REAL*4 VX(NLINES),VY(NLINES)
      COMMON/TPOLY WXYZ/ILAST,ILIST,VX,VY,IPNT
C
C     all returned vertices will be scaled to this region:
      PARAMETER (XMIN = -1.0)
      PARAMETER (XMAX = 1.0)
      PARAMETER (YMIN = -1.0)
      PARAMETER (YMAX = 1.0)
C
C Note:  The input grid is assumed to be in row-major form (C-style).
C        Since FORTRAN uses column-major 2-D arrays, this subroutine
C        associates rows with the Y-axis and columns with the X-axis.
C
      ILIST=0
      ILAST=NLINES+1
      NRM=NR-1
      NCM=NC-1
C
      XD=(XMAX-XMIN)/(NR-1)
      YD=(YMAX-YMIN)/(NC-1)
C
C     CONSTRUCT SQUARE BOUNDARY LINES
      VX(1)=XMAX
      VY(1)=YMIN
      VX(2)=XMIN
      VY(2)=YMIN
      VX(3)=XMIN
      VY(3)=YMIN
      VX(4)=XMIN
      VY(4)=YMAX
      VX(5)=XMIN
      VY(5)=YMAX
      VX(6)=XMAX
      VY(6)=YMAX
      VX(7)=XMAX
      VY(7)=YMAX
      VX(8)=XMAX
      VY(8)=YMIN
      ILIST=8
C
C SET UP MARK ARRAY
C MARK= 0 IF AVAIL FOR LABEL CENTER,
C       2 IF IN LABEL, AND
C       1 IF NOT AVAILABLE AND NOT IN LABEL
C
C LR AND LC GIVE LABEL SIZE IN GRID BOXES
C LRR AND LCC GIVE UNAVAILABLE RADIUS
      LR=1+(NR-2)/50
      LC=1+(NC-2)/10
      LC2=LC/2
      LR2=LR/2
      LRR=1+(NR-2)/8
      LCC=1+(NC-2)/8
C
      DO 10 IR=1,NR
         DO 10 IC=1,NC
            MARK(IR,IC)=0
10    CONTINUE
      DO 12 IC=1,NC
         DO 12 IR=1,LR
           MARK(IR,IC)=1
           MARK(NR-IR,IC)=1
12    CONTINUE
      DO 14 IR=1,NR
         DO 14 IC=1,LC
            MARK(IR,IC)=1
            MARK(IR,NC-IC)=1
14    CONTINUE
      NRNCP=NR*NC+1
      DO 16 I=1,NRNCP
         IPNT(I)=-1
16    CONTINUE
C
C NOW DRAW THE CONTOURS
      IP=0
      DO 2000 IR=1,NRM
         XX=XD*(IR-1.0)+XMIN
         DO 2000 IC=1,NCM
            YY=YD*(IC-1.0)+YMIN
            IP=IP+1
            IPNT(IP)=ILIST
            GA=G(IR,IC)
            IF(GA .GT. 1.E30) GO TO 2000
            GB=G(IR+1,IC)
            IF(GB .GT. 1.E30) GO TO 2000
            GC=G(IR,IC+1)
            IF(GC .GT. 1.E30) GO TO 2000
            GD=G(IR+1,IC+1)
            IF(GD .GT. 1.E30) GO TO 2000
            GV=(GA+GB+GC+GD)/4.0
            GN=AMIN1(GA,GB,GC,GD)
            GX=AMAX1(GA,GB,GC,GD)
C
            CLOW=CINT*(NINT(GN/CINT)-1)
1           IF(CLOW .GE. GN) GO TO 2
            CLOW=CLOW+CINT
            GO TO 1
C
2           CHI=CINT*(NINT(GX/CINT)+1)
3           IF(CHI .LE. GX) GO TO 4
            CHI=CHI-CINT
            GO TO 3
C
4           NUMC=1+NINT((CHI-CLOW)/CINT)
C
            GG=CLOW
            DO 1800 IL=1,NUMC
               IF(GG .LT. GN) GO TO 1800
               IF(GG .GT. GX) GO TO 1900
               II=0
               IF(GG .GT. GA) II=1
               IF(GG .GT. GB) II=II+2
               IF(GG .GT. GC) II=II+4
               IF(GG .GT. GD) II=II+8
               IF(II .GT. 7) II=15-II
               IF(II .LE. 0) GO TO 1800
C           DO LABEL HERE
               IF(MARK(IR,IC) .NE. 0) GO TO 90
C           BOX TO AVOID
               KC=IC-LC2-LCC
               KR=IR-LR2-LRR
               MC=KC+2*LCC+LC-1
               MR=KR+2*LRR+LR-1
               DO 30 JC=KC,MC
                  IF(JC .LT. 1 .OR. JC .GT. NC) GO TO 30
                  DO 29 JR=KR,MR
                     IF (JR .LT. 1 .OR. JR .GT. NR) GO TO 29
                     IF (MARK(JR,JC) .NE. 2) MARK(JR,JC)=1
29                CONTINUE
30             CONTINUE
C           BOX TO HOLD LABEL
               KC=IC-LC2
               KR=IR-LR2
               MC=KC+LC-1
               MR=KR+LR-1
               DO 35 JC=KC,MC
                  IF(JC .LT. 1 .OR. JC .GT. NC) GO TO 35
                  DO 34 JR=KR,MR
                     IF(JR .LT. 1 .OR. JR .GT. NR) GO TO 34
                     MARK(JR,JC)=2
34                CONTINUE
35             CONTINUE
               XK=XD*(KR-1.0)+XMIN
               YK=YD*(KC-1.0)+YMIN
               XM=XD*MR+XMIN
               YM=YD*MC+YMIN
               HH=CLOW+CINT*FLOAT((NUMC-1)/2)
C
               CALL PLTDG WXYZ (HH,XK,YK,XM,YM)
C
90             IF(ILIST+100 .GE. ILAST) GO TO 900
               GOTO (100,200,300,400,500,600,700) II
C
100            GBA=GB-GA
               GCA=GC-GA
               ILIST=ILIST+1
               IF(ABS(GBA) .LT. 0.0000001) THEN
                  VX(ILIST)=XX
               ELSE
                  VX(ILIST)=XX+XD*(GG-GA)/GBA
               ENDIF
               VY(ILIST)=YY
               ILIST=ILIST+1
               IF(ABS(GCA) .LT. 0.0000001) THEN
                  VY(ILIST)=YY
               ELSE
                  VY(ILIST)=YY+YD*(GG-GA)/GCA
               ENDIF
               VX(ILIST)=XX
               GO TO 1800
C
200            GBA=GB-GA
               GDB=GD-GB
               ILIST=ILIST+1
               IF(ABS(GBA) .LT. 0.0000001) THEN
                 VX(ILIST)=XX
               ELSE
                 VX(ILIST)=XX+XD*(GG-GA)/GBA
               ENDIF
               VY(ILIST)=YY
               ILIST=ILIST+1
               IF(ABS(GDB) .LT. 0.0000001) THEN
                 VY(ILIST)=YY
               ELSE
                 VY(ILIST)=YY+YD*(GG-GB)/GDB
               ENDIF
               VX(ILIST)=XX+XD
               GO TO 1800
C
300            GCA=GC-GA
               GDB=GD-GB
               ILIST=ILIST+1
               IF(ABS(GCA) .LT. 0.0000001) THEN
                 VY(ILIST)=YY
               ELSE
                 VY(ILIST)=YY+YD*(GG-GA)/GCA
               ENDIF
               VX(ILIST)=XX
               ILIST=ILIST+1
               IF(ABS(GDB) .LT. 0.0000001) THEN
                 VY(ILIST)=YY
               ELSE
                 VY(ILIST)=YY+YD*(GG-GB)/GDB
               ENDIF
               VX(ILIST)=XX+XD
               GO TO 1800
C
400            GCA=GC-GA
               GDC=GD-GC
               ILIST=ILIST+1
               IF(ABS(GCA) .LT. 0.0000001) THEN
                 VY(ILIST)=YY
               ELSE
                 VY(ILIST)=YY+YD*(GG-GA)/GCA
               ENDIF
               VX(ILIST)=XX
               ILIST=ILIST+1
               IF(ABS(GDC) .LT. 0.0000001) THEN
                 VX(ILIST)=XX
               ELSE
                 VX(ILIST)=XX+XD*(GG-GC)/GDC
               ENDIF
               VY(ILIST)=YY+YD
               GO TO 1800
C
500            GBA=GB-GA
               GDC=GD-GC
               ILIST=ILIST+1
               IF(ABS(GBA) .LT. 0.0000001) THEN
                 VX(ILIST)=XX
               ELSE
                 VX(ILIST)=XX+XD*(GG-GA)/GBA
               ENDIF
               VY(ILIST)=YY
               ILIST=ILIST+1
               IF(ABS(GDC) .LT. 0.0000001) THEN
                 VX(ILIST)=XX
               ELSE
                 VX(ILIST)=XX+XD*(GG-GC)/GDC
               ENDIF
               VY(ILIST)=YY+YD
               GO TO 1800
C
600            GBA=GB-GA
               GDC=GD-GC
               GCA=GC-GA
               GDB=GD-GB
               ILIST=ILIST+1
               IF(ABS(GBA) .LT. 0.0000001) THEN
                 VX(ILIST)=XX
               ELSE
                 VX(ILIST)=XX+XD*(GG-GA)/GBA
               ENDIF
               VY(ILIST)=YY
C HERE'S A BRAIN TEASER
               IF(GG .GT. GV .XOR. GA .LT. GB) THEN
                  ILIST=ILIST+1
                  IF(ABS(GCA) .LT. 0.0000001) THEN
                     VY(ILIST)=YY
                  ELSE
                     VY(ILIST)=YY+YD*(GG-GA)/GCA
                  ENDIF
                  VX(ILIST)=XX
                  ILIST=ILIST+1
                  IF(ABS(GDB) .LT. 0.0000001) THEN
                     VY(ILIST)=YY
                  ELSE
                     VY(ILIST)=YY+YD*(GG-GB)/GDB
                  ENDIF
                  VX(ILIST)=XX+XD
               ELSE
                  ILIST=ILIST+1
                  IF(ABS(GDB) .LT. 0.0000001) THEN
                     VY(ILIST)=YY
                  ELSE
                     VY(ILIST)=YY+YD*(GG-GB)/GDB
                  ENDIF
                  VX(ILIST)=XX+XD
                  ILIST=ILIST+1
                  IF(ABS(GCA) .LT. 0.0000001) THEN
                     VY(ILIST)=YY
                  ELSE
                     VY(ILIST)=YY+YD*(GG-GA)/GCA
                  ENDIF
                  VX(ILIST)=XX
               ENDIF
               ILIST=ILIST+1
               IF(ABS(GDC) .LT. 0.0000001) THEN
                  VX(ILIST)=XX
               ELSE
                  VX(ILIST)=XX+XD*(GG-GC)/GDC
               ENDIF
               VY(ILIST)=YY+YD
               GO TO 1800
C
700            GDB=GD-GB
               GDC=GD-GC
               ILIST=ILIST+1
               IF(ABS(GDB) .LT. 0.0000001) THEN
                  VY(ILIST)=YY
               ELSE
                  VY(ILIST)=YY+YD*(GG-GB)/GDB
               ENDIF
               VX(ILIST)=XX+XD
               ILIST=ILIST+1
               IF(ABS(GDC) .LT. 0.0000001) THEN
                  VX(ILIST)=XX
               ELSE
                  VX(ILIST)=XX+XD*(GG-GC)/GDC
               ENDIF
               VY(ILIST)=YY+YD
               GO TO 1800
C
1800           GG=GG+CINT
1900     CONTINUE
2000  CONTINUE
C
      IP=IP+1
      IPNT(IP)=ILIST
C
C     'unravel' the VX,VY, IPNT arrays into the VX1,VY1,VX2,VY2,VX3,VY3
C     vertex arrays.
C
C     copy vertices of square to VX1,VY1 arrays
      DO 3015 I=1,8
         VX1(I) = VX(I)
         VY1(I) = VY(I)
3015  CONTINUE
C
C     Copy contour and 'hidden' line segment vertices to VX1,VY1,
C     and VX2,VY2 arrays.
      NUMV1=8
      NUMV2=0
      IP=0
      DO 3040 IR=1,NRM
         DO 3040 IC=1,NCM
            IP=IP+1
            IQ=IPNT(IP)
            LEN=IPNT(IP+1)-IQ
            IF (LEN .LE. 0) GO TO 3040
            IF (MARK(IR,IC) .EQ. 2) THEN
               DO 3020 I=1,LEN
                  IF (NUMV2 .EQ. MAXV2) GOTO 3040
                  NUMV2 = NUMV2 + 1
                  VX2(NUMV2) = VX(IQ+I)
                  VY2(NUMV2) = VY(IQ+I)
3020           CONTINUE
            ELSE
               DO 3030 I=1,LEN
                  IF (NUMV1 .EQ. MAXV1) GOTO 3040
                  NUMV1 = NUMV1 + 1
                  VX1(NUMV1) = VX(IQ+I)
                  VY1(NUMV1) = VY(IQ+I)
3030           CONTINUE
            ENDIF
3040  CONTINUE
C
C     Copy label line segment vertices to VX3,VY3 arrays.
      NP=NLINES+1
      NUMV3=NLINES-ILAST+1
      IF (NUMV3 .GT. MAXV3) NUMV3 = MAXV3
      DO 3050 I=1,NUMV3
         VX3(I) = VX(NP-I)
         VY3(I) = VY(NP-I)
3050  CONTINUE
      RETURN
C
C     In case of error:
900   NUMV1 = 0
      NUMV2 = 0
      NUMV3 = 0
      RETURN
      END
C
C
C
C
      SUBROUTINE PLTDG WXYZ (GG,XK,YK,XM,YM)
C     Draw the value 'GG' as a set of line segments in the region
C     bounded by (XK,YK) to (XM,YM).  The segment vertices will
C     by stored in the VX,VY arrays starting at the end.
C     Note: ROWS = XK TO XM, COLS = YK TO YM
C
      PARAMETER (NLINES=20000)
      INTEGER*4 LB(90),LT(11)
      INTEGER*4 IPNT(NLINES)
      REAL*4 VX(NLINES),VY(NLINES)
      COMMON/TPOLY WXYZ/ILAST,ILIST,VX,VY,IPNT
      DATA LB/105,102,80,20,02,05,27,87,105,85,103,3,1,5,87,105,102,80,
     *  60,7,0,87,105,102,80,70,52,54,52,30,20,2,5,27,104,57,50,100,0,
     *  100,107,67,62,40,20,2,5,27,80,102,105,87,27,5,2,20,30,52,57,
     *  107,100,4,105,102,80,70,52,55,37,27,5,2,20,30,52,55,77,87,105,
     *  27, 5,2,20,80,102,105,87,77,55,50/
      DATA LT/1,10,15,22,35,40,49,60,63,80,91/
C
C
C     extract digits from GG:
      CALL EXDIG( GG, DIG, JG, J1,J2,J3, ISIGN )
C
      XMK=ABS(XM-XK)
      YMK=ABS(YM-YK)
C
      HGT=XMK/1.2
      H=YMK/(DIG+0.2)
      IF(H .LT. HGT) HGT=H
      ROW=AMIN1(XM,XK)+0.5*(XMK-HGT)
      COL=AMAX1(YM,YK)-0.5*(YMK-DIG*HGT)
      H=HGT/10.0
C
C PLOT 1000THS
      IF (J3 .EQ. 0) GO TO 5
      IB=LT(J3+1)
      IE=LT(J3+2)-1
      DO 2 I=IB,IE
         LLIN=LB(I)/10
         LLEL=LB(I)-LLIN*10
         HL=H*LLIN
         HE=H*LLEL
         IF(I .EQ. IB) GO TO 1 
         ILAST=ILAST-1
         VX(ILAST)=RS
         VY(ILAST)=CS
         ILAST=ILAST-1
         VX(ILAST)=ROW+HL
         VY(ILAST)=COL-HE
1        RS=ROW+HL
         CS=COL-HE
2     CONTINUE
C SPACE FOR COLUMN OF DIGIT
      COL=COL-HGT
C PLOT 100THS
5     IF (J2 .EQ. 0 .AND. J3 .EQ. 0) GO TO 8
      IB=LT(J2+1)
      IE=LT(J2+2)-1
      DO 7 I=IB,IE
         LLIN=LB(I)/10
         LLEL=LB(I)-LLIN*10
         HL=H*LLIN
         HE=H*LLEL
         IF(I .EQ. IB) GO TO 6
         ILAST=ILAST-1
         VX(ILAST)=RS
         VY(ILAST)=CS
         ILAST=ILAST-1
         VX(ILAST)=ROW+HL
         VY(ILAST)=COL-HE
6        RS=ROW+HL
         CS=COL-HE
7     CONTINUE
C SPACE FOR COLUMN OF DIGIT
      COL=COL-HGT
C PLOT 10THS
8     IF(J1 .EQ. 0 .AND. J2 .EQ. 0 .AND. J3 .EQ. 0) GO TO 100
C PLOT DIGIT RIGHT OF DECIMAL
      IB=LT(J1+1)
      IE=LT(J1+2)-1
      DO 20 I=IB,IE
         LLIN=LB(I)/10
         LLEL=LB(I)-LLIN*10
         HL=H*LLIN
         HE=H*LLEL
         IF(I .EQ. IB) GO TO 10
         ILAST=ILAST-1
         VX(ILAST)=RS
         VY(ILAST)=CS
         ILAST=ILAST-1
         VX(ILAST)=ROW+HL
         VY(ILAST)=COL-HE
10       RS=ROW+HL
         CS=COL-HE
20    CONTINUE
C SPACE FOR COLUMN OF DIGIT
      COL=COL-HGT
C PLOT DECIMAL CROSS
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.1*HGT
      VY(ILAST)=COL-0.2*HGT
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.2*HGT
      VY(ILAST)=COL-0.3*HGT
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.2*HGT
      VY(ILAST)=COL-0.2*HGT
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.1*HGT
      VY(ILAST)=COL-0.3*HGT
C HALF SPACE FOR COLUMN OF DECIMAL CROSS
      COL=COL-0.5*HGT
C PLOT DIGITS LEFT OF DECIMAL
100   M=JG-(JG/10)*10
      IB=LT(M+1)
      IE=LT(M+2)-1
      DO 120 I=IB,IE
         LLIN=LB(I)/10
         LLEL=LB(I)-LLIN*10
         HL=H*LLIN
         HE=H*LLEL
         IF(I .EQ. IB) GO TO 110
         ILAST=ILAST-1
         VX(ILAST)=RS
         VY(ILAST)=CS
         ILAST=ILAST-1
         VX(ILAST)=ROW+HL
         VY(ILAST)=COL-HE
110      RS=ROW+HL
         CS=COL-HE
120   CONTINUE
      JG=JG/10
C SPACE FOR COLUMN OF DIGIT
      COL=COL-HGT
      IF(JG .NE. 0) GO TO 100
      IF(ISIGN .GT. 0) RETURN
C PLOT MINUS SIGN
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.5*HGT
      VY(ILAST)=COL-0.4*HGT
      ILAST=ILAST-1
      VX(ILAST)=ROW+0.5*HGT
      VY(ILAST)=COL
      RETURN
      END
