subroutine anacon(xarr,nx,ny,xlev,numlev,itvord, & xl,yl,xsc,ysc,cut) c c xarr is data array to contour (real*4) c nx,ny is size of xarr c xlev is a table of contour levels real*4 c numlev are the number of levels to contour c itvord is the desired flip of the array - see inccon below for details c xsc,ysc are the coords. of the upper rhs c xl,yl are the coordinates of the lower left corner for the contour c real*4 xarr(nx*ny),aa(4),bb(4) real*4 xlev(numlev),xa(4),ya(4) integer*4 ltvcol(8) real*4 cell(5),xmax,xmin c type *,'numlev, nx,ny =',numlev, nx,ny c do i=1,numlev c type *,i,xlev(i) c enddo c-- compute internal values for spacing c type *,'xl, yl, xsc, ysc', xl, yl, xsc, ysc xsc=xsc-xl ysc=ysc-yl if(mod(itvord,2).eq.0)then !match tektronix x,y with data x,y ysc=ysc/ny xsc=xsc/nx else ysc=ysc/nx xsc=xsc/ny endif c must reset ix,iy to center the contour on the raster xz=xl+xsc/2. yz=yl+ysc/2. c type *,'xz, yz, xsc, ysc', xz, yz, xsc, ysc c c set up increments and limits for array processing according to itvord c call inccon(itvord,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4) minc=ny-1 if(ninc.eq.minc)minc=nx-1 do irow=1,minc !row counter for position on grinnell do icol=1,ninc !col counter for position of grinnell cell(1)=xarr(ibase) cell(2)=xarr(ibase+n2) cell(3)=xarr(ibase+n3) cell(4)=xarr(ibase+n4) cell(5)=cell(1) c set ibase for next cell ibase=ibase+ibindx xmax=cell(1) xmin=cell(1) imax=1 do j=2,4 !find max,min of cell xmax=amax1(xmax,cell(j)) xmin=amin1(xmin,cell(j)) if(cell(j).eq.xmax)imax=j !keep track of max enddo c c now check each cell for contours c do k=1,numlev if(xlev(k).le.xmax.and.xlev(k).ge.xmin)then c level goes through cell ndsh=1 if(xlev(k).lt.cut)ndsh=2 ns=0 !number of sides do l=1,4 if(xlev(k).ge.cell(l).and.xlev(k).lt. & cell(l+1).or.xlev(k).lt.cell(l).and. & xlev(k).ge.cell(l+1))then c c calculate endpoints of line segment on the side of the cell c ns=ns+1 if(l.eq.1)then !top xa(ns)=float(icol-1)+(xlev(k)-cell(l)) & /(cell(l+1)-cell(l)) ya(ns)=float(irow-1) else if(l.eq.2)then !right side xa(ns)=float(icol) ya(ns)=float(irow-1)+(xlev(k)-cell(l)) & /(cell(l+1)-cell(l)) else if (l.eq.3)then !bottom side xa(ns)=float(icol-1)+(xlev(k)- & cell(l+1))/(cell(l)-cell(l+1)) ya(ns)=float(irow) else !left side xa(ns)=float(icol-1) ya(ns)=float(irow-1)+(xlev(k)- & cell(l+1))/(cell(l)-cell(l+1)) endif endif enddo c if(ns.gt.2)then ntim=2 xmean=cellavg(cell,4) c if the present level is greater than the mean, then the contour c is drawn between the max point and the center. (imax gives max c point) c first calculate all screen coordinates in=imax do ii=1,4 aa(ii)=xz+xa(in)*xsc bb(ii)=yz+ya(in)*ysc in=in-1 if(in.lt.1)in=4 enddo if(xlev(k).le.xmean)then c switch 2 and 4 so that point 1 connects to point 2 and point 3 to 4. tmp=aa(4) aa(4)=aa(2) aa(2)=tmp tmp=bb(4) bb(4)=bb(2) bb(2)=tmp endif else if(ns.eq.2)then ntim=1 do ii=1,2 aa(ii)=xz+xa(ii)*xsc bb(ii)=yz+ya(ii)*ysc enddo endif !ns.gt.2 c now draw the contours if (ns.ge.2) then call tkdash(aa,bb,ndsh,ntim) endif endif !xlev(k).gt.... enddo !k=1,numlev enddo !icol=1,nx-1 ibase=ibase+incr !reset for next row or col enddo !irow=1,ny-1 999 return end c function cellavg(cell,nav) real*4 cell(*) cellavg=0.0 do ii=1,nav cellavg=cellavg+cell(ii) enddo cellavg=cellavg/float(nav) 999 return end c subroutine inccon(itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4) c get limits and increments for processing the raw array c inputs: c itv is the order to process : c 0 left to right, bottom to top c 1 bottom to top, left to right c 2 left to right, top to bottom c 3 top to bottom, left to right c 4 right to left, top to bottom c 5 top to bottom, right to left c 6 right to left, bottom to top c 7 bottom to top, right to left c nx,ny size of raw array (columns,rows) c c outputs: c ibase base index to begin processing c ibindx index used to increment base address during processing c incr used at end of a row or column to reset base address c n2,n3,n4 offsets in base address for 4 point cell c c note that even though itv may indicate array is to be processed c vertically in a certain direction (either top to bottom or bottom to c top), the outputs are such that the array is processed in the c opposite direction. this is because the contours are actually drawn c from the bottom of the screen to the top. c c-- 7/7/86 orders 1 and 5 seem to be reversed, patch to fix c the comments and logic in parts of this seem wrong, should be cleaned c up someday, beware ! if(mod(itv,2).eq.0)then !process horizontal before vertical ibindx=1 if(itv.gt.3)ibindx=-1 !right to left incr=1 if(itv.eq.2.or.itv.eq.6)incr=2*nx-1 if(itv.eq.2.or.itv.eq.4)incr=-incr ninc=nx-1 !process to end of row before reset else !process vertical before horizontal ibindx=nx !processing top to bottom if(itv.eq.5.or.itv.eq.7)ibindx=-nx !nope, bottom to top incr=(ny-1)*nx+1 !reset a col to begin processing next col if(itv.eq.1.or.itv.eq.5)incr=incr-2 if(itv.eq.3.or.itv.eq.1)incr=-incr ninc=ny-1 !process to end of column before reset endif c now get initial base address if(itv.eq.2.or.itv.eq.7)then ibase=(ny-1)*nx+1 !start at lower left else if(itv.eq.4.or.itv.eq.5)then ibase=nx*ny !lower right else if(itv.eq.0.or.itv.eq.1)then ibase=1 !upper left else ibase=nx !upper right endif c now need to fill the cell properly - figure offsets from base address if(itv.eq.0)then n2=1 n3=nx+1 n4=nx else if(itv.eq.3)then n2=nx n3=nx-1 n4=-1 else if(itv.eq.2)then n2=1 n3=-nx+1 n4=-nx else if(itv.eq.5)then n2=-nx n3=-nx-1 n4=-1 else if(itv.eq.4)then n2=-1 n3=-nx-1 n4=-nx else if(itv.eq.7)then n2=-nx n3=-nx+1 n4=1 else if(itv.eq.6)then n2=-1 n3=nx-1 n4=nx else if(itv.eq.1)then n2=nx n3=nx+1 n4=1 endif c type *,'variables in inccon:' c type 1001,itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4 c1001 format(1x,'itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4',/,1x,10i6) return end