REM =============================================================== REM Documentation REM ============= REM CONREC performs contouring of a rectangular grid of data points. REM Any number of contour levels are possible but they must REM be provided in increasing order.. REM The routine calls a user defined subroutine VECOUT with REM the start and stop coordinates of each line segment and REM the contour level for that line segment. REM History REM Modified from CONREC routine written in FORTRAN-77 REM by Dr M.D.Johns for a PDP-11 at Auckland University Physics REM Department. (25-Nov-83) REM Input variables to CONREC REM d(0:iub,0:jub) ' matrix for the data surface REM iub, jub ' index bounds of the data array REM x(0:iub) ' data array for column coordinates REM y(0:jub) ' data array for row coordinates REM nc ' number of contour levels REM z(0:nc-1) ' contour levels in increasing order REM VECOUT ' An external subroutine to plot the contour lines REM False and true boolean values conrec: REM Local declarations for CONREC REM =========================== DIM h(4) ' relative heights of the box above contour DIM ish(4) ' sign of h() DIM xh(4) ' x coordinates of box DIM yh(4) ' y coordinates of box DIM im(3) ' mapping from vertex numbers to x offsets im(0)=0 : im(1)=1 : im(2)=1 : im(3)=0 DIM jm(3) ' mapping from vertex numbers to y offsets jm(0)=0 : jm(1)=0 : jm(2)=1 : jm(3)=1 DIM castab(2,2,2) ' case switch table DATA 0,0,8,0,2,5,7,6,9,0,3,4,1,3,1,4,3,0,9 ,6,7,5,2,0,8,0,0 FOR k=0 TO 2 : FOR j=0 TO 2 : FOR i=0 TO 2 READ castab(k,j,i) NEXT i : NEXT j : NEXT k REM Check the input parameters for validity REM =================================== prmerr=false IF (iub<=0 OR jub<=0) THEN prmerr=true IF (nc<=0) THEN prmerr=true FOR k=1 TO nc-1 : IF (z(k)<=z(k-1)) THEN prmerr=true : NEXT k IF (prmerr) THEN msg\$="Error in input parameters" : RETURN REM Scan the array, top down, left to right REM ================================== FOR j=jub-1 TO 0 STEP -1 FOR i=0 TO iub-1 REM Find the lowest vertex IF (d(i,j)d(i,j+1)) THEN dmax=d(i,j) ELSE dmax=d(i,j+1) IF (d(i+1,j)>dmax) THEN dmax=d(i+1,j) IF (d(i+1,j+1)>dmax) THEN dmax=d(i+1,j+1) IF (dmaxz(nc-1)) THEN GOTO noneinbox REM Draw each contour within this box FOR k=0 TO nc-1 IF ((z(k)dmax)) THEN GOTO noneintri FOR m=4 TO 0 STEP -1 IF m>0 THEN h(m)=d(i+im(m-1),j+jm(m-1))-z(k) : xh(m)=x(i+im(m-1)) : yh(m)=y(j+jm(m-1)) IF m=0 THEN h(0)=(h(1)+h(2)+h(3)+h(4))/4 : xh(0)=(x(i)+x(i+1))/2 : yh(0)=(y(j)+y(j+1))/2 IF h(m)>0 THEN ish(m)=2 : ELSE IF (h(m)<0) THEN ish(m)=0 : ELSE ish(m)=1 NEXT m REM Scan each triangle in the box FOR m=1 TO 4 m1=m : m2=0 : m3=m+1 : IF (m3=5) THEN m3=1 case=CINT(castab(ish(m1),ish(m2),ish(m3))) IF (case=0) THEN GOTO case0 ON case GOTO case1,case2,case3,case4,case5,case6,case7,case8,case9 REM Line between vertices m1 and m2 case1: x1=xh(m1) : y1=yh(m1) : x2=xh(m2) : y2=yh(m2) GOTO drawit REM Line between vertices m2 and m3 case2: x1=xh(m2) : y1=yh(m2) : x2=xh(m3) :y2=yh(m3) GOTO drawit REM Line between vertices m3 and m1 case3: x1=xh(m3) : y1=yh(m3) : x2=xh(m1) : y2=yh(m1) GOTO drawit REM Line between vertex m1 and side m2-m3 case4: x1=xh(m1) : y1=yh(m1) x2=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) y2=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) GOTO drawit REM Line between vertex m2 and side m3-m1 case5: x1=xh(m2) : y1=yh(m2) x2=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) y2=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) GOTO drawit REM Line between vertex m3 and side m1-m2 case6: x1=xh(m3) : y1=yh(m3) x2=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) y2=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) GOTO drawit REM Line between sides m1-m2 and m2-m3 case7: x1=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) y1=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) x2=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) y2=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) GOTO drawit REM Line between sides m2-m3 and m3-m1 case8: x1=(h(m3)*xh(m2)-h(m2)*xh(m3))/(h(m3)-h(m2)) y1=(h(m3)*yh(m2)-h(m2)*yh(m3))/(h(m3)-h(m2)) x2=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) y2=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) GOTO drawit REM Line between sides m3-m1 and m1-m2 case9: x1=(h(m1)*xh(m3)-h(m3)*xh(m1))/(h(m1)-h(m3)) y1=(h(m1)*yh(m3)-h(m3)*yh(m1))/(h(m1)-h(m3)) x2=(h(m2)*xh(m1)-h(m1)*xh(m2))/(h(m2)-h(m1)) y2=(h(m2)*yh(m1)-h(m1)*yh(m2))/(h(m2)-h(m1)) drawit: CALL vecout(x1,y1,x2,y2,z(k)) case0: NEXT m noneintri: NEXT k noneinbox: NEXT i : NEXT j RETURN REM ===============================================================