REM CONTOUR  18 Jan 1987               

REM  This is a simplified contour plotting program
REM  for microcomputer--see BYTE Nov 1983 p487
REM  Plots from file where row 0 & col 0 have centre
REM  values of array elements.
 
DIM SHARED f(50,50)
DIM SHARED c(20)
INPUT "Name of file ";N$
IF LEN(N$)<=1 GOTO fin

m0=-999

REM file read

getfile:
REM el painful way to read in data, but it doesn't make too many assumptions!
OPEN N$ FOR INPUT AS #1
iMax=0:jMax=0
FOR j=0 TO 50
  p0=1:i=0
  LINE INPUT #1,a$
  a$=a$+" "
  i1=LEN(a$)
  IF i1<=1 GOTO done
    fl=0
    FOR p=1 TO i1
      IF MID$(a$,p,1)<>" " THEN 
        fl=1
      ELSEIF fl<>0 THEN
        f(i,j)=VAL(MID$(a$,p0,p-p0))
        i=i+1:p0=p:fl=0
      END IF
    NEXT p
    IF i>iMax THEN iMax=i-1
NEXT j
done:
jMax=j-1
CLOSE #1
 

printarray:            
FOR j=0 TO jMax         'J corresp to y-axis
  FOR i=0 TO iMax       'I corresp to x-axis
   x=f(i,j)
   IF x=m0 THEN
     PRINT "    ";
   ELSE
     PRINT USING "####";x;
   END IF
  NEXT i
  PRINT
NEXT j  

getsize:
PRINT
fMax=-1E+10:fMin=1E+10
FOR j=1 TO jMax
  FOR i=1 TO iMax
    x=f(i,j)
    IF fMax<x THEN fMax=x
    IF fMin>x THEN fMin=x
  NEXT i
NEXT j  

setup:
PRINT "highest=";fMax;"lowest=";fMin
PRINT
LINE INPUT "Name of X-scale " ; x$
LINE INPUT "Name of Y-scale " ; y$
INPUT "Min, Max X-scale to plot "; xMin, xMax
INPUT "Min, Max Y-scale to plot "; yMin, yMax
INPUT "How many contours " ; C1 : C1=C1+1
c(1)=fMin-1 : x=(fMax-fMin)/C1
FOR i=2 TO C1 : c(i)=c(i-1)+x : NEXT i

pr1:
PRINT "Contour values :- ";
FOR i=1 TO C1 : PRINT c(i); : NEXT i
PRINT:PRINT "OK ? ";
LINE INPUT a$ : a$=MID$(a$,1,1): IF a$<>"N" AND a$<>"n" GOTO pr2
FOR i=1 TO C1
  LINE INPUT "Contour?";a$
  IF LEN(a$)>=1 THEN c(i)=VAL(a$)
NEXT i
GOTO pr1

pr2:  
INPUT "X, Y Courseness [0..1]";s2,s3


con:
HCelSiz= (f(iMax,0) - f(1,0))/(iMax)
VCelSiz=-(f(0,jMax) - f(0,1))/(jMax)
xLo=INT(1+(xMin-f(1,0))/HCelSiz)
xHi=INT(1+(xMax-f(1,0))/HCelSiz)
yLo=jMax-INT((yMax-f(0,jMax))/VCelSiz)+1
yHi=jMax-INT((yMin-f(0,jMax))/VCelSiz)+1

axis:   

CLS

Y0=(0-f(0,1))/(f(0,jMax)-f(0,1))
IF xLo>=xHi THEN SWAP xLo,xHi
IF yLo>=yHi THEN SWAP yLo,yHi
sx= 550/(xHi-xLo)
sy= 150/(yHi-yLo)
xb=50 : yb= 16

LINE (50,170)-(50+(xHi-xLo)*sx,170),1
LINE (50,170)-(50, 170-(yHi-yLo)*sy),1
FOR i=1 TO (xHi-xLo)
  LINE (50+i*sx,170) - (50+i*sx,165),1
NEXT i
FOR i=1 TO (yHi-yLo)
  LINE (50,170-i*sy) - (55,170-i*sy),1
NEXT i

LOCATE 20,30 : PRINT x$
LOCATE 12,1 : PRINT y$
FOR i=xLo TO xHi : LOCATE 23,5+(i-xLo)*sx/8 
  IF i>0 AND i<=iMax THEN PRINT f(i,0);
NEXT i
FOR i=yLo TO yHi : LOCATE 22-(yHi-i)*sy/8, 3
  IF i>0 AND i<=jMax THEN PRINT f(0,i); 
NEXT i
LOCATE 1,1

draw:
CALL plot(xLo,xHi,yLo,yHi)
END




SUB plot(xLo,xHi,yLo,yHi) STATIC
  SHARED iMax,jMax,sx,sy,s2,s3,m0,xb,yb,C1
  c0 = c(1)
  FOR i=xLo TO xHi
    IF i<0 OR i>=iMax GOTO confin1
    FOR j=yLo TO yHi
      IF j<0 OR j>=jMax GOTO confin
   
      REM define 4 corners of a cell

      x1=f(i,j)
      x2=f(i,j+1)  
      x3=f(i+1,j)
      x4=f(i+1,j+1)

      REM if all 4 corners are less than lowest contour or
      REM any corner has a missing value,go to next cell

      IF x1<c0 AND x2<c0 AND x3<c0 AND x4<c0 GOTO confin
      IF x1=m0 OR x2=m0 OR x3=m0 OR x4=m0 GOTO confin

      REM  I-direction interpolate over cell

      FOR k=0 TO 1-s2 STEP s2
        z1=x1-k*(x1-x3)
        z2=x2-k*(x2-x4)
        IF z1<c0 AND z2<c0 GOTO con1
        GOSUB cross

        REM CROSS gets subscripts of contours crossed

        IF c4<c3 GOTO con1
        m=z2-z1
        r2=i+k
        b=z1-m*j

        REM find R1, the j-direction crossing coordinate

        FOR c5=c3 TO c4
          r1=(c(c5)-b)/m
          PSET ( sx*(r2-xLo)+xb,sy*(r1-yLo)+yb),(c5 MOD 3)+1 'corresp TO move+draw
        NEXT c5
con1:
      NEXT k
      REM      J-direction interpolate
      FOR k=0 TO 1-s3 STEP s3
        z1=x1-k*(x1-x2)
        z2=x3-k*(x3-x4)
        IF z1<c(1) AND z2<c(1) GOTO con2
        GOSUB cross
        IF c3>c4 GOTO con2
        m=z2-z1
        b=z1-m*i
        r1=j+k
        FOR c5=c3 TO c4
          r2=(c(c5)-b)/m
          PSET ( sx*(r2-xLo)+xb,sy*(r1-yLo)+yb),(c5 MOD 3)+1
        NEXT c5
con2: NEXT k
confin:
    NEXT j
confin1:
  NEXT i
  PRINT "done"
  GOTO fin
 

cross:  'check FOR contour crossing between Z1 AND Z2
  IF z1>z2 GOTO cross1
  Y1=z1
  Y2=z2
  GOTO cross2

cross1:
  Y1=z2
  Y2=z1

cross2:
  FOR c3=1 TO C1
    IF Y1<c(c3) GOTO cross3
  NEXT c3
  c4=0
RETURN

cross3:
  FOR c4=c3 TO C1
    IF Y2<=c(c4) GOTO cross4
  NEXT c4

cross4:
  c4=c4-1
RETURN

fin:
END SUB

END