#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'README' <<'END_OF_FILE' XDavid Scott's ASH routines X X Xash1 1D ash (density plot) Xbin1 forms bin counts X Xbin2 2D bin counts used by Xash2 2D ash X Xashn multi-D ash and 3D visualization program END_OF_FILE if test 170 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test -f 'Instructions' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Instructions'\" else echo shar: Extracting \"'Instructions'\" \(5571 characters\) sed "s/^X//" >'Instructions' <<'END_OF_FILE' XThis suite contains three separate sets of programs for 1D, 2D and 3D Xdensity estimation. The last, ashn, is combined with a visualization program. X X1. Put all the files in a safe place, say SHOME/library/ash X X2. Edit Makefile as needed, especially SHOME and S. Experiment may show that X other libraries are needed in the ld command for loadmod.o: that supplied X works on a Sun 4. X X3. make X or make ash for 1D & 2D only, X make ashn for 3D only X X X1D and 2D X========= X Xdyn.load("ash.o") Xsuntools() # or openlook(), motif(), .... Xx <- rnorm(100) # data Xf <- ash1(bin1(x,nbin=50),5) # compute ash estimate Xplot( f , type="l" ) # line plot of estimate Xab <- c(-5,5) # bin interval Xbins <- bin1(x,ab,10) # bin x into 10 bins over ab X X Xx <- matrix( rnorm(200), 100 , 2) # bivariate normal n=100 Xab <- matrix( c(-5,-5,5,5), 2, 2) # interval [-5,5) x [-5,5) Xnbin <- c( 20, 20) # 400 bins Xbins <- bin2(x, ab, nbin) # bin counts,ab,nskip Xm <- c(5,5) Xf <- ash2(bins,m) Xpersp( f$z / max(f$z) ) Xcontour( f ) X# note: it is a good idea to not plot contours at level = 0 X X X3D X== X X4. To run graphics on a colour machine under OpenWindows or Motif, within X S-plus invoke openlook(OL.col) or motif(OL.col). Note: this will not X work if ~/.sgraphrc sets the colours, as it will if you have saved X them. In that case set up a new colour scheme with these X details in OL.col (but you don't need \"). X X5. To use the suntools() driver add the following line to ~/.defaults: X X/Splus/ColorScheme11 "Scott: (1,1,1) (255,0,0) 48 (255,255,255) (0,0,255) 48 (0,255,255)" X X and select this colour scheme for the window. X X6. If this is a library and your S/S-Plus is late enough, .First.lib will X automatically load in "loadmod.o". Otherwise your code will have to do it. X7. Actually, I have an automatic mode that you X can use to see if everything is working. X Exit Splus if you are in it. X % chmod +x norm.1 and edit it as appropriate X % norm.1 X a. this fires Splus up (the local version) X b. starts the graphics window X c. generates 5 views of the 999 trivariate X normal points in "xyz" using the X "action option" choices in "stream.s" X X8. The automatic version just zips along, primarily to X produce files with the *.quad extension. These X can be feed into the MinneView program. X X9. Otherwise, here is a typical session that you might try. X X > openlook(OL.col) # a nice color window X > dyn.load("loadmod.o") # if necessary X > pie( rep(1,100), col=1:100 ) # look at the nice colors X > ashn(xyz) # it does some automatic things X # you can ignore for now X 6 X 2 5 1 # this changes the "eye" location X 2 X 4 4 4 # this changes the smoothing parameters X 3 X .5 # select the shell at alpha-level X # 0.50 * f(mode) 0 < alpha < 1 X 0 # plot it X 3 X .2 # plot the alpha = 0.20 shell X 13 X 0 # don't erase window for next plot X 3 X .5 X 0 # adds to the plot [this will work X # great on the SGI if the X # shells have been made transparent X 13 X 1 # back to erasing screeen X 3 X .1 X 8 X 2 1 1 # leave some of the shell off X 11 X 1 # turn on *.quad output X 14 X 1 0 0 .5 # RBG color + transparency fraction X # for shell [all between 0 and 1] X 0 # plot it X 99 # quit X > q() X X X Action options: X 0 = plot with current parameter settings [recomputes ash if necessary] X 1 = t slice; value of unseed variable f(x,y,z,t); X an integer between 1 and nbin[4] -- bin center; X if 5-D f(x,y,z,t1,t2), then specify t1 and t2 as integers X 2 = m vector of smoothing parameters (length nvar) X 3 = f alpha-shell level; fraction between 0 and 1 X 4 = axes show/don't show axes on plot X 5 = cube show/don't show box on plot X 6 = eye vector of length 3 giving eye coordinates X [ the program views the world as living in X the cube [-1,1]^3 ] X 7 = slicelimit limits region where shells are computed to X hyperrectangle defined by these 6 integers X 8 = sliceincrement normally 1; >1 skips bins X 9 = iris can specific coloroffset=999 on an SGI Iris X 10 = orth only plot triangles oriented towards viewer X 11 = print 1 = produce *.quad files X 12 = coltbl lookup from old S days X 13 = blank screen yes/no erase screen before plot X 14 = alpha RGBA or {red,green,blue,alpha} all between 0 and 1 X this alpha refers to the transparency of the shell X 15 = fmaxin if more than 3 variables, allow the use to specify X the value of f(x,y,z,t) at the global mode X 16 = kernel 0=ash 1,2,3,... use kernels of the form (1-x^2)^kernel X 99 = exit ashn X X Note that ashn attempts to plot an ash shell based on the calling X sequence in the S function ashn. By default, f=4 (which is obviously X greater than 1), so that no shell will be found unless "ff" is X specified (note that "m" is called "mv" in "ashn"): X X ashn(xyz) is equivalent to ashn(xyz, ff=4, mv=c(1,1,1)) X which does nothing very interesting until you X using the action options X ashn(xyz, ff=.5, mv=c(4,4,4)) would at least do something X initially, and would not have to recompute the X ash if you changed the "eye", for example X X X10. This version of ashn actually will look at shells of X slices of the first 3 variables of a multivariate X ASH of up to 10 variables. I really haven't tried X it beyond 5, but the loops are there. Also there X is a hidden regression surface option I can document X if you are interested. X X11. Thanks for your interest and patience. Let me know what X else would be helpful. X X------------------------------------ X END_OF_FILE if test 5571 -ne `wc -c <'Instructions'`; then echo shar: \"'Instructions'\" unpacked with wrong size! fi # end of 'Instructions' fi if test -f 'Makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Makefile'\" else echo shar: Extracting \"'Makefile'\" \(849 characters\) sed "s/^X//" >'Makefile' <<'END_OF_FILE' XS = Splus XSHOME = /home/splus3.1 XCFLAGS = -O -I$(SHOME)/include XFFLAGS = -O X XOBJS = tmp2.o sort.o dumpit.o writeit.o cubeint.o simpint.o myinput.o X XHELP = ash1 ash2 bin1 bin2 X XHLP = $(HELP:%=.Data/.Help/%) X X.SUFFIXES: X.SUFFIXES: .o .c .f .r X Xall: ash ashn X Xash: .Data .Data/.Help ash.o .Data/ash1 $(HLP) Xashn: .Data .Data/.Help loadmod.o .Data/compress .Data/ct .Data/ashn .Data/.Help/ashn X X.Data: X mkdir .Data X X.Data/.Help: X mkdir .Data/.Help X X.Data/.Help/%: %.d X cp $*.d .Data/.Help/$* X Xloadmod.o: $(OBJS) X ld -r -o loadmod.o $(OBJS) /lib/libm.a X X.Data/ashn: ashn.s X $S < ashn.s X.Data/ash1: ash.s X $S < ash.s X.Data/compress: compress.s X $S < compress.s X.Data/ct: ct mymaster init.s X $S < init.s Xclean: X rm -rf .Data *.o X Xshar: X shar -o scott.shar README Instructions Makefile *.s *.d *c tmp2.r sort.f \ X ash.f ct mymaster stream.s.norm.1 norm.1 END_OF_FILE if test 849 -ne `wc -c <'Makefile'`; then echo shar: \"'Makefile'\" unpacked with wrong size! fi # end of 'Makefile' fi if test -f 'ash.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash.s'\" else echo shar: Extracting \"'ash.s'\" \(2289 characters\) sed "s/^X//" >'ash.s' <<'END_OF_FILE' Xbin1 <- function(x,ab=nicerange(x),nbin=50) { X n <- length(x) X if(ab[1]>=ab[2]) X fatal("Interval vector has negative orientation") X if(nbin<=0) X fatal("Number of bin intervals nonpositive") X r<-.Fortran("bin1", X as.single(x), X as.integer(n), X as.single(ab), X as.integer(nbin), X nc=integer(nbin), X nskip=integer(1)) X list(nc=r$nc,ab=ab,nskip=r$nskip) X} X Xash1 <- function(bins,m=5,kopt=c(2,2)){ X nc <- bins$nc X ab <- bins$ab X nbin <- length(nc) X r <- .Fortran("ash1", X as.integer(m), X as.integer(nc), X as.integer(nbin), X as.single(ab), X as.integer(kopt), X t=single(nbin), X f=single(nbin), X single(m), X ier=integer(1)) X if(r$ier==1) print("ash estimate nonzero outside interval ab") X list(x=r$t,y=r$f,m=m,ab=ab,kopt=kopt,ier=r$ier) X} X Xnicerange <- function(x, beta = .1) { X ab <- range(x) # interval surrounding data X del <- ((ab[2] - ab[1]) * beta)/2 X return(c(ab + c( - del, del))) X} X Xbin2 <- function(x,ab,nbin=c(20,20)){ X if(missing(ab)){ X ab <- t(array(c(nicerange(x[,1]),nicerange(x[,2])),c(2,2))) X } X n <- nrow(x) X r <- .Fortran("bin2", X as.single(x), X as.integer(n), X as.single(ab), X as.integer(nbin), X nc=integer(nbin[1]*nbin[2]), X nskip=integer(1)) X list(nc=matrix(r$nc,nbin[1],nbin[2]),ab=ab,nskip=r$nskip) X} X Xash2 <- function(bins,m=c(5,5),kopt=c(2,2)){ X nc <- bins$nc; if(!is.matrix(nc)) stop("bins does not contain bin count matrix") X ab <- bins$ab; if(!is.matrix(ab)) stop("ab not a matrix - should be 2 by 2") X nbin <- dim(nc) X r <- .Fortran("ash2", X as.integer(m), X as.integer(nc), X as.integer(nbin[1]), X as.integer(nbin[2]), X as.single(ab), X as.integer(kopt), X f=single(nbin[1]*nbin[2]), X single(m[1]*m[2]), X ier=integer(1)) X if(r$ier==1) print(" estimate nonzero outside ab rectangle") X list(z=matrix(r$f,nbin[1],nbin[2]), X x=center(ab[1,],nbin[1])[[1]],y=center(ab[2,],nbin[2])[[1]], X ab=ab,m=m,kopt=kopt,ier=r$ier) X} X Xcenter <- function(ab,k) { X h <- diff(ab)/k X list( seq(ab[1]+h/2,by=h,length=k) ) } X X END_OF_FILE if test 2289 -ne `wc -c <'ash.s'`; then echo shar: \"'ash.s'\" unpacked with wrong size! fi # end of 'ash.s' fi if test -f 'ashn.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ashn.s'\" else echo shar: Extracting \"'ashn.s'\" \(3497 characters\) sed "s/^X//" >'ashn.s' <<'END_OF_FILE' Xashn <- function( ix, ioauto=0, master=mymaster, ctable=ct, ff=4, X mv=rep(1,ncol(ix$nc)), square=c(0,450), sun50=0, mmax=10, ireg=0 ) { X X##### S function to compute and plot n-d ASH slice over cube [-1,1]^3 X##### Copyright 1988,1992 David W. Scott X X# ix ----------- data structure containing: X# nc - integer data matrix (n x nvar) - ignores skipped data X# ab - bin range matrix (nvar x 2) X# [ ab[1,1] , ab[1,2] ] contains x-axis range X# [ ab[2,1] , ab[2,2] ] contains y-axis range X# [ ab[3,1] , ab[3,2] ] contains z-axis range X# nbin - vector length nvar; number of bins along each axis X# mmax - maximum allowed value of m; dimensioned (m+1)^3 X# iexp - alternate kernel (0==triangle 1=Epan 2=biweight 3=triweight) X# regress - optional array of regression estimates (in "ix") X Xnvar <- ncol(ix$nc); if(nvar<3) stop("less than 3 data variables") Xif( !is.matrix(ix$nc) | !is.integer(ix$nc) ) stop("ix$nc should be integer matrix") Xif( any(dim(ix$ab) != c(nvar,2)) ) stop("dim of ix$ab should be nvarx2") Xif( !is.integer(ix$nbin) | length(ix$nbin)!=nvar ) stop("ix$nbin should be integer vector nvar") Xif( !is.matrix(ctable) | !is.single(ctable) ) stop("ctable not an single matrix") X X# tri - 220 x 4 matrix X# a. sides connected for triangle X# b. corner that indicates decreasing/increasing X# cuts - list of 500 triangles for the 254 cases X# contour - 254 x 2 (256 cases: 0 and 255 imply no triangles) X# a. pointer to first triangle X# b. number of triangles in this case X# orth? - default 0==> all triangles 1==> only "outside" triangles X Xif( any(dim(master$tri) != c(220,4)) ) stop("master$tri not 220 x 4") Xif( length(master$cuts) != 500 ) stop("master$cuts not 500 long") Xif( any(dim(master$contour) != c(254,2)) ) stop("master$contour not 254 x 2") X X# data have been compressed into integer matrix X# mmax is maximum m allowed for kernel array (default 6) X X mp1 <- rep( mmax+1, 3) # dimensioned (0:mmax) X # WEIGHT ARRAY REALLY has 3 in dimension here X X# ctable: 1 row == background color ###of rows should be odd X# rows 2-(nr+1)/2 are the OUTER colors -- edge to middle X# last rows are INNER colors -- edge to middle X# S refers to these colors as 0 to nr-1 X Xcolor <- integer(4) Xcolor[1] <- 1; color[4] <- nrow(ctable)-1 Xcolor[2] <- max(color[4]/2,color[1]); color[3] <- min(color[2]+1,color[4]) X # color[1:2] towards eye color[3:4] away from eye X Xwk <- single( prod(mp1) ) # STRUCTURE( wk/REAL,ARRAY,3,mp1/ ) # ASH weight array Xxfmax <- single( ix$nbin[1] ) Xyfmax <- single( ix$nbin[2] ) Xzfmax <- single( ix$nbin[3] ) Xtx <- single( ix$nbin[1] ) Xty <- single( ix$nbin[2] ) Xtz <- single( ix$nbin[3] ) Xif ( ireg != 1 ) {ash <- single( prod(ix$nbin[1:3]) ) } # ARRAY X else { if( !is.array(ix$regress) | is.integer(ix$regress) ) { X stop("regress not a real array") } X if( any( dim(ix$regress) != ix$nbin) ) X stop("ix$nbin and ix$regress do not agree") X ash <- ix$regress } Xn <- nrow(ix$nc) X Xz <- .Fortran("ashn", as.integer(ix$nc), as.integer(n),as.integer(nvar), X as.single(ix$ab), as.integer(ix$nbin), ash=as.single(ash), X as.single(xfmax), as.single(yfmax), as.single(zfmax), X as.single(tx), as.single(ty), as.single(tz), X as.single(wk), as.integer(mmax), X as.integer(master$contour), as.integer(master$cuts), X as.integer(master$tri), as.integer(color), X as.integer(square), as.integer(sun50), X as.single(ctable), as.integer(nrow(ctable)), X as.integer(mv), as.single(ff), as.integer(ioauto), as.integer(ireg) ) X Xinvisible( array(z$ash,ix$nbin[1:3]) ) X X} END_OF_FILE if test 3497 -ne `wc -c <'ashn.s'`; then echo shar: \"'ashn.s'\" unpacked with wrong size! fi # end of 'ashn.s' fi if test -f 'compress.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'compress.s'\" else echo shar: Extracting \"'compress.s'\" \(1575 characters\) sed "s/^X//" >'compress.s' <<'END_OF_FILE' Xcompress <- function ( x, ab=t(apply(x,2,nicerange)), nbin=rep(30,ncol(x)) ) { X n <- nrow(x); nvar <- ncol(x); if( ncol(ab) != 2 ) stop("ab should have 2 columns") X if( nrow(ab)!=nvar||length(nbin)!=nvar) stop("dimensions of x,ab,nbin do not agree") X if( min(nbin) <= 0 ) stop("nbin must be positive") X X##### S function to compress data into binned integer form X##### Copyright 1988,1992 David W. Scott X X# x - input trivariate data matrix (n x nvar) X# ab - input bin range matrix (nvar x 2) X# x binned over [ ab[1,1], ab[1,2] ) X# y binned over [ ab[2,1], ab[2,2] ) X# z binned over [ ab[3,1], ab[3,2] ) etc. X# nbin - vector length nvar - number of bins along each axis X# nout - subset of n data points in bin range X# nc - output matrix (nout x nvar) containing trivariate bin numbers X Xdel <- (ab[,2]-ab[,1])/nbin # bin widths Xiskip <- 0 # number of points not in trivariate bins Xnout <- 0 # number of points in output array X X # loop on data Xnc <- matrix(0,n,nvar); k <- NULL Xfor( i in 1:nvar ) { if( ab[i,1] >= ab[i,2]) stop("ab array incorrectly ordered") X nc[,i] <- floor( (x[,i]-ab[i,1])/del[i] + 1.0 ) X k <- unique( c( k, seq(n)[ nc[,i]<1 | nc[,i]>nbin[i] ] ) ) } X Xiskip <- length(k); if(iskip>0) {nc <- nc[-k,]; X print(paste(as.character(iskip),"point(s) not in hyper-rectangle")) } X Xstorage.mode(nc) <- "integer" Xstorage.mode(ab) <- "single" Xstorage.mode(nbin) <- "integer" Xreturn(nc=nc,ab=ab,nbin=nbin) } X Xnicerange <- function(x,beta=.1) { ab <- range(x) # interval surrounding data X del <- (ab[2]-ab[1])*beta/2; return( c( ab+c(-del,del) ) ) } X END_OF_FILE if test 1575 -ne `wc -c <'compress.s'`; then echo shar: \"'compress.s'\" unpacked with wrong size! fi # end of 'compress.s' fi if test -f 'init.s' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'init.s'\" else echo shar: Extracting \"'init.s'\" \(874 characters\) sed "s/^X//" >'init.s' <<'END_OF_FILE' Xmymaster <- dget("mymaster") Xct <- dget("ct") X XX11.ash2 <- c("-xrm 'splus*colors: black red 48 white blue 48 cyan'", X "-xrm splus@nHalftones:0", "-xrm splus*halftoneImage:No") XOL.col <- "-xrm 'sgraphOpenLook.colorSchemes : name: \"Scott\"; background : black; lines : white; text : white; polygons : red 48 white blue 48 cyan; images : red 48 white blue 48 cyan'" X Xxyz <- compress( matrix( rnorm(3000),1000,3 ) ) # an example X Xstorage.mode(ct) <- "single" # dput/dget not perfect Xstorage.mode(mymaster$contour) <- "integer" Xstorage.mode(mymaster$cuts) <- "integer" Xstorage.mode(mymaster$tri) <- "integer" X X.First.lib <- function(lib.loc, section) { Xpath <- paste(lib.loc, section,"loadmod.o", sep="/") Xif(!is.loaded(symbol.C("exvalue"))) dyn.load(path) Xpath <- paste(lib.loc, section,"ash.o", sep="/") Xif(!is.loaded(symbol.For("ash1"))) dyn.load(path) X} END_OF_FILE if test 874 -ne `wc -c <'init.s'`; then echo shar: \"'init.s'\" unpacked with wrong size! fi # end of 'init.s' fi if test -f 'ash1.d' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash1.d'\" else echo shar: Extracting \"'ash1.d'\" \(821 characters\) sed "s/^X//" >'ash1.d' <<'END_OF_FILE' X.BG X.FN ash1 X.TL Xash1: Computes univariate averaged shifted histogram (polynomial kernel) X.CS Xash1(bins, m, kopt) X.PP X.AG bins X(input list) $nc=integer vector of bin counts and $ab=bin interval X.AG m X(input) optional integer smoothing parameter; default=5. X.AG kopt X(input) vector of length 2 specifying the kernel, which Xis proportional to ( 1 - abs(i/m)**kopt(1) )**kopt(2); X(2,2)=biweight (default); (0,0)=uniform; (1,0)=triangle; X(2,1)=Epanechnikov; (2,3)=triweight. X.RT Xreturns structure suitable for input to "plot" X.RC x=t Xvector of bin center locations X.RC y=f Xvector of ash estimates X.RC ier X0=normal exit; 1=estimate nonzero outside interval ab X.EX Xx <- rnorm(100) # data Xf <- ash1(bin1(x,nbin=50),5) # compute ash estimate Xplot( f , type="l" ) # line plot of estimate X.KW density estimation X.WR END_OF_FILE if test 821 -ne `wc -c <'ash1.d'`; then echo shar: \"'ash1.d'\" unpacked with wrong size! fi # end of 'ash1.d' fi if test -f 'ash2.d' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash2.d'\" else echo shar: Extracting \"'ash2.d'\" \(722 characters\) sed "s/^X//" >'ash2.d' <<'END_OF_FILE' X.BG X.FN ash2 X.TL Xash2: Compute bivariate ASH estimate (product polynomial kernel) X.CS Xash2(bins, m, kopt) X.PP X.AG bins X(input list) bin count matrix nc and interval matrix ab from "bin2" X.AG m X(input integer vector of length 2) x and y direction smoothing Xparameters. Default is 5 by 5. X.RT XMatrix of ASH estimates returned. XComponents x,y,z can be given to the contour function directly. XOther input variables returned in list for record keeping. X.EX XContinuing example from help("bin2") X X1. To use with "persp" X Xm <- c(5,5) Xf <- ash2(bins,m) Xpersp( f$z / max(f$z) ) X X2. To use with "contour" (continuing part 1) X Xcontour( f ) X X# note: it is a good idea to not plot contours at level = 0 X X.KW density estimation X.WR END_OF_FILE if test 722 -ne `wc -c <'ash2.d'`; then echo shar: \"'ash2.d'\" unpacked with wrong size! fi # end of 'ash2.d' fi if test -f 'ashn.d' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ashn.d'\" else echo shar: Extracting \"'ashn.d'\" \(725 characters\) sed "s/^X//" >'ashn.d' <<'END_OF_FILE' X.BG X.FN ashn X.TL Xvisualization of multi-dimensional averaged shifted histograms X.CS Xashn(ix, ioauto=0, master=mymaster, ctable=ct, ff=4, mv=rep(1, ncol(ix$nc)), square=c(0, 450), sun50=0, mmax=10, ireg=0) X.RA X.AG ix Xdata structure containing: Xnc Xinteger data matrix (n x nvar) - ignores skipped data Xab Xbin range matrix (nvar x 2) X[ ab[1,1] , ab[1,2] ] contains x-axis range X[ ab[2,1] , ab[2,2] ] contains y-axis range X[ ab[3,1] , ab[3,2] ] contains z-axis range Xnbin Xvector length nvar; number of bins along each axis Xmmax Xmaximum allowed value of m; dimensioned (m+1)^3 Xiexp Xalternate kernel (0==triangle 1=Epan 2=biweight 3=triweight) Xregress Xoptional array of regression estimates (in "ix") X.KW density estimation X.WR END_OF_FILE if test 725 -ne `wc -c <'ashn.d'`; then echo shar: \"'ashn.d'\" unpacked with wrong size! fi # end of 'ashn.d' fi if test -f 'bin1.d' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'bin1.d'\" else echo shar: Extracting \"'bin1.d'\" \(741 characters\) sed "s/^X//" >'bin1.d' <<'END_OF_FILE' X.BG X.FN bin1 X.TL Xbin1: Function to compute array of bin counts for a data vector X.CS Xbin1(x, ab, nbin=50) X.PP X.AG x X(input) data vector X.AG ab X(input vector of length 2): half-open interval for bins [a,b). If no Xvalue is specified, the range of x is stretched by 5% at each end and Xused the interval. X.AG nbin X(input integer): number of bins desired. Default 50. X.AG opt X(input OPTIONAL logical): to suppress output messages, set opt=TRUE X.RT Xbin1 returns a list including the vector of integer bin counts and Xthe ab vector and the number of points outside the ab interval. X.EX Xx <- rnorm(100) # data vector Xab <- c(-5,5) # bin interval Xbins <- bin1(x,ab,10) # bin x into 10 bins over ab X.KW density estimation X.WR END_OF_FILE if test 741 -ne `wc -c <'bin1.d'`; then echo shar: \"'bin1.d'\" unpacked with wrong size! fi # end of 'bin1.d' fi if test -f 'bin2.d' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'bin2.d'\" else echo shar: Extracting \"'bin2.d'\" \(809 characters\) sed "s/^X//" >'bin2.d' <<'END_OF_FILE' X.BG X.FN bin2 X.TL Xbin2: Bin bivariate data x X.CS Xbin2(x, ab, nbin) X.PP X.AG x X(input matrix with 2 columns) data sample X.AG ab X(input 2 x 2 matrix) rows 1 and 2 contain x and y axis bin intervals, Xrespectively. If not specified, the ranges are stretched by 5% Xat each end for each dimension. X.AG nbin X(input vector of length 2) number of bins along x and y axes. Default Xis 20 by 20. X.RT Xbin2 returns a list including the bivariate bin matrix (see "persp" Xfunction for orientation) and the number of points outside the ab rectangle. X.EX Xx <- matrix( rnorm(200), 100 , 2) # bivariate normal n=100 Xab <- matrix( c(-5,-5,5,5), 2, 2) # interval [-5,5) x [-5,5) Xnbin <- c( 20, 20) # 400 bins Xbins <- bin2(x, ab, nbin) # bin counts,ab,nskip X.KW density estimation X.WR END_OF_FILE if test 809 -ne `wc -c <'bin2.d'`; then echo shar: \"'bin2.d'\" unpacked with wrong size! fi # end of 'bin2.d' fi if test -f 'cubeint.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'cubeint.c'\" else echo shar: Extracting \"'cubeint.c'\" \(905 characters\) sed "s/^X//" >'cubeint.c' <<'END_OF_FILE' Xstruct { X int a, b; X} inner[6] = {{4,6}, {2,6}, {2,3}, {1,3}, {1,5}, {4,5}}; X Xcubeint(origin, basis, F, X, size) Xfloat *origin, **basis, *F, **X; Xint *size; X{ X int i, npatch; X float c[24], *corner[8], *vertex[4], value[4], **x; X X for(i = 0; i < 8; i++) X corner[i] = c + 3*i; X make_corners(origin, basis, corner); X vertex[0] = corner[0]; value[0] = F[0]; X vertex[3] = corner[7]; value[3] = F[7]; X x = X; X npatch = 0; X for(i = 0; i < 6; i++) { X vertex[1] = corner[inner[i].a]; value[1] = F[inner[i].a]; X vertex[2] = corner[inner[i].b]; value[2] = F[inner[i].b]; X size[i] = simpint(vertex, value, x); X x += size[i]; X npatch += size[i] > 0; X } X return(npatch); X} X Xmake_corners(o, b, c) Xfloat *o, **b, **c; X{ X int i, j, k, l, n; X X n = 0; X for(k = 0; k < 2; k++) for(j = 0; j < 2; j++) for(i = 0; i < 2; i++) { X for(l = 0; l < 3; l++) X c[n][l] = o[l] + i*b[0][l] + j*b[1][l] + k*b[2][l]; X n++; X } X} X X END_OF_FILE if test 905 -ne `wc -c <'cubeint.c'`; then echo shar: \"'cubeint.c'\" unpacked with wrong size! fi # end of 'cubeint.c' fi if test -f 'dumpit.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dumpit.c'\" else echo shar: Extracting \"'dumpit.c'\" \(902 characters\) sed "s/^X//" >'dumpit.c' <<'END_OF_FILE' X#include X#include "S.h" X Xstatic FILE *f = NULL; X XF77_SUB(startdump) (n) Xlong *n; X{ X char filename[20]; X X sprintf(filename, "figure.%d.quad", *n); X f = fopen(filename, "w"); X if(f == NULL) X fprintf(stderr, "cannot create %s for writing\n", filename); X fprintf(f, "CNQUAD\n"); X} X XF77_SUB(dumpit) (p, q, r, n, a) Xfloat *p, *q, *r, *n, *a; X{ X if(f == NULL) X return; X fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", X p[0], p[1], p[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); X fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", X q[0], q[1], q[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); X fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", X r[0], r[1], r[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); X fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", X p[0], p[1], p[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); X fprintf(f, "\n"); X} X XF77_SUB(enddump) () X{ X if(f != NULL) X fclose(f); X} END_OF_FILE if test 902 -ne `wc -c <'dumpit.c'`; then echo shar: \"'dumpit.c'\" unpacked with wrong size! fi # end of 'dumpit.c' fi if test -f 'myinput.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'myinput.c'\" else echo shar: Extracting \"'myinput.c'\" \(4705 characters\) sed "s/^X//" >'myinput.c' <<'END_OF_FILE' X#include X#include "S.h" Xstatic FILE *fin; XF77_SUB(myio) (iopt, intx, realx, n) Xint *iopt; int *intx; float *realx; int *n; X{ X int i; X switch (*iopt) X { X case 0: X fin = fopen("stream.s", "r"); X break; X case 1: X for (i=0; i < *n; i++) { X fscanf(fin, " %d", intx+i); X } X break; X case 2: X for (i=0; i < *n; i++) { X fscanf(fin, " %f", realx+i); X } X break; X case 3: X fclose(fin); X break; X } X} X X XF77_SUB(mess) (iopt, intx, realx, n) Xint *iopt; int *intx; float *realx; int *n; X{ X int i; X Xswitch (*iopt) { X case 0: X printf("\nACTION OPTION") ; X printf("\n 0=plot 1=t 2=m 3=f 4=axes 5=cube 6=eye 7=slicelimit 8=sliceincr"); X printf("\n 9=iris 10=orth 11=print 12=coltbl lookup 13=blank screen") ; X printf("\n 14=alpha 15=fmaxin 16=kernel 99=quit > ") ; X fflush(stdout); break; X case 1: X printf("\n z slice:"); fflush(stdout); break; X case 2: X printf("%d",*intx); fflush(stdout); break; X case 3: X printf("\ninside ASH check some calling parameters") ; X printf("\n nvar = %d",*intx) ; printf("\n n = %d",*n); fflush(stdout); break; X case 4: X printf("\n mv ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; X printf("\n hv ="); for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; X fflush(stdout); break; X case 5: X printf(" echoing: ichoice = %d",*intx); break; X case 6: X printf("\n number of triangle patches found = %d",*intx); fflush(stdout); break; X case 7: X printf("\n sorting...."); fflush(stdout); break; X case 8: X printf("sort finished....beginning plot"); fflush(stdout); break; X case 9: X printf("\n"); fflush(stdout); break; X case 10: X printf("\n enter slice [ivar slice] = %d %d > ",*(intx+0),*(intx+1)); fflush(stdout); break; X case 11: X printf("\n mv ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; X printf("; enter m > "); fflush(stdout); break; X case 12: X printf("\nError: smoothing parameter mmax exceeded for m[1,2,3], mmax = %d",*intx); X fflush(stdout); break; X case 13: X printf("\nError: smoothing parameters not positive"); fflush(stdout); break; X case 14: X printf("\n ff = %g enter f level > ",*realx); fflush(stdout); break; X case 15: X printf("\n draw axes? 1=yes 0=no > "); fflush(stdout); break; X case 16: X printf("\n draw cube? 1=yes 0=no > "); fflush(stdout); break; X case 17: X printf("\n eye ="); for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; X printf("; enter eye > "); fflush(stdout); break; X case 18: X printf("\nError: eye=0 try again"); fflush(stdout); break; X case 19: X printf("\n ix iy iz limits ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; X printf("; enter > "); fflush(stdout); break; X case 20: X printf("\nError: out of range: ix iy iz maxrange ="); X for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; fflush(stdout); break; X case 21: X printf(" enter ix,iy,iz increments > "); fflush(stdout); break; X case 22: X printf("\nError: increments less than 1;"); fflush(stdout); break; X case 23: X printf("\n color offset 999 for an iris; 1=yes > "); fflush(stdout); break; X case 24: X printf("\n show orthogonal triangles only? 1=yes > "); fflush(stdout); break; X case 25: X printf("\n do you want to print 3D vertices 0-1-2-3 for 00 01 10 11 > "); X fflush(stdout); break; X case 26: X printf("\n color table lookup ictopt = %d; enter > ",*intx); X fflush(stdout); break; X case 27: X printf("\n blank screen? 0=no 1=yes; enter > "); fflush(stdout); break; X case 28: X printf("\n texture=a[1,2,3] alpha=a[4]: "); X for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; X printf("; enter > "); fflush(stdout); break; X case 29: X printf("\n fmax = %f; fmaxin = %f; enter fmaxin (neg ==> use fmax) > ", X *(realx+0),*(realx+1)); fflush(stdout); break; X case 30: X printf("\nWarning: only 3 dimensions"); fflush(stdout); break; X case 31: X printf("\n working on point "); fflush(stdout); break; X case 32: X printf("\n fmax = %f; fmaxin = %f",*(realx+0),*(realx+1)); fflush(stdout); break; X case 33: X printf("\n n show nskip = %d %d %d",*(intx+0),*(intx+1),*(intx+2)); fflush(stdout); break; X case 34: X printf("\n ff = %g",*realx); fflush(stdout); break; X case 35: X printf("; c = %f",*realx); fflush(stdout); break; X case 36: X printf("\n t slice(s) ="); X for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; fflush(stdout); break; X case 37: X printf("\n iexp = %d 0=tri 1=Epan 2=biwt 3=triwt",*intx); X printf("; enter iexp > "); fflush(stdout); break; X case 38: X printf("\nWarning: iexp>0 kernel shape correct but not normalized"); X fflush(stdout); break; X case 39: X printf("."); fflush(stdout); break; X case 40: X printf(" %d",*intx); fflush(stdout); break; X } X} X END_OF_FILE if test 4705 -ne `wc -c <'myinput.c'`; then echo shar: \"'myinput.c'\" unpacked with wrong size! fi # end of 'myinput.c' fi if test -f 'simpint.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'simpint.c'\" else echo shar: Extracting \"'simpint.c'\" \(2088 characters\) sed "s/^X//" >'simpint.c' <<'END_OF_FILE' X#include Xstruct connection { X int n, start[5], end[5]; X} C[5][5] = { X {0},{3,{0,0,0},{1,2,3}},{4,{0,2,1,3},{2,1,3,0}},{3,{3,3,3},{0,1,2}},{0}, X {0},{3,{0,0,1},{2,3,1}},{3,{0,1,2},{3,3,2}},{0},{0}, X {0},{3,{0,1,2},{3,1,2}},{0},{0},{0}, X {3,{0,1,2},{0,1,2}},{3,{1,2,3},{1,2,3}},{0},{0},{0}, X {0},{0},{0},{0},{0} X}; X Xsimpint(P, F, X) Xfloat **P, *F, **X; X{ X float *p[4], f[4]; X int ord[4], nneg, nzero, nvertex, orientation, i, j; X struct connection c; X X sort4(F, ord); X nneg = nzero = 0; X for(i = 0; i < 4; i++) { X p[i] = P[ord[i]]; X f[i] = F[ord[i]]; X nneg += f[i] < 0; X nzero += f[i] == 0; X } X c = C[nzero][nneg]; X nvertex = c.n; X if(nvertex == 0) X return(0); X for(i = 0; i < nvertex; i++) { X int u = c.start[i], v = c.end[i], w; X double a; X X if(u == v) X for(j = 0; j < 3; j++) X X[i][j] = p[u][j]; X else { X if(f[u] < 0) {w = u; u = v; v = w;} X a = f[u] / (f[u] - f[v]); X for(j = 0; j < 3; j++) X X[i][j] = (1 - a) * p[u][j] + a * p[v][j]; X } X } X if(nneg) X orientation = which_side(X[0], X[1], X[2], p[0]); X else X orientation = -which_side(X[0], X[1], X[2], p[3]); X if(orientation < 0) { X float *tmp[4]; X for(i = 0; i < nvertex; i++) X tmp[i] = X[i]; X for(i = 0; i < nvertex; i++) X X[i] = tmp[nvertex-i-1]; X } X return(nvertex); X} X X#define ORDER(p,q,r) {o[0]=p; o[1]=q; o[2]=r;} X Xsort4(a, o) Xfloat *a; Xint *o; X{ X int i, j; X X if(a[0] < a[1]) X if(a[2] < a[0]) ORDER(2,0,1) X else if(a[2] > a[1]) ORDER(0,1,2) X else ORDER(0,2,1) X else X if(a[2] < a[1]) ORDER(2,1,0) X else if(a[2] > a[0]) ORDER(1,0,2) X else ORDER(1,2,0) X for(i = 0; i < 3; i++) X if(a[3] < a[o[i]]) X break; X for(j = 3; j > i; j--) X o[j] = o[j-1]; X o[i] = 3; X} X Xwhich_side(x0, x1, x2, p) Xfloat *x0, *x1, *x2, *p; X{ X int i; X float v1[3], v2[3], v[3], cross[3], dot; X X for(i = 0; i < 3; i++) { X v1[i] = x2[i] - x1[i]; X v2[i] = x0[i] - x1[i]; X v[i] = p[i] - x1[i]; X } X cross[0] = v1[1]*v2[2] - v1[2]*v2[1]; X cross[1] = v1[2]*v2[0] - v1[0]*v2[2]; X cross[2] = v1[0]*v2[1] - v1[1]*v2[0]; X dot = 0; X for(i = 0; i < 3; i++) X dot += cross[i] * v[i]; X return(dot > 0 ? 1 : -1); X} X END_OF_FILE if test 2088 -ne `wc -c <'simpint.c'`; then echo shar: \"'simpint.c'\" unpacked with wrong size! fi # end of 'simpint.c' fi if test -f 'writeit.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'writeit.c'\" else echo shar: Extracting \"'writeit.c'\" \(1254 characters\) sed "s/^X//" >'writeit.c' <<'END_OF_FILE' X#include X#include "S.h" X Xstatic FILE *f = NULL; X XF77_SUB(startwrite) (n) Xlong *n; X{ X char filename[20]; X X sprintf(filename, "fig.%d.quad", *n); X f = fopen(filename, "w"); X if(f == NULL) X fprintf(stderr, "cannot create %s for writing\n", filename); X fprintf(f, "CQUAD\n"); X} X X#define VEC(x, m) \ X{ \ X for(k = 0; k < m; k++) \ X fprintf(f, " %g", x[k]); \ X for(k = 0; k < 4; k++) \ X fprintf(f, " %g", color[k]); \ X fprintf(f, "\n"); \ X} X Xstatic float b0[3] = {1,0,0}, b1[3] = {0,1,0}, b2[3] = {0,0,1}; Xstatic float *basis[3] = {b0, b1, b2}; X XF77_SUB(docube) (origin, cube_size, values, contour, color) Xfloat *origin, *cube_size, *values, *contour, *color; X{ X int i, j, k, size[6], npatch; X float x[72], *points[24], **patch, F[8]; X X if(f == NULL) X return; X for(i = 0; i < 3; i++) X basis[i][i] = cube_size[i]; X for(i = 0; i < 8; i++) X F[i] = values[i] - *contour; X for(i = 0; i < 24; i++) X points[i] = x + 3*i; X npatch = cubeint(origin, basis, F, points, size); X patch = points; X for(i = 0; i < 6; i++) { X if(size[i] == 0) X continue; X for(j = 0; j < size[i]; j++) X VEC(patch[j], 3) X if(size[i] == 3) X VEC(patch[0], 3) X fprintf(f, "\n"); X patch += size[i]; X } X} X XF77_SUB(endwrite) () X{ X if(f != NULL) X fclose(f); X} X END_OF_FILE if test 1254 -ne `wc -c <'writeit.c'`; then echo shar: \"'writeit.c'\" unpacked with wrong size! fi # end of 'writeit.c' fi if test -f 'tmp2.r' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'tmp2.r'\" else echo shar: Extracting \"'tmp2.r'\" \(25958 characters\) sed "s/^X//" >'tmp2.r' <<'END_OF_FILE' X#-July 20, 1988 X# plot triangularization September 13, 1990 X#-August 1, 1992 automatic I/O option X Xsubroutine ashn ( nc, n, nvar, ab, nbin, ash, xfmax, yfmax, zfmax, X tx, ty, tz, wk, mmax, contour, cuts, tri, color, X square, sun50, ictable, ict, mv, ff, ioauto, ireg ) X X dimension hv(10),nbin(nvar),mv(nvar),itbin(10) # nvar <=10 assumed X dimension nc(n,nvar),ab(nvar,2),ash(nbin(1),nbin(2),nbin(3)) X dimension v1(3),v2(3),v3(3),icolor(5),eye(3) X dimension ixyz(3),ixyzl(3),ixyzr(3) X dimension ileft(3),irght(3),idel(3) X dimension xfmax(nbin(1)),yfmax(nbin(2)),zfmax(nbin(3)) X dimension tx(nbin(1)),ty(nbin(2)),tz(nbin(3)) X dimension tdel(3),txyz(3),xyz(3,12),ee(3,3),ttxyz(3),tnormal(3),talpha(4) X dimension px(12),py(12),pz(12) X dimension wk(0:mmax,0:mmax,0:mmax) # kernel weights X dimension cb(16,3) # cube X X integer ictable(ict,1) # color table lookup X integer orth, contour(254,2), cuts(500), tri(220,4) X dimension vout(3),vout2(3) # poutward normal vector X integer ic(3,11) # integer coordinates of corner but [0,1]^3 + 3 faces X dimension as(0:1,0:1,0:1),rc(3,11) X integer color(4), square(2), sun50 X X X dimension pxy(3,2,25000),ztri(25000) # big long list here X real angle(25000) # pos==>normal towards eye; neg==>away X integer iout(10) X dimension rout(10) X integer iperm(25000) # sorted ztri pointers (could combine with ishow X # if rewrote sort to handle subsets other X # than 1,2,3,4,.....actually since X # orth is input, shouldn't need this X X# hv--bin width in data units X# nbin--number bins each dimension X# mv--number shifts in averaging X X Xreal am(200) Xcommon/bgrp/am X Xcommon/encbfc/encbf0 Xcommon/encbfl/iqqbfl,iqqfil Xcharacter*256 encbf0 ; integer iqqbfl,iqqfil X Xcommon/decbfc/decbf0 Xcharacter*512 decbf0 Xcommon/decbfl/jqqbfl,jqqinf,jqqflp,jqqfle,jqqlsc,jqdchk Xinteger jqqbfl,jqqinf,jqqflp,jqqlsc; logical jqqfle,jqdchk X X X external qproject, qtriout, qoutward X X data ntrimax/25000/ X# cube X data cb/ 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1, 1,-1,-1,-1,-1, X 1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, X -1,-1,-1,-1, 1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1/ X## data ff/.50/ # default slice level --- 9/11/90 IN CALL X X# corners X data ic/ 0,0,0, 1,0,0, 1,1,0, 0,1,0, X 0,0,1, 1,0,1, 1,1,1, 0,1,1, X 0,0,0, 0,0,0, 0,0,0/ # fills columns...last three refer to corner 1 X data rc/ -1.,-1.,-1., 1.,-1.,-1., 1., 1.,-1., -1., 1.,-1., X -1.,-1., 1., 1.,-1., 1., 1., 1., 1., -1., 1., 1., X 0., 0.,-1., -1., 0., 0., 0.,-1., 0./ # fills columns ... last 3 are faces X X data tnormal/1.,1.,1./ # dummy normal for iris X data talpha/1.,0.,0.,1./ # texture plus alpha channel X X if (ioauto.eq.1) call myio(0,iout,rout,0) X### call beginz # a null call to put iblank right X { Xrxxxx1 = 0.1; rxxxx2 = 0.1; rxxxx3 = 0.1; rxxxx4 = 0.1 Xcall spmarz(rxxxx1,rxxxx2,rxxxx3,rxxxx4) X; } X s3 = sqrt(3.0) X { Xam(61) = square(1) Xam(62) = square(2) Xam(63) = square(1) Xam(64) = square(2) Xcall zscalz X; } X do j=1,8{do i=1,3{rc(i,j)=rc(i,j)/s3}} # direction vectors of corners (j=1,8 not 11) X # 3 face vectors already length 1 X call mess(3,nvar,0,n) # echoing nvar,n X X# default values X X npict = 0 # number of pictures X iprint = 0 # print triangle vertices CHANGED X iris=1; icoloroffset = 999 # color offset (0 for Sun, 999 for Iris) 9/10/90 CHANGED X inew = 1 # new parameters, recompute ash X iblank = 1 # blank out screen X orth = 0 # default plot all triangles X iaxis = 1 # slice on x-axis X ix1=1; iy1=1; iz1=1; ix2=nbin(1)-1; iy2=nbin(2)-1; iz2=nbin(3)-1 # slice limits X ix3=1; iy3=1; iz3=1 # slice increments X icolor(1)=4; icolor(2)=1; icolor(3)=0; icolor(4)=3; icolor(5)=2 X iexp=0; iexp1=1; iexp2=1; if(iexp>0) {iexp1=2; iexp2=iexp } # kernel option X wk(0,0,0)=1.0 # kernel array is 1 x 1 x 1 X##INCALL## do i=1,nvar { mv(i)=1} # default smoothing parameters CHANGED X do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) X do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) X do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} X do i=1,nvar { hv(i)=(ab(i,2)-ab(i,1))/nbin(i) } X X call mess(4,mv,hv,nvar) # echoing mv,hv X X eye(1)=1; eye(2)=1; eye(3)=0 X v1(1)=-sqrt(2.)/2; v1(2)=-v1(1); v1(3)=0 X v2(1)=0; v2(2)=0; v2(3)=1 X v3(1)=v1(2); v3(2)=v3(1); v3(3)=0 X fmax=0.0; fmaxin=-1000.0 # fmaxin max of f(.) X iskip=0; ishow=0 # numbers points skipped and shown X do i=1,3{ileft(i)=1; irght(i)=nbin(i); idel(i)=1}# select slices along axes to display X do i=4,nvar{itbin(i)=1} # default t-slice X ifleft = 1; ifright = 1; ifdel = 1 # select density levels to display X iqaxis = 1 # draw coordinate axes at center X icube = 1 # draw unit cube X ictopt = 0 # use color table lookup values (right for postscript!) CHANGED X X goto 6663 X5555 continue X X# change default values X Xcall mess(0,iout,rout,0) # ACTION OPTIONS X Xif(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(ichoice,.true.) X}} else {call myio(1,ichoice,rout,1)} Xif(ichoice!=99) call mess(5,ichoice,rout,0) # echoing ichoice Xif(ichoice==0) goto 6663 Xswitch(ichoice){ Xcase 1: # change t slice X if(nvar==3) { call mess(30,0,0,0); goto 5555} # warning: only 3 dimensions X inew=1 # recompute ash X do i=4,nvar {iout(1)=i;iout(2)=itbin(i) X call mess(10,iout,0,2) # ivar slice X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(itbin(i),.true.) X}} X else {call myio(1,itbin(i),rout,1)} } Xcase 2: # change smoothing parameters X 6671 inew = 1 # recompute ash X call mess(11,mv,0,nvar) # current m values X if(nvar==3) { X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(mv(1),.true.) Xcall deci00(mv(2),.true.) Xcall deci00(mv(3),.true.) X}} X else {call myio(1,mv,rout,3)} } X else { X if(nvar==4) { if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(mv(1),.true.) Xcall deci00(mv(2),.true.) Xcall deci00(mv(3),.true.) Xcall deci00(mv(4),.true.) X}} X else {call myio(1,mv,rout,4)} } X else { do i=1,nvar{ iout(1)=i;iout(2)=mv(i) X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(mv(i),.true.) X}} X else {call myio(1,mv(i),rout,1)} } } } X if(max0(mv(1),mv(2),mv(3)).gt.mmax) { call mess(12,mmax,0,0); goto 6671} # m too big X mmin=mmax; do i=1,nvar{mmin=min0(mmin,mv(i))}; if(mmin<1) { X call mess(13,0,0,0) ; goto 6671} # m too small X do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) # change ASH wts X do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) X do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} Xcase 3: # change contour levels X call mess(14,0,ff,1) # new f level X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall decr00(ff) X}} else {call myio(2,iout,ff,1)} Xcase 4: # draw coordinate axes? X call mess(15,0,0,0) # draw axes? X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(iqaxis,.true.) X}} else {call myio(1,iqaxis,rout,1)} Xcase 5: # draw unit cube? X call mess(16,0,0,0) # draw cube? X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(icube,.true.) X}} else {call myio(1,icube,rout,1)} Xcase 6: # change eye X 6662 call mess(17,0,eye,3) # enter eye X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall decr00(eye(1)) Xcall decr00(eye(2)) Xcall decr00(eye(3)) X}} X else {call myio(2,iout,eye,3)} X e1 = eye(1); e2 = eye(2); e3 = eye(3) X rho = e1**2 + e2**2 + e3**2; rho=sqrt(rho) X if(rho.le.0.) {call mess(18,0,0,0); goto 6662} # eye is zero X phi = acos(e3/rho); pi = 3.14159265 X if( e1!=0. | e2!=0. ) theta = atan2(eye(2),eye(1)) # atan2 handles special cases X if( e1==0. & e2==0. & e3>0. ) theta=-pi/2 X if( e1==0. & e2==0. & e3<0. ) theta= pi/2 # don't have RATFOR manual X sp=sin(phi); cp=cos(phi); st=sin(theta); ct=cos(theta) X v1(1)=-st; v1(2)=ct; v1(3)=0.0 X v2(1)=-ct*cp; v2(2)=-st*cp; v2(3)=sp X v3(1)=ct*sp; v3(2)=st*sp; v3(3)=cp #right hand system (not graphics LHS) Xcase 7: # axis contour slices shown X 6677 iout(1)=ix1;iout(2)=ix2;iout(3)=iy1;iout(4)=iy2;iout(5)=iz1;iout(6)=iz2 X call mess(19,iout,0,6) # ix iy iz limits X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(ix1,.true.) Xcall deci00(ix2,.true.) Xcall deci00(iy1,.true.) Xcall deci00(iy2,.true.) Xcall deci00(iz1,.true.) Xcall deci00(iz2,.true.) X}} X else { call myio(1,ix1,rout,1); call myio(1,ix2,rout,1); call myio(1,iy1,rout,1); X call myio(1,iy2,rout,1); call myio(1,iz1,rout,1); call myio(1,iz2,rout,1) } X if(ix1<1|iy1<1|iz1<1|ix2>=nbin(1)|iy2>=nbin(2)|iz2>=nbin(3)) { X iout(1)=1;iout(2)=nbin(1)-1;iout(3)=1;iout(4)=nbin(2)-1;iout(5)=1;iout(6)=nbin(3)-1 X call mess(20,iout,0,6); goto6677 } # ix,iy,iz out of limits Xcase 8: # contour slice increments X 6678 call mess(21,0,0,0) X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(ix3,.true.) Xcall deci00(iy3,.true.) Xcall deci00(iz3,.true.) X}} X else { call myio(1,ix3,rout,1); call myio(1,iy3,rout,1); call myio(1,iz3,rout,1) } X if(min0(ix3,iy3,iz3)<1) {call mess(22,0,0,0); goto 6678} Xcase 9: #-which contour levels to plot? X call mess(23,0,0,0) # color offset 999 for an iris X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(iris,.true.) X}} else {call myio(1,iris,rout,1)} X if(iris.ne.1) icoloroffset=0 # 9/10/90 X if(iris.eq.1) icoloroffset=999 # 9/10/90 Xcase 10: # show only "outside" triangles X call mess(24,0,0,0) #show orthogonal triangles only? X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(orth,.true.) X}} else {call myio(1,orth,rout,1)} Xcase 11: # print triangle vertices X call mess(25,0,0,0) # do you want to print 3D vertices 0-1-2-3 for 00 01 10 11 X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(iprint,.true.) X}} else {call myio(1,iprint,rout,1)} X if(iprint==0) { call enddump; call endwrite } Xcase 12: # use color table lookup or not! X call mess(26,ictopt,0,0) # ictopt= X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(ictopt,.true.) X}} else {call myio(1,ictopt,rout,1)} Xcase 13: # blank screen option X call mess(27,0,0,0) # blank screen 0-1? X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(iblank,.true.) X}} else {call myio(1,iblank,rout,1)} Xcase 14: # read alpha channel stuff X call mess(28,0,talpha,4) # texture=a[1,2,3] alpha=a[4] X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall decr00(talpha(1)) Xcall decr00(talpha(2)) Xcall decr00(talpha(3)) Xcall decr00(talpha(4)) X}} X else {call myio(2,iout,talpha,4)} Xcase 15: # input fmaxin X rout(1)=fmax;rout(2)=fmaxin X call mess(29,0,rout,2) # enter fmaxin: neg==> use fmax X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall decr00(fmaxin) X}} else {call myio(2,iout,fmaxin,1)} X if(fmaxin.le.0.0) fmaxin=fmax Xcase 16: # change kernel iexp 0==triangle 1==Epan 2==Biweight 3==Triwt X 6971 inew = 1 # recompute ash X call mess(37,iexp,0,1) # current "iexp" value X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(iexp,.true.) X}} X else {call myio(1,iexp,rout,1)} X if(iexp<0) { goto 6971} # iexp negative X if(iexp>0) call mess(38,0,0,0) X iexp1=1; iexp2=1; if(iexp>0) {iexp1=2; iexp2=iexp } # kernel option X do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) # change ASH wts X do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) X do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} Xcase 99: #quit X 9000 if (ioauto.eq.1) call myio(3,iout,rout,0) X call enddump; call endwrite X call mess(9,0,0,0) X goto 9005 X } X Xgoto 5555 X X6663 continue X X if (inew.eq.0) goto 6611 # do not have to recompute ash X X#-zero ASH array -- computed at centers of trivariate bins X X6666 call mess(31,0,0,0) # initializing ASH array; working on point X Xfmax = 0. # ASH value at mode Xiskip =0 ; ishow =0 X Xdo iz = 1,nbin(3) { do iy = 1,nbin(2) { do ix = 1,nbin(1) { X if(ireg==1) {fmax=amax1(fmax,ash(ix,iy,iz))} else {ash(ix,iy,iz) = 0.0} }}} Xif (ireg==1) goto 8642 # visualize regression data X Xdo i = 1,n { X if((i/500)*500.eq.i) {call mess(2,i,0,0) } # looping thru data, on point i X else{ if((i/100)*100.eq.i) {call mess(39,0,0,0)} } X tfact=1.0 # part of product kernel in reference dimensions (4 to nvar) X X do j = 4,nvar { # null loop if nvar=3 X idt = iabs(nc(i,j)-itbin(j)) X if( idt >= mv(j) ){ # this point not covered by hyperrectangle X iskip=iskip+1; goto 1110 } X tfact = tfact * xkern(idt,mv(j),iexp1,iexp2) X } X X ishow = ishow+1 X X do j = 1,3 { X ixyz(j) = nc(i,j) X ixyzl(j) = max0(1,ixyz(j)-(mv(j)-1)) X ixyzr(j) = min0(nbin(j),ixyz(j)+(mv(j)-1)) X } X X#-update ASH array in (x,y,z|t) space X X do iz = ixyzl(3),ixyzr(3) { X idz = iabs(iz-ixyz(3)) X do iy = ixyzl(2),ixyzr(2) { X idy = iabs(iy-ixyz(2)) X do ix = ixyzl(1),ixyzr(1) { X idx = iabs(ix-ixyz(1)) X ash(ix,iy,iz)=ash(ix,iy,iz) + wk(idx,idy,idz)*tfact X } } } X1110 continue # next data point X} X X#-scale by n*hx*hy*hz; find max, etc X8642 continue # regression visualization: skip density estimation above Xdo i = 1,nbin(1) { xfmax(i) = 0 } Xdo i = 1,nbin(2) { yfmax(i) = 0 } Xdo i = 1,nbin(3) { zfmax(i) = 0 } # max of "CELL" in that plane X Xc = n; do i = 1,nvar { c = c * (mv(i)*hv(i)) }; if(ireg==1) {c=1} # leave regression "ash" values alone X Xdo iz = 1,nbin(3) { X do iy = 1,nbin(2) { X do ix = 1,nbin(1) { X a = ash(ix,iy,iz) / c X ash(ix,iy,iz) = a X xfmax(ix) = amax1(a,xfmax(ix)) X yfmax(iy) = amax1(a,yfmax(iy)) X zfmax(iz) = amax1(a,zfmax(iz)) X} } } X X X fmax = 0; do i = 1,nbin(1) { fmax = amax1(fmax,xfmax(i)) } X if( nvar==3) fmaxin=fmax # don't worry about higher dimensions X X475 continue X rout(1)=fmax;rout(2)=fmaxin; call mess(32,0,rout,2) # print fmax,fmaxin X iout(1)=n;iout(2)=ishow;iout(3)=iskip; call mess(33,iout,0,3) # n show skip X X6611 continue X X#-draw triangularization X Xif (iblank.eq.1) { # clear and redraw screen if iblank=1 X#??? call finisz X X if (iprint.ne.0) {call enddump; call endwrite X npict = npict + 1 X if(iprint==1|iprint==3) call startdump(npict) X if(iprint.ge.2) call startwrite(npict) } X X call beginz; call boxz X} X X do i=1,3 { tdel(i) = 2.0/nbin(i) } X c=-1.-tdel(1)/2; do i=1,nbin(1) { tx(i)=c+i*tdel(1) } # make into 3x3 array X c=-1.-tdel(2)/2; do i=1,nbin(2) { ty(i)=c+i*tdel(2) } #####!!!!! X c=-1.-tdel(3)/2; do i=1,nbin(3) { tz(i)=c+i*tdel(3) } X Xntri=0 #### number of triangles found X Xif( fmaxin>0.0 ) c = fmaxin * ff # contour level X else c = fmax * ff Xcall mess(34,0,ff,1) # print ff Xcall mess(35,0, c,1) # print contour level c X Xif (iprint==1 | iprint==3) { X ttxyz(1)=tx(ix1);ttxyz(2)=ty(iy1);ttxyz(3)=tz(iz1) X call dumpit(ttxyz,ttxyz,ttxyz,tnormal,talpha) X ttxyz(1)=tx(ix2+1);ttxyz(2)=ty(iy2+1);ttxyz(3)=tz(iz2+1) X call dumpit(ttxyz,ttxyz,ttxyz,tnormal,talpha)} X X#----------------------------------------------- Xif(nvar>3) call mess(36,itbin(4),0,nvar-3) # print t slice(s) Xcall mess(1,0,0,0) # z slice Xdo iz = iz1,iz2,iz3{ X if(zfmax(iz)<=c & zfmax(iz+1)<=c) goto 603 X call mess(40,iz,0,0) # z slice value is iz X txyz(3)=tz(iz) X do iy = iy1,iy2,iy3{ X if(yfmax(iy)<=c & yfmax(iy+1)<=c) goto 602 X txyz(2)=ty(iy) X do ix = ix1,ix2,ix3{ X if(xfmax(ix)<=c & xfmax(ix+1)<=c) goto 601 X txyz(1)=tx(ix) X do idz=0,1{do idy=0,1{do idx=0,1{ X a=ash(ix+idx,iy+idy,iz+idz); as(idx,idy,idz)=a X if(a==c) {call intpr("opps contour level equality problem===adjusting ff and c", X 56,0,0); ff=1.00001*ff; c=1.00001*c }}}} # Fatal: adjust ff by a bit X X if( iprint.ge.2 ) call docube( txyz, tdel, as, c, talpha) X X ncase=0 X if(as(0,1,1)>c)ncase=ncase+1 # corner 8 X if(as(1,1,1)>c)ncase=ncase+2 # corner 7 X if(as(1,0,1)>c)ncase=ncase+4 # corner 6 X if(as(0,0,1)>c)ncase=ncase+8 # corner 5 X if(as(0,1,0)>c)ncase=ncase+16 # corner 4 X if(as(1,1,0)>c)ncase=ncase+32 # corner 3 X if(as(1,0,0)>c)ncase=ncase+64 # corner 2 X if(as(0,0,0)>c)ncase=ncase+128 # corner 1 Xif(ncase==0 | ncase==255) go to 601 X Xif(ix2-ix1+iy2-iy1+iz2-iz1.eq.0) call realpr("f xyz",5,as,8) X Xdo i=1,12{px(i)=0.;py(i)=0.} # not really necessary if working X call qproject(c,0,0,0,as(0,0,0),as(1,0,0),1,txyz,tdel,1, X v1,v2,v3,px,py,pz,xyz(1,1)) X call qproject(c,1,0,0,as(1,0,0),as(1,1,0),2,txyz,tdel,2, X v1,v2,v3,px,py,pz,xyz(1,2)) X call qproject(c,0,1,0,as(0,1,0),as(1,1,0),3,txyz,tdel,1, X v1,v2,v3,px,py,pz,xyz(1,3)) X call qproject(c,0,0,0,as(0,0,0),as(0,1,0),4,txyz,tdel,2, X v1,v2,v3,px,py,pz,xyz(1,4)) X call qproject(c,0,0,1,as(0,0,1),as(1,0,1),5,txyz,tdel,1, X v1,v2,v3,px,py,pz,xyz(1,5)) X call qproject(c,1,0,1,as(1,0,1),as(1,1,1),6,txyz,tdel,2, X v1,v2,v3,px,py,pz,xyz(1,6)) X call qproject(c,0,1,1,as(0,1,1),as(1,1,1),7,txyz,tdel,1, X v1,v2,v3,px,py,pz,xyz(1,7)) X call qproject(c,0,0,1,as(0,0,1),as(0,1,1),8,txyz,tdel,2, X v1,v2,v3,px,py,pz,xyz(1,8)) X call qproject(c,0,0,0,as(0,0,0),as(0,0,1),9,txyz,tdel,3, X v1,v2,v3,px,py,pz,xyz(1,9)) X call qproject(c,1,0,0,as(1,0,0),as(1,0,1),10,txyz,tdel,3, X v1,v2,v3,px,py,pz,xyz(1,10)) X call qproject(c,1,1,0,as(1,1,0),as(1,1,1),11,txyz,tdel,3, X v1,v2,v3,px,py,pz,xyz(1,11)) X call qproject(c,0,1,0,as(0,1,0),as(0,1,1),12,txyz,tdel,3, X v1,v2,v3,px,py,pz,xyz(1,12)) X X do k=1,contour(ncase,2){ # number of triangles X ktri = cuts( contour(ncase,1)+k-1 ) # triangle number X is1=tri(ktri,1); is2=tri(ktri,2); is3=tri(ktri,3) X X iplot=1 # plot this triangle X X#### always checks now if(orth==1) { # check qoutward normal vs. eye X X do ii=1,3{ee(ii,1)=xyz(ii,is2)-xyz(ii,is1)} # vector side 1 X do ii=1,3{ee(ii,2)=xyz(ii,is3)-xyz(ii,is2)} X do ii=1,3{ee(ii,3)=xyz(ii,is1)-xyz(ii,is3)} X call qoutward(ee,vout) # a normal vector X jj=tri(ktri,4) # a distance corner from triangle X if(jj<=0){call intpr("-----mistake: used triangle",27,0,0); return} X pp=0.; do ii=1,3{pp=pp+vout(ii)*rc(ii,jj)} X if(pp<0.) do ii=1,3{vout(ii)=-vout(ii)} # direction towards corner X if(as(ic(1,jj),ic(2,jj),ic(3,jj))>c) { X do ii=1,3{vout(ii)=-vout(ii)}} # direction now towards out X pp=0.; do ii=1,3{pp=pp+vout(ii)*v3(ii)}; if(pp<=0.)iplot=0 X X## vout2(1) = ee(2,1)*ee(3,2)-ee(3,1)*ee(2,2) X## vout2(2) = -(ee(1,1)*ee(3,2)-ee(3,1)*ee(1,2)) X## vout2(3) = ee(1,1)*ee(2,2)-ee(2,1)*ee(1,2) X##flip=0.; do ii=1,3{flip=flip+vout(ii)*vout2(ii)} X Xflip=0.; do ii=1,3{flip=flip+vout(ii)**2}; flip=sqrt(flip); do ii=1,3{vout2(ii)=vout(ii)/flip} Xif(iprint==1|iprint==3) call dumpit(xyz(1,is1),xyz(1,is2),xyz(1,is3),vout2,talpha) X### call realpr("cosine eye",11,pp,1) X############### } X X if( !(orth==1 & iplot==0)) { # plot this triangle X### call segmtz(px(is1),py(is1),px(is2),py(is2),1) X### call segmtz(px(is2),py(is2),px(is3),py(is3),1) X### call segmtz(px(is3),py(is3),px(is1),py(is1),1) X X ntri = ntri + 1 X### call intpr("------------ntr--------",23,ntri,1) X if(ntri>ntrimax) {ntri=ntri-1000 X call intpr("----ntrimax----deleted 1000 triangles----",41,0,0)} X angle(ntri)=pp X pxy(1,1,ntri)=px(is1); pxy(1,2,ntri)=py(is1) X pxy(2,1,ntri)=px(is2); pxy(2,2,ntri)=py(is2) X pxy(3,1,ntri)=px(is3); pxy(3,2,ntri)=py(is3) X ztri(ntri)=amin1(pz(is1),pz(is2),pz(is3)) X } X Xdo i=1,3{ii=tri(ktri,i);if(px(ii)==0. & py(ii)==0.) call intpr("=======error=======",19,0,0)} X } X X601 continue # next ix X } # end ix loop X602 continue # next iy X } # end iy loop X603 continue # next iz X} # end iz loop X X# ------------------------------------ X Xcall mess(6,ntri,rout,0) # number of triangle patches found Xif(ntri>0) {call mess(7,0,0,0); do i=1,ntri{iperm(i)=i} X call sort (ztri,iperm,1,ntri); call mess(8,0,0,0) } # sorting Xistop = 2500000 # SHOULD BE 0 Xdo ii=1,ntri{ Xif(ii>istop) {call intpr("enter stopping point",20,0,0) X if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer Xcall deci00(istop,.true.) X}} else {call myio(1,istop,rout,1)} } X X i=iperm(ii) Xif(istop<0) i=iperm(-istop) ### for debugging purposes X z=angle(i) X if( ictopt.ne.1 ) { ### not postscript X if(z>=0.) { k=color(1)+z*(color(2)-color(1)+.99)} X else { k=color(3)+abs(z)*(color(4)-color(3)+.99)} # abs(z) 7-28 Xif(k>color(4)) call intpr("================color too large===============",46,0,0) } X if( ictopt.eq.1 ) { ### for postscript 9-25-89 X if(z>=0.) { k=color(1)+z*(color(2)-color(1)+.99)} X else { k=color(3)+abs(z)*(color(4)-color(3)+.99)} X k=ictable(k,1)} X X if(sun50==0) {{ Xam(10) = k+icoloroffset X; }} X else {{ Xam(8) = k X; }} X X call qtriout(pxy(1,1,i),square) X X## call segmtz(pxy(1,1,i),pxy(1,2,i),pxy(2,1,i),pxy(2,2,i),1) X## call segmtz(pxy(2,1,i),pxy(2,2,i),pxy(3,1,i),pxy(3,2,i),1) X## call segmtz(pxy(3,1,i),pxy(3,2,i),pxy(1,1,i),pxy(1,2,i),1) X} X X# print axes X X if(sun50==0) {{ Xam(10) = 1+icoloroffset X; }} X else {{ Xam(8) = 1 X; }} X Xif (iqaxis==1) { X p1x=.5*v1(1); p1y=.5*v2(1) # x-axis X p2x=.5*v1(2); p2y=.5*v2(2) # y-axis X p3x=.5*v1(3); p3y=.5*v2(3) # z-axis X p4x=0.; p4y=0. # origin X Xp1x=square(1) + (square(2)-square(1))*(p1x+s3)/(2*s3) + .999 Xp2x=square(1) + (square(2)-square(1))*(p2x+s3)/(2*s3) + .999 Xp3x=square(1) + (square(2)-square(1))*(p3x+s3)/(2*s3) + .999 Xp4x=square(1) + (square(2)-square(1))*(p4x+s3)/(2*s3) + .999 Xp1y=square(1) + (square(2)-square(1))*(p1y+s3)/(2*s3) + .999 Xp2y=square(1) + (square(2)-square(1))*(p2y+s3)/(2*s3) + .999 Xp3y=square(1) + (square(2)-square(1))*(p3y+s3)/(2*s3) + .999 Xp4y=square(1) + (square(2)-square(1))*(p4y+s3)/(2*s3) + .999 X X call segmtz(p1x,p1y,p4x,p4y,1) X call segmtz(p2x,p2y,p4x,p4y,1) X call segmtz(p3x,p3y,p4x,p4y,1) X call textz(p1x,p1y,"x"//char(0)) X call textz(p2x,p2y,"y"//char(0)) X call textz(p3x,p3y,"z"//char(0)) X} X X# draw the cube !!! Xif (icube==1) { X p1x=0.; do j=1,3{p1x=p1x+v1(j)*cb(1,j)} X p1y=0.; do j=1,3{p1y=p1y+v2(j)*cb(1,j)} X do i = 2,16{ X p2x=0.; do j=1,3{p2x=p2x+v1(j)*cb(i,j)} X p2y=0.; do j=1,3{p2y=p2y+v2(j)*cb(i,j)} X X p1xx=square(1) + (square(2)-square(1))*(p1x+s3)/(2*s3) + .999 X p2xx=square(1) + (square(2)-square(1))*(p2x+s3)/(2*s3) + .999 X p1yy=square(1) + (square(2)-square(1))*(p1y+s3)/(2*s3) + .999 X p2yy=square(1) + (square(2)-square(1))*(p2y+s3)/(2*s3) + .999 X X call segmtz(p1xx,p1yy,p2xx,p2yy,1) X p1x=p2x; p1y=p2y X } X} X Xinew = 0 # can use current ash for many parameter changes X Xgo to 5555 X X9005 return X end X Xsubroutine qproject(c,ix,iy,iz,f1,f2,ns,txyz,tdel,iptr,v1,v2,v3,px,py,pz,xyz) Xdimension v1(3),v2(3),v3(3),txyz(3),tdel(3),px(12),py(12),pz(12),xyz(3),ixyz(3) X X# computes intersection of (f1,f2) at level c X# txyz contains point of intersection X# iptr indicates subscript in middle of a side X# ixyz indicates beginning of line segment (== side) X# txyz contains the (xyz) coordinates of origin X# ns contains number of side X# xyz should be the "ns-th" column of xyz above X Xif( amax1(f1,f2)<=c | amin1(f1,f2)>=c ) return X Xixyz(1)=ix; ixyz(2)=iy;ixyz(3)=iz Xdenom = f2-f1; if(denom==0.) {call intpr("denom========0",14,0,0);return} Xfc = (c-f1)/denom Xif (fc>0. & fc<1.) { # found contour X do i=1,3{xyz(i)=txyz(i)+ixyz(i)*tdel(i)} X xyz(iptr) = txyz(iptr) + fc*tdel(iptr) X px(ns)=v1(1)*xyz(1) + v1(2)*xyz(2) + v1(3)*xyz(3) X py(ns)=v2(1)*xyz(1) + v2(2)*xyz(2) + v2(3)*xyz(3) X pz(ns)=v3(1)*xyz(1) + v3(2)*xyz(2) + v3(3)*xyz(3) X###call intpr(" intersection found on side",31,ns,1) X } Xreturn Xend X Xsubroutine qoutward(e,vout) # compute a normal vector Xdimension e(3,3),vout(3), det(3) X Xdet(1) = e(2,1)*e(3,2) - e(3,1)*e(2,2) Xdet(2) = e(1,1)*e(3,2) - e(3,1)*e(1,2) Xdet(3) = e(1,1)*e(2,2) - e(2,1)*e(1,2) X Xid=1; detm=abs(det(1)) Xdo i=2,3{if(abs(det(i))>detm){id=i; detm=abs(det(i))}} Xif( detm==0.) { call intpr("determinant is zero!!!",22,0,0);return } X X## id = 1 ####---------------TEMP FIX----------##### X Xswitch(id){ Xcase 1: X vout(1)=1. X vout(2) = ( -e(1,1)*e(3,2) + e(1,2)*e(3,1) ) / det(1) X vout(3) = ( -e(2,1)*e(1,2) + e(2,2)*e(1,1) ) / det(1) Xcase 2: X vout(1) = ( -e(2,1)*e(3,2) + e(2,2)*e(3,1) ) / det(2) X vout(2)=1. X vout(3) = ( -e(1,1)*e(2,2) + e(1,2)*e(2,1) ) / det(2) Xcase 3: X vout(1) = ( -e(3,1)*e(2,2) + e(3,2)*e(2,1) ) / det(3) X vout(2) = ( -e(1,1)*e(3,2) + e(1,2)*e(3,1) ) / det(3) X vout(3)=1. X} X Xpp=0.; do i=1,3{pp=pp+vout(i)**2}; pp=sqrt(pp) Xdo i=1,3{vout(i)=vout(i)/pp} X Xreturn Xend X Xsubroutine qtriout(xyorig,square) Xreal xyorig(3,2), xy(3,2), y(3) # y(3) rather than y(2) if an error Xinteger square(2) X X# find vertical lines intersecting triangle X X Xreal am(200) Xcommon/bgrp/am X Xcommon/encbfc/encbf0 Xcommon/encbfl/iqqbfl,iqqfil Xcharacter*256 encbf0 ; integer iqqbfl,iqqfil X Xcommon/decbfc/decbf0 Xcharacter*512 decbf0 Xcommon/decbfl/jqqbfl,jqqinf,jqqflp,jqqfle,jqqlsc,jqdchk Xinteger jqqbfl,jqqinf,jqqflp,jqqlsc; logical jqqfle,jqdchk X X X Xdo j=1,2{ do i=1,3{ X xy(i,j)=square(1) + (xyorig(i,j)+sqrt(3.))*(square(2)-square(1))/(2*sqrt(3.)) + .999}} X X call zpolyz( xy(1,1), xy(1,2), 3 ) # 7-11/89 shaded polygon?? Xreturn X################-------------------replaced last part-----------------############# X X987 continue Xx1=amin1(xy(1,1),xy(2,1),xy(3,1)); x2=amax1(xy(1,1),xy(2,1),xy(3,1)) Xk1 = x1 + .999 Xk2 = x2 X X##############work on this so rounds rather than truncates xy(*,*) X###### k1 and k2 as well X X###print(tstring(x1 x2 y1 y2 k1 k2),R(x1),R(x2),R(y1),R(y2),I(k1),I(k2)) Xdo k=k1,k2 { X iptr=0 X x=k X X#### could add case where denom=0 (ie x1=x2) but x=x1=x2 (a vertical line, get 2 points) X X denom=(xy(2,1)-xy(1,1)) X if(denom!=0.) { X a = (x-xy(1,1))/denom X if(a>=0.&a<=1.) { X yy = (xy(2,2)-xy(1,2))*a + xy(1,2) X iptr=iptr+1;y(iptr)=yy ####;print(I(1),R(x),R(yy)) X } X } X denom=(xy(3,1)-xy(2,1)) X if(denom!=0.) { X a = (x-xy(2,1))/denom X if(a>=0.&a<=1.) { X yy = (xy(3,2)-xy(2,2))*a + xy(2,2) X iptr=iptr+1;y(iptr)=yy ###;print(I(2),R(x),R(yy)) X } X } X denom=(xy(1,1)-xy(3,1)) X if(denom!=0.) { X a = (x-xy(3,1))/denom X if(a>=0.&a<=1.) { X yy = (xy(1,2)-xy(3,2))*a + xy(3,2) X iptr=iptr+1;y(iptr)=yy ###;print(I(3),R(x),R(yy)) X } X } X X if(iptr==2) { call segmtz(x,y(1),x,y(2),1) } X else {call intpr("------error iptr==============",31,0,0)} X} Xreturn Xend X Xreal function xkern(k,mv,iexp1,iexp2) Xxkern = (1.0 - (k/float(mv))**iexp1)**iexp2 Xreturn Xend X END_OF_FILE if test 25958 -ne `wc -c <'tmp2.r'`; then echo shar: \"'tmp2.r'\" unpacked with wrong size! fi # end of 'tmp2.r' fi if test -f 'sort.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sort.f'\" else echo shar: Extracting \"'sort.f'\" \(2050 characters\) sed "s/^X//" >'sort.f' <<'END_OF_FILE' Xc dimension v(100) Xc integer a(100) Xc1 print *,'enter nv plus values' Xc read *,nv,(v(i),i=1,nv) Xc if (nv.eq.0) goto 999 Xc print *,'enter na plus integer pointers' Xc read *,na,(a(i),i=1,na) Xc call sort(v,a,1,na) Xc n=max0(na,nv) Xc print *,' perm=',(a(i),i=1,n) Xc print *,' sort=',(v(i),i=1,n) Xc goto1 Xc999 stop Xc end X X subroutine sort (v,a,ii,jj) X Xc puts into a the permutation vector which sorts v into Xc increasing order. only elements from ii to jj are considered. Xc arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements X Xc the vector v is sorted upon return X Xc this is a modification of cacm algorithm #347 by r. c. singleton, Xc which is a modified hoare quicksort. X X dimension a(jj),v(1),iu(20),il(20) X integer t,tt X integer a X real v X m=1 X i=ii X j=jj X 10 if (i.ge.j) go to 80 X 20 k=i X ij=(j+i)/2 X t=a(ij) X vt=v(ij) X if (v(i).le.vt) go to 30 X a(ij)=a(i) X a(i)=t X t=a(ij) X v(ij)=v(i) X v(i)=vt X vt=v(ij) X 30 l=j X if (v(j).ge.vt) go to 50 X a(ij)=a(j) X a(j)=t X t=a(ij) X v(ij)=v(j) X v(j)=vt X vt=v(ij) X if (v(i).le.vt) go to 50 X a(ij)=a(i) X a(i)=t X t=a(ij) X v(ij)=v(i) X v(i)=vt X vt=v(ij) X go to 50 X 40 a(l)=a(k) X a(k)=tt X v(l)=v(k) X v(k)=vtt X 50 l=l-1 X if (v(l).gt.vt) go to 50 X tt=a(l) X vtt=v(l) X 60 k=k+1 X if (v(k).lt.vt) go to 60 X if (k.le.l) go to 40 X if (l-i.le.j-k) go to 70 X il(m)=i X iu(m)=l X i=k X m=m+1 X go to 90 X 70 il(m)=k X iu(m)=j X j=l X m=m+1 X go to 90 X 80 m=m-1 X if (m.eq.0) return X i=il(m) X j=iu(m) X 90 if (j-i.gt.10) go to 20 X if (i.eq.ii) go to 10 X i=i-1 X 100 i=i+1 X if (i.eq.j) go to 80 X t=a(i+1) X vt=v(i+1) X if (v(i).le.vt) go to 100 X k=i X 110 a(k+1)=a(k) X v(k+1)=v(k) X k=k-1 X if (vt.lt.v(k)) go to 110 X a(k+1)=t X v(k+1)=vt X go to 100 X end X END_OF_FILE if test 2050 -ne `wc -c <'sort.f'`; then echo shar: \"'sort.f'\" unpacked with wrong size! fi # end of 'sort.f' fi if test -f 'ash.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ash.f'\" else echo shar: Extracting \"'ash.f'\" \(7135 characters\) sed "s/^X//" >'ash.f' <<'END_OF_FILE' Xc ash.f---------------------------------------------------------------------- Xc April 8, 1986 Xc Xc Find bin counts of data array "x(n)" for ASH estimator Xc Xc "nbin" bins are formed over the interval [a,b) Xc Xc bin counts returned in array "nc" - # pts outside [a,b) = "nskip" X Xc ##### Copyright 1986 David W. Scott X X subroutine bin1 ( x , n , ab , nbin , nc , nskip ) X dimension x(n) , nc(nbin) , ab(2) X X nskip = 0 X a = ab(1) X b = ab(2) X X do 5 i = 1,nbin X nc(i) = 0 X5 continue X X d = (b-a) / nbin X X do 10 i = 1,n X k = (x(i)-a) / d + 1.0 X if (k.ge.1 .and. k.le.nbin) then X nc(k) = nc(k) + 1 X else X nskip = nskip + 1 X end if X10 continue X X return X end X X Xc April 8, 1986 Xc Xc Computer ASH density estimate; Quartic (biweight) kernel Xc Xc Average of "m" shifted histograms Xc Xc Bin counts in array "nc(nbin)" - from routine "bin1" Xc Xc "nbin" bins are formed over the interval [a,b) Xc Xc ASH estimates returned in array "f(nbin)" Xc Xc FP-ASH plotted at a+d/2 ... b-d/2 where d = (b-a)/nbin Xc Xc Note: If "nskip" was nonzero, ASH estimates incorrect near boundary Xc Note: Should leave "m" empty bins on each end of array "nc" so f OK X Xc ##### Copyright 1986 David W. Scott X X subroutine ash1 ( m, nc, nbin, ab, kopt, t, f, w, ier ) X dimension nc(nbin), ab(2), t(nbin), f(nbin), w(m), kopt(2) X X ier = 0 X a = ab(1) X b = ab(2) X n = 0 X Xc-compute weights cons * ( 1-abs((i/m))^kopt1)^kopt2 Xc -- should sum to "m" 5-8-91 Xc w-array shifted by 1 X X mm1 = m-1 X xm = m Xc cons = sum of weights from -(m-1) to (m-1) = 1 + 2 (sum from 1 to m-1) X X w(1) = 1.0 X cons = 1.0 X do 5 i = 1,mm1 X w(i+1) = ( 1.0 - abs(i/xm)**kopt(1) )**kopt(2) X cons = cons + 2*w(i+1) X5 continue X cons = float(m)/cons X do 6 i=1,m X w(i) = cons*w(i) X6 continue X Xc-check if estimate extends beyond mesh X X do 7 i = 1,mm1 X if( nc(i)+nc(nbin+1-i) .gt. 0) ier = 1 X7 continue X Xc-compute ash(m) estimate X X delta = (b-a) / nbin X h = m*delta X X do 10 i = 1,nbin X t(i) = a + (i-0.5)*delta X f(i) = 0.0 X n = n + nc(i) X10 continue X X do 20 i = 1,nbin X if (nc(i).eq.0) goto 20 X X c = nc(i) / (n*h) X do 15 k = max0(1,i-mm1) , min0(nbin,i+mm1) X f(k) = f(k) + c * w( iabs(k-i)+1 ) X15 continue X20 continue X X return X end X X Xc April 12, 1986 bin2.f Xc Xc Find bin counts of data array "x(n,2)" for ASH estimator Xc Xc "nbin(1)" by "nbin(2)" bins are formed Xc Xc x:axis [ ab(1,1) , ab(1,2) ) Xc y:axis [ ab(2,1) , ab(2,2) ) half-open Xc Xc bin counts returned in array "nc" - # pts outside [a,b) = "nskip" X Xc ##### Copyright 1986 David W. Scott X X subroutine bin2 ( x , n , ab , nbin , nc , nskip ) X dimension x(n,2) , nbin(2) , nc(nbin(1),nbin(2)) , ab(2,2) X X nskip = 0 X ax = ab(1,1) X bx = ab(1,2) X ay = ab(2,1) X by = ab(2,2) X X do 5 j = 1,nbin(2) X do 4 i = 1,nbin(1) X nc(i,j) = 0 X4 continue X5 continue X X dx = (bx-ax) / nbin(1) X dy = (by-ay) / nbin(2) X X do 10 i = 1,n X kx = (x(i,1)-ax) / dx + 1.0 X ky = (x(i,2)-ay) / dy + 1.0 X if (kx.ge.1 .and. kx.le.nbin(1) .and. X 1 ky.ge.1 .and. ky.le.nbin(2)) then X nc(kx,ky) = nc(kx,ky) + 1 X else X nskip = nskip + 1 X end if X10 continue X X return X end X X Xc April 12, 1986 ash2.f Xc Xc Computer ASH density estimate; Product Quartic (biweight) kernel Xc Xc Average of "m[1] by m[2]" shifted histograms Xc Xc Bin counts in matrix "nc" - from routine "nbin2" Xc Xc ASH estimates returned in matrix "f" Xc Xc FP-ASH plotted at a+d/2 ... b-d/2 where d = (b-a)/nbin Xc Xc Note: If "nskip" was nonzero, ASH estimates incorrect near boundary Xc Note: Should leave "m" empty bins on each end of array "nc" so f OK X Xc #### Copyright 1986 David W. Scott Xc #### added kernel option kopt 5-8-91 X X subroutine ash2 ( m, nc, nbinx, nbiny, ab, kopt, f, w, ier ) X dimension m(2) , nc(nbinx,nbiny) , ab(2,2) , f(nbinx,nbiny) X dimension w( m(1) , m(2) ), kopt(2) X X ier = 0 X ax = ab(1,1) X bx = ab(1,2) X ay = ab(2,1) X by = ab(2,2) X Xc-compute weights cons * ( 1-abs(i/m)^kopt1)^kopt2 Xc -- should sum to "m" 5-8-91 Xc w-array shifted by 1 X X mx = m(1) X my = m(2) X mxm1 = mx-1 X mym1 = my-1 X xm = mx X ym = my X Xc-----ASSUMES f dimensioned larger than m(1) or m(2) X Xc-put marginal weights in f array as work array X X f(1,1) = 1.0 X f(2,1) = 1.0 Xc consx = sum of weights from -(m-1) to (m-1) = 1 + 2 (sum from 1 to m-1) X consx = 1.0 X consy = 1.0 X X do 5 i = 1,mxm1 X f(1,i+1) = ( 1.0 - abs(i/xm)**kopt(1) )**kopt(2) X consx = consx + 2*f(1,i+1) X5 continue X consx = float(mx)/consx X do 6 i = 1,mym1 X f(2,i+1) = ( 1.0 - abs(i/ym)**kopt(1) )**kopt(2) X consy = consy + 2*f(2,i+1) X6 continue X consy = float(my)/consy X Xc-computer product weight array (avoids later multiplications) X X do 3 j = 1,my X do 2 i = 1,mx X w(i,j) = (consx*f(1,i)) * (consy*f(2,j)) X2 continue X3 continue X Xc-compute ash(m) estimate X X n = 0 X do 10 j = 1,nbiny X do 9 i = 1,nbinx X f(i,j) = 0.0 X n = n + nc(i,j) X9 continue X10 continue Xc-check if estimate extends beyond mesh X X ncheck = 0 X do 12 j = my , nbiny+1-my X do 11 i = mx , nbinx+1-mx X ncheck = ncheck + nc(i,j) X11 continue X12 continue X if (ncheck .ne. n) ier = 1 X X dx = (bx-ax) / nbinx X dy = (by-ay) / nbiny X X hx = mx*dx X hy = my*dy X X do 20 j = 1,nbiny X do 19 i = 1,nbinx X X if (nc(i,j).eq.0) goto 19 X X c = nc(i,j) / (n*hx*hy) X X do 18 ky = max0(1,j-mym1) , min0(nbiny,j+mym1) X do 17 kx = max0(1,i-mxm1) , min0(nbinx,i+mxm1) X f(kx,ky) = f(kx,ky) + c * X 1 w( iabs(kx-i)+1 , iabs(ky-j)+1 ) X17 continue X18 continue X X19 continue X20 continue X X return X end X X END_OF_FILE if test 7135 -ne `wc -c <'ash.f'`; then echo shar: \"'ash.f'\" unpacked with wrong size! fi # end of 'ash.f' fi if test -f 'ct' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'ct'\" else echo shar: Extracting \"'ct'\" \(1749 characters\) sed "s/^X//" >'ct' <<'END_OF_FILE' Xstructure(.Dim = c(101, 3), .Data = c(128, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 155, 157.0408, 159.0816, X 161.1225, 163.1633, 165.2041, 167.2449, 169.2857, 171.3265, 173.3673, X 175.4082, 177.449, 179.4898, 181.5306, 183.5714, 185.6122, 187.6531, X 189.6939, 191.7347, 193.7755, 195.8163, 197.8571, 199.898, 201.9388, X 203.9796, 206.0204, 208.0612, 210.102, 212.1429, 214.1837, 216.2245, X 218.2653, 220.3061, 222.3469, 224.3878, 226.4286, 228.4694, 230.5102, X 232.551, 234.5918, 236.6327, 238.6735, 240.7143, 242.7551, 244.7959, X 246.8367, 248.8775, 250.9184, 252.9592, 255, 128, 0, 5.204082, 10.40816, X 15.61224, 20.81633, 26.02041, 31.22449, 36.42857, 41.63265, 46.83673, X 52.04082, 57.2449, 62.44898, 67.65306, 72.85714, 78.06123, 83.2653, X 88.46938, 93.67347, 98.87755, 104.0816, 109.2857, 114.4898, 119.6939, X 124.898, 130.102, 135.3061, 140.5102, 145.7143, 150.9184, 156.1225, X 161.3265, 166.5306, 171.7347, 176.9388, 182.1429, 187.3469, 192.551, X 197.7551, 202.9592, 208.1633, 213.3673, 218.5714, 223.7755, 228.9796, X 234.1837, 239.3878, 244.5918, 249.7959, 255, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 128, 255, X 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, X 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, X 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, X 255, 255, 255, 255, 255, 255, 255, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, X 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) X X END_OF_FILE if test 1749 -ne `wc -c <'ct'`; then echo shar: \"'ct'\" unpacked with wrong size! fi # end of 'ct' fi if test -f 'mymaster' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'mymaster'\" else echo shar: Extracting \"'mymaster'\" \(7205 characters\) sed "s/^X//" >'mymaster' <<'END_OF_FILE' Xlist(contour = structure(.Dim = c(254, 2), .Data = c(8, 7, 63, 6, 61, 59, 230, X 5, 57, 55, 227, 53, 224, 221, 499, 4, 51, 49, 218, 47, 215, 212, 495, X 45, 209, 206, 491, 203, 487, 483, 65, 3, 43, 41, 200, 39, 197, 194, X 479, 37, 191, 188, 475, 185, 471, 467, 68, 35, 182, 179, 465, 176, X 461, 457, 71, 173, 453, 449, 74, 445, 77, 80, 9, 2, 33, 31, 170, 29, X 167, 164, 441, 27, 161, 158, 437, 155, 433, 429, 83, 25, 152, 149, X 425, 146, 421, 417, 86, 143, 413, 409, 89, 405, 92, 95, 11, 23, 140, X 137, 401, 134, 397, 395, 98, 131, 391, 387, 101, 383, 104, 107, 13, X 128, 379, 375, 110, 371, 113, 116, 15, 367, 119, 122, 17, 125, 19, X 21, 1, 1, 21, 19, 125, 17, 122, 119, 363, 15, 116, 113, 359, 110, 355, X 351, 128, 13, 107, 104, 347, 101, 343, 339, 131, 98, 337, 333, 134, X 329, 137, 140, 23, 11, 95, 92, 325, 89, 321, 317, 143, 86, 313, 309, X 146, 305, 149, 152, 25, 83, 301, 297, 155, 293, 158, 161, 27, 289, X 164, 167, 29, 170, 31, 33, 2, 9, 80, 77, 285, 74, 281, 277, 173, 71, X 273, 269, 176, 267, 179, 182, 35, 68, 263, 259, 185, 255, 188, 191, X 37, 251, 194, 197, 39, 200, 41, 43, 3, 65, 247, 243, 203, 239, 206, X 209, 45, 235, 212, 215, 47, 218, 49, 51, 4, 233, 221, 224, 53, 227, X 55, 57, 5, 230, 59, 61, 6, 63, 7, 8, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, X 3, 2, 3, 3, 2, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 1, 2, X 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 2, 3, 4, 4, 3, 3, X 4, 4, 3, 4, 3, 3, 2, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, X 2, 3, 3, 4, 3, 4, 4, 3, 3, 4, 4, 3, 4, 3, 3, 2, 2, 3, 3, 4, 3, 4, 2, X 3, 3, 4, 4, 3, 4, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, X 2, 1, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 4, 3, X 4, 4, 3, 3, 2, 4, 3, 4, 3, 3, 2, 2, 3, 3, 4, 3, 4, 4, 3, 3, 4, 4, 3, X 4, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 2, 3, 3, X 4, 3, 4, 4, 3, 3, 4, 4, 3, 2, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, X 3, 2, 3, 2, 2, 1, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 2, X 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1)), cuts = c(24, 8, 63, 108, X 176, 168, 189, 204, 70, 159, 24, 63, 16, 133, 30, 23, 24, 168, 24, X 189, 24, 204, 18, 53, 8, 108, 8, 176, 73, 3, 8, 189, 8, 204, 72, 100, X 63, 176, 63, 168, 116, 58, 63, 204, 108, 176, 108, 168, 108, 189, 150, X 103, 195, 191, 176, 189, 175, 182, 208, 173, 168, 204, 216, 193, 159, X 106, 134, 95, 61, 133, 139, 141, 70, 69, 76, 73, 70, 159, 189, 70, X 159, 204, 52, 10, 100, 30, 23, 63, 24, 63, 168, 116, 58, 24, 24, 63, X 204, 19, 49, 30, 16, 133, 168, 16, 133, 189, 122, 127, 16, 192, 47, X 23, 30, 23, 189, 20, 143, 175, 208, 173, 24, 24, 168, 204, 216, 193, X 24, 164, 26, 53, 18, 53, 176, 28, 39, 18, 17, 119, 116, 18, 53, 204, X 8, 108, 176, 73, 3, 108, 8, 108, 189, 150, 103, 8, 4, 37, 191, 8, 176, X 189, 175, 182, 8, 174, 78, 3, 73, 3, 204, 216, 193, 8, 72, 100, 176, X 72, 100, 168, 80, 90, 72, 71, 153, 150, 195, 191, 63, 116, 58, 176, X 175, 182, 63, 59, 88, 173, 63, 168, 204, 194, 121, 58, 195, 191, 108, X 108, 176, 189, 172, 151, 103, 208, 173, 108, 150, 103, 168, 104, 129, X 193, 211, 202, 208, 207, 187, 195, 200, 170, 182, 183, 178, 216, 218, X 220, 104, 129, 214, 177, 120, 109, 101, 140, 103, 151, 205, 188, 159, X 106, 134, 204, 77, 75, 60, 130, 58, 121, 197, 167, 95, 61, 133, 189, X 59, 88, 205, 201, 66, 145, 139, 141, 70, 189, 173, 88, 67, 154, 71, X 153, 138, 140, 69, 76, 73, 204, 70, 159, 216, 193, 3, 78, 185, 179, X 52, 10, 100, 168, 4, 37, 197, 190, 46, 41, 5, 89, 192, 47, 23, 63, X 30, 23, 116, 58, 20, 143, 175, 63, 59, 88, 173, 24, 24, 63, 168, 204, X 194, 121, 58, 24, 17, 119, 117, 130, 19, 49, 30, 189, 12, 110, 16, X 133, 208, 173, 122, 127, 16, 168, 191, 37, 13, 120, 53, 26, 157, 203, X 188, 209, 54, 27, 28, 39, 55, 27, 183, 178, 216, 24, 164, 26, 53, 176, X 20, 143, 185, 169, 154, 144, 21, 38, 23, 47, 214, 203, 16, 127, 117, X 120, 17, 119, 116, 176, 18, 53, 175, 182, 110, 12, 28, 39, 18, 204, X 19, 49, 36, 38, 4, 37, 191, 108, 8, 108, 176, 189, 172, 151, 103, 8, X 174, 78, 3, 108, 73, 3, 150, 103, 104, 129, 193, 8, 203, 212, 96, 7, X 52, 10, 85, 190, 200, 170, 182, 8, 100, 10, 34, 179, 72, 100, 195, X 191, 80, 90, 72, 176, 174, 78, 65, 140, 70, 141, 138, 154, 71, 153, X 150, 168, 145, 66, 95, 61, 123, 201, 207, 187, 195, 63, 133, 61, 82, X 167, 77, 57, 111, 130, 211, 202, 208, 108, 159, 106, 124, 188, 120, X 102, 137, 140, 134, 106, 156, 177, 218, 220), tri = structure(.Dim = c( X 220, 4), .Data = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, X 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, X 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, X 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, X 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, X 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, X 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, X 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, X 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, X 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 10, 2, 2, X 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, X 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, X 8, 9, 9, 9, 10, 10, 11, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, X 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, X 8, 9, 9, 9, 10, 10, 11, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, X 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 5, X 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, X 9, 10, 10, 11, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, X 10, 10, 11, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 8, 8, 8, X 8, 9, 9, 9, 10, 10, 11, 9, 9, 9, 10, 10, 11, 10, 10, 11, 11, 3, 4, X 5, 6, 7, 8, 9, 10, 11, 12, 4, 5, 6, 7, 8, 9, 10, 11, 12, 5, 6, 7, 8, X 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, X 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 4, 5, 6, 7, 8, 9, 10, X 11, 12, 5, 6, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, X 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, X 5, 6, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, X 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 6, 7, X 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, X 12, 10, 11, 12, 11, 12, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, X 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 8, 9, 10, 11, 12, 9, 10, 11, X 12, 10, 11, 12, 11, 12, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, X 10, 11, 12, 11, 12, 12, 11, 12, 12, 12, 1, 1, 6, 6, 2, 2, 3, 2, 1, X 4, 1, 4, 2, 2, 4, 4, 3, 2, 1, 5, 1, 1, 5, 1, 4, 3, 2, 2, 2, 1, 1, 1, X 6, 5, 1, 5, 6, 1, 2, 2, 2, 1, 2, 1, 2, 2, 5, 2, 1, 1, 1, 4, 3, 2, 3, X 1, 3, 7, 7, 3, 1, 4, 3, 2, 1, 1, 3, 3, 2, 1, 4, 3, 2, 6, 3, 2, 3, 6, X 4, 3, 3, 6, 2, 2, 7, 2, 3, 7, 2, 3, 3, 2, 3, 2, 1, 3, 4, 2, 2, 4, 4, X 4, 8, 8, 3, 2, 1, 4, 4, 4, 3, 3, 4, 3, 4, 3, 7, 1, 3, 4, 7, 4, 8, 7, X 3, 3, 4, 4, 8, 3, 1, 3, 1, 2, 4, 3, 4, 8, 1, 4, 1, 1, 5, 1, 1, 1, 4, X 1, 4, 4, 8, 2, 4, 1, 1, 5, 8, 1, 2, 1, 1, 2, 4, 3, 5, 5, 1, 6, 3, 8, X 5, 8, 7, 6, 5, 5, 2, 7, 4, 1, 3, 8, 7, 4, 7, 5, 5, 2, 7, 4, 6, 5, 8, X 7, 5, 1, 8, 2, 4, 8, 1, 6, 3, 8, 5, 1, 5, 6, 2, 3, 6, 3, 1, 6, 2, 7, X 1, 2, 1, 2))) X END_OF_FILE if test 7205 -ne `wc -c <'mymaster'`; then echo shar: \"'mymaster'\" unpacked with wrong size! fi # end of 'mymaster' fi if test -f 'stream.s.norm.1' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'stream.s.norm.1'\" else echo shar: Extracting \"'stream.s.norm.1'\" \(263 characters\) sed "s/^X//" >'stream.s.norm.1' <<'END_OF_FILE' X16 1 X2 4 4 4 X6 2 5 1 X11 1 X14 0 0 1 1 X3 .100 0 X8 2 1 1 0 X8 2 2 1 0 X8 2 2 2 0 X8 1 1 1 X 3 .1 14 0 0 1 .3 0 13 0 X 3 .4 14 1 0 0 .5 0 X 3 .8 14 0 1 0 .7 0 X99 X Xlevel f=.1 full contour/then with stripes/escher X X3 solid nested contours/transparent X END_OF_FILE if test 263 -ne `wc -c <'stream.s.norm.1'`; then echo shar: \"'stream.s.norm.1'\" unpacked with wrong size! fi # end of 'stream.s.norm.1' fi if test -f 'norm.1' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'norm.1'\" else echo shar: Extracting \"'norm.1'\" \(333 characters\) sed "s/^X//" >'norm.1' <<'END_OF_FILE' X#!/bin/sh X\rm -f stream.s Xcp stream.s.norm.1 stream.s XSplus <