C                  INTERNATIONAL AVS CENTER
C        (This disclaimer must remain at the top of all files)
C
C WARRANTY DISCLAIMER
C
C This module and the files associated with it are distributed free of charge.
C It is placed in the public domain and permission is granted for anyone to use,
C duplicate, modify, and redistribute it unless otherwise noted.  Some modules
C may be copyrighted.  You agree to abide by the conditions also included in
C the AVS Licensing Agreement, version 1.0, located in the main module
C directory located at the International AVS Center ftp site and to include
C the AVS Licensing Agreement when you distribute any files downloaded from
C that site.
C
C The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module provide absolutely
C NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
C the quality and performance of this software is with the user.  IN NO EVENT
C WILL The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module BE LIABLE TO
C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.
C
C This AVS module and associated files are public domain software unless
C otherwise noted.  Permission is hereby granted to do whatever you like with
C it, subject to the conditions that may exist in copyrighted materials. Should
C you wish to make a contribution toward the improvement, modification, or
C general performance of this module, please send us your comments:  why you
C liked or disliked it, how you use it, and most important, how it helps your
C work. We will receive your comments at avs@ncsc.org.
C
C Please send AVS module bug reports to avs@ncsc.org.
C
C  NOTE:  THIS MODULE AND SOURCE CODE IS FOR USE 
C  WITH THE AVS SOFTWARE ENVIRONMENT ONLY 
      SUBROUTINE CELL_LINES( FLEVEL, FLAG3D, X,Y,Z, F, N, XPT,YPT,ZPT )
*
*--   Function to find contour lines in a grid cell 
*--   using a marching-squares algorithm, i.e. put a recipe for every ca
     $se
*---  in a look up table then calculate an index into the table
*
*--   Arguments
*
*--     input contour level, x and y coordinates, vertex field values
        REAL     FLEVEL, X(4), Y(4), Z(4), F(4)
        LOGICAL  FLAG3D
*
*--     number of lines returned, either 0,1,2
        INTEGER  N
*
*--     arrays for the line output indexed for ( line number, vertex num
     $ber )
        REAL     XPT(2,2), YPT(2,2), ZPT(2,2)
*
*--   Local variables
*
        INTEGER  K, J
        REAL     AVF, Z1, Z3, Z5, Z7
        REAL     F1,X1,Y1, F3,X3,Y3, F5,X5,Y5, F7,X7,Y7, FACTOR
*
*--   Look up table
*
*--     devise a recipe encoding
*--     each point is an interpolation between two vertices
*--     vertex labelling starts at the bottom left 1 (xmin,ymin) then
*--     proceeds anti-clockwise around to 4 (xmin,ymax)
*--     the recipe will be a number of lines then twice that number of p
     $airs of
*--     vertices for interpolation
*
        INTEGER LUT(0:8,0:15)
*
        DATA  (LUT(I,0),  I=0,0)   / 0  /              
        DATA  (LUT(I,1),  I=0,4)   / 1, 3,4, 4,1  /
        DATA  (LUT(I,2),  I=0,4)   / 1, 2,3, 3,4  /
        DATA  (LUT(I,3),  I=0,4)   / 1, 2,3, 4,1  /
        DATA  (LUT(I,4),  I=0,4)   / 1, 1,2, 2,3  /
        DATA  (LUT(I,5),  I=0,8)   / 2, 1,2, 2,3, 3,4, 4,1 /
        DATA  (LUT(I,6),  I=0,4)   / 1, 1,2, 3,4  /
        DATA  (LUT(I,7),  I=0,4)   / 1, 1,2, 4,1  /
        DATA  (LUT(I,8),  I=0,4)   / 1, 1,2, 4,1  /
        DATA  (LUT(I,9),  I=0,4)   / 1, 1,2, 3,4  /
        DATA  (LUT(I,10), I=0,8)   / 2, 1,2, 4,1, 2,3, 3,4 /
        DATA  (LUT(I,11), I=0,4)   / 1, 1,2, 2,3  /
        DATA  (LUT(I,12), I=0,4)   / 1, 2,3, 4,1  /
        DATA  (LUT(I,13), I=0,4)   / 1, 2,3, 3,4  /
        DATA  (LUT(I,14), I=0,4)   / 1, 3,4, 4,1  /
        DATA  (LUT(I,15), I=0,0)   / 0  /
*
*
*--  Calulate index
*
      J = 0
      DO 100 K=1,4
        IF (F(K).GT.FLEVEL)  J = J + 2**(4-K)
  100 CONTINUE
*
*--   test for saddles
      IF (J.EQ.5) THEN
        AVF = (F(1)+F(2)+F(3)+F(4)) / 4.0
        IF (AVF.GT.FLEVEL) J = 10
      ELSE IF (J.EQ.10) THEN
        AVF = (F(1)+F(2)+F(3)+F(4)) / 4.0
        IF (AVF.GT.FLEVEL) J = 5
      END IF
*
*--   Draw the lines
*
*--   number of lines
      N = LUT(0,J)
*
*--   first line
      IF (N.GT.0) THEN
*
*--     first interpolation 
        F1 = F(LUT(1,J))
        FACTOR = ( FLEVEL - F1 ) / ( F(LUT(2,J)) - F1 )
        X1 = X(LUT(1,J))
        XPT(1,1) = X1 + FACTOR * ( X(LUT(2,J)) - X1 )
        Y1 = Y(LUT(1,J))
        YPT(1,1) = Y1 + FACTOR * ( Y(LUT(2,J)) - Y1 )
        IF (FLAG3D) THEN
          Z1 = Z(LUT(1,J))
          ZPT(1,1) = Z1 + FACTOR * ( Z(LUT(2,J)) - Z1 )
        END IF
* 
*--     second interpolation 
        F3 = F(LUT(3,J))
        FACTOR = ( FLEVEL - F3 ) / ( F(LUT(4,J)) - F3 )
        X3 = X(LUT(3,J))
        XPT(1,2) = X3 + FACTOR * ( X(LUT(4,J)) - X3 )
        Y3 = Y(LUT(3,J))
        YPT(1,2) = Y3 + FACTOR * ( Y(LUT(4,J)) - Y3 )
        IF (FLAG3D) THEN
          Z3 = Z(LUT(3,J))
          ZPT(1,2) = Z3 + FACTOR * ( Z(LUT(4,J)) - Z3 )
        END IF
*
      END IF
*
*--   second line
      IF (N.EQ.2) THEN
*
*--     first interpolation 
        F5 = F(LUT(5,J))
        FACTOR = ( FLEVEL - F5 ) / ( F(LUT(6,J)) - F5 )
        X5 = X(LUT(5,J))
        XPT(2,1) = X5 + FACTOR * ( X(LUT(6,J)) - X5 )
        Y5 = Y(LUT(5,J))
        YPT(2,1) = Y5 + FACTOR * ( Y(LUT(6,J)) - Y5 )
        IF (FLAG3D) THEN
          Z5 = Z(LUT(5,J))
          ZPT(2,1) = Z5 + FACTOR * ( Z(LUT(6,J)) - Z5 )
        END IF
* 
*--     second interpolation
        F7 = F(LUT(7,J))
        FACTOR = ( FLEVEL - F7 ) / ( F(LUT(8,J)) - F7 )
        X7 = X(LUT(7,J))
        XPT(2,2) = X7 + FACTOR * ( X(LUT(8,J)) - X7 )
        Y7 = Y(LUT(7,J))
        YPT(2,2) = Y7 + FACTOR * ( Y(LUT(8,J)) - Y7 )
        IF (FLAG3D) THEN
          Z7 = Z(LUT(7,J))
          ZPT(2,2) = Z7 + FACTOR * ( Z(LUT(8,J)) - Z7 )
        END IF
*
      END IF
*
      RETURN
      END