path: root/sys/gio/ncarutil/isosrf.f

SUBROUTINE ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,SLAB,TISO,IFLAG)
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
C DIMENSION OF           T(LU,LV,MW),EYE(3),SLAB(MUVWP2,MUVWP2)
C ARGUMENTS
C
C LATEST REVISION        DECEMBER 1984
C
C PURPOSE                ISOSRF DRAWS AN APPROXIMATION OF AN ISO-VALUED
C                        SURFACE FROM A THREE-DIMENSIONAL ARRAY WITH
C                        HIDDEN LINES REMOVED.
C
C USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
C
C                            CALL EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
C
C                               ASSUMPTIONS:
C                                 -- ALL OF THE T ARRAY IS TO BE USED.
C                                 -- IFLAG IS CHOSEN INTERNALLY.
C                                 -- FRAME IS CALLED BY EZISOS.
C
C                        IF THE ASSUMPTIONS ARE NOT MET, USE
C
C                            CALL ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,
C                                         SLAB,TISO,IFLAG)
C
C ARGUMENTS
C
C ON INPUT               T
C                          THREE DIMENSIONAL ARRAY OF DATA THAT DEFINES
C                          THE ISO-VALUED SURFACE.
C
C                        LU
C                          FIRST DIMENSION OF T IN THE CALLING PROGRAM.
C
C                        MU
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE U DIRECTION (THE FIRST
C                          SUBSCRIPT DIRECTION).  WHEN PROCESSING THE
C                          ENTIRE ARRAY, LU = MU (AND LV = MV).
C
C                        LV
C                          SECOND DIMENSION OF T IN THE CALLING PROGRAM.
C
C                        MV
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE V DIRECTION (THE SECOND
C                          SUBSCRIPT DIRECTION).
C
C                        MV
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE W DIRECTION (THE THIRD
C                          SUBSCRIPT DIRECTION).
C
C                        EYE
C                          THE POSITION OF THE EYE IN THREE-SPACE.  T IS
C                          CONSIDERED TO BE IN A BOX WITH OPPOSITE
C                          CORNERS (1,1,1) AND (MU,MV,MW).  THE EYE IS
C                          AT (EYE(1),EYE(2),EYE(3)), WHICH MUST BE
C                          OUTSIDE THE BOX THAT CONTAINS T.  WHILE GAINING
C                          EXPERIENCE WITH THE ROUTINE, A GOOD CHOICE
C                          FOR EYE MIGHT BE (5.0*MU,3.5*MV,2.0*MW).
C
C                        MUVWP2
C                          THE MAXIMUM OF (MU,MV,MW)+2; THAT IS,
C                          MUVWP2 = MAX(MU,MV,MW)+2).
C
C                        SLAB
C                          A WORK SPACE USED FOR INTERNAL STORAGE.  SLAB
C                          MUST BE AT LEAST MUVWP2*MUVWP2 WORDS LONG.
C
C                        TISO
C                          THE ISO-VALUE USED TO DEFINE THE SURFACE.  THE
C                          SURFACE DRAWN WILL SEPARATE VOLUMES OF T THAT
C                          HAVE VALUES GREATER THAN OR EQUAL TO TISO FROM
C                          VOLUMES OF T THAT HAVE VALUES LESS THAN TISO.
C
C                        IFLAG
C                          THIS FLAG SERVES TWO PURPOSES.
C                          .  FIRST, THE ABSOLUTE VALUE OF IFLAG
C                             DETERMINES WHICH TYPES OF LINES ARE DRAWN
C                             TO APPROXIMATE THE SURFACE.  THREE TYPES
C                             OF LINES ARE CONSIDERED:  LINES OF
C                             CONSTANT U, LINES OF CONSTANT V AND LINES
C                             OF CONSTANT W.  THE FOLLOWING TABLE LISTS
C                             THE TYPES OF LINES DRAWN.
C
C                                                LINES OF CONSTANT
C                                                -----------------
C                             IABS(IFLAG)           U    V    W
C                             -----------          ---  ---  ---
C                                  1               NO   NO   YES
C                                  2               NO   YES  NO
C                                  3               NO   YES  YES
C                                  4               YES  NO   NO
C                                  5               YES  NO   YES
C                                  6               YES  YES  NO
C                             0, 7 OR MORE         YES  YES  YES
C
C                          .  SECOND, THE SIGN OF IFLAG DETERMINES WHAT
C                             IS INSIDE AND WHAT IS OUTSIDE, HENCE,
C                             WHICH LINES ARE VISIBLE AND WHAT IS DONE
C                             AT THE BOUNDARY OF T.  FOR IFLAG:
C
C                             POSITIVE   T VALUES GREATER THAN TISO ARE
C                                        ASSUMED TO BE INSIDE THE SOLID
C                                        FORMED BY THE DRAWN SURFACE.
C                             NEGATIVE   T VALUES LESS THAN TISO ARE
C                                        ASSUMED TO BE INSIDE THE SOLID
C                                        FORMED BY THE DRAWN SURFACE.
C                             IF THE ALGORITHM DRAWS A CUBE, REVERSE THE
C                             SIGN OF IFLAG.
C
C ON OUTPUT              T,LU,MU,LV,MV,MW,EYE,MUVWP2,TISO AND IFLAG ARE
C                        UNCHANGED.  SLAB HAS BEEN WRITTEN IN.
C
C NOTE                   .   THIS ROUTINE IS FOR LOWER RESOLUTION ARRAYS
C                            THAN ISOSRFHR.  40 BY 40 BY 40 IS A
C                            PRACTICAL MAXIMUM.
C                        .   TRANSFORMATIONS CAN BE ACHIEVED BY
C                            ADJUSTING SCALING STATEMENT FUNCTIONS IN
C                            ISOSRF, SET3D AND TR32.
C                        .   THE HIDDEN-LINE ALGORITHM IS NOT EXACT, SO
C                            VISIBILITY ERRORS CAN OCCUR.
C                        .   THREE-DIMENSIONAL PERSPECTIVE CHARACTER
C                            LABELING OF ISOSRF IS POSSIBLE BY USING
C                            THE UTILITY PWRZI.  FOR A DESCRIPTION OF
C                            THE USAGE, SEE THE PWRZI DOCUMENTATION.
C
C ENTRY POINTS           ISOSRF, EZISOS, SET3D, TRN32I, ZEROSC,
C                        STCNTR, DRCNTR, TR32, FRSTS, KURV1S, KURV2S,
C                        FRSTC, FILLIN, DRAWI, ISOSRB, MMASK
C
C COMMON BLOCKS          ISOSR1, ISOSR2, ISOSR3, ISOSR4, ISOSR5,
C                        ISOSR6, ISOSR7, ISOSR8, ISOSR9, TEMPR,
C                        PWRZ1I
C
C REQUIRED LIBRARY       THE ERPRT77 PACKAGE AND THE SPPS.
C ROUTINES
C
C I/O                    PLOTS SURFACE
C
C PRECISION              SINGLE
C
C LANGUAGE               FORTRAN 77
C
C HISTORY                DEVELOPED FOR USERS OF ISOSRFHR WITH SMALLER
C                        ARRAYS.
C
C ALGORITHM              CUTS THROUGH THE THREE-DIMENSIONAL ARRAY ARE
C                        CONTOURED WITH A SMOOTHING CONTOURER WHICH ALSO
C                        MARKS A MODEL OF THE PLOTTING PLANE.  INTERIORS
C                        OF BOUNDARIES ARE FILLED IN AND THE RESULT IS
C                        .OR.ED INTO ANOTHER MODEL OF THE PLOTTING PLANE
C                        WHICH IS USED TO TEST SUBSEQUENT CONTOUR LINES
C                        FOR VISIBILITY.
C
C TIMING                 VARIES WIDELY WITH SIZE OF T AND THE VOLUME OF
C                        THE SPACE ENCLOSED BY THE SURFACE DRAWN.
C
C  **NOTE**                  SPACE REQUIREMENTS CAN BE REDUCED BY
C                            CHANGING THE SIZE OF THE ARRAYS ISCR, ISCA
C                            (FOUND IN COMMON ISOSR2), MASK(FOUND IN
C                            COMMON ISOSR5) AND THE VARIABLE NBPW
C                            (COMMON ISOSR5).
C                            ISCR AND ISCA NEED 128X128 BITS.  SO ON A
C                            64 BIT MACHINE ISCR, ISCA CAN BE
C                            DIMENSIONED TO (2,128). NBPW SET IN
C                            SUBROUTINE MMASK SHOULD CONTAIN THE
C                            NUMBER OF BITS PER WORD YOU WISH TO
C                            UTILIZE.
C                            THE DIMENSION OF MASK AND NMASK SHOULD
C                            EQUAL THE VALUE OF NBPW.
C                            LS SHOULD BE SET TO THE FIRST DIMENSION
C                            OF ISCA AND ISCR.
C
C               EXAMPLES:
C                            ON A 60 BIT MACHINE:
C                                 DIMENSION ISCA(4,128), ISCR(4,128)
C                                 DIMENSION MASK(32)
C                                 NBPW = 32
C                            ON A 64 BIT MACHINE:
C                                 DIMENSION ISCA(2,128), ISCR(2,128)
C                                 DIMENSION MASK(64)
C                                 NBPW = 64
C
C INTERNAL PARAMETERS    NAME  DEFAULT           FUNCTION
C                        ----  -------           --------
C                        IREF     1     FLAG TO CONTROL DRAWING OF AXES.
C                                       .IREF=NONZERO   DRAW AXES.
C                                       .IREF=ZERO      DO NOT DRAW AXES.
C
C
      SAVE
      DIMENSION       T(LU,LV,MW),EYE(3)     ,SLAB(MUVWP2,MUVWP2)
C
      COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
      COMMON /ISOSR4/ RX         ,RY
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR9/ BIG        ,IXBIT
C
      LOGICAL         GENDON
      DATA IREF/1/
C
      AVE(A,B) = (A+B)*.5
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
C +NOAO - Blockdata ISOSRB rewritten as run time initialization
C     EXTERNAL ISOSRB
      call isosrb
C -NOAO
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('NSSL','ISOSRF','ISOSRF','VERSION 12')
      NERR = 0
C
C 3-SPACE  U,V,W,IU,IV,IW,ETC
C 2-SPACE  X,Y,IX,IY,ETC
C
C  INITIALIZE MASKS
C
      IF (.NOT.GENDON) CALL MMASK
C
C  SET SHIFT VALUE FOR X,Y PACKING
C
C  IF YOUR MACHINE HAS MORE THAN 16 BITS PER WORD THIS CHECK MAY BE
C  MODIFIED
C
      IF (LU .LE. 256) GO TO  10
      NERR = NERR + 1
      CALL SETER('DIMENSION OF CUBE EXCEEDS 256',NERR,2)
      RETURN
   10 DO  20 J=1,30
         IF (LU .LE. 2**(J-1)) GO TO  30
   20 CONTINUE
   30 IXBIT = J
      NU = MU
      NUP2 = NU+2
      NV = MV
      NVP2 = NV+2
      NW = MW
      NWP2 = NW+2
      FNU = NU
      FNV = NV
      FNW = NW
      SU1 = SU(1.)
      SV1 = SV(1.)
      SW1 = SW(1.)
      SUNU = SU(FNU)
      SVNV = SV(FNV)
      SWNW = SW(FNW)
      AVEU = AVE(SU1,SUNU)
      AVEV = AVE(SV1,SVNV)
      AVEW = AVE(SW1,SWNW)
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
      NUVWP2 = MUVWP2
      TVAL = TISO
      NFLAG = IABS(IFLAG)
      IF (NFLAG.EQ.0 .OR. NFLAG.GE.8) NFLAG = 7
C
C SET UP SCALING
C
      FACT = -ISIGN(1,IFLAG)
      CALL SET3D (EYE,1.,FNU,1.,FNV,1.,FNW)
C
C BOUND LOWER AND LEFT EDGE OF SLAB
C
      EDGE = SIGN(BIG,FACT)
      DO  40 IUVW=1,NUVWP2
         SLAB(IUVW,1) = EDGE
         SLAB(1,IUVW) = EDGE
   40 CONTINUE
C
C SLICES PERPENDICULAR TO U. THAT IS, V W SLICES. T OF CONSTANT U.
C
      IF (NFLAG .LT. 4) GO TO 100
      CALL ZEROSC
      ISLBT = -1
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO  50 IV=2,NVP2
         SLAB(IV,NWP2) = EDGE
   50 CONTINUE
      DO  60 IW=2,NWP2
         SLAB(NVP2,IW) = EDGE
   60 CONTINUE
C
C GO THRU 3-D ARRAY IN U DIRECTION.  IUEW=IU EITHER WAY.
C PICK IU BASED ON EYEU.
C
      DO  90 IUEW=1,NU
         IU = IUEW
         IF (EYEU .GT. AVEU) IU = NU+1-IUEW
         U = IU
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO  80 IV=1,NV
            DO  70 IW=1,NW
               SLAB(IV+1,IW+1) = T(IU,IV,IW)
   70       CONTINUE
   80    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NVP2,NWP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
   90 CONTINUE
C
C SLICES PERPENDICULAR TO V. U W SLICES. T OF CONSTANT V.
C
  100 IF (MOD(NFLAG/2,2) .EQ. 0) GO TO 160
      CALL ZEROSC
      ISLBT = 0
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO 110 IU=2,NUP2
         SLAB(IU,NWP2) = EDGE
  110 CONTINUE
      DO 120 IW=2,NWP2
         SLAB(NUP2,IW) = EDGE
  120 CONTINUE
C
C GO THRU T IN V DIRECTION.  IVEW=IV EITHER WAY.
C
      DO 150 IVEW=1,NV
         IV = IVEW
         IF (EYEV .GT. AVEV) IV = NV+1-IVEW
         V = IV
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO 140 IU=1,NU
            DO 130 IW=1,NW
               SLAB(IU+1,IW+1) = T(IU,IV,IW)
  130       CONTINUE
  140    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NUP2,NWP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
  150 CONTINUE
C
C SLICES PERPENDICULAR TO W. U V SLICES. T OF CONSTANT W.
C
  160 IF (MOD(NFLAG,2) .EQ. 0) GO TO 220
      CALL ZEROSC
C
      ISLBT = 1
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO 170 IU=2,NUP2
         SLAB(IU,NVP2) = EDGE
  170 CONTINUE
      DO 180 IV=2,NVP2
         SLAB(NUP2,IV) = EDGE
  180 CONTINUE
C
C GO THRU T IN W DIRECTION.
C
      DO 210 IWEW=1,NW
         IW = IWEW
         IF (EYEW .GT. AVEW) IW = NW+1-IWEW
         W = IW
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO 200 IU=1,NU
            DO 190 IV=1,NV
               SLAB(IU+1,IV+1) = T(IU,IV,IW)
  190       CONTINUE
  200    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NUP2,NVP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
  210 CONTINUE
C
C DRAW REFERENCE PLANE EDGES AND W AXIS.
C
  220 IF (IREF .EQ. 0) RETURN
      CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
      IF (EYEV .LT. SV1) GO TO 240
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 230 IU=2,NU
         CALL TRN32I (SU(FLOAT(IU)),SV1,SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  230 CONTINUE
      GO TO 250
  240 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SUNU,SV1,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  250 IF (EYEU .GT. SUNU) GO TO 270
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 260 IV=2,NV
         CALL TRN32I (SUNU,SV(FLOAT(IV)),SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  260 CONTINUE
      GO TO 280
  270 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SUNU,SVNV,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  280 IF (EYEV .GT. SVNV) GO TO 300
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 290 IUOW=2,NU
         CALL TRN32I (SU(FLOAT(NU-IUOW+1)),SVNV,SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  290 CONTINUE
      GO TO 310
  300 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SVNV,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  310 IF (EYEU .LT. SU1) GO TO 330
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 320 IVOW=2,NV
         CALL TRN32I (SU1,SV(FLOAT(NV-IVOW+1)),SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  320 CONTINUE
      GO TO 340
  330 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  340 IF (EYEU.LE.SU1 .OR. EYEV.LE.SV1) GO TO 360
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 350 IW=2,NW
         CALL TRN32I (SU1,SV1,SW(FLOAT(IW)),XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  350 CONTINUE
C +NOAO - Plotit buffer needs to be flushed before returning.
      call plotit (0, 0, 2)
C -NOAO
      RETURN
  360 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SV1,SWNW,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
C +NOAO - Plotit buffer needs to be flushed before returning.
      call plotit (0, 0, 2)
C -NOAO
      RETURN
      END
      SUBROUTINE EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
C
      SAVE
      DIMENSION       T(MU,MV,MW),EYE(3)
C
        DATA ANG,PI/.35,3.141592/
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('NSSL','ISOSRF','EZISOS','VERSION 12')
C
C ARGUMENTS DESCRIBED IN ISOSRF
C
C PICK TYPES OF LINES TO DRAW
C
      NU = MU
      NV = MV
      NW = MW
      TVAL = TISO
      MAX = MAX0(NU,NV,NW)+2
      ATU = NU/2
      ATV = NV/2
      ATW = NW/2
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
      RU = EYEU-ATU
      RV = EYEV-ATV
      RW = EYEW-ATW
      RU2 = RU*RU
      RV2 = RV*RV
      RW2 = RW*RW
      DU = SQRT(RV2+RW2)
      DV = SQRT(RU2+RW2)
      DW = SQRT(RU2+RV2)
      DR = 1./SQRT(RU2+RV2+RW2)
C
C  COMPUTE THE ARCCOSINE
C
      TU = DU*DR
      ANGU = ATAN(ABS(SQRT(1.-TU*TU)/TU))
      IF (TU .LE. 0.) ANGU = PI-ANGU
      TV = DV*DR
      ANGV = ATAN(ABS(SQRT(1.-TV*TV)/TV))
      IF (TV .LE. 0.) ANGV = PI-ANGV
      TW = DW*DR
      ANGW = ATAN(ABS(SQRT(1.-TW*TW)/TW))
      IF (TW .LE. 0.) ANGW = PI-ANGW
C
C BREAK POINT IS ABOUT 20 DEGREES OR ABOUT .35 RADIANS
C
      IFLAG = 0
      IF (ANGU .GT. ANG) IFLAG = IFLAG+4
      IF (ANGV .GT. ANG) IFLAG = IFLAG+2
      IF (ANGW .GT. ANG) IFLAG = IFLAG+1
C
C FIND SIGN OF IFLAG
C
      ICNT = 0
      IF (ABS(RU) .LE. ATU) GO TO  30
      IU = 1
      IF (EYEU .GT. ATU) IU = NU
      DO  20 IW=1,NW
         DO  10 IV=1,NV
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   10    CONTINUE
   20 CONTINUE
   30 IF (ABS(RV) .LE. ATV) GO TO  60
      IV = 1
      IF (EYEV .GT. ATV) IV = NV
      DO  50 IW=1,NW
         DO  40 IU=1,NU
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   40    CONTINUE
   50 CONTINUE
   60 IF (ABS(RW) .LE. ATW) GO TO  90
      IW = 1
      IF (EYEW .GT. ATW) IW = NW
      DO  80 IV=1,NV
         DO  70 IU=1,NU
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   70    CONTINUE
   80 CONTINUE
   90 IFLAG = ISIGN(IFLAG,ICNT)
      CALL ISOSRF (T,NU,NU,NV,NV,NW,EYE,MAX,SLAB,TVAL,IFLAG)
C +NOAO - Call to frame is suppressed.
C     CALL FRAME
C -NOAO
      RETURN
      END
      SUBROUTINE SET3D (EYE,ULO,UHI,VLO,VHI,WLO,WHI)
      SAVE
      COMMON /TEMPR/  RZERO
C
      DIMENSION       EYE(3)
C
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
      COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
     1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
     2                EYEV       ,EYEW
C
C
      AVE(A,B) = (A+B)*.5
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
C CONSTANTS FOR PWRZ
C
      UUMIN = ULO
      UUMAX = UHI
      VVMIN = VLO
      VVMAX = VHI
      WWMIN = WLO
      WWMAX = WHI
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
C
C FIND CORNERS IN 2-SPACE FOR 3-SPACE BOX CONTAINING OBJECT
C
      ISCALE = 0
      ATU = AVE(SU(UUMIN),SU(UUMAX))
      ATV = AVE(SV(VVMIN),SV(VVMAX))
      ATW = AVE(SW(WWMIN),SW(WWMAX))
      BIGD = 0.
      IF (RZERO .LE. 0.) GO TO  10
C
C RELETIVE SIZE FEATURE IN USE.
C GENERATE EYE POSITION THAT MAKES BOX HAVE MAXIMUM PROJECTED SIZE.
C
      ALPHA = -(VVMIN-ATV)/(UUMIN-ATU)
      VVEYE = -RZERO/SQRT(1.+ALPHA*ALPHA)
      UUEYE = VVEYE*ALPHA
      VVEYE = VVEYE+ATV
      UUEYE = UUEYE+ATU
      WWEYE = ATW
      CALL TRN32I (ATU,ATV,ATW,UUEYE,VVEYE,WWEYE,1)
      CALL TRN32I (UUMIN,VVMIN,ATW,XMIN,DUMM,DUMM,2)
      CALL TRN32I (UUMAX,VVMIN,WWMIN,DUMM,YMIN,DUMM,2)
      CALL TRN32I (UUMAX,VVMAX,ATW,XMAX,DUMM,DUMM,2)
      CALL TRN32I (UUMAX,VVMIN,WWMAX,DUMM,YMAX,DUMM,2)
      BIGD = SQRT((UUMAX-UUMIN)**2+(VVMAX-VVMIN)**2+(WWMAX-WWMIN)**2)*.5
      R0 = RZERO
      GO TO  20
   10 CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
      CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMIN),X1,Y1,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMAX),X2,Y2,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMIN),X3,Y3,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMAX),X4,Y4,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMIN),X5,Y5,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMAX),X6,Y6,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMIN),X7,Y7,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMAX),X8,Y8,DUM,2)
      XMIN = AMIN1(X1,X2,X3,X4,X5,X6,X7,X8)
      XMAX = AMAX1(X1,X2,X3,X4,X5,X6,X7,X8)
      YMIN = AMIN1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
      YMAX = AMAX1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
C
C ADD RIGHT AMOUNT TO KEEP PICTURE SQUARE
C
   20 WIDTH = XMAX-XMIN
      HIGHT = YMAX-YMIN
      DIF = .5*(WIDTH-HIGHT)
      IF (DIF)  30, 50, 40
   30 XMIN = XMIN+DIF
      XMAX = XMAX-DIF
      GO TO  50
   40 YMIN = YMIN-DIF
      YMAX = YMAX+DIF
   50 ISCALE = 1
      CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
      RETURN
      END
      SUBROUTINE TRN32I (U,V,W,XT,YT,ZT,IENT)
C
C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR-
C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE
C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15,
C 2, 193-204,1968.
C ARGUMENTS FOR SET
C U,V,W    ARE THE 3-SPACE COORDINATES OF THE INTERSECTION
C          OF THE LINE OF SIGHT AND THE IMAGE PLANE.  THIS
C          POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT.
C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION.
C
C TRN32 ARGUMENTS
C U,V,W    ARE THE 3-SPACE COORDINATES OF A POINT TO BE
C          TRANSFORMED.
C XT,YT    THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR-
C          MATION.  WHEN ISCALE=0, XT AND YT ANR IN THE SAME
C          UNITS AS U,V, AND W.  WHEN ISCALE'0, XT AND YT
C          ARE IN PLOTTER COORDINATES.
C ZT       NOT USED.
C
      SAVE
      COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
     1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
     2                EYEV       ,EYEW
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
C
C RANGE OF PLOTTER COORDINATES
C
C
C  WARNING
C  IF PLOTTER MAXIMUM VALUE RANGES (IN X OR Y DIRECTION) FALL BELOW
C  101, THEN CHANGES MUST BE MADE IN SUBROUTINE FRSTC. THE REQUIRED
C  CHANGES ARE MARKED BY WARNING COMMENTS IN FRSTC.
      DATA NLX,NBY,NRX,NTY/10,10,32760,32760/
        DATA PI/3.1411592/
C
C STORE THE PARAMETERS OF THE SET CALL FOR USE
C WITH THE TRANSLATE CALL
C
C DECIDE IF SET OR TRANSLATE CALL
C
      IF (IENT .NE. 1) GO TO  50
      AU = U
      AV = V
      AW = W
      EU = XT
      EV = YT
      EW = ZT
C
C
C
C
C
      DU = AU-EU
      DV = AV-EV
      DW = AW-EW
      D = SQRT(DU*DU+DV*DV+DW*DW)
      COSAL = DU/D
      COSBE = DV/D
      COSGA = DW/D
C
C  COMPUTE THE ARCCOSINE
C
      AL = ATAN(ABS(SQRT(1.-COSAL*COSAL)/COSAL))
      IF (COSAL .LE. 0.) AL = PI-AL
      BE = ATAN(ABS(SQRT(1.-COSBE*COSBE)/COSBE))
      IF (COSBE .LE. 0.) BE = PI-BE
      GA = ATAN(ABS(SQRT(1.-COSGA*COSGA)/COSGA))
      IF (COSGA .LE. 0.) GA = PI-GA
      SINGA = SIN(GA)
C
C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF
C THE 2-SPACE.  THE 3-SPACE W AXIS IS TRANSFORMED INTO THE
C 2-SPACE Y AXIS.  IF THE LINE OF SIGHT IS CLOSE TO PARALLEL
C TO THE 3-SPACE W AXIS, THE 3-SPACE V AXIS IS CHOSEN (IN-
C STEAD OF THE 3-SPACE W AXIS) TO BE TRANSFORMED INTO THE
C 2-SPACE Y AXIS.
C
      ASSIGN  90 TO JDONE
      IF (ISCALE)  10, 30, 10
   10 X0 = XMIN
      Y0 = YMIN
      X1 = NLX
      Y1 = NBY
      X2 = NRX-NLX
      Y2 = NTY-NBY
      X3 = X2/(XMAX-XMIN)
      Y3 = Y2/(YMAX-YMIN)
      X4 = NRX
      Y4 = NTY
      FACT = 1.
      IF (BIGD .LE. 0.) GO TO  20
      X0 = -BIGD
      Y0 = -BIGD
      X3 = X2/(2.*BIGD)
      Y3 = Y2/(2.*BIGD)
      FACT = R0/D
   20 DELCRT = X2
      ASSIGN  80 TO JDONE
   30 IF (SINGA .LT. 0.0001) GO TO  40
      R = 1./SINGA
      ASSIGN  70 TO JUMP
      RETURN
   40 SINBE = SIN(BE)
      R = 1./SINBE
      ASSIGN  60 TO JUMP
      RETURN
C
C********************  ENTRY TRN32  ************************
C     ENTRY TRN32 (U,V,W,XT,YT,ZT)
C
   50 UU = U
      VV = V
      WW = W
      Q = D/((UU-EU)*COSAL+(VV-EV)*COSBE+(WW-EW)*COSGA)
      GO TO JUMP,( 60, 70)
   60 UU = ((EW+Q*(WW-EW)-AW)*COSAL-(EU+Q*(UU-EU)-AU)*COSGA)*R
      VV = (EV+Q*(VV-EV)-AV)*R
      GO TO JDONE,( 80, 90)
   70 UU = ((EU+Q*(UU-EU)-AU)*COSBE-(EV+Q*(VV-EV)-AV)*COSAL)*R
      VV = (EW+Q*(WW-EW)-AW)*R
      GO TO JDONE,( 80, 90)
   80 XT = AMIN1(X4,AMAX1(X1,X1+X3*(FACT*UU-X0)))
      YT = AMIN1(Y4,AMAX1(Y1,Y1+Y3*(FACT*VV-Y0)))
      RETURN
   90 XT = UU
      YT = VV
      RETURN
      END
      SUBROUTINE ZEROSC
      SAVE
C
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
C
C ZERO BOTH SCRENE MODELS.
C
      DO  20 I=1,LX
         DO  10 J=1,NY
            ISCR(I,J) = 0
            ISCA(I,J) = 0
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE STCNTR (Z,L,M,N,CONV)
C
      SAVE
      DIMENSION       Z(L,N)
C
C THIS ROUTINE FINDS THE BEGINNINGS OF ALL CONTOUR LINES AT LEVEL CONV.
C FIRST THE EDGES ARE SEARCHED FOR LINES INTERSECTING THE EDGE (OPEN
C LINES) THEN THE INTERIOR IS SEARCHED FOR LINES WHICH DO NOT INTERSECT
C THE EDGE (CLOSED LINES).  BEGINNINGS ARE STORED IN IR TO PREVENT RE-
C TRACING OF LINES.  IF IR IS FILLED, THE SEARCH IS STOPPED FOR THIS
C CONV.
C
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR9/ BIG        ,IXBIT
C
C  PACK X AND Y
C
      IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
C
      IENTRY = 0
      NP = 0
      CV = CONV
C
C THE FOLLOWING CODE SHOULD BE RE-ENABLED IF THIS ROUTINE IS USED FOR
C GENERAL CONTOURING
C
C     ISS=0
C     DO 2 IP1=2,M
C     I=IP1-1
C     IF(Z(I,1).GE.CV.OR.Z(IP1,1).LT.CV) GO TO 1
C     IX=IP1
C     IY=1
C     IDX=-1
C     IDY=0
C     IS=1
C     CALL DRLINE(Z,L,M,N)
C   1 IF(Z(IP1,N).GE.CV.OR.Z(I,N).LT.CV) GO TO 2
C     IX=I
C     IY=N
C     IDX=1
C     IDY=0
C     IS=5
C     CALL DRLINE(Z,L,M,N)
C   2 CONTINUE
C     DO 4 JP1=2,N
C     J=JP1-1
C     IF(Z(M,J).GE.CV.OR.Z(M,JP1).LT.CV) GO TO 3
C     IX=M
C     IY=JP1
C     IDX=0
C     IDY=-1
C     IS=7
C     CALL DRLINE(Z,L,M,N)
C   3 IF(Z(1,JP1).GE.CV.OR.Z(1,J).LT.CV) GO TO 4
C     IX=1
C     IY=J
C     IDX=0
C     IDY=1
C     IS=3
C     CALL DRLINE(Z,L,M,N)
C   4 CONTINUE
C
      ISS = 1
      DO  40 JP1=3,N
         J = JP1-1
         DO  30 IP1=2,M
            I = IP1-1
            IF (Z(I,J).GE.CV .OR. Z(IP1,J).LT.CV) GO TO  30
            IXY = IPXY(IP1,J)
            IF (NP .EQ. 0) GO TO  20
            DO  10 K=1,NP
               IF (IR(K) .EQ. IXY) GO TO  30
   10       CONTINUE
   20       NP = NP+1
            IF (NP .GT. NR) RETURN
            IR(NP) = IXY
            IX = IP1
            IY = J
            IDX = -1
            IDY = 0
            IS = 1
            CALL DRCNTR (Z,L,M,N)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE DRCNTR (Z,L,MM,NN)
      SAVE
C
      DIMENSION       Z(L,NN)
C
C THIS ROUTINE TRACES A CONTOUR LINE WHEN GIVEN THE BEGINNING BY STLINE.
C TRANSFORMATIONS CAN BE ADDED BY DELETING THE STATEMENT FUNCTIONS FOR
C FX AND FY IN DRLINE AND MINMAX AND ADDING EXTERNAL FUNCTIONS.
C X=1. AT Z(1,J), X=FLOAT(M) AT Z(M,J). X TAKES ON NON-INTEGER VALUES.
C Y=1. AT Z(I,1), Y=FLOAT(N) AT Z(I,N). Y TAKES ON NON-INTEGER VALUES.
C
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR9/ BIG        ,IXBIT
C
      LOGICAL         IPEN       ,IPENO
C
      DATA IOFFP,SPVAL/0,0./
      DATA IPEN,IPENO/.TRUE.,.TRUE./
C
C  PACK X AND Y
C
      IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
      FX(X1,Y1) = X1
      FY(X1,Y1) = Y1
      C(P11,P21) = (P11-CV)/(P11-P21)
C
      M = MM
      N = NN
      IF (IOFFP .EQ. 0) GO TO  10
      ASSIGN 100 TO JUMP1
      ASSIGN 150 TO JUMP2
      GO TO  20
   10 ASSIGN 120 TO JUMP1
      ASSIGN 160 TO JUMP2
   20 IX0 = IX
      IY0 = IY
      IS0 = IS
      IF (IOFFP .EQ. 0) GO TO  30
      IX2 = IX+INX(IS)
      IY2 = IY+INY(IS)
      IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL
      IPENO = IPEN
   30 IF (IDX .EQ. 0) GO TO  40
      Y = IY
      ISUB = IX+IDX
      X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
      GO TO  50
   40 X = IX
      ISUB = IY+IDY
      Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
   50 IF (IPEN) CALL FRSTS (FX(X,Y),FY(X,Y),1)
   60 IS = IS+1
      IF (IS .GT. 8) IS = IS-8
      IDX = INX(IS)
      IDY = INY(IS)
      IX2 = IX+IDX
      IY2 = IY+IDY
      IF (ISS .NE. 0) GO TO  70
      IF (IX2.GT.M .OR. IY2.GT.N .OR. IX2.LT.1 .OR. IY2.LT.1) GO TO 190
   70 IF (CV-Z(IX2,IY2))  80, 80, 90
   80 IS = IS+4
      IX = IX2
      IY = IY2
      GO TO  60
   90 IF (IS/2*2 .EQ. IS) GO TO  60
      GO TO JUMP1,(100,120)
  100 ISBIG = IS+(8-IS)/6*8
      IX3 = IX+INX(ISBIG-1)
      IY3 = IY+INY(ISBIG-1)
      IX4 = IX+INX(ISBIG-2)
      IY4 = IY+INY(ISBIG-2)
      IPENO = IPEN
      IF (ISS .NE. 0) GO TO 110
      IF (IX3.GT.M .OR. IY3.GT.N .OR. IX3.LT.1 .OR. IY3.LT.1) GO TO 190
      IF (IX4.GT.M .OR. IY4.GT.N .OR. IX4.LT.1 .OR. IY4.LT.1) GO TO 190
  110 IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL .AND.
     1       Z(IX3,IY3).NE.SPVAL .AND. Z(IX4,IY4).NE.SPVAL
  120 IF (IDX .EQ. 0) GO TO 130
      Y = IY
      ISUB = IX+IDX
      X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
      GO TO 140
  130 X = IX
      ISUB = IY+IDY
      Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
  140 GO TO JUMP2,(150,160)
  150 IF (.NOT.IPEN) GO TO 170
      IF (IPENO) GO TO 160
C
C END OF LINE SEGMENT
C
      CALL FRSTS (D1,D2,3)
      CALL FRSTS (FX(XOLD,YOLD),FY(XOLD,YOLD),1)
C
C CONTINUE LINE SEGMENT
C
  160 CALL FRSTS (FX(X,Y),FY(X,Y),2)
  170 XOLD = X
      YOLD = Y
      IF (IS .NE. 1) GO TO 180
      NP = NP+1
      IF (NP .GT. NR) GO TO 190
      IR(NP) = IPXY(IX,IY)
  180 IF (ISS .EQ. 0) GO TO  60
      IF (IX.NE.IX0 .OR. IY.NE.IY0 .OR. IS.NE.IS0) GO TO  60
C
C END OF LINE
C
  190 CALL FRSTS (D1,D2,3)
      RETURN
      END
      SUBROUTINE TR32 (X,Y,MX,MY)
      SAVE
C
      COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
      XX = X
      YY = Y
      IF (ISLBT)  10, 20, 30
   10 CALL TRN32I (SU(U),SV(XX-1.),SW(YY-1.),XT,YT,DUM,2)
      GO TO  40
   20 CALL TRN32I (SU(XX-1.),SV(V),SW(YY-1.),XT,YT,DUM,2)
      GO TO  40
   30 CALL TRN32I (SU(XX-1.),SV(YY-1.),SW(W),XT,YT,DUM,2)
   40 MX = XT
      MY = YT
      RETURN
      END
      SUBROUTINE FRSTS (XX,YY,IENT)
C
C THIS IS A SPECIAL VERSION OF THE SMOOTHING DASHED LINE PACKAGE.  LINES
C ARE SMOOTHED IN THE SAME WAY, BUT NO SOFTFARE DASHED LINES ARE USED.
C CONDITIONAL PLOTTING ROUTINES ARE CALL WHICH DETERMINE THE VISIBILITY
C OF A LINE SEGMENT BEFORE PLOTTING.
C
      SAVE
      DIMENSION       XSAVE(70)  ,YSAVE(70)  ,XP(70)     ,YP(70)     ,
     1                TEMP(70)
C
      COMMON /ISOSR7/ IENTRY     ,IONES
C
      DATA NP/150/
      DATA L1/70/
      DATA TENSN/2.5/
      DATA PI/3.14159265358/
      DATA SMALL/128./
C
      AVE(A,B) = .5*(A+B)
C
C  DECIDE IF FRSTS,VECTS,LASTS CALL
C
      GO TO ( 10, 20, 40),IENT
   10 DEG = 180./PI
      X = XX
      Y = YY
      LASTFL = 0
      SSLP1 = 0.0
      SSLPN = 0.0
      XSVN = 0.0
      YSVN = 0.0
C
C INITIALIZE THE POINT AND SEGMENT COUNTER
C N COUNTS THE NUMBER OF POINTS/SEGMENT
C
      N = 0
C
C NSEG = 0       FIRST SEGMENT
C NSEG = 1       MORE THAN ONE SEGMENT
C
      NSEG = 0
      CALL TR32 (X,Y,MX,MY)
C
C SAVE THE X,Y COORDINATES OF THE FIRST POINT
C XSV1           CONTAINS THE X COORDINATE OF THE FIRST POINT
C                OF A LINE
C YSV1           CONTAINS THE Y COORDINATE OF THE FIRST POINT
C                OF A LINE
C
      XSV1 = MX
      YSV1 = MY
      GO TO  30
C
C ************************* ENTRY VECTS *************************
C     ENTRY VECTS (XX,YY)
C
   20 X = XX
      Y = YY
C
C VECTS          SAVES THE X,Y COORDINATES OF THE ACCEPTED
C                POINTS ON A LINE SEGMENT
C
      CALL TR32 (X,Y,MX,MY)
C
CIF THE NEW POINT IS TOO CLOSE TO THE PREVIOUS POINT, IGNORE IT
C
      IF (ABS(FLOAT(IFIX(XSVN)-MX))+ABS(FLOAT(IFIX(YSVN)-MY)) .LT.
     1    SMALL) RETURN
      IFLAG = 0
   30 N = N+1
C
C SAVE THE X,Y COORDINATES OF EACH POINT OF THE SEGMENT
C XSAVE          THE ARRAY OF X COORDINATES OF LINE SEGMENT
C YSAVE          THE ARRAY OF Y COORDINATES OF LINE SEGMENT
C
      XSAVE(N) = MX
      YSAVE(N) = MY
      XSVN = XSAVE(N)
      YSVN = YSAVE(N)
      IF (N .GE. L1-1) GO TO  50
      RETURN
C
C ************************* ENTRY LASTS *************************
C     ENTRY LASTS
C
   40 LASTFL = 1
C
C LASTS          CHECKS FOR PERIODIC LINES AND SETS UP
C                  THE CALLS TO KURV1S AND KURV2S
C
C IFLAG = 0      OK TO CALL LASTS DIRECTLY
C IFLAG = 1      LASTS WAS JUST CALLED FROM BY VECTS
C                IGNORE CALL TO LASTS
C
      IF (IFLAG .EQ. 1) RETURN
C
C COMPARE THE LAST POINT OF SEGMENT WITH FIRST POINT OF LINE
C
   50 IFLAG = 1
C
C IPRD = 0       PERIODIC LINE
C IPRD = 1       NON-PERIODIC LINE
C
      IPRD = 1
      IF (ABS(XSV1-XSVN)+ABS(YSV1-YSVN) .LT. SMALL) IPRD = 0
C
C TAKE CARE OF THE CASE OF ONLY TWO DISTINCT P0INTS ON A LINE
C
      IF (NSEG .GE. 1) GO TO  70
      IF (N-2) 160,150, 60
   60 IF (N .GE. 4) GO TO  70
      DX = XSAVE(2)-XSAVE(1)
      DY = YSAVE(2)-YSAVE(1)
      SLOPE = ATAN2(DY,DX)*DEG+90.
      IF (SLOPE .GE. 360.) SLOPE = SLOPE-360.
      IF (SLOPE .LE. 0.) SLOPE = SLOPE+360.
      SLP1 = SLOPE
      SLPN = SLOPE
      ISLPSW = 0
      SIGMA = TENSN
      GO TO 110
   70 SIGMA = TENSN
      IF (IPRD .GE. 1) GO TO  90
      IF (NSEG .GE. 1) GO TO  80
C
C SET UP FLAGS FOR A  1  SEGMENT, PERIODIC LINE
C
      ISLPSW = 4
      XSAVE(N) = XSV1
      YSAVE(N) = YSV1
      GO TO 110
C
C SET UP FLAGS FOR AN N-SEGMENT, PERIODIC LINE
C
   80 SLP1 = SSLPN
      SLPN = SSLP1
      ISLPSW = 0
      GO TO 110
   90 IF (NSEG .GE. 1) GO TO 100
C
C SET UP FLAGS FOR THE 1ST SEGMENT OF A NON-PERIODIC LINE
C
      ISLPSW = 3
      GO TO 110
C
C SET UP FLAGS FOR THE NTH SEGMENT OF A NON-PERIODIC LINE
C
  100 SLP1 = SSLPN
      ISLPSW = 1
C
C CALL THE SMOOTHING ROUTINES
C
  110 CALL KURV1S (N,XSAVE,YSAVE,SLP1,SLPN,XP,YP,TEMP,S,SIGMA,ISLPSW)
      IF (IPRD.EQ.0 .AND. NSEG.EQ.0 .AND. S.LT.70.) GO TO 170
      IENTRY = 1
C
C DETERMINE THE NUMBER OF POINTS TO INTERPOLATE FOR EACH SEGMENT
C
      IF (NSEG.GE.1 .AND. N.LT.L1-1) GO TO 120
      NPRIME = FLOAT(NP)-(S*FLOAT(NP))/(2.*32768.)
      IF (S .GE. 32768.) NPRIME = .5*FLOAT(NP)
      NPL = FLOAT(NPRIME)*S/32768.
      IF (NPL .LT. 2) NPL = 2
  120 DT = 1./FLOAT(NPL)
      IF (NSEG .LE. 0) CALL FRSTC (IFIX(XSAVE(1)),IFIX(YSAVE(1)),1)
      T = 0.0
      NSLPSW = 1
      IF (NSEG .GE. 1) NSLPSW = 0
      NSEG = 1
      CALL KURV2S (T,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
C
C SAVE SLOPE AT THE FIRST POINT OF THE LINE
C
      IF (NSLPSW .GE. 1) SSLP1 = SLP
      NSLPSW = 0
      XSOLD = XSAVE(1)
      YSOLD = YSAVE(1)
      DO 130 I=1,NPL
         T = T+DT
         TT = -T
         IF (I .EQ. NPL) NSLPSW = 1
         CALL KURV2S (TT,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
C
C SAVE THE LAST SLOPE OF THIS LINE SEGMENT
C
         IF (NSLPSW .GE. 1) SSLPN = SLP
C
C DRAW EACH PART OF THE LINE SEGMENT
C
         CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
         CALL FRSTC (IFIX(XS),IFIX(YS),2)
         XSOLD = XS
         YSOLD = YS
  130 CONTINUE
      IF (IPRD .NE. 0) GO TO 140
C
C CONNECT THE LAST POINT WITH THE FIRST POINT OF A PERIODIC LINE
C
      CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
      CALL FRSTC (IFIX(XSV1),IFIX(YSV1),2)
C
C BEGIN THE NEXT LINE SEGMENT WITH THE LAST POINT OF THIS SEGMENT
C
  140 XSAVE(1) = XS
      YSAVE(1) = YS
      N = 1
  150 CONTINUE
  160 RETURN
  170 N = 0
      RETURN
      END
      SUBROUTINE FRSTC (MX,MY,IENT)
      SAVE
C
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR4/ RX         ,RY
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
     1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
     2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
     3                IHB        ,IHS        ,IHV        ,IVOLD      ,
     4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
     5                IHF
      LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
C
C
C  DRAW LINE TO THE POINT MX,MY
C
C  ENTER THE POINT INTO THE CURRENT SCREEN, ISCR, IF THE POINT CONFORMS
C  TO THE SHADING ALGORITHM.
C  THE POINT IS NOT ENTERED WHEN;
C  1. IT IS THE SAME POINT USED IN THE LAST CALL, RESOLUTION PROBLEM
C  2. IT IS PART OF A HORIZONTAL LINE BUT NOT AN END POINT
C  3.  THE ENTIRE CONTOUR RESTS ON A HORIZONTAL PLANE
C
C  WHEN DRAWING A HORIZONTAL LINE THREE CONDITIONS EXIST;
C  1. WHEN THE LINE IS A HORIZONTAL STEP ENTER ONLY THE OUTSIDE POINT.
C     A HORIZONTAL STEP IS DEFINED BY THE ENTERING AND EXITING Y
C     DIRECTION THAT IS THE SAME.
C  2. ENTER BOTH END POINTS OF A HORIZONTAL TURNING POINT. A HORIZONTAL
C     TURNING POINT IS A LINE WITH GREATER THAN 1 HORIZONTAL BITS
C     AND THE ENTERING AND EXITING Y DIRECTION IS DIFFIRENT.
C  3.  WHEN THE ENTIRE CONTOUR IS A HORIZONTAL LINE NO POINTS ARE
C      ENTERED. THIS CONDITION IS DETECTED BY THE STATUS OF YCHANG.
C      IF IT IS TRUE THEN THE CONTOUR IS NOT A SINGLE HORIZONTAL LINE.
C
C  THE PREVIOUS POINT IS ERASED IF IT IS A VERTICAL TURNING POINT.
C  A VERTICAL TURNING POINT IS A HORIZONTAL LINE WITH ONLY 1 POINT
C  AND THE ENTERING AND EXITING Y DIRECTION DIFFERS.THIS DATA IS
C  IN THE VARIABLES IOSLSN-OLD SLOPE AND ISLSGN-NEW SLOPE.
C  THE CHANGE IN SLOPE MUST BE -1 TO 1 OR 1 TO -1.
C
C  OTHERWISE THE POINT IS ENTERED INTO ISCR.
C
C  THE TWO ENTRY POINTS ARE REQUIRED BY THE HARDWARE DRAWING ROUTINES.
C  FIRSTC IS USED FOR THE FIRST POINT ON THE CONTOUR. THE REMAINING
C  POINTS ON THE SAME CONTOUR ARE ENTERED VIA VECTC.
C
      DATA IONE/1/
      AVE(A,B) = (A+B)*.5
C
C  COMPUTE VISIBILITY OF THIS POINT
C
C  WARNING
C  IF X OR Y PLOTTER MAXIMUM VALUE RANGES FALL BELOW 101 THEN THE
C  FOLLOWING TWO STATEMENTS WHICH SET IX AND IY MUST BE CHANGED.
C  REPLACE THE CONSTANT 1.0 BY 0.5 IN THE STATEMENTS WHERE THE
C  MAXIMUM PLOTTER VALUE IS LESS THAN 101 FOR THAT DIRECTION.  THE
C  PLOTTER CORDINATE RANGES ARE SET IN SET32.
C
      IX = FLOAT(MX-1)*RX+1.0
      NRLX = IX
      IY = FLOAT(MY-1)*RY+1.0
      IBIT = NBPW-MOD(IX,NBPW)
      IX = IX/NBPW+1
      IVNOW = IAND(ISHIFT(ISCA(IX,IY),1-IBIT),IONE)
C
C  DECIDE IF FRSTC OR VECTC CALL
C
      IF (IENT .NE. 1) GO TO  10
C
      XOLD = MX
      YOLD = MY
C
C
C  SET INITIAL VALUES
C
      IHF = .FALSE.
      IYDIR = 0
      ITPD = 0
      IVAL = 0
      IOSLSN = 0
      IFSX = NRLX
      IFSY = IY
      LASTV = IVNOW
      HBFLAG = .FALSE.
      YCHANG = .FALSE.
      CALL PLOTIT (IFIX(XOLD),IFIX(YOLD),0)
      GO TO 180
C
C****************************  ENTRY VECTC  ****************************
C     ENTRY VECTC (MX,MY)
C
   10 XNOW = MX
      YNOW = MY
      JUMP = IVNOW*2+LASTV+1
      GO TO ( 20, 30, 40, 50),JUMP
C
C BOTH VISIBLE
C
   20 CALL PLOTIT (IFIX(XNOW),IFIX(YNOW),1)
      GO TO  50
C
C JUST TURNED VISIBLE
C
   30 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),0)
      GO TO  50
C
C JUST TURNED INVISIBLE
C
   40 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),1)
C
C BOTH INVISIBLE
C
   50 XOLD = XNOW
      YOLD = YNOW
      LASTV = IVNOW
C
C  TEST FOR RESOLUTION PROBLEM
C
      IF (NRLX.EQ.LRLX .AND. IY.EQ.IYOLD) RETURN
C
C  TEST FOR HORIZONTAL BITS
C
      IF (IYOLD .NE. IY) GO TO  70
C
C  HORIZONTAL BITS DETECTED. SET FLAG AND EXIT.
C  THIS AND THE NEXT HORIZONTAL BIT TEST IS NECESSARY FOR ISCR TO
C  CONFORM TO THE SHADING ALGORITHM IN SUBROUTINE FILLIN
C
C
C  IF HORIZONTAL LINE PREVIOUSLY DETECTED EXIT
C
      IF (.NOT.HBFLAG) GO TO  60
C
C  IF END OF CONTOUR ON A HORIZONTAL LINE BRANCH FOR SPECIAL PROCESSING.
C
      IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
      GO TO 200
C
C  SAVE SLOPE PRIOR TO HORIZONTAL LINE
C
   60 IHX = IXOLD
      IHB = IBTOLD
      IHS = IOSLSN
      IOSLSN = 0
      HBFLAG = .TRUE.
      IHRX = LRLX
      IHV = IVOLD
      IF (LRLX.EQ.IFSX .AND. IYOLD.EQ.IFSY) IHF = .TRUE.
C
C  THIS IS THE SECOND TRAP FOR END OF CONTOUR ON A HORIZONTAL LINE.
C
      IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
      GO TO 200
C
C  COMPUTE THE SLOPE TO THIS POINT
C
   70 IF (IY-IYOLD)  80, 90,100
   80 ISLSGN = 1
      GO TO 110
   90 ISLSGN = 0
      GO TO 120
  100 ISLSGN = -1
  110 IF (IYDIR .EQ. 0) IYDIR = ISLSGN
  120 CONTINUE
C
C  IF PROCESS REACHES THIS CODE THE CONTOUR IS NOT CONTAINED ON A SINGLE
C  HORIZONTAL  PLANE, SO RECORD THIS FACT BY SETTING Y CHANGE FLAG.
C
      YCHANG = .TRUE.
C
C  TEST FOR END OF HORIZONTAL LINE
C
      IF (.NOT.HBFLAG) GO TO 160
      HBFLAG = .FALSE.
C
C  HORIZONTAL LINE JUST ENDED
C
C  TEST FOR REDRAW
C
      ITEMP = IAND(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
      IF ((IHV .EQ. 0) .AND. (ITEMP .EQ. 0)) GO TO 130
C
C  REDRAWING ERASE THIS POINT
C
      ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
      ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
      GO TO 170
C
C  TEST FOR STEP PROBLEM
C
  130 IF (IHS .NE. ISLSGN) GO TO 140
C
C  STEP PROBLEM
C
      GO TO 170
C
C  TURNING PROBLEM HORIZONTAL LINE IS A TURNING POINT
C
  140 CONTINUE
C
C  ENTER THE TURNING POINT ONLY IF IT IS NOT THE SECOND SUCCEEDING
C  EVENT IN A ROW
C
      ICTPD = 1
      IF (IHRX .GT. NRLX) ICTPD = -1
      IF (ICTPD .NE. ITPD) GO TO 150
      ITPD = 0
C
C  ERASE THE FIRST POINT
C
      ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
      GO TO 170
C
C  ENTER THE TURNING POINT
C
  150 CONTINUE
      ITPD = ICTPD
C
C  ENTER THE SECOND POINT
C
      ISCR(IXOLD,IYOLD) = IOR(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
      GO TO 170
C
C  CHECK IF PREVIOUS ENTRY WAS A VERTICAL TURNING POINT.
C  IF SO ERASE IT.
C
  160 IF (ISLSGN.EQ.IOSLSN .OR. (IOSLSN.EQ.0 .OR. ISLSGN.EQ.0))
     1    GO TO 170
      ITPD = 0
      ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
C
  170 IOSLSN = ISLSGN
C
C  CHECK IF THIS GRID POINT PREVIOUSLY ACTIVATED
C
      IVAL = IAND(ISCR(IX,IY),MASK(IBIT))
C
C  IF GRID POINTS ACTIVATED BRANCH
C
      IF (IVAL .NE. 0) GO TO 190
C
C  GRID POINT NOT ACTIVATED   SET AND EXIT
C
  180 CONTINUE
      ISCR(IX,IY) = IOR(ISCR(IX,IY),MASK(IBIT))
      GO TO 200
C
C  THIS POINT IS BEING REDRAWN SO ERASE IT.
C  (THIS IS TO CONFORM WITH THE SHADING ALGORITHM, FILLIN.
C  HOWEVER IF BACK TO STARTING POINT DO NOT ERASE
C
  190 IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) RETURN
      ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
C
C
  200 IXOLD = IX
      LRLX = NRLX
      IYOLD = IY
      IBTOLD = IBIT
      IVOLD = IVAL
      RETURN
C
C  PERFORM THIS OPERATION WHEN A CONTOUR STARTS OR ENDS ON A HORIZONTAL
C  LINE.
C
  210 CONTINUE
C
C  ERASE THE FIRST POINT OF A CONTOUR WHEN IT IS PART OF A HORIZONTAL
C  LINE SEGMENT AND IS NOT THE ENDPOINT OF THE SEGMENT
C
      IF (.NOT.IHF) GO TO 220
      ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
  220 CONTINUE
C
C  ERASE THE FIRST POINT OF A HORIZONTAL LINE SEGMENT WHEN IT ENDS
C  THE CONTOUR AND IS NOT THE HIGHEST LINE SEG ON THS SIDE.
C
      IF (.NOT.YCHANG) GO TO 230
      IF (IYDIR .NE. IHS) GO TO 200
  230 ISCR(IHX,IY) = IAND(ISCR(IHX,IY),NMASK(IHB))
      GO TO 200
      END
      SUBROUTINE FILLIN
C
      SAVE
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR7/ IENTRY     ,IONES
C
      IF (IENTRY .EQ. 0) RETURN
C
C  THIS IS A SHADING ALGORITHM IT IS USED TO DETERMINE CONTOUR LINES
C  THAT ARE HIDDEN BY THE PRESENT LINE. THE ALGORITHM PROCESSES
C  HORIZONTAL ROWS. IT ASSUMES THAT THE BIT PATTERN PASSED TO IT
C  HAS ONLY BITS SET TO MARK THE START AND END OF SHADING. THE
C  ALGORITHM ALSO ASSUMES THAT WHEN AN ON BIT IS ENCOUNTERED THAT A
C  CORRESPONDING OFF BIT IS INCLUDED IN THE SAME ROW.
C
C
C  PULL OUT ROWS OF THE CONTOUR PATTERN
C
      IBVAL = 0
      DO  80 IYNOW=1,NY
         DO  40 IXNOW=1,LX
C
C  IF NO ACTIVATED BITS BRANCH
C
            ICRWD = ISCR(IXNOW,IYNOW)
            IF (ICRWD .EQ. 0) GO TO  30
C
C  ACTIVATED BITS IN WORD SET SHADING FLAG
C
C  CHECK BIT BY BIT FOR ON/OFF FLAGS
C
            DO  20 IB=1,NBPW
               IBIT = (NBPW+1)-IB
C
C
C  PULL OUT THE CURRENT GRID POINT VALUE
C
               IVAL = IAND(ICRWD,MASK(IBIT))
C
C  IF IVAL SET, THIS IS AN ON/OFF FLAG
C
               IF (IVAL .EQ. 0) GO TO  10
C
C  FLAG BIT, ALWAYS SET
C
               IBVAL = MOD(IBVAL+1,2)
               GO TO  20
C
C  SHADE THE SCREEN ACCORDING TO THE STATUS OF IBVAL
C
   10          IF (IBVAL .NE. 0) ICRWD = IOR(ICRWD,MASK(IBIT))
C
   20       CONTINUE
C
C  ZERO OUT THE SCREEN
C
            ISCR(IXNOW,IYNOW) = 0
            ISCA(IXNOW,IYNOW) = IOR(ICRWD,ISCA(IXNOW,IYNOW))
            GO TO  40
C
   30       IF (IBVAL .NE. 0) ISCA(IXNOW,IYNOW) = IONES
   40    CONTINUE
C
C  FIX FOR NONCORRECTABLE RUNAWAYS
C
         IF (IBVAL .EQ. 0) GO TO  80
         IBVAL = 0
         DO  70 K=1,LX
            ITEST = 0
            IF (IYNOW .EQ. 1) GO TO  50
            ITEST = ISCA(K,IYNOW-1)
            IF (IYNOW .EQ. NY) GO TO  60
   50       ITEST = IOR(ITEST,ISCA(K,IYNOW+1))
   60       ISCA(K,IYNOW) = ITEST
   70    CONTINUE
C
   80 CONTINUE
      RETURN
      END
      SUBROUTINE DRAWI (IXA,IYA,IXB,IYB)
C
C INCLUDED FOR USE BY PWRZ
C
      SAVE
      CALL FRSTC (IXA,IYA,1)
      CALL FRSTC (IXB,IYB,2)
      RETURN
      END
      SUBROUTINE MMASK
C
C  MAKE THE MACHINE DEPENDENT MASKS USED IN THE CONTOUR DRAWING
C  AND SHADING ALGORITHMS
C
      SAVE
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
     1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
     2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
     3                IHB        ,IHS        ,IHV        ,IVOLD      ,
     4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
     5                IHF
      COMMON /ISOSR9/ BIG        ,IXBIT
      LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
      GENDON = .TRUE.
      NBPW = 16
C
C  GET BIGGEST REAL NUMBER
C
      BIG = R1MACH(2)
C
C  MASKS TO SELECT A SPECIFIC BIT
C
      DO  10 K=1,NBPW
         MASK(K) = ISHIFT(1,K-1)
   10 CONTINUE
C
C  GENERATE  THE BIT PATTERN 177777 OCTAL
C
      ITEMP1 = 0
      ITEMP = MASK(NBPW)
      IST = NBPW-1
      DO  20 K=1,IST
         ITEMP1 = IOR(ITEMP,ISHIFT(ITEMP1,-1))
   20 CONTINUE
      MFIX = IOR(ITEMP1,1)
C
C  MASKS TO CLEAR A SPECIFIC BIT
C
      DO  30 K=1,NBPW
         NMASK(K) = IAND(ITEMP1,MFIX)
         ITEMP1 = IOR(ISHIFT(ITEMP1,1),1)
   30 CONTINUE
      IONES = MFIX
      RETURN
C
C REVISION HISTORY---
C
C JANUARY 1978     DELETED REFERENCES TO THE  *COSY  CARDS AND
C                  ADDED REVISION HISTORY
C JANUARY 1979     NEW SHADING ALGORITHM
C MARCH 1979       MADE CODE MACHINE INDEPENDENT AND CONFORM
C                  TO 66 FORTRAN STANDARD
C JUNE 1979        THIS VERSION PLACED ON ULIB.
C SEPTEMBER 1979   FIXED PROBLEM IN EZISOS DEALING WITH
C                  DETERMINATION OF VISIBILITY OF W PLANE.
C DECEMBER 1979    FIXED PROBLEM WITH PEN DOWN ON CONTOUR
C                  INITIALIZATION IN SUBROUTINE FRSTC
C MARCH            CHANGED ROUTINE NAMES  TRN32I  AND  DRAW  TO
C                  TRN32I  AND  DRAWI  TO BE CONSISTENT WITH THE
C                  USAGE OF THE NEW ROUTINE  PWRZI.
C JUNE  1980       FIXED PROBLEM WITH ZERO INDEX COMPUTATION IN
C                  SUBROUTINE FRSTC.  ADDED INPUT PARAMETER
C                  DIMENSION STATEMENT MISSING IN EZISOS.
C                  FIXED ERROR IN COMPUTATION OF ARCCOSINE
C                  IN EZISOS AND TRN32I.
C DECEMBER 1984    CONVERTED TO GKS LEVEL 0A AND STANDARD FORTRAN 77
C-----------------------------------------------------------------------
C
      END