allbery@ncoast.UUCP (12/04/87)
NetLanders;
The contour routine is not big. It is actually quite small.
The key to using the routines is as follows:
The routines expect the data to be upon some plane of coordinates x
and y or whatever is being used. The random points are located within this
2d mesh of x and y points. The routine requires that 3 2d arrays should be
present the first array is the location of the x point in mesh coordinates,
and so on. You can read about it in the Man page supplied. This could be
used to do random point contouring by masking the unused coordinates within
the arrays with some value that you know the random values will not touch.
Either way above the max or way below the minimun value of the function to
be plotted. The source and man page follow. It uses unix type plot calls
from a package developed here at Purdue (CRC Graphics). I am sure your
system has a similiar package....Ed
---------------8<------------------------8<--------------------8<-------------
# This is a shell archive. Save this into a file, edit it
# and delete all lines above this comment. Then give this
# file to sh by executing the command "sh file". The files
# will be extracted into the current directory owned by
# you with default permissions.
#
# The files contained herein are:
# UContour.f UContour.3G
#
echo 'x - UContour.f'
sed 's/^X//' <<'________This_Is_The_END________' >>UContour.f
Xc
Xc UContour was originally written by Dr. W.J. Usab for the
Xc PEPL graphics package. The routines have been modified to
Xc use the CRC graphics package by Edward L. Haletky and Tom
Xc Jentink.
Xc
Xc The following modifications were made: 11/13/87
Xc Use of output files and not direct printing (gf)
Xc gf options: filename - plot to file filename
Xc scr - plot to standard output
Xc lpr - plot to Printronix
Xc
Xc added ml option - specifies length of column in 2D matrix
Xc Call to plots,fname were added to routines
Xc Added Label option
Xc Call to plot(0,0,999) added
Xc Number call changed to reflect C syntax on formats
Xc Transformed to lower case
Xc
Xc To Compile:
Xc f77 -c UContour.f -lG
Xc
Xc
Xc subroutine contour(x,y,f,n,m,cmin,cmax,inc)
X subroutine ucontour(x,y,f,ml,in,jmax,cmin,cmax,inc,gf,
X +lbl,bd)
X real xmin,xmax,ymin,ymax,gscale
X common /climit/ xmin,xmax,ymin,ymax,gscale
X character *8 option
X character *20 lbl
X character *3 dev1,dev2,gf
X integer dev,blank
X real x(ml,*),y(ml,*),f(ml,*),bd(4),jmax(ml)
Xc define constants
X option = ' '
X dev1 = "scr"
X dev2 = "lpr"
X dev = 0
X if (gf.eq.dev1) dev = 6
X if (gf.eq.dev2) dev = 16
X blank = 0
X xfac = 10.0
X yfac = 10.0
X
Xc define min and max ranges for x,y,f
X
X call plots(dev,blank,option)
X if ((dev.ne.16).and.(dev.ne.6)) call fname(gf)
X call symbol(2,9.5,0.075,lbl,0.0)
Xc
Xc define min and max ranges for x,y,f
Xc
X xmin = x(1,1)
X xmax = x(1,1)
X ymin = y(1,1)
X ymax = y(1,1)
X fmin = f(1,1)
X fmax = f(1,1)
Xc
Xc do 10 j = 1,m
Xc do 10 i = 1,n
X do 10 i=1,in
X m=jmax(i)
X do 10 j=1,m
X xmin = min(x(i,j),xmin)
X xmax = max(x(i,j),xmax)
X ymin = min(y(i,j),ymin)
X ymax = max(y(i,j),ymax)
X fmin = min(f(i,j),fmin)
X fmax = max(f(i,j),fmax)
X 10 continue
Xc
X xscale = 1.0*xfac/(xmax-xmin)
X yscale = 1.0*yfac/(ymax-ymin)
X gscale = min(xscale,yscale)
Xc
Xc do 20 j = 1,m-1
Xc do 20 i = 1,n-1
X do 20 i=1,in-1
X m=jmax(i)
X do 20 j=1,m-1
X x1 = gscale*(x(i,j)-xmin)
X y1 = gscale*(y(i,j)-ymin)
X f1 = f(i,j)
X x2 = gscale*(x(i,j+1)-xmin)
X y2 = gscale*(y(i,j+1)-ymin)
X f2 = f(i,j+1)
Xc******************************
X if (j.eq.(m-1)) then
X x3=gscale*x(i+1,j)
X y3=gscale*y(i+1,j)
X f3=f(i+1,j)
X else
X x3 = gscale*(x(i+1,j+1)-xmin)
X y3 = gscale*(y(i+1,j+1)-ymin)
X f3 = f(i+1,j+1)
X endif
Xc******************************
X x4 = gscale*(x(i+1,j)-xmin)
X y4 = gscale*(y(i+1,j)-ymin)
X f4 = f(i+1,j)
X diag1 = (x3-x1)**2+(y3-y1)**2
X diag2 = (x4-x2)**2+(y4-y2)**2
X if(diag1.le.diag2) then
X call tricon(x1,y1,f1,x2,y2,f2,x3,y3,f3,
X 1 cmin,cmax,inc)
X call tricon(x1,y1,f1,x3,y3,f3,x4,y4,f4,
X 1 cmin,cmax,inc)
X else
X call tricon(x1,y1,f1,x2,y2,f2,x4,y4,f4,
X 1 cmin,cmax,inc)
X call tricon(x2,y2,f2,x3,y3,f3,x4,y4,f4,
X 1 cmin,cmax,inc)
X endif
Xc
X 20 continue
Xc
Xc plot contour region boundary (clockwise from 0,0)
Xc
X x1 = gscale*(x(1,1)-xmin)
X y1 = gscale*(y(1,1)-ymin)
X call plot(x1,y1,3)
Xc
X m=jmax(1)
X x1=gscale*x(1,m)
X y1=gscale*y(1,m)
X call plot(x1,y1,2)
Xc
X
X x1=gscale*bd(1)
X y1=gscale*bd(2)
X call plot(x1,y1,2)
Xc
X x1=gscale*bd(3)
X y1=gscale*bd(4)
X call plot(x1,y1,2)
Xc
X if (bd(3).lt.x(in,1)) then
X x1=gscale*x(in,1)
X y1=gscale*bd(4)
X call plot(x1,y1,2)
Xc
X x1 = gscale*x(in,1)
X y1 = gscale*y(in,1)
X call plot(x1,y1,2)
Xc
X x1 = gscale*(x(1,1)-xmin)
X y1 = gscale*(y(1,1)-ymin)
X call plot(x1,y1,2)
X endif
X call plot(0,0,999)
Xc
Xc
X return
X end
Xc
X subroutine tricon(xa,ya,fa,xb,yb,fb,xc,yc,fc,cmin,cmax,inc)
X real x(3),y(3),f(3)
Xc
X x(1) = xa
X x(2) = xb
X x(3) = xc
X y(1) = ya
X y(2) = yb
X y(3) = yc
X f(1) = fa
X f(2) = fb
X f(3) = fc
Xc
Xc sort points in ascending order
Xc
X do 10 n=1,2
X do 10 i=1,3-n
X if(f(i).gt.f(i+1))then
X xt = x(i)
X yt = y(i)
X ft = f(i)
X x(i) = x(i+1)
X y(i) = y(i+1)
X f(i) = f(i+1)
X x(i+1) = xt
X y(i+1) = yt
X f(i+1) = ft
X endif
X 10 continue
Xc
Xc locate contour range within triangle
Xc
X dc = (cmax-cmin)/float(inc)
X icmin = int((f(1)-cmin)/dc)
X icmin = max(icmin,0)
X icmax = int((f(3)-cmin)/dc)
X icmax = min(icmax,inc)
X if(icmax.lt.icmin) return
Xc
Xc plot contours within range
Xc
X do 20 i = icmin,icmax
Xc do 20 i=0,inc
X clev = cmin+dc*float(i)
Xc if(clev.lt.f(1)) go to 20
X if (clev.le.f(1)) goto 20
Xc if(clev.gt.f(3)) go to 20
X if (clev.ge.f(3)) goto 20
X if(f(2).gt.clev) then
X fac = (clev-f(1))/(f(2)-f(1))
X x1 = x(1)+fac*(x(2)-x(1))
X y1 = y(1)+fac*(y(2)-y(1))
X else
X fac = (clev-f(2))/(f(3)-f(2))
X x1 = x(2)+fac*(x(3)-x(2))
X y1 = y(2)+fac*(y(3)-y(2))
X endif
X fac = (clev-f(1))/(f(3)-f(1))
X x2 = x(1)+fac*(x(3)-x(1))
X y2 = y(1)+fac*(y(3)-y(1))
Xc
X call plot(x1,y1,3)
X call plot(x2,y2,2)
Xc call labinc(i,clev,x1,y1,x2,y2)
X 20 continue
Xc
X return
X end
Xc
X subroutine labinc(i,clev,x1,y1,x2,y2)
Xc common /clabinc/ icount
Xc data icount/0/
X tdist = abs(x2-x1)+abs(y2-y1)
X if(tdist.lt.1e-20) then
X write(6,*)' ** error in labinc: tdist=',tdist
X return
X endif
X interv = 21
X icount = icount+1
X if(icount.lt.interv) return
X icount = 0
X pi = acos(-1.0)
X hpi = 0.5*pi
X theta = atan2((y2-y1),(x2-x1))
X if(theta.ge.hpi) then
X theta = theta-pi
X else if(theta.lt.-hpi) then
X theta = theta+pi
X endif
X thetan = theta+hpi
X x = 0.5*(x1+x2)+0.05*cos(thetan)-0.225*cos(theta)
X y = 0.5*(y1+y2)+0.05*sin(thetan)-0.225*sin(theta)
X thetad = 180.0*theta/pi
Xcc call number(x,y,.1,i,0.0,'i2')
Xc call number(x,y,.075,clev,thetad,'f6.3')
Xc call number(x,y,.055,thetad,"%6.3f",clev)
X return
X end
________This_Is_The_END________
echo 'x - UContour.3G'
sed 's/^X//' <<'________This_Is_The_END________' >>UContour.3G
X.\" troff -man
X.EN
X.TH UContour 3G 11/13/87
X.SH NAME
XUContour - EGraphics Contour Routine.
X.SH AUTHORS
XDr. W.J.Usab Jr.
X.br
XEdward L. Haletky
X.br
XTom Jentink
X.SH SYNOPSIS
Xcall ucontour(x,y,f,ml,in,jmax,cmin,cmax,inc,gf,lbl,bd)
X.SH DESCRIPTION
X.fi
XUContour will plot a family of curves of number
X.I inc
Xwithin the range
X.I [cmin:cmax].
XAll the curves will be of a constant value of the function,
X.I f.
X.sp
X.RS
X.fi
X.B x
X- two dimensional array for x-coordinate at grid point (n,m) at which
Xthe function,
X.I f,
Xwas computed.
X.sp 1
X.B y
X- two dimensional array for y-coordinate at grid point (n,m) at which
Xthe function,
X.I f,
Xwas computed.
X.sp 1
X.B f
X- two dimensional array for the function at grid point (n,m).
X.sp 1
X.B ml
X- maximun length of the arrays
X.I x,
X.I y,
X.I f.
XThe arrays are declared in the following manor:
X.br
X.ce 1
Xreal x(ml,*),y(ml,*),f(ml,*)
X.br
X.sp 1
X.B in
X- maximun number of grid points in n-direction.
X.sp 1
X.B jmax
X- an array that gives maximun number of grid points in the m-direction,
Xfor each,
X.I n
Xstation.
X.sp 1
X.B cmin
X- minimun value of the constant-value contour of the function,
X.I f.
X.sp 1
X.B cmax
X- maximun value of the constant-value contour of the function,
X.I f.
X.sp 1
X.B inc
X- integer number of contours between the range
X.I [cmin:cmax].
X.sp 1
X.B gf
X- device which to plot
X.RE
X.in +10
X.B scr
X- plot to stdout.
X.br
X.B lpr
X- plot to Printronix Printer.
X.br
X.B gf
X- plot to file. File can be printed with the
X.I gplp(1G)
Xcommand. The
X.I filename
Xcan not be more than 3 characters long.
X.in -10
X.RS
X.sp 1
X.B lbl
X- label the contour plot.
X.B lbl
Xcan be a maximun of 20 characters long.
X.br
X.sp 1
X.B bd
X- an array that holds the bounding area values.
X.RE
X.in 10
X.B bd(1)
X- xb1
X.br
X.B bd(2)
X- yb1
X.br
X.B bd(3)
X- xb2
X.br
X.B bd(4)
X- yb2
X.br
X.in -10
X.RS
X.sp 1
X.B "To Compile"
Xf77 prog.f -L/x/edward/lib -lE -lG
X.RE
X.SH FILES
X.RS
X/x/edward/lib/libE.a
X.br
X/x/edward/src/UContour.f
X.RE
X.SH DIAGNOSTICS
X.LP
XUContour will blank the page before plotting. No CRC graphics options are used.
XThe Border of the Channel will be printed, it is a crude representation. The
Xlabeling of the lines will be implemented at a later date. There are a few bugs
Xstill.
X.SH "SEE ALSO"
X.B "CRC Graphics Manaul ,"
X.B gplp(1G)
X.SH BUGS
X.RS
XPlease Report all bugs to:
X.br
Xedward@ga.ecn.purdue.edu
X.RE
________This_Is_The_END________
exit
---------------8<------------------------8<--------------------8<-------------
===================================== ^ =================================
Edward L. Haletky |E |o| U| "To race against the sky..."
Usenet: ~!pur-ee!edward |L ^/| |\^ S| "The world of the
Arpa: edward@ga.ecn.purdue.edu |H / | | \ A| immagination is boundless."
=====================================/__-|-__\=================================