V10/cmd/view2d/co.e

#  simple contour plotting
#    connect-the-dots, rather than approximate hyperbolas.
#    contour (implicitly) perturbed downward by infinitesimal.
#  extensively revised 12 Apr 1985
#  added NaN  9 Jul 1985

procedure g2co(w,d,n1,n2,f,nc,cl,NaN)
integer n1, n2, j0, j1, j2, nc
real w(6), d(2,2), f(n1,n2), cl(nc)
real c, dx(2), ll(2), NaN
integer j, jj, xp, nt
real fl, fh, t(2,4)
real swap
real xx(2,60)
dx(1)=(d(1,2)-d(1,1))/(n1-1)
dx(2)=(d(2,2)-d(2,1))/(n2-1)
xp=0
do j0 = 1,nc{
   c=cl(j0)
   do j1 = 1,n1-1{
      ll(1)=d(1,1)+(j1-1)*dx(1)
      do j2 = 1,n2-1{
         ll(2)=d(2,1)+(j2-1)*dx(2)
         #write(,j1:i(4),j2:i(4))
         #write(,f(j1,j2+1):f(7,1),f(j1+1,j2+1):f(7,1))
         #write(,f(j1,j2  ):f(7,1),f(j1+1,j2  ):f(7,1))

         #   find intersections of contours with sides
         nt=0
         #   ... west ...
         fl=f(j1,j2)
         fh=f(j1,j2+1)
         if(fl<=NaN) next
         if(fh<=NaN) next
         if( (c-fl)*(c-fh)<0 |
               c==fl & fh<c  |
               c==fh & fl<c  ){
            nt=nt+1
            t(1,nt)=0
            t(2,nt)=(c-fl)/(fh-fl)
         }
         #   ... north ...
         fl=fh
         fh=f(j1+1,j2+1)
         if(fh<=NaN) next
         if( (c-fl)*(c-fh)<0 |
               c==fl & fh<c  |
               c==fh & fl<c  ){
            nt=nt+1
            t(1,nt)=(c-fl)/(fh-fl)
            t(2,nt)=1
         }
         #   ... south ...
         fl=f(j1,j2)
         fh=f(j1+1,j2)
         if(fh<=NaN) next
         if( (c-fl)*(c-fh)<0 |
               c==fl & fh<c  |
               c==fh & fl<c  ){
            nt=nt+1
            t(1,nt)=(c-fl)/(fh-fl)
            t(2,nt)=0
         }
         #   ... east ...
         fl=fh
         fh=f(j1+1,j2+1)
         if( (c-fl)*(c-fh)<0 |
               c==fl & fh<c  |
               c==fh & fl<c  ){
            nt=nt+1
            t(1,nt)=1
            t(2,nt)=(c-fl)/(fh-fl)
         }

         #  swap north and south to correctly match with east and west
         if(nt==4){
            if( t(1,2)>t(1,3) |
                (t(1,2)==t(1,3) & f(j1,j2)<f(j1+1,j2))  ){
               swap=t(1,2)
               t(1,2)=t(1,3)
               t(1,3)=swap
               swap=t(2,2)
               t(2,2)=t(2,3)
               t(2,3)=swap
            }
         }

         if ( nt==2 | nt==4 ){
           do jj=1,nt/2{
             j=2*jj-1
             if( t(1,j+1)!=t(1,j) | t(2,j+1)!=t(2,j) ){
               xp=xp+2
               xx(1,xp-1)=ll(1)+dx(1)*(t(1,j))
               xx(2,xp-1)=ll(2)+dx(2)*(t(2,j))
               xx(1,xp)=ll(1)+dx(1)*(t(1,j+1))
               xx(2,xp)=ll(2)+dx(2)*(t(2,j+1))
             }
           }
         }else if( nt!=0 ){
           write(,"can't happen!  nt=",nt,j1,j2)
           write(,f(j1,j2+1),f(j1+1,j2+1))
           write(,f(j1,j2  ),f(j1+1,j2  ))
           stop
         }

         if(xp>=50){
            g2la(w,xx,xp)
            xp=0
            }

         }  # do j2
      } # do j1

   if(xp>0){
      g2la(w,xx,xp)
      xp=0
      }
   } # do j0
end