define FOGR,0  rem  section dealing with graph as a whole
define FSGR,1
define FAGR,2
define TCGR,3
define LCGR,4
define LWGR,5
define PDGR,6
define PSGR,7
define BGGR,8
define PWGR,9
define PHGR,10
define BWGR,11
define SCGR,12
define CIGR,13
define GRGR,14
define ROWS,15
define COLUMNS,16
define TDGR,17
define GRAPHTYPE,18
define XLABEL,19
define YLABELS,20
define THICKNESSGR,26
define FRGR,27
define OLGR,28
define OTGR,29
define ORGR,30
define OBGR,31
define ILGR,32
define ITGR,33
define IRGR,34
define IBGR,35
define MTGR,36
define MSGR,37
define PIGR,38
define PSIZE,39
define FOA,0   rem  Section dealing with axis
define FSA,1
define FAA,2
define TCA,3
define LCA,4
define LWA,5
define STA,6
define LGA,7
define INA,8
define UPA,9
define LOA,10
define NSA,11
define AA,12
define BA,13
define WIDTHA,14
define HEIGHTA,15
define HOTXA,16
define HOTYA,17
define FORMATA,18
define LENGTHA,19
define PAXIS,20
define FOD,0    rem Section dealing with data set
define FSD,1
define FAD,2
define TCD,3
define LWD,4
define LCD,5
define PDD,6
define PSD,7
define FCD,8
define MTD,9
define MSD,10
define MCD,11
define FRD,12
define TDD,13
define PDATA, 14
define FOM,0    rem Section dealing with marker keys - text properties
define FSM,1
define FAM,2
define TCM,3
define LWM,4    rem line properties
define LCM,5
define PDM,6
define PSM,7
define MTM,8    rem marker properties- repeated 5 times
define MSM,9
define MCM,10
define XSM,28
define YSM,29
define XEM,30
define YEM,31
define NSM,32
define MDATA,33
define FOC,0    rem Section dealing with colour keys - text properties
define FSC,1
define FAC,2
define TCC,3
define LWC,4    rem line properties
define LCC,5
define PDC,6
define PSC,7
define FCC,8    rem colour (repeated 16 times)
define XSC,20
define YSC,21
define XEC,22
define YEC,23
define CDATA,24
define LINES_ONLY,0   rem types of 2D graph
define POINTS_ONLY,1
define LINES_AND_POINTS,2
define VERTICAL_HISTOGRAM,3
define HORIZONTAL_HISTOGRAM,4
define STACKED_VERTICAL_HISTOGRAM,5
define STACKED_HORIZONTAL_HISTOGRAM,6
define PICTOGRAM,7
define PIE_CHART,8
define TWO_LINES,9
define BOX_AND_TAILS,10
define DISTIBUTION,11
define GENERAL,12
define AT,0
define BT,1
define CT,2
define DT,3
define PP,4
define X,5
define Y,6
define D,7
define S,8
define CVLENGTH,8
define DRAG_ALLOWED, 0x1
define SELECT_ALLOWED, 0x2
define TEXT_EFFECTS, 0x4
define FILL_EFFECTS, 0x8
define LINE_EFFECTS, 0x10
define GRAPH_COLOUR_ALLOWED, 0x20
define THREE_DEE_ALLOWED, 0x40
define GRID_ALLOWED, 0x80
define EXPLODE_ALLOWED, 0x100
define COLOURS_ALLOWED, 0x200
define BAR_WIDTH_ALLOWED, 0x400
define POINT_DETAILS_ALLOWED, 0x800
define TICKS_ALLOWED, 0x1000
define DIMENSIONS_ALLOWED, 0x2000
define TOP_ALLOWED, 0x4000
define DELETE_ALLOWED, 0x8000
define ALL_SELECTED, 0x10000
define NONE_SELECTED, 0x20000
define REDRAW_NEEDED, 0x40000


rem - Contour lines

rem - This module is given an array of heights at equally spaced points
rem - on a grid.  It plots a contour map of the surface implied by the 
rem - heights.

rem - Method :

rem - 1.  The data is scanned and the extreme values found.  Then the
rem - standard 'scale' function is used to find about 10 'round' values
rem - between the extremes to plot as contours.

rem - Each contour line is plotted in a different colour. The process is
rem - as follows:

rem -     The surface is deemed to consist of plane triangles, arranged
rem - thus:

rem -        * ----- * ----- * ----- *
rem -        | \   / | \   / | \   / |
rem -        |  \ /  |  \ /  |  \ /  |
rem -        |   x   |   x   |   x   |
rem -        |  / \  |  / \  |  / \  |
rem -        | /   \ | /   \ | /   \ |
rem -        * ----- * ----- * ----- *
rem -        | \   / | \   / | \   / |
rem -        |  \ /  |  \ /  |  \ /  |
rem -        |   x   |   x   |   x   |
rem -        |  / \  |  / \  |  / \  |
rem -        | /   \ | /   \ | /   \ |
rem -        * ----- * ----- * ----- *

rem - where the *'s are the grid points, and the x's are centre points
rem - whose values are the average of the four sorrounding grid points
rem - There are two types of triangle:
rem - top,right,bottom and left.

rem - We use symbols for eight directions of movement, as follows.
rem -   
rem -       N for north
rem -       A for northeast
rem -       E for east
rem -       B for southeast
rem -       S for south
rem -       C for southwest
rem -       W for west    
rem -       D for northwest

rem - a top triangle can be entered moving S, A,D and left moving N,B,C
rem - a right triangle can be entered moving W,A,B and left moving E,C,D
rem - a bottom triangle can be entered moving N,B,C and left moving S,A,D
rem - a left triangle can be entered moving E,C,D and left moving W,A,B 


rem - A given contour line may or may not cross a triangle. If it does, it 
rem - will have a point and a side of entry, and a point and a side of 
rem - exit. The way these are allocated is arbitrary.

rem - It is easy to see that any contour line must either be a closed
rem - loop, or join two points on the edge of the picture.  We regard
rem - a line which intersects the edge as part of a closed loop which has
rem - 1-4 straight sections round the edge; the direction being such as to
rem - ensure that the loop encloses higher ground.

rem - To plot a given contour, we proceed as follows:

rem - a)   We set up an array of marker words, one for each square cell.
rem - The contents are bit-significant, with one bit for each of the four
rem - triangles in the cell : 
rem - Bit 0 - top; bit 1 - right; bit 2 - bottom; bit 3 - left
rem - Each marker will record whether contour passes through that triangle.
rem - Initially all markers are cleared to false.

rem - b)   We repeat the following sequence:
rem -      
rem -      Find an unvisited triangle through which the contour passes.  
rem - If such a triangle cannot be found, the plot of this contour is 
rem - complete. If a triangle is found, then the finder function returns 
rem - various interesting facts:
rem -  
rem -    The identity of the triangle (X and Y coordinates)
rem -    The position of an exit point
rem -    The side of the exit point.

rem -    If the exit point is at an edge of the picture, the finder 
rem -    indicates whether to turn Left or Right to search for the 
rem -    re-entry point.

rem -    When the first triangle is found, the system opens up a path
rem -    record, plants the exit point in it, and marks the triangle as 
rem -    visited.

rem -    The next task is to track the contour line until the loop closes.

rem -    If the exit point is not a picture edge, it immediately identifies
rem -    the next triangle.  Given the entry point, a function delivers the
rem -    corresponding exit point, and so on.

rem -    If the exit point is a picture edge, the edge of the picture is
rem -    scanned until a re-entry point is found. It may be necessary to 
rem -    scan round several edges of the picture.

rem -    When the scan returns to the original triangle, the system closes 
rem -    the polygon and writes it to the output file.

rem -    Several of the functions described below return Exit Point 
rem -    Descriptors.   This is a string, arranged as follows:
rem -
rem -    First one of the six letters NEWSABCD to indicate the direction 
rem -    of exit
rem -    Second a letter : 
rem -       L means : You have crossed the edge of the graph, Go left
rem -       R means : you have crossed the edge of te graph. Go right
rem -       C means : You have not crossed the edge of the graph
rem -    Next, the sequence X=nnn, where nnn is the x coordinate of the 
rem -    crossing point (assuming unit spacing of points)
rem -    Next, the sequence Y=nnn, where nnn is the y-coordinate of the 
rem -    crossing point (assuming unit spacing of points)
rem -
rem -    A null exit point is coded as "0"    
 

rem -    We use a special structure called c to pass round various values 
rem -    to dowith contour lines.   The elements of the structure are as 
rem -    follows:
rem -        AT, BT, maps unit x-xoordinates on to the picture, using the
rem -        transformation  x' = x*c(AT) + c(BT)

rem -        CT,DT are used similarly for the y coordinates.  Note that 
rem -        CT is negative because screen coordinates go up but array 
rem -        coordinates go down

rem -        X and Y are the coordinates of the current triangle (in array
rem -         coordinates)

rem -        D is 0 if the current triangle is upper; 1 if right;
rem -        2 if bottom amd 3 if left
rem -        PP is a pointer for the current path

        

rem - ******************************************************************

rem - macro top-cross examines a given upper triangle to see if it is
rem - crossed by a contour with the stated value. If so it returns an
rem - Exit Point descriptor.

rem - parameters are :   d = data matrix (r rows by c columns)
rem -       m - intermediate matrix (r-1 rows by c-1 columns
rem - coordinates of current square (c(Y),c(X))
rem - contour value - val
rem - a code showing direction of entry as S,A,D or X (unknown) - ep

macro top_cross (d(),m(),c(),rows,cols,val,ep)
local p,q,r,ex,ey,edge,yy,xx,delta
   yy = c(Y)
   xx = c(X)
   p = d(yy,xx)-val
   q = d(yy,xx+1)-val
   r = m(yy,xx)-val
   if (p*q < 0) and (ep <> "S") then
                   rem North exit point
       ex = xx + p/(p-q)
       ey  = yy
       edge = if yy >0 then "NCX=" else if p>q then "NLX=" else "NRX=" endif endif
       = edge cat ex cat "Y=" cat ey
   endif
   if (p*r < 0) and (ep <>"A") then
                   rem SW (C) exit point
       delta= 0.5*p/(p-r)
       ex = xx + delta
       ey  = yy + delta
       = "CCX=" cat ex cat "Y=" cat ey
   endif
   if ( q*r < 0) then
                    rem Southeast exit point
      delta = 0.5*q/(q-r)
      ex = xx+1-delta
      ey = yy + delta
      = "BCX=" cat ex cat "Y=" cat ey
   endif
   ="0"
endmacro

rem - ****************************************************************

rem - macro right-cross examines a given right triangle to see if it is
rem - crossed by a contour with the stated value. If so it returns an
rem - Exit Point descriptor.

rem - parameters are :   d = data matrix (r rows by c columns)
rem -       m - intermediate matrix (r-1 rows by c-1 columns
rem - coordinates of current square (c(Y),c(X))
rem - contour value - val
rem - a code showing direction of entry as W,A,B or X (unknown) - ep

macro right_cross (d(),m(),c(),rows,cols,val,ep)
local p,q,r,ex,ey,edge,yy,xx,delta
   yy = c(Y)
   xx = c(X)
   p = d(yy,xx+1)-val
   q = d(yy+1,xx+1)-val
   r = m(yy,xx)-val
   if (p*q < 0) and (ep <> "W") then
                   rem East exit point
       ex = xx + 1
       ey  = yy + p/(p-q)
       edge = if xx < cols-2 then "ECX=" else if p>q then "ELX=" else "ERX=" endif endif
       = edge cat ex cat "Y=" cat ey
   endif
   if (p*r < 0) and (ep <>"B") then
                   rem NW (D) exit point
       delta= 0.5*p/(p-r)
       ex = xx + 1-delta
       ey  = yy + delta
       = "DCX=" cat ex cat "Y=" cat ey
   endif
   if ( q*r < 0) then
                    rem SouthWest exit point
      delta = 0.5*q/(q-r)
      ex = xx+1-delta
      ey = yy +1- delta
      = "CCX=" cat ex cat "Y=" cat ey
   endif
   ="0"
endmacro
         


rem - ****************************************************************

rem - macro bottom-cross examines a given bottom triangle to see if it is
rem - crossed by a contour with the stated value. If so it returns an
rem - Exit Point descriptor.

rem - parameters are :   d = data matrix (r rows by c columns)
rem -       m - intermediate matrix (r-1 rows by c-1 columns
rem - coordinates of current square (c(Y),c(X))
rem - contour value - val
rem - a code showing direction of entry as N,B, C or X (unknown) - ep

macro bottom_cross (d(),m(),c(),rows,cols,val,ep)
local p,q,r,ex,ey,edge,yy,xx,delta
   yy = c(Y)
   xx = c(X)
   p = d(yy+1,xx)-val
   q = d(yy+1,xx+1)-val
   r = m(yy,xx)-val
   if (p*q < 0) and (ep <> "N") then
                   rem South exit point
       ex = xx + p/(p-q)
       ey  = yy + 1
       edge = if yy < rows-2 then "SCX=" else if p>q then "SRX=" else "SLX=" endif endif
       = edge cat ex cat "Y=" cat ey
   endif
   if (p*r < 0) and (ep <>"B") then
                   rem NW (D) exit point
       delta= 0.5*p/(p-r)
       ex = xx + delta
       ey  = yy + 1-delta
       = "DCX=" cat ex cat "Y=" cat ey
   endif
   if ( q*r < 0) then
                    rem Northeast exit point
      delta = 0.5*q/(q-r)
      ex = xx+1-delta
      ey = yy +1- delta
      = "ACX=" cat ex cat "Y=" cat ey
   endif
   ="0"
endmacro


rem - ****************************************************************

rem - macro left-cross examines a given bottom triangle to see if it is
rem - crossed by a contour with the stated value. If so it returns an
rem - Exit Point descriptor.

rem - parameters are :   d = data matrix (r rows by c columns)
rem -       m - intermediate matrix (r-1 rows by c-1 columns
rem - coordinates of current square (c(Y),c(X))
rem - contour value - val
rem - a code showing direction of entry as E, C D or X (unknown) - ep

macro left_cross (d(),m(),c(),rows,cols,val,ep)
local p,q,r,ex,ey,edge,yy,xx,delta
   yy = c(Y)
   xx = c(X)
   p = d(yy,xx)-val
   q = d(yy+1,xx)-val
   r = m(yy,xx)-val
   if (p*q < 0) and (ep <> "E") then
                   rem West exit point
       ex = xx
       ey  = yy + p/(p-q)
       edge = if xx > 0 then "WCX=" else if p>q then "WRX=" else "WLX=" endif endif
       = edge cat ex cat "Y=" cat ey
   endif
   if (p*r < 0) and (ep <>"C") then
                   rem NE (A) exit point
       delta= 0.5*p/(p-r)
       ex = xx + delta
       ey  = yy + delta
       = "ACX=" cat ex cat "Y=" cat ey
   endif
   if ( q*r < 0) then
                    rem Southeast exit point
      delta = 0.5*q/(q-r)
      ex = xx+delta
      ey = yy +1- delta
      = "BCX=" cat ex cat "Y=" cat ey
   endif
   ="0"
endmacro


rem - ****************************************************************
rem - This macro looks for a triangle crossed by the contour at val.
rem - The parameters are :
rem -           d()  The data block
rem -           m()  The intermediate matrix
rem -           c()  The plot parameters
rem -           rows, cols:  The size of the data block
rem -           up()  Registration array for triangles
rem -           val : The contour being sought
rem - If the function finds a triangle which has not already been crossed
rem - it puts its details into c(X, c(Y) and c(D) and returns an exit
rem - descriptor.  
rem - If not found, it returns "0" (with c(X), c(Y) and c(D) undefined)

macro find_triangle (d(),m(),c(),rows,cols, up(),val)
local x,y,q
   for y = 0 to rows-2
       c(Y) = y
       for x = 0 to cols-2
           c(X)=x
           if (up(y,x) & 1) = 0 then
              c(D) = 1
              q =  top_cross(d,m,c,rows,cols,val,"X")
              if q <> "0" then =q
           endif
           if (up (y,x) & 2) = 0 then
               c(D)=2
               q = right_cross(d,m,c,rows,cols,val,"X")
               if q <> "0" then =q
           endif
           if (up (y,x) & 4) = 0 then
               c(D)=4
               q = bottom_cross(d,m,c,rows,cols,val,"X")
               if q <> "0" then =q
           endif
           if (up (y,x) & 8) = 0 then
               c(D)=8
               q = left_cross(d,m,c,rows,cols,val,"X")
               if q <> "0" then =q
           endif
       next x
   next y
="0"
endmacro 

rem : *************************************************************
rem  This macro adds a straight-line segment to the current path
rem  to got to position xx,yy (which are given in array coordinates)

macro grow_path (path(),c(),xx,yy)
   path(c(PP)) = 8
   path(c(PP)+1) = xx*c(AT)+c(BT)
   path(c(PP)+2) = yy*c(CT)+c(DT)
   c(PP) = c(PP)+3
   =0
endmacro


rem : **********************************************************
rem :  This function is entered with an exit point descriptor
rem :  which crosses the edge of the plot.  The function searches
rem :  round the edge until it finds a re-entry point and returns
rem :  a dummy exit descriptor for this point. It also adds between
rem :  1 and 4 line segments to the current output path, and returns
rem :  the new value of the path pointer
rem : 
rem :  Also planted in c(X), c(Y) and c(D) is the identity of the 
rem : triangle just found
  
macro track_edge(d(),c(),rows,cols,path(),val,ed)
local en ,xx,yy,j,k,res,w
   en = left(ed,2)
   xx = floor(extract(ed,"X="))
   yy = floor(extract(ed,"Y="))
   res = grow_path(path,c,xx,yy)
   k = 2*cols+2*rows
   for j = 0 to k
             rem set maximum length of search
NL:
             rem search left along top edge
      if en= "NL" then
      xx =xx-1
      if xx < 0 then
            en="WL"     :      rem turn corner to left
            res = grow_path(path,c,0,0)
            yy=-1
            xx = 0
            goto WL
         endif
         if d(0,xx) < val then
            w = xx+(val-d(0,xx))/(d(0,xx+1)-d(0,xx))
            res = grow_path(path,c,w,0)
            c(Y)=0 : c(X)=xx: c(D)=1
            = "SCX=" cat w cat "Y=" cat 0
         endif
      endif
WL:
           rem search down along left edge
      if en= "WL" then
         yy=yy+1
         if yy = rows-1 then
            en="SL"           :rem turn corner
            res = grow_path(path,c,0,yy)
            yy=rows-1
            xx=-1
            goto SL
         endif
         if d(yy+1,0) < val then
            w = yy+1-(val-d(yy+1,0))/(d(yy,0)-d(yy+1,0))
            res = grow_path(path,c,0,w)
            c(Y)=yy: c(X) = 0: c(D) = 8
            = "ECX=" cat 0 cat "Y=" cat w
         endif
      endif
SL:
   rem search right along bottom edge
   if en= "SL" then
         xx=xx+1
         if xx = cols-1 then
            en="EL"
            res = grow_path(path,c,cols-1,rows-1)
            yy = rows 
            goto EL
         endif
         if d(yy,xx+1) < val then
            w = xx+1-(val-d(yy,xx+1))/(d(yy,xx)-d(yy,xx+1))
            res = grow_path(path,c,w,yy)
            c(Y)=yy-1: c(X) = xx: c(D) = 4
            = "NCX=" cat w cat "Y=" cat yy
         endif
      endif
EL:
   rem search up along right edge
      if en= "EL" then
         if yy = 0 then
            xx = cols-1
            en="NL"
            res = grow_path(path,c,xx,0)
            goto NL
         endif
         yy=yy-1
         if d(yy,xx) < val then
            w = yy+ (val-d(yy,xx))/(d(yy+1,xx)-d(yy,xx))
            res = grow_path(path,c,xx,w)
            c(Y)=yy: c(X) = xx-1: c(D) = 2
            = "WCX=" cat xx cat "Y=" cat w
         endif
      endif
NR:
    rem search right along top
      if en= "NR" then
         xx=xx+1
         if xx = cols-1 then
            en="ER"
            res = grow_path(path,c,xx,0)
            yy = -1
            goto ER
         endif
         if d(0,xx+1) < val then
            w = xx+1-(val-d(0,xx+1))/(d(0,xx)-d(0,xx+1))
            res = grow_path(path,c,w,0)
            c(Y)=0: c(X) = xx: c(D) = 1
            = "SCX=" cat w cat "Y=" cat 0
         endif
      endif
ER:
  rem search down along right side
      if en= "ER" then
         yy=yy+1
         if yy = rows-1 then
            en="SR"
            res = grow_path(path,c,xx,yy)
            xx=cols
            goto SR
         endif
         if d(yy+1,xx) < val then
            w = yy+1-(val-d(yy+1,xx))/(d(yy,xx)-d(yy+1,xx))
            res = grow_path(path,c,xx,w)
            c(Y)=yy: c(X) = xx-1: c(D) = 2
            = "WCX=" cat xx cat "Y=" cat w
         endif
      endif
SR:
  rem search left along bottom
      if en= "SR" then
         if xx <  1 then
            xx=0
            en="WR"
            res = grow_path(path,c,0,rows-1)
            goto WR
         endif
         xx=xx-1
         if d(yy,xx) < val then
            w = xx+(val-d(yy,xx))/(d(yy,xx+1)-d(yy,xx))
            res = grow_path(path,c,w,yy)
            c(Y)=yy-1: c(X) = xx: c(D) = 4
            = "NCX=" cat w cat "Y=" cat yy
         endif
      endif
WR:
      rem search up along left side
      if en= "WR" then
         if yy = 0 then
            en="NR"
            res = grow_path(path,c,0,0)
            xx = -1
            goto NR
         endif
         yy=yy-1
         if d(yy,xx) < val then
            w = yy+ (val-d(yy,xx))/(d(yy+1,xx)-d(yy,xx))
            res = grow_path(path,c,xx,w)
            c(Y)=yy: c(X) = xx: c(D) = 8
            = "ECX=" cat xx cat "Y=" cat w
         endif
      endif
   next j
   = "track-edge failure"
endmacro


rem : ***************************************************************
rem : The main contour-plotting function 

macro contour_plot(a(),s,v,p(),height,width)
local triangles(height-1,width-1), mid_points(height-1,width-1)
local path(5*height*width+50) 
local val,minval,maxval
local res, x,y,w,az
local mincon, maxcon, stepcon
local c(CVLENGTH)
local gr_tag,z,ww
local tstart, xstart, ystart
local data_up, data_right, data_width, data_height     
            rem scan the data and reject eveything if any data points
            rem are missing
    az = 0
    minval = a(0,0)
    maxval = a(0,0)
    for y = 0 to height-1
       for x = 0 to width-1
          if type(a(y,x)) > 2 then  ="Invalid data"
          if a(y,x) > maxval then maxval = a(y,x)
          if a(y,x) < minval then minval = a(y,x)
       next x
    next y
    for y = 0 to height-2
       for x = 0 to width-2
          mid_points(y,x) = (a(y,x)+a(y,x+1)+a(y+1,x)+a(y+1,x+1))/4
       next x
    next y
               rem temporary
    mincon = 80.01
    maxcon =120.01
    stepcon = 10
    
              rem emd of temporary
              rem now open window and set transform data :
              rem        xx = atx+bt,
              rem        yy = ctx + dt

   gr_tag = psexsub(v)
                                        rem open window
   res = settag(gr_tag,77)
   w = gstart(p(PWGR),p(PHGR))
   data_right = p(OLGR) + p(ILGR)
   data_up = p(OBGR) + p(IBGR)
   data_width = p(PWGR) - data_right - p(ORGR)-p(IRGR)
   data_height = p(PHGR) - data_up - p(OTGR)-p(ITGR)
   c(AT) = data_width/width
   c(BT) = data_right
   c(CT) = -data_height/height
   c(DT) = data_up + data_height
   res = lineattributes(w,1,1,0,0)
   res = line(w,0*c(AT)+c(BT),0*c(CT)+c(DT),(width-1)*c(AT),0)
   res = line(w,(width-1)*c(AT)+c(BT),0*c(CT)+c(DT),0,(height-1)*c(CT))
   res = line(w,(width-1)*c(AT)+c(BT),(height-1)*c(CT)+c(DT),-(width-1)*c(AT),0)
   res = line(w,0*c(AT)+c(BT),(height-1)*c(CT)+c(DT),0,-(height-1)*c(CT))
   res = lineattributes(w,1,1,0xE0E0E000,0)
   for x=1 to width-2
      res = line(w,x*c(AT)+c(BT),c(DT),0,(height-1)*c(CT))
   next x
   for y=1 to height-2
      res = line(w, c(BT), y*c(CT)+c(DT),(width-1)*c(AT),0)
   next y

   for val = mincon to maxcon step stepcon
      c(PP) = 0       rem set path pointer
                      rem clear upper and lower triangles
      for y = 0 to height-2
         for x = 0 to width-2
           triangles(y,x) = 0
         next x
      next y

      while (1)
      c(PP)=0
       z = find_triangle (a,mid_points,c,height,width,triangles,val)
              if z = "0" then goto loop
                     rem a bit of contour located. start path
       res = grow_path(path,c,extract(z,"X="),extract(z,"Y="))
       path(0) = 2
       xstart = c(X)
       ystart = c(Y)
       tstart = c(D)
           triangles(ystart,xstart) =  triangles(ystart,xstart) or c(D)
                  rem now we find the next triangle
inner:
          ww = left(z,1)
          if (find(z,"L")>=0) or (find(z,"R")>=0) then
                 z = track_edge(a,c,height,width,path,val,z)
          ww = left(z,1)

                goto horrible
          endif
              rem syspect if then else doesn't work here 

             if ww = "A" then  c(D) = if c(D)=8 then 1 else 2 endif
             if ww = "B" then  c(D) = if c(D)=1 then 2 else 4 endif
             if ww = "C" then  c(D) = if c(D)=2 then 4 else 8 endif
             if ww = "D" then  c(D) = if c(D)=4 then 8 else 1 endif

             if ( ww = "N") then c(Y) = c(Y)-1: c(D) = 4
             if ( ww = "W") then c(X) = c(X)-1: c(D) = 2
             if ( ww = "S") then c(Y) = c(Y)+1: c(D) = 1
             if ( ww = "E") then c(X) = c(X)+1: c(D) = 8
          
                 rem now find exit point to next triangle
horrible:

          
          triangles(c(Y),c(X)) = triangles(c(Y),c(X)) or c(D)
          if c(D) = 1 then
          z = top_cross(a,mid_points,c,height,width,val,ww)
          endif
          if c(D) = 2 then
             z = right_cross(a,mid_points,c,height,width,val,ww)
          endif
          if c(D) = 4 then
             z = bottom_cross(a,mid_points,c,height,width,val,ww)
          endif
          if c(D) = 8 then
             z = left_cross(a,mid_points,c,height,width,val,ww)
          endif
   
          res=grow_path(path,c,extract(z,"X="),extract(z,"Y="))
        if (xstart <> c(X) or ystart <> c(Y) or tstart <> c(D)) then goto inner

          path(c(PP)) = 5
          path(c(PP)+1) = 0
          c(PP) = c(PP) +2
         res = lineattributes(w,1,1,0,0)
         res = drawobject(w,path, c(PP), 0x01010100 * val)         
     endwhile
loop:           
   next val
   =w   
endmacro

rem : *****************************************************************
rem : This is the high-level macro called to make a contour plot.
rem : The data consists entirely of numbers; no labels allowed
rem : anywhere
rem : The graph may of course be post-edited to label it.
macro contours(a(),s)
graphmacro "Contour",0x280,""
   local v,q,n,res,h,w,j
   local p(PSIZE)

   v = psregister(s)
   q = psstep(v)
   if iserror(q) then = q
   if q<> code("A") then = "Missing parameter substring"
   res = get_gr_params(p,v)
   p(GRAPHTYPE) = 0
   h = first(a)
   w = second(a)
   res = contour_plot(a,s,v,p,h,w)
   q= close_string(v)
   = res
endmacro
