August 6, 1992 January 21, 1993 To: Anonymous From: David W. Scott Re: ashn code Reference: Scott, David W. (1992). "Multivariate Density Estimation: Theory, Practice, and Visualization." John Wiley, New York. Here are the steps as I see them. 1. Create a safe directory somewhere, and cut and paste the files as named in the =========== header ========== % mkdir .Data [local Splus data directory] 2. chmod +x make.ashn then run it to make a local version of Splusn 3. must run S in the directory that contains the "local.Sqpe" file 4. run "S" then source("ashn.s"); source("compress.s") note: compress doesn't really do much compression these days, but it is a preprocessor that replaces the raw data by bin locations and outputs a list suitable for running ashn.... 5. some S files required can be read by a. mymaster <- dget("mymaster") b. ct <- dget("ct") c. X11.ash2 <- dget("X11.ash2") d. xyz <- compress( matrix( rnorm(3000),1000,3 ) ) # an example e. storage.mode(ct ) <- "single" # dput/dget not perfect storage.mode(mymaster$contour) <- "integer" storage.mode(mymaster$cuts ) <- "integer" storage.mode(mymaster$tri ) <- "integer" 6. I assume you are running X-windows on a color Sun. In the console, run xrdb .Xresources to get the color information used in X11.ash2 to work.... fire up an X window inside Splus by typing X11(X11.ash2) 7. Actually, I have an automatic mode that you can use to see if everything is working. Exit Splus if you are in it. % chmod +x norm.1 % norm.1 a. this fires Splus up (the local version) b. starts the X11 window c. generates 5 views of the 999 trivariate normal points in "xyz" using the "action option" choices in "stream.s" 8. The automatic version just zips along, primarily to produce files with the *.quad extension. These can be feed into the MinneView program, which I'll get to you when ftp is up. I only have a binary version. 9. Otherwise, here is a typical session that you might try. > x <- matrix( rnorm(3000),1000,3 ) # trivariate normal data > xyz <- compress(x) # bins the data > X11(X11.ash2) # a nice color window > pie( rep(1,100), col=1:100 ) # look at the nice colors > ashn(xyz) # it does some automatic things # you can ignore for now 6 2 5 1 # this changes the "eye" location 2 4 4 4 # this changes the smoothing parameters 3 .5 # select the shell at alpha-level # 0.50 * f(mode) 0 < alpha < 1 0 # plot it 3 .2 # plot the alpha = 0.20 shell 13 0 # don't erase window for next plot 3 .5 0 # adds to the plot [this will work # great on the SGI if the # shells have been made transparent 13 1 # back to erasing screeen 3 .1 8 2 1 1 # leave some of the shell off 11 1 # turn on *.quad output 14 1 0 0 .5 # RBG color + transparency fraction # for shell [all between 0 and 1] 99 # quit > q() Action options: 0 = plot with current parameter settings [recomputes ash if necessary] 1 = t slice; value of unseed variable f(x,y,z,t); an integer between 1 and nbin[4] -- bin center; if 5-D f(x,y,z,t1,t2), then specify t1 and t2 as integers 2 = m vector of smoothing parameters (length nvar) 3 = f alpha-shell level; fraction between 0 and 1 4 = axes show/don't show axes on plot 5 = cube show/don't show box on plot 6 = eye vector of length 3 giving eye coordinates [ the program views the world as living in the cube [-1,1]^3 ] 7 = slicelimit limits region where shells are computed to hyperrectangle defined by these 6 integers 8 = sliceincrement normally 1; >1 skips bins 9 = iris can specific coloroffset=999 on an SGI Iris 10 = orth only plot triangles oriented towards viewer 11 = print 1 = produce *.quad files 12 = coltbl lookup from old S days 13 = blank screen yes/no erase screen before plot 14 = alpha RGBA or {red,green,blue,alpha} all between 0 and 1 this alpha refers to the transparency of the shell 15 = fmaxin if more than 3 variables, allow the use to specify the value of f(x,y,z,t) at the global mode 16 = kernel 0=ash 1,2,3,... use kernels of the form (1-x^2)^kernel 99 = exit ashn Note that ashn attempts to plot an ash shell based on the calling sequence in the S function ashn. By default, f=4 (which is obviously greater than 1), so that no shell will be found unless "ff" is specified (note that "m" is called "mv" in "ashn"): ashn(xyz) is equivalent to ashn(xyz, ff=4, mv=c(1,1,1)) which does nothing very interesting until you using the action options ashn(xyz, ff=.5, mv=c(4,4,4)) would at least do something initially, and would not have to recompute the ash if you changed the "eye", for example 10. This version of ashn actually will look at shells of slices of the first 3 variables of a multivariate ASH of up to 10 variables. I really haven't tried it beyond 5, but the loops are there. Also there is a hidden regression surface option I can document if you are interested. 11. Thanks for your interest and patience. Let me know what else would be helpful. ----------------------FILES BELOW------------------- ==================== norm.1 ======================= 8/6/92 #!/bin/sh \rm -f stream.s cp stream.s.norm.1 stream.s Splus<= 8 /* COLOR MONITORS WITH MANY COLOR PLANES */ /* The following settings are for color monitors with large colormaps (64 or more entries). The image scheme is identical to that specified by the S-PLUS data set "xcm.heatA". This colormap uses the first 7 colors for regular plotting, and rest of the color specification create a ramp of smoothly blending colors The splus*firstImageColor tells the image command to start with the 7th color when creating images. The first color listed is used as the background color, and isn't counted in the first image color. */ splus*colors : black yellow cyan magenta green MediumBlue red \ black 11 blue 12 magenta 12 red 6 yellow 6 white splus*firstImageColor : 7 splus*nHalftones : 0 splus*halftonePolygon : No splus*halftoneImage : Yes /* Set the colors used in the S-PLUS X11 control panel */ splus*background : blue splus*Text*background : #a0f splus*Text*foreground : yellow splus*Label*background : cyan splus*Label*foreground : red splus*Command*background : red splus*Command*foreground : cyan #else /* COLOR MONITORS WITH ONLY A FEW COLOR PLANES */ /* Don't use as many colors for monitors with fewer color planes. Use a moderate number of halftones to generate enough shade for images */ splus*colors : black yellow cyan magenta green MediumBlue red black yellow splus*nHalftones : 16 splus*firstImageColor : 103 splus*halftonePolygon : No splus*halftoneImage : Yes #endif #else /* BLACK AND WHITE MONITORS */ splus*colors : black white splus*firstImageColor : 0 splus*nHalftones : 16 splus*halftonePolygon : Yes splus*halftoneImage : Yes #endif /* End of COLOR RESOURCES */ /* OTHER RESOURCES */ /* These resources are set by your system administrator, you should not change uncomment these resources unless you know what you are doing. */ /* The printMethod resource determines the starting printing method in the window A value of 0 sets up the graphics menu with the PostScript method, 1 start the graphics window menu with HP LaserJet as the print method. This overrides the value that your system administrator set. */ /* splus*printMethod*state : 0 */ /* Since the printing commands should get their values from the remote machine, you should not declare the printing commands in the resources database with the xrdb command */ /* splus*postScriptPrintCommand*value : lpr -h -r */ /* splus*laserJetPrintCommand*value : lpr -Php -h -r */ /* The next resource sets the initial resolution of the HP LaserJet plot in Dots Per Inch. The values 0,1,2,3 correspond to 75, 100, 150, or 300 dots per inch. This overrides the value your system administrator set. */ /* splus*printDPI*state : 2 */ ==================== mymaster ========contouring structure=== 8/6/92 list(contour = structure(.Dim = c(254, 2), .Data = c(8, 7, 63, 6, 61, 59, 230, 5, 57, 55, 227, 53, 224, 221, 499, 4, 51, 49, 218, 47, 215, 212, 495, 45, 209, 206, 491, 203, 487, 483, 65, 3, 43, 41, 200, 39, 197, 194, 479, 37, 191, 188, 475, 185, 471, 467, 68, 35, 182, 179, 465, 176, 461, 457, 71, 173, 453, 449, 74, 445, 77, 80, 9, 2, 33, 31, 170, 29, 167, 164, 441, 27, 161, 158, 437, 155, 433, 429, 83, 25, 152, 149, 425, 146, 421, 417, 86, 143, 413, 409, 89, 405, 92, 95, 11, 23, 140, 137, 401, 134, 397, 395, 98, 131, 391, 387, 101, 383, 104, 107, 13, 128, 379, 375, 110, 371, 113, 116, 15, 367, 119, 122, 17, 125, 19, 21, 1, 1, 21, 19, 125, 17, 122, 119, 363, 15, 116, 113, 359, 110, 355, 351, 128, 13, 107, 104, 347, 101, 343, 339, 131, 98, 337, 333, 134, 329, 137, 140, 23, 11, 95, 92, 325, 89, 321, 317, 143, 86, 313, 309, 146, 305, 149, 152, 25, 83, 301, 297, 155, 293, 158, 161, 27, 289, 164, 167, 29, 170, 31, 33, 2, 9, 80, 77, 285, 74, 281, 277, 173, 71, 273, 269, 176, 267, 179, 182, 35, 68, 263, 259, 185, 255, 188, 191, 37, 251, 194, 197, 39, 200, 41, 43, 3, 65, 247, 243, 203, 239, 206, 209, 45, 235, 212, 215, 47, 218, 49, 51, 4, 233, 221, 224, 53, 227, 55, 57, 5, 230, 59, 61, 6, 63, 7, 8, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 2, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 2, 3, 4, 4, 3, 3, 4, 4, 3, 4, 3, 3, 2, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 4, 3, 4, 4, 3, 3, 4, 4, 3, 4, 3, 3, 2, 2, 3, 3, 4, 3, 4, 2, 3, 3, 4, 4, 3, 4, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 4, 3, 4, 4, 3, 3, 2, 4, 3, 4, 3, 3, 2, 2, 3, 3, 4, 3, 4, 4, 3, 3, 4, 4, 3, 4, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 2, 3, 3, 4, 3, 4, 4, 3, 3, 4, 4, 3, 2, 3, 3, 2, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 3, 4, 4, 3, 4, 3, 3, 2, 4, 3, 3, 2, 3, 2, 2, 1, 2, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1)), cuts = c(24, 8, 63, 108, 176, 168, 189, 204, 70, 159, 24, 63, 16, 133, 30, 23, 24, 168, 24, 189, 24, 204, 18, 53, 8, 108, 8, 176, 73, 3, 8, 189, 8, 204, 72, 100, 63, 176, 63, 168, 116, 58, 63, 204, 108, 176, 108, 168, 108, 189, 150, 103, 195, 191, 176, 189, 175, 182, 208, 173, 168, 204, 216, 193, 159, 106, 134, 95, 61, 133, 139, 141, 70, 69, 76, 73, 70, 159, 189, 70, 159, 204, 52, 10, 100, 30, 23, 63, 24, 63, 168, 116, 58, 24, 24, 63, 204, 19, 49, 30, 16, 133, 168, 16, 133, 189, 122, 127, 16, 192, 47, 23, 30, 23, 189, 20, 143, 175, 208, 173, 24, 24, 168, 204, 216, 193, 24, 164, 26, 53, 18, 53, 176, 28, 39, 18, 17, 119, 116, 18, 53, 204, 8, 108, 176, 73, 3, 108, 8, 108, 189, 150, 103, 8, 4, 37, 191, 8, 176, 189, 175, 182, 8, 174, 78, 3, 73, 3, 204, 216, 193, 8, 72, 100, 176, 72, 100, 168, 80, 90, 72, 71, 153, 150, 195, 191, 63, 116, 58, 176, 175, 182, 63, 59, 88, 173, 63, 168, 204, 194, 121, 58, 195, 191, 108, 108, 176, 189, 172, 151, 103, 208, 173, 108, 150, 103, 168, 104, 129, 193, 211, 202, 208, 207, 187, 195, 200, 170, 182, 183, 178, 216, 218, 220, 104, 129, 214, 177, 120, 109, 101, 140, 103, 151, 205, 188, 159, 106, 134, 204, 77, 75, 60, 130, 58, 121, 197, 167, 95, 61, 133, 189, 59, 88, 205, 201, 66, 145, 139, 141, 70, 189, 173, 88, 67, 154, 71, 153, 138, 140, 69, 76, 73, 204, 70, 159, 216, 193, 3, 78, 185, 179, 52, 10, 100, 168, 4, 37, 197, 190, 46, 41, 5, 89, 192, 47, 23, 63, 30, 23, 116, 58, 20, 143, 175, 63, 59, 88, 173, 24, 24, 63, 168, 204, 194, 121, 58, 24, 17, 119, 117, 130, 19, 49, 30, 189, 12, 110, 16, 133, 208, 173, 122, 127, 16, 168, 191, 37, 13, 120, 53, 26, 157, 203, 188, 209, 54, 27, 28, 39, 55, 27, 183, 178, 216, 24, 164, 26, 53, 176, 20, 143, 185, 169, 154, 144, 21, 38, 23, 47, 214, 203, 16, 127, 117, 120, 17, 119, 116, 176, 18, 53, 175, 182, 110, 12, 28, 39, 18, 204, 19, 49, 36, 38, 4, 37, 191, 108, 8, 108, 176, 189, 172, 151, 103, 8, 174, 78, 3, 108, 73, 3, 150, 103, 104, 129, 193, 8, 203, 212, 96, 7, 52, 10, 85, 190, 200, 170, 182, 8, 100, 10, 34, 179, 72, 100, 195, 191, 80, 90, 72, 176, 174, 78, 65, 140, 70, 141, 138, 154, 71, 153, 150, 168, 145, 66, 95, 61, 123, 201, 207, 187, 195, 63, 133, 61, 82, 167, 77, 57, 111, 130, 211, 202, 208, 108, 159, 106, 124, 188, 120, 102, 137, 140, 134, 106, 156, 177, 218, 220), tri = structure(.Dim = c( 220, 4), .Data = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 10, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 8, 8, 8, 8, 9, 9, 9, 10, 10, 11, 9, 9, 9, 10, 10, 11, 10, 10, 11, 11, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 4, 5, 6, 7, 8, 9, 10, 11, 12, 5, 6, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 4, 5, 6, 7, 8, 9, 10, 11, 12, 5, 6, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 5, 6, 7, 8, 9, 10, 11, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 6, 7, 8, 9, 10, 11, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 7, 8, 9, 10, 11, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 8, 9, 10, 11, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 9, 10, 11, 12, 10, 11, 12, 11, 12, 12, 10, 11, 12, 11, 12, 12, 11, 12, 12, 12, 1, 1, 6, 6, 2, 2, 3, 2, 1, 4, 1, 4, 2, 2, 4, 4, 3, 2, 1, 5, 1, 1, 5, 1, 4, 3, 2, 2, 2, 1, 1, 1, 6, 5, 1, 5, 6, 1, 2, 2, 2, 1, 2, 1, 2, 2, 5, 2, 1, 1, 1, 4, 3, 2, 3, 1, 3, 7, 7, 3, 1, 4, 3, 2, 1, 1, 3, 3, 2, 1, 4, 3, 2, 6, 3, 2, 3, 6, 4, 3, 3, 6, 2, 2, 7, 2, 3, 7, 2, 3, 3, 2, 3, 2, 1, 3, 4, 2, 2, 4, 4, 4, 8, 8, 3, 2, 1, 4, 4, 4, 3, 3, 4, 3, 4, 3, 7, 1, 3, 4, 7, 4, 8, 7, 3, 3, 4, 4, 8, 3, 1, 3, 1, 2, 4, 3, 4, 8, 1, 4, 1, 1, 5, 1, 1, 1, 4, 1, 4, 4, 8, 2, 4, 1, 1, 5, 8, 1, 2, 1, 1, 2, 4, 3, 5, 5, 1, 6, 3, 8, 5, 8, 7, 6, 5, 5, 2, 7, 4, 1, 3, 8, 7, 4, 7, 5, 5, 2, 7, 4, 6, 5, 8, 7, 5, 1, 8, 2, 4, 8, 1, 6, 3, 8, 5, 1, 5, 6, 2, 3, 6, 3, 1, 6, 2, 7, 1, 2, 1, 2))) ==================== X11.ash2 ======================= 8/6/92 c("-xrm 'splus*colors: black red 48 white blue 48 cyan'", "-xrm splus@nHalftones:0", "-xrm splus*halftoneImage:No") ==================== ct ==============color table======== 8/6/92 structure(.Dim = c(101, 3), .Data = c(128, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 155, 157.0408, 159.0816, 161.1225, 163.1633, 165.2041, 167.2449, 169.2857, 171.3265, 173.3673, 175.4082, 177.449, 179.4898, 181.5306, 183.5714, 185.6122, 187.6531, 189.6939, 191.7347, 193.7755, 195.8163, 197.8571, 199.898, 201.9388, 203.9796, 206.0204, 208.0612, 210.102, 212.1429, 214.1837, 216.2245, 218.2653, 220.3061, 222.3469, 224.3878, 226.4286, 228.4694, 230.5102, 232.551, 234.5918, 236.6327, 238.6735, 240.7143, 242.7551, 244.7959, 246.8367, 248.8775, 250.9184, 252.9592, 255, 128, 0, 5.204082, 10.40816, 15.61224, 20.81633, 26.02041, 31.22449, 36.42857, 41.63265, 46.83673, 52.04082, 57.2449, 62.44898, 67.65306, 72.85714, 78.06123, 83.2653, 88.46938, 93.67347, 98.87755, 104.0816, 109.2857, 114.4898, 119.6939, 124.898, 130.102, 135.3061, 140.5102, 145.7143, 150.9184, 156.1225, 161.3265, 166.5306, 171.7347, 176.9388, 182.1429, 187.3469, 192.551, 197.7551, 202.9592, 208.1633, 213.3673, 218.5714, 223.7755, 228.9796, 234.1837, 239.3878, 244.5918, 249.7959, 255, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 128, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) =================== ashn.s ======================= 8/6/92 ashn <- function( 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 ) { ##### S function to compute and plot n-d ASH slice over cube [-1,1]^3 ##### Copyright 1988,1992 David W. Scott # ix ----------- data structure containing: # nc - integer data matrix (n x nvar) - ignores skipped data # ab - bin range matrix (nvar x 2) # [ ab[1,1] , ab[1,2] ] contains x-axis range # [ ab[2,1] , ab[2,2] ] contains y-axis range # [ ab[3,1] , ab[3,2] ] contains z-axis range # nbin - vector length nvar; number of bins along each axis # mmax - maximum allowed value of m; dimensioned (m+1)^3 # iexp - alternate kernel (0==triangle 1=Epan 2=biweight 3=triweight) # regress - optional array of regression estimates (in "ix") nvar <- ncol(ix$nc); if(nvar<3) stop("less than 3 data variables") if( !is.matrix(ix$nc) | !is.integer(ix$nc) ) stop("ix$nc should be integer matrix") if( any(dim(ix$ab) != c(nvar,2)) ) stop("dim of ix$ab should be nvarx2") if( !is.integer(ix$nbin) | length(ix$nbin)!=nvar ) stop("ix$nbin should be integer vector nvar") if( !is.matrix(ctable) | !is.single(ctable) ) stop("ctable not an single matrix") # tri - 220 x 4 matrix # a. sides connected for triangle # b. corner that indicates decreasing/increasing # cuts - list of 500 triangles for the 254 cases # contour - 254 x 2 (256 cases: 0 and 255 imply no triangles) # a. pointer to first triangle # b. number of triangles in this case # orth? - default 0==> all triangles 1==> only "outside" triangles if( any(dim(master$tri) != c(220,4)) ) stop("master$tri not 220 x 4") if( length(master$cuts) != 500 ) stop("master$cuts not 500 long") if( any(dim(master$contour) != c(254,2)) ) stop("master$contour not 254 x 2") # data have been compressed into integer matrix # mmax is maximum m allowed for kernel array (default 6) mp1 <- rep( mmax+1, 3) # dimensioned (0:mmax) # WEIGHT ARRAY REALLY has 3 in dimension here # ctable: 1 row == background color ###of rows should be odd # rows 2-(nr+1)/2 are the OUTER colors -- edge to middle # last rows are INNER colors -- edge to middle # S refers to these colors as 0 to nr-1 color <- integer(4) color[1] <- 1; color[4] <- nrow(ctable)-1 color[2] <- max(color[4]/2,color[1]); color[3] <- min(color[2]+1,color[4]) # color[1:2] towards eye color[3:4] away from eye wk <- single( prod(mp1) ) # STRUCTURE( wk/REAL,ARRAY,3,mp1/ ) # ASH weight array xfmax <- single( ix$nbin[1] ) yfmax <- single( ix$nbin[2] ) zfmax <- single( ix$nbin[3] ) tx <- single( ix$nbin[1] ) ty <- single( ix$nbin[2] ) tz <- single( ix$nbin[3] ) if ( ireg != 1 ) {ash <- single( prod(ix$nbin[1:3]) ) } # ARRAY else { if( !is.array(ix$regress) | is.integer(ix$regress) ) { stop("regress not a real array") } if( any( dim(ix$regress) != ix$nbin) ) stop("ix$nbin and ix$regress do not agree") ash <- ix$regress } n <- nrow(ix$nc) z <- .Fortran("ashn", as.integer(ix$nc), as.integer(n),as.integer(nvar), as.single(ix$ab), as.integer(ix$nbin), ash=as.single(ash), as.single(xfmax), as.single(yfmax), as.single(zfmax), as.single(tx), as.single(ty), as.single(tz), as.single(wk), as.integer(mmax), as.integer(master$contour), as.integer(master$cuts), as.integer(master$tri), as.integer(color), as.integer(square), as.integer(sun50), as.single(ctable), as.integer(nrow(ctable)), as.integer(mv), as.single(ff), as.integer(ioauto), as.integer(ireg) ) invisible( array(z$ash,ix$nbin[1:3]) ) } ================== compress.s ======================= 8/6/92 compress <- function ( x, ab=t(apply(x,2,nicerange)), nbin=rep(30,ncol(x)) ) { n <- nrow(x); nvar <- ncol(x); if( ncol(ab) != 2 ) stop("ab should have 2 columns") if( nrow(ab)!=nvar||length(nbin)!=nvar) stop("dimensions of x,ab,nbin do not agree") if( min(nbin) <= 0 ) stop("nbin must be positive") ##### S function to compress data into binned integer form ##### Copyright 1988,1992 David W. Scott # x - input trivariate data matrix (n x nvar) # ab - input bin range matrix (nvar x 2) # x binned over [ ab[1,1], ab[1,2] ) # y binned over [ ab[2,1], ab[2,2] ) # z binned over [ ab[3,1], ab[3,2] ) etc. # nbin - vector length nvar - number of bins along each axis # nout - subset of n data points in bin range # nc - output matrix (nout x nvar) containing trivariate bin numbers del <- (ab[,2]-ab[,1])/nbin # bin widths iskip <- 0 # number of points not in trivariate bins nout <- 0 # number of points in output array # loop on data nc <- matrix(0,n,nvar); k <- NULL for( i in 1:nvar ) { if( ab[i,1] >= ab[i,2]) stop("ab array incorrectly ordered") nc[,i] <- floor( (x[,i]-ab[i,1])/del[i] + 1.0 ) k <- unique( c( k, seq(n)[ nc[,i]<1 | nc[,i]>nbin[i] ] ) ) } iskip <- length(k); if(iskip>0) {nc <- nc[-k,]; print(paste(as.character(iskip),"point(s) not in hyper-rectangle")) } storage.mode(nc) <- "integer" storage.mode(ab) <- "single" storage.mode(nbin) <- "integer" return(nc=nc,ab=ab,nbin=nbin) } nicerange <- function(x,beta=.1) { ab <- range(x) # interval surrounding data del <- (ab[2]-ab[1])*beta/2; return( c( ab+c(-del,del) ) ) } =================== make.ashn ====in local.Sqpe========= 8/6/92 #!/bin/csh f77 -c tmp2.r f77 -c sort.f Splus LOAD tmp2.o sort.o dumpit.c writeit.c cubeint.c simpint.c myinput.c =================== simpint.c ============== 8/6/92 #include struct connection { int n, start[5], end[5]; } C[5][5] = { {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}, {0},{3,{0,0,1},{2,3,1}},{3,{0,1,2},{3,3,2}},{0},{0}, {0},{3,{0,1,2},{3,1,2}},{0},{0},{0}, {3,{0,1,2},{0,1,2}},{3,{1,2,3},{1,2,3}},{0},{0},{0}, {0},{0},{0},{0},{0} }; simpint(P, F, X) float **P, *F, **X; { float *p[4], f[4]; int ord[4], nneg, nzero, nvertex, orientation, i, j; struct connection c; sort4(F, ord); nneg = nzero = 0; for(i = 0; i < 4; i++) { p[i] = P[ord[i]]; f[i] = F[ord[i]]; nneg += f[i] < 0; nzero += f[i] == 0; } c = C[nzero][nneg]; nvertex = c.n; if(nvertex == 0) return(0); for(i = 0; i < nvertex; i++) { int u = c.start[i], v = c.end[i], w; double a; if(u == v) for(j = 0; j < 3; j++) X[i][j] = p[u][j]; else { if(f[u] < 0) {w = u; u = v; v = w;} a = f[u] / (f[u] - f[v]); for(j = 0; j < 3; j++) X[i][j] = (1 - a) * p[u][j] + a * p[v][j]; } } if(nneg) orientation = which_side(X[0], X[1], X[2], p[0]); else orientation = -which_side(X[0], X[1], X[2], p[3]); if(orientation < 0) { float *tmp[4]; for(i = 0; i < nvertex; i++) tmp[i] = X[i]; for(i = 0; i < nvertex; i++) X[i] = tmp[nvertex-i-1]; } return(nvertex); } #define ORDER(p,q,r) {o[0]=p; o[1]=q; o[2]=r;} sort4(a, o) float *a; int *o; { int i, j; if(a[0] < a[1]) if(a[2] < a[0]) ORDER(2,0,1) else if(a[2] > a[1]) ORDER(0,1,2) else ORDER(0,2,1) else if(a[2] < a[1]) ORDER(2,1,0) else if(a[2] > a[0]) ORDER(1,0,2) else ORDER(1,2,0) for(i = 0; i < 3; i++) if(a[3] < a[o[i]]) break; for(j = 3; j > i; j--) o[j] = o[j-1]; o[i] = 3; } which_side(x0, x1, x2, p) float *x0, *x1, *x2, *p; { int i; float v1[3], v2[3], v[3], cross[3], dot; for(i = 0; i < 3; i++) { v1[i] = x2[i] - x1[i]; v2[i] = x0[i] - x1[i]; v[i] = p[i] - x1[i]; } cross[0] = v1[1]*v2[2] - v1[2]*v2[1]; cross[1] = v1[2]*v2[0] - v1[0]*v2[2]; cross[2] = v1[0]*v2[1] - v1[1]*v2[0]; dot = 0; for(i = 0; i < 3; i++) dot += cross[i] * v[i]; return(dot > 0 ? 1 : -1); } =================== cubeint.c =============== 8/6/92 struct { int a, b; } inner[6] = {{4,6}, {2,6}, {2,3}, {1,3}, {1,5}, {4,5}}; cubeint(origin, basis, F, X, size) float *origin, **basis, *F, **X; int *size; { int i, npatch; float c[24], *corner[8], *vertex[4], value[4], **x; for(i = 0; i < 8; i++) corner[i] = c + 3*i; make_corners(origin, basis, corner); vertex[0] = corner[0]; value[0] = F[0]; vertex[3] = corner[7]; value[3] = F[7]; x = X; npatch = 0; for(i = 0; i < 6; i++) { vertex[1] = corner[inner[i].a]; value[1] = F[inner[i].a]; vertex[2] = corner[inner[i].b]; value[2] = F[inner[i].b]; size[i] = simpint(vertex, value, x); x += size[i]; npatch += size[i] > 0; } return(npatch); } make_corners(o, b, c) float *o, **b, **c; { int i, j, k, l, n; n = 0; for(k = 0; k < 2; k++) for(j = 0; j < 2; j++) for(i = 0; i < 2; i++) { for(l = 0; l < 3; l++) c[n][l] = o[l] + i*b[0][l] + j*b[1][l] + k*b[2][l]; n++; } } =================== writeit.c ================= 8/6/92 #include #include "S.h" static FILE *f = NULL; F77_SUB(startwrite) (n) long *n; { char filename[20]; sprintf(filename, "fig.%d.quad", *n); f = fopen(filename, "w"); if(f == NULL) fprintf(stderr, "cannot create %s for writing\n", filename); fprintf(f, "CQUAD\n"); } #define VEC(x, m) \ { \ for(k = 0; k < m; k++) \ fprintf(f, " %g", x[k]); \ for(k = 0; k < 4; k++) \ fprintf(f, " %g", color[k]); \ fprintf(f, "\n"); \ } static float b0[3] = {1,0,0}, b1[3] = {0,1,0}, b2[3] = {0,0,1}; static float *basis[3] = {b0, b1, b2}; F77_SUB(docube) (origin, cube_size, values, contour, color) float *origin, *cube_size, *values, *contour, *color; { int i, j, k, size[6], npatch; float x[72], *points[24], **patch, F[8]; if(f == NULL) return; for(i = 0; i < 3; i++) basis[i][i] = cube_size[i]; for(i = 0; i < 8; i++) F[i] = values[i] - *contour; for(i = 0; i < 24; i++) points[i] = x + 3*i; npatch = cubeint(origin, basis, F, points, size); patch = points; for(i = 0; i < 6; i++) { if(size[i] == 0) continue; for(j = 0; j < size[i]; j++) VEC(patch[j], 3) if(size[i] == 3) VEC(patch[0], 3) fprintf(f, "\n"); patch += size[i]; } } F77_SUB(endwrite) () { if(f != NULL) fclose(f); } =================== dumpit.c ==================== 8/6/92 #include #include "S.h" static FILE *f = NULL; F77_SUB(startdump) (n) long *n; { char filename[20]; sprintf(filename, "figure.%d.quad", *n); f = fopen(filename, "w"); if(f == NULL) fprintf(stderr, "cannot create %s for writing\n", filename); fprintf(f, "CNQUAD\n"); } F77_SUB(dumpit) (p, q, r, n, a) float *p, *q, *r, *n, *a; { if(f == NULL) return; fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", p[0], p[1], p[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", q[0], q[1], q[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", r[0], r[1], r[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); fprintf(f, "%g %g %g %g %g %g %g %g %g %g\n", p[0], p[1], p[2], n[0], n[1], n[2], a[0], a[1], a[2], a[3]); fprintf(f, "\n"); } F77_SUB(enddump) () { if(f != NULL) fclose(f); } =================== sort.f ============================ 8/6/92 c dimension v(100) c integer a(100) c1 print *,'enter nv plus values' c read *,nv,(v(i),i=1,nv) c if (nv.eq.0) goto 999 c print *,'enter na plus integer pointers' c read *,na,(a(i),i=1,na) c call sort(v,a,1,na) c n=max0(na,nv) c print *,' perm=',(a(i),i=1,n) c print *,' sort=',(v(i),i=1,n) c goto1 c999 stop c end subroutine sort (v,a,ii,jj) c puts into a the permutation vector which sorts v into c increasing order. only elements from ii to jj are considered. c arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements c the vector v is sorted upon return c this is a modification of cacm algorithm #347 by r. c. singleton, c which is a modified hoare quicksort. dimension a(jj),v(1),iu(20),il(20) integer t,tt integer a real v m=1 i=ii j=jj 10 if (i.ge.j) go to 80 20 k=i ij=(j+i)/2 t=a(ij) vt=v(ij) if (v(i).le.vt) go to 30 a(ij)=a(i) a(i)=t t=a(ij) v(ij)=v(i) v(i)=vt vt=v(ij) 30 l=j if (v(j).ge.vt) go to 50 a(ij)=a(j) a(j)=t t=a(ij) v(ij)=v(j) v(j)=vt vt=v(ij) if (v(i).le.vt) go to 50 a(ij)=a(i) a(i)=t t=a(ij) v(ij)=v(i) v(i)=vt vt=v(ij) go to 50 40 a(l)=a(k) a(k)=tt v(l)=v(k) v(k)=vtt 50 l=l-1 if (v(l).gt.vt) go to 50 tt=a(l) vtt=v(l) 60 k=k+1 if (v(k).lt.vt) go to 60 if (k.le.l) go to 40 if (l-i.le.j-k) go to 70 il(m)=i iu(m)=l i=k m=m+1 go to 90 70 il(m)=k iu(m)=j j=l m=m+1 go to 90 80 m=m-1 if (m.eq.0) return i=il(m) j=iu(m) 90 if (j-i.gt.10) go to 20 if (i.eq.ii) go to 10 i=i-1 100 i=i+1 if (i.eq.j) go to 80 t=a(i+1) vt=v(i+1) if (v(i).le.vt) go to 100 k=i 110 a(k+1)=a(k) v(k+1)=v(k) k=k-1 if (vt.lt.v(k)) go to 110 a(k+1)=t v(k+1)=vt go to 100 end ==================== myinput.c =========================== 8/6/92 #include #include "S.h" static FILE *fin; F77_SUB(myio) (iopt, intx, realx, n) int *iopt; int *intx; float *realx; int *n; { int i; switch (*iopt) { case 0: fin = fopen("stream.s", "r"); break; case 1: for (i=0; i < *n; i++) { fscanf(fin, " %d", intx+i); } break; case 2: for (i=0; i < *n; i++) { fscanf(fin, " %f", realx+i); } break; case 3: fclose(fin); break; } } F77_SUB(mess) (iopt, intx, realx, n) int *iopt; int *intx; float *realx; int *n; { int i; switch (*iopt) { case 0: printf("\nACTION OPTION") ; printf("\n 0=plot 1=t 2=m 3=f 4=axes 5=cube 6=eye 7=slicelimit 8=sliceincr"); printf("\n 9=iris 10=orth 11=print 12=coltbl lookup 13=blank screen") ; printf("\n 14=alpha 15=fmaxin 16=kernel 99=quit > ") ; fflush(stdout); break; case 1: printf("\n z slice:"); fflush(stdout); break; case 2: printf("%d",*intx); fflush(stdout); break; case 3: printf("\ninside ASH check some calling parameters") ; printf("\n nvar = %d",*intx) ; printf("\n n = %d",*n); fflush(stdout); break; case 4: printf("\n mv ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; printf("\n hv ="); for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; fflush(stdout); break; case 5: printf(" echoing: ichoice = %d",*intx); break; case 6: printf("\n number of triangle patches found = %d",*intx); fflush(stdout); break; case 7: printf("\n sorting...."); fflush(stdout); break; case 8: printf("sort finished....beginning plot"); fflush(stdout); break; case 9: printf("\n"); fflush(stdout); break; case 10: printf("\n enter slice [ivar slice] = %d %d > ",*(intx+0),*(intx+1)); fflush(stdout); break; case 11: printf("\n mv ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; printf("; enter m > "); fflush(stdout); break; case 12: printf("\nError: smoothing parameter mmax exceeded for m[1,2,3], mmax = %d",*intx); fflush(stdout); break; case 13: printf("\nError: smoothing parameters not positive"); fflush(stdout); break; case 14: printf("\n ff = %g enter f level > ",*realx); fflush(stdout); break; case 15: printf("\n draw axes? 1=yes 0=no > "); fflush(stdout); break; case 16: printf("\n draw cube? 1=yes 0=no > "); fflush(stdout); break; case 17: printf("\n eye ="); for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; printf("; enter eye > "); fflush(stdout); break; case 18: printf("\nError: eye=0 try again"); fflush(stdout); break; case 19: printf("\n ix iy iz limits ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; printf("; enter > "); fflush(stdout); break; case 20: printf("\nError: out of range: ix iy iz maxrange ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; fflush(stdout); break; case 21: printf(" enter ix,iy,iz increments > "); fflush(stdout); break; case 22: printf("\nError: increments less than 1;"); fflush(stdout); break; case 23: printf("\n color offset 999 for an iris; 1=yes > "); fflush(stdout); break; case 24: printf("\n show orthogonal triangles only? 1=yes > "); fflush(stdout); break; case 25: printf("\n do you want to print 3D vertices 0-1-2-3 for 00 01 10 11 > "); fflush(stdout); break; case 26: printf("\n color table lookup ictopt = %d; enter > ",*intx); fflush(stdout); break; case 27: printf("\n blank screen? 0=no 1=yes; enter > "); fflush(stdout); break; case 28: printf("\n texture=a[1,2,3] alpha=a[4]: "); for (i=0; i < *n; i++) { printf(" %g",*(realx+i));}; printf("; enter > "); fflush(stdout); break; case 29: printf("\n fmax = %f; fmaxin = %f; enter fmaxin (neg ==> use fmax) > ", *(realx+0),*(realx+1)); fflush(stdout); break; case 30: printf("\nWarning: only 3 dimensions"); fflush(stdout); break; case 31: printf("\n working on point "); fflush(stdout); break; case 32: printf("\n fmax = %f; fmaxin = %f",*(realx+0),*(realx+1)); fflush(stdout); break; case 33: printf("\n n show nskip = %d %d %d",*(intx+0),*(intx+1),*(intx+2)); fflush(stdout); break; case 34: printf("\n ff = %g",*realx); fflush(stdout); break; case 35: printf("; c = %f",*realx); fflush(stdout); break; case 36: printf("\n t slice(s) ="); for (i=0; i < *n; i++) { printf(" %d",*(intx+i));}; fflush(stdout); break; case 37: printf("\n iexp = %d 0=tri 1=Epan 2=biwt 3=triwt",*intx); printf("; enter iexp > "); fflush(stdout); break; case 38: printf("\nWarning: iexp>0 kernel shape correct but not normalized"); fflush(stdout); break; case 39: printf("."); fflush(stdout); break; case 40: printf(" %d",*intx); fflush(stdout); break; } } ==================== tmp2.r ====this is ashn.r run through m4======== 8/6/92 #-July 20, 1988 # plot triangularization September 13, 1990 #-August 1, 1992 automatic I/O option subroutine ashn ( nc, n, nvar, ab, nbin, ash, xfmax, yfmax, zfmax, tx, ty, tz, wk, mmax, contour, cuts, tri, color, square, sun50, ictable, ict, mv, ff, ioauto, ireg ) dimension hv(10),nbin(nvar),mv(nvar),itbin(10) # nvar <=10 assumed dimension nc(n,nvar),ab(nvar,2),ash(nbin(1),nbin(2),nbin(3)) dimension v1(3),v2(3),v3(3),icolor(5),eye(3) dimension ixyz(3),ixyzl(3),ixyzr(3) dimension ileft(3),irght(3),idel(3) dimension xfmax(nbin(1)),yfmax(nbin(2)),zfmax(nbin(3)) dimension tx(nbin(1)),ty(nbin(2)),tz(nbin(3)) dimension tdel(3),txyz(3),xyz(3,12),ee(3,3),ttxyz(3),tnormal(3),talpha(4) dimension px(12),py(12),pz(12) dimension wk(0:mmax,0:mmax,0:mmax) # kernel weights dimension cb(16,3) # cube integer ictable(ict,1) # color table lookup integer orth, contour(254,2), cuts(500), tri(220,4) dimension vout(3),vout2(3) # poutward normal vector integer ic(3,11) # integer coordinates of corner but [0,1]^3 + 3 faces dimension as(0:1,0:1,0:1),rc(3,11) integer color(4), square(2), sun50 dimension pxy(3,2,25000),ztri(25000) # big long list here real angle(25000) # pos==>normal towards eye; neg==>away integer iout(10) dimension rout(10) integer iperm(25000) # sorted ztri pointers (could combine with ishow # if rewrote sort to handle subsets other # than 1,2,3,4,.....actually since # orth is input, shouldn't need this # hv--bin width in data units # nbin--number bins each dimension # mv--number shifts in averaging real am(200) common/bgrp/am common/encbfc/encbf0 common/encbfl/iqqbfl,iqqfil character*256 encbf0 ; integer iqqbfl,iqqfil common/decbfc/decbf0 character*512 decbf0 common/decbfl/jqqbfl,jqqinf,jqqflp,jqqfle,jqqlsc,jqdchk integer jqqbfl,jqqinf,jqqflp,jqqlsc; logical jqqfle,jqdchk external qproject, qtriout, qoutward data ntrimax/25000/ # cube data cb/ 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, -1,-1,-1,-1, 1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1/ ## data ff/.50/ # default slice level --- 9/11/90 IN CALL # corners data ic/ 0,0,0, 1,0,0, 1,1,0, 0,1,0, 0,0,1, 1,0,1, 1,1,1, 0,1,1, 0,0,0, 0,0,0, 0,0,0/ # fills columns...last three refer to corner 1 data rc/ -1.,-1.,-1., 1.,-1.,-1., 1., 1.,-1., -1., 1.,-1., -1.,-1., 1., 1.,-1., 1., 1., 1., 1., -1., 1., 1., 0., 0.,-1., -1., 0., 0., 0.,-1., 0./ # fills columns ... last 3 are faces data tnormal/1.,1.,1./ # dummy normal for iris data talpha/1.,0.,0.,1./ # texture plus alpha channel if (ioauto.eq.1) call myio(0,iout,rout,0) ### call beginz # a null call to put iblank right { rxxxx1 = 0.1; rxxxx2 = 0.1; rxxxx3 = 0.1; rxxxx4 = 0.1 call spmarz(rxxxx1,rxxxx2,rxxxx3,rxxxx4) ; } s3 = sqrt(3.0) { am(61) = square(1) am(62) = square(2) am(63) = square(1) am(64) = square(2) call zscalz ; } do j=1,8{do i=1,3{rc(i,j)=rc(i,j)/s3}} # direction vectors of corners (j=1,8 not 11) # 3 face vectors already length 1 call mess(3,nvar,0,n) # echoing nvar,n # default values npict = 0 # number of pictures iprint = 0 # print triangle vertices CHANGED iris=1; icoloroffset = 999 # color offset (0 for Sun, 999 for Iris) 9/10/90 CHANGED inew = 1 # new parameters, recompute ash iblank = 1 # blank out screen orth = 0 # default plot all triangles iaxis = 1 # slice on x-axis ix1=1; iy1=1; iz1=1; ix2=nbin(1)-1; iy2=nbin(2)-1; iz2=nbin(3)-1 # slice limits ix3=1; iy3=1; iz3=1 # slice increments icolor(1)=4; icolor(2)=1; icolor(3)=0; icolor(4)=3; icolor(5)=2 iexp=0; iexp1=1; iexp2=1; if(iexp>0) {iexp1=2; iexp2=iexp } # kernel option wk(0,0,0)=1.0 # kernel array is 1 x 1 x 1 ##INCALL## do i=1,nvar { mv(i)=1} # default smoothing parameters CHANGED do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} do i=1,nvar { hv(i)=(ab(i,2)-ab(i,1))/nbin(i) } call mess(4,mv,hv,nvar) # echoing mv,hv eye(1)=1; eye(2)=1; eye(3)=0 v1(1)=-sqrt(2.)/2; v1(2)=-v1(1); v1(3)=0 v2(1)=0; v2(2)=0; v2(3)=1 v3(1)=v1(2); v3(2)=v3(1); v3(3)=0 fmax=0.0; fmaxin=-1000.0 # fmaxin max of f(.) iskip=0; ishow=0 # numbers points skipped and shown do i=1,3{ileft(i)=1; irght(i)=nbin(i); idel(i)=1}# select slices along axes to display do i=4,nvar{itbin(i)=1} # default t-slice ifleft = 1; ifright = 1; ifdel = 1 # select density levels to display iqaxis = 1 # draw coordinate axes at center icube = 1 # draw unit cube ictopt = 0 # use color table lookup values (right for postscript!) CHANGED goto 6663 5555 continue # change default values call mess(0,iout,rout,0) # ACTION OPTIONS if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(ichoice,.true.) }} else {call myio(1,ichoice,rout,1)} if(ichoice!=99) call mess(5,ichoice,rout,0) # echoing ichoice if(ichoice==0) goto 6663 switch(ichoice){ case 1: # change t slice if(nvar==3) { call mess(30,0,0,0); goto 5555} # warning: only 3 dimensions inew=1 # recompute ash do i=4,nvar {iout(1)=i;iout(2)=itbin(i) call mess(10,iout,0,2) # ivar slice if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(itbin(i),.true.) }} else {call myio(1,itbin(i),rout,1)} } case 2: # change smoothing parameters 6671 inew = 1 # recompute ash call mess(11,mv,0,nvar) # current m values if(nvar==3) { if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(mv(1),.true.) call deci00(mv(2),.true.) call deci00(mv(3),.true.) }} else {call myio(1,mv,rout,3)} } else { if(nvar==4) { if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(mv(1),.true.) call deci00(mv(2),.true.) call deci00(mv(3),.true.) call deci00(mv(4),.true.) }} else {call myio(1,mv,rout,4)} } else { do i=1,nvar{ iout(1)=i;iout(2)=mv(i) if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(mv(i),.true.) }} else {call myio(1,mv(i),rout,1)} } } } if(max0(mv(1),mv(2),mv(3)).gt.mmax) { call mess(12,mmax,0,0); goto 6671} # m too big mmin=mmax; do i=1,nvar{mmin=min0(mmin,mv(i))}; if(mmin<1) { call mess(13,0,0,0) ; goto 6671} # m too small do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) # change ASH wts do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} case 3: # change contour levels call mess(14,0,ff,1) # new f level if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call decr00(ff) }} else {call myio(2,iout,ff,1)} case 4: # draw coordinate axes? call mess(15,0,0,0) # draw axes? if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(iqaxis,.true.) }} else {call myio(1,iqaxis,rout,1)} case 5: # draw unit cube? call mess(16,0,0,0) # draw cube? if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(icube,.true.) }} else {call myio(1,icube,rout,1)} case 6: # change eye 6662 call mess(17,0,eye,3) # enter eye if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call decr00(eye(1)) call decr00(eye(2)) call decr00(eye(3)) }} else {call myio(2,iout,eye,3)} e1 = eye(1); e2 = eye(2); e3 = eye(3) rho = e1**2 + e2**2 + e3**2; rho=sqrt(rho) if(rho.le.0.) {call mess(18,0,0,0); goto 6662} # eye is zero phi = acos(e3/rho); pi = 3.14159265 if( e1!=0. | e2!=0. ) theta = atan2(eye(2),eye(1)) # atan2 handles special cases if( e1==0. & e2==0. & e3>0. ) theta=-pi/2 if( e1==0. & e2==0. & e3<0. ) theta= pi/2 # don't have RATFOR manual sp=sin(phi); cp=cos(phi); st=sin(theta); ct=cos(theta) v1(1)=-st; v1(2)=ct; v1(3)=0.0 v2(1)=-ct*cp; v2(2)=-st*cp; v2(3)=sp v3(1)=ct*sp; v3(2)=st*sp; v3(3)=cp #right hand system (not graphics LHS) case 7: # axis contour slices shown 6677 iout(1)=ix1;iout(2)=ix2;iout(3)=iy1;iout(4)=iy2;iout(5)=iz1;iout(6)=iz2 call mess(19,iout,0,6) # ix iy iz limits if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(ix1,.true.) call deci00(ix2,.true.) call deci00(iy1,.true.) call deci00(iy2,.true.) call deci00(iz1,.true.) call deci00(iz2,.true.) }} else { call myio(1,ix1,rout,1); call myio(1,ix2,rout,1); call myio(1,iy1,rout,1); call myio(1,iy2,rout,1); call myio(1,iz1,rout,1); call myio(1,iz2,rout,1) } if(ix1<1|iy1<1|iz1<1|ix2>=nbin(1)|iy2>=nbin(2)|iz2>=nbin(3)) { iout(1)=1;iout(2)=nbin(1)-1;iout(3)=1;iout(4)=nbin(2)-1;iout(5)=1;iout(6)=nbin(3)-1 call mess(20,iout,0,6); goto6677 } # ix,iy,iz out of limits case 8: # contour slice increments 6678 call mess(21,0,0,0) if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(ix3,.true.) call deci00(iy3,.true.) call deci00(iz3,.true.) }} else { call myio(1,ix3,rout,1); call myio(1,iy3,rout,1); call myio(1,iz3,rout,1) } if(min0(ix3,iy3,iz3)<1) {call mess(22,0,0,0); goto 6678} case 9: #-which contour levels to plot? call mess(23,0,0,0) # color offset 999 for an iris if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(iris,.true.) }} else {call myio(1,iris,rout,1)} if(iris.ne.1) icoloroffset=0 # 9/10/90 if(iris.eq.1) icoloroffset=999 # 9/10/90 case 10: # show only "outside" triangles call mess(24,0,0,0) #show orthogonal triangles only? if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(orth,.true.) }} else {call myio(1,orth,rout,1)} case 11: # print triangle vertices call mess(25,0,0,0) # do you want to print 3D vertices 0-1-2-3 for 00 01 10 11 if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(iprint,.true.) }} else {call myio(1,iprint,rout,1)} if(iprint==0) { call enddump; call endwrite } case 12: # use color table lookup or not! call mess(26,ictopt,0,0) # ictopt= if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(ictopt,.true.) }} else {call myio(1,ictopt,rout,1)} case 13: # blank screen option call mess(27,0,0,0) # blank screen 0-1? if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(iblank,.true.) }} else {call myio(1,iblank,rout,1)} case 14: # read alpha channel stuff call mess(28,0,talpha,4) # texture=a[1,2,3] alpha=a[4] if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call decr00(talpha(1)) call decr00(talpha(2)) call decr00(talpha(3)) call decr00(talpha(4)) }} else {call myio(2,iout,talpha,4)} case 15: # input fmaxin rout(1)=fmax;rout(2)=fmaxin call mess(29,0,rout,2) # enter fmaxin: neg==> use fmax if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call decr00(fmaxin) }} else {call myio(2,iout,fmaxin,1)} if(fmaxin.le.0.0) fmaxin=fmax case 16: # change kernel iexp 0==triangle 1==Epan 2==Biweight 3==Triwt 6971 inew = 1 # recompute ash call mess(37,iexp,0,1) # current "iexp" value if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(iexp,.true.) }} else {call myio(1,iexp,rout,1)} if(iexp<0) { goto 6971} # iexp negative if(iexp>0) call mess(38,0,0,0) iexp1=1; iexp2=1; if(iexp>0) {iexp1=2; iexp2=iexp } # kernel option do k=0,mv(3)-1 { zfact = xkern(k,mv(3),iexp1,iexp2) # change ASH wts do j=0,mv(2)-1 { yfact = xkern(j,mv(2),iexp1,iexp2) do i=0,mv(1)-1 {wk(i,j,k) = xkern(i,mv(1),iexp1,iexp2)*yfact*zfact}}} case 99: #quit 9000 if (ioauto.eq.1) call myio(3,iout,rout,0) call enddump; call endwrite call mess(9,0,0,0) goto 9005 } goto 5555 6663 continue if (inew.eq.0) goto 6611 # do not have to recompute ash #-zero ASH array -- computed at centers of trivariate bins 6666 call mess(31,0,0,0) # initializing ASH array; working on point fmax = 0. # ASH value at mode iskip =0 ; ishow =0 do iz = 1,nbin(3) { do iy = 1,nbin(2) { do ix = 1,nbin(1) { if(ireg==1) {fmax=amax1(fmax,ash(ix,iy,iz))} else {ash(ix,iy,iz) = 0.0} }}} if (ireg==1) goto 8642 # visualize regression data do i = 1,n { if((i/500)*500.eq.i) {call mess(2,i,0,0) } # looping thru data, on point i else{ if((i/100)*100.eq.i) {call mess(39,0,0,0)} } tfact=1.0 # part of product kernel in reference dimensions (4 to nvar) do j = 4,nvar { # null loop if nvar=3 idt = iabs(nc(i,j)-itbin(j)) if( idt >= mv(j) ){ # this point not covered by hyperrectangle iskip=iskip+1; goto 1110 } tfact = tfact * xkern(idt,mv(j),iexp1,iexp2) } ishow = ishow+1 do j = 1,3 { ixyz(j) = nc(i,j) ixyzl(j) = max0(1,ixyz(j)-(mv(j)-1)) ixyzr(j) = min0(nbin(j),ixyz(j)+(mv(j)-1)) } #-update ASH array in (x,y,z|t) space do iz = ixyzl(3),ixyzr(3) { idz = iabs(iz-ixyz(3)) do iy = ixyzl(2),ixyzr(2) { idy = iabs(iy-ixyz(2)) do ix = ixyzl(1),ixyzr(1) { idx = iabs(ix-ixyz(1)) ash(ix,iy,iz)=ash(ix,iy,iz) + wk(idx,idy,idz)*tfact } } } 1110 continue # next data point } #-scale by n*hx*hy*hz; find max, etc 8642 continue # regression visualization: skip density estimation above do i = 1,nbin(1) { xfmax(i) = 0 } do i = 1,nbin(2) { yfmax(i) = 0 } do i = 1,nbin(3) { zfmax(i) = 0 } # max of "CELL" in that plane c = n; do i = 1,nvar { c = c * (mv(i)*hv(i)) }; if(ireg==1) {c=1} # leave regression "ash" values alone do iz = 1,nbin(3) { do iy = 1,nbin(2) { do ix = 1,nbin(1) { a = ash(ix,iy,iz) / c ash(ix,iy,iz) = a xfmax(ix) = amax1(a,xfmax(ix)) yfmax(iy) = amax1(a,yfmax(iy)) zfmax(iz) = amax1(a,zfmax(iz)) } } } fmax = 0; do i = 1,nbin(1) { fmax = amax1(fmax,xfmax(i)) } if( nvar==3) fmaxin=fmax # don't worry about higher dimensions 475 continue rout(1)=fmax;rout(2)=fmaxin; call mess(32,0,rout,2) # print fmax,fmaxin iout(1)=n;iout(2)=ishow;iout(3)=iskip; call mess(33,iout,0,3) # n show skip 6611 continue #-draw triangularization if (iblank.eq.1) { # clear and redraw screen if iblank=1 #??? call finisz if (iprint.ne.0) {call enddump; call endwrite npict = npict + 1 if(iprint==1|iprint==3) call startdump(npict) if(iprint.ge.2) call startwrite(npict) } call beginz; call boxz } do i=1,3 { tdel(i) = 2.0/nbin(i) } c=-1.-tdel(1)/2; do i=1,nbin(1) { tx(i)=c+i*tdel(1) } # make into 3x3 array c=-1.-tdel(2)/2; do i=1,nbin(2) { ty(i)=c+i*tdel(2) } #####!!!!! c=-1.-tdel(3)/2; do i=1,nbin(3) { tz(i)=c+i*tdel(3) } ntri=0 #### number of triangles found if( fmaxin>0.0 ) c = fmaxin * ff # contour level else c = fmax * ff call mess(34,0,ff,1) # print ff call mess(35,0, c,1) # print contour level c if (iprint==1 | iprint==3) { ttxyz(1)=tx(ix1);ttxyz(2)=ty(iy1);ttxyz(3)=tz(iz1) call dumpit(ttxyz,ttxyz,ttxyz,tnormal,talpha) ttxyz(1)=tx(ix2+1);ttxyz(2)=ty(iy2+1);ttxyz(3)=tz(iz2+1) call dumpit(ttxyz,ttxyz,ttxyz,tnormal,talpha)} #----------------------------------------------- if(nvar>3) call mess(36,itbin(4),0,nvar-3) # print t slice(s) call mess(1,0,0,0) # z slice do iz = iz1,iz2,iz3{ if(zfmax(iz)<=c & zfmax(iz+1)<=c) goto 603 call mess(40,iz,0,0) # z slice value is iz txyz(3)=tz(iz) do iy = iy1,iy2,iy3{ if(yfmax(iy)<=c & yfmax(iy+1)<=c) goto 602 txyz(2)=ty(iy) do ix = ix1,ix2,ix3{ if(xfmax(ix)<=c & xfmax(ix+1)<=c) goto 601 txyz(1)=tx(ix) do idz=0,1{do idy=0,1{do idx=0,1{ a=ash(ix+idx,iy+idy,iz+idz); as(idx,idy,idz)=a if(a==c) {call intpr("opps contour level equality problem===adjusting ff and c", 56,0,0); ff=1.00001*ff; c=1.00001*c }}}} # Fatal: adjust ff by a bit if( iprint.ge.2 ) call docube( txyz, tdel, as, c, talpha) ncase=0 if(as(0,1,1)>c)ncase=ncase+1 # corner 8 if(as(1,1,1)>c)ncase=ncase+2 # corner 7 if(as(1,0,1)>c)ncase=ncase+4 # corner 6 if(as(0,0,1)>c)ncase=ncase+8 # corner 5 if(as(0,1,0)>c)ncase=ncase+16 # corner 4 if(as(1,1,0)>c)ncase=ncase+32 # corner 3 if(as(1,0,0)>c)ncase=ncase+64 # corner 2 if(as(0,0,0)>c)ncase=ncase+128 # corner 1 if(ncase==0 | ncase==255) go to 601 if(ix2-ix1+iy2-iy1+iz2-iz1.eq.0) call realpr("f xyz",5,as,8) do i=1,12{px(i)=0.;py(i)=0.} # not really necessary if working call qproject(c,0,0,0,as(0,0,0),as(1,0,0),1,txyz,tdel,1, v1,v2,v3,px,py,pz,xyz(1,1)) call qproject(c,1,0,0,as(1,0,0),as(1,1,0),2,txyz,tdel,2, v1,v2,v3,px,py,pz,xyz(1,2)) call qproject(c,0,1,0,as(0,1,0),as(1,1,0),3,txyz,tdel,1, v1,v2,v3,px,py,pz,xyz(1,3)) call qproject(c,0,0,0,as(0,0,0),as(0,1,0),4,txyz,tdel,2, v1,v2,v3,px,py,pz,xyz(1,4)) call qproject(c,0,0,1,as(0,0,1),as(1,0,1),5,txyz,tdel,1, v1,v2,v3,px,py,pz,xyz(1,5)) call qproject(c,1,0,1,as(1,0,1),as(1,1,1),6,txyz,tdel,2, v1,v2,v3,px,py,pz,xyz(1,6)) call qproject(c,0,1,1,as(0,1,1),as(1,1,1),7,txyz,tdel,1, v1,v2,v3,px,py,pz,xyz(1,7)) call qproject(c,0,0,1,as(0,0,1),as(0,1,1),8,txyz,tdel,2, v1,v2,v3,px,py,pz,xyz(1,8)) call qproject(c,0,0,0,as(0,0,0),as(0,0,1),9,txyz,tdel,3, v1,v2,v3,px,py,pz,xyz(1,9)) call qproject(c,1,0,0,as(1,0,0),as(1,0,1),10,txyz,tdel,3, v1,v2,v3,px,py,pz,xyz(1,10)) call qproject(c,1,1,0,as(1,1,0),as(1,1,1),11,txyz,tdel,3, v1,v2,v3,px,py,pz,xyz(1,11)) call qproject(c,0,1,0,as(0,1,0),as(0,1,1),12,txyz,tdel,3, v1,v2,v3,px,py,pz,xyz(1,12)) do k=1,contour(ncase,2){ # number of triangles ktri = cuts( contour(ncase,1)+k-1 ) # triangle number is1=tri(ktri,1); is2=tri(ktri,2); is3=tri(ktri,3) iplot=1 # plot this triangle #### always checks now if(orth==1) { # check qoutward normal vs. eye do ii=1,3{ee(ii,1)=xyz(ii,is2)-xyz(ii,is1)} # vector side 1 do ii=1,3{ee(ii,2)=xyz(ii,is3)-xyz(ii,is2)} do ii=1,3{ee(ii,3)=xyz(ii,is1)-xyz(ii,is3)} call qoutward(ee,vout) # a normal vector jj=tri(ktri,4) # a distance corner from triangle if(jj<=0){call intpr("-----mistake: used triangle",27,0,0); return} pp=0.; do ii=1,3{pp=pp+vout(ii)*rc(ii,jj)} if(pp<0.) do ii=1,3{vout(ii)=-vout(ii)} # direction towards corner if(as(ic(1,jj),ic(2,jj),ic(3,jj))>c) { do ii=1,3{vout(ii)=-vout(ii)}} # direction now towards out pp=0.; do ii=1,3{pp=pp+vout(ii)*v3(ii)}; if(pp<=0.)iplot=0 ## vout2(1) = ee(2,1)*ee(3,2)-ee(3,1)*ee(2,2) ## vout2(2) = -(ee(1,1)*ee(3,2)-ee(3,1)*ee(1,2)) ## vout2(3) = ee(1,1)*ee(2,2)-ee(2,1)*ee(1,2) ##flip=0.; do ii=1,3{flip=flip+vout(ii)*vout2(ii)} flip=0.; do ii=1,3{flip=flip+vout(ii)**2}; flip=sqrt(flip); do ii=1,3{vout2(ii)=vout(ii)/flip} if(iprint==1|iprint==3) call dumpit(xyz(1,is1),xyz(1,is2),xyz(1,is3),vout2,talpha) ### call realpr("cosine eye",11,pp,1) ############### } if( !(orth==1 & iplot==0)) { # plot this triangle ### call segmtz(px(is1),py(is1),px(is2),py(is2),1) ### call segmtz(px(is2),py(is2),px(is3),py(is3),1) ### call segmtz(px(is3),py(is3),px(is1),py(is1),1) ntri = ntri + 1 ### call intpr("------------ntr--------",23,ntri,1) if(ntri>ntrimax) {ntri=ntri-1000 call intpr("----ntrimax----deleted 1000 triangles----",41,0,0)} angle(ntri)=pp pxy(1,1,ntri)=px(is1); pxy(1,2,ntri)=py(is1) pxy(2,1,ntri)=px(is2); pxy(2,2,ntri)=py(is2) pxy(3,1,ntri)=px(is3); pxy(3,2,ntri)=py(is3) ztri(ntri)=amin1(pz(is1),pz(is2),pz(is3)) } do i=1,3{ii=tri(ktri,i);if(px(ii)==0. & py(ii)==0.) call intpr("=======error=======",19,0,0)} } 601 continue # next ix } # end ix loop 602 continue # next iy } # end iy loop 603 continue # next iz } # end iz loop # ------------------------------------ call mess(6,ntri,rout,0) # number of triangle patches found if(ntri>0) {call mess(7,0,0,0); do i=1,ntri{iperm(i)=i} call sort (ztri,iperm,1,ntri); call mess(8,0,0,0) } # sorting istop = 2500000 # SHOULD BE 0 do ii=1,ntri{ if(ii>istop) {call intpr("enter stopping point",20,0,0) if(ioauto!=1) {{jqqinf=0; jqdchk=.true.; call gtinbf #initialize buffer call deci00(istop,.true.) }} else {call myio(1,istop,rout,1)} } i=iperm(ii) if(istop<0) i=iperm(-istop) ### for debugging purposes z=angle(i) if( ictopt.ne.1 ) { ### not postscript if(z>=0.) { k=color(1)+z*(color(2)-color(1)+.99)} else { k=color(3)+abs(z)*(color(4)-color(3)+.99)} # abs(z) 7-28 if(k>color(4)) call intpr("================color too large===============",46,0,0) } if( ictopt.eq.1 ) { ### for postscript 9-25-89 if(z>=0.) { k=color(1)+z*(color(2)-color(1)+.99)} else { k=color(3)+abs(z)*(color(4)-color(3)+.99)} k=ictable(k,1)} if(sun50==0) {{ am(10) = k+icoloroffset ; }} else {{ am(8) = k ; }} call qtriout(pxy(1,1,i),square) ## call segmtz(pxy(1,1,i),pxy(1,2,i),pxy(2,1,i),pxy(2,2,i),1) ## call segmtz(pxy(2,1,i),pxy(2,2,i),pxy(3,1,i),pxy(3,2,i),1) ## call segmtz(pxy(3,1,i),pxy(3,2,i),pxy(1,1,i),pxy(1,2,i),1) } # print axes if(sun50==0) {{ am(10) = 1+icoloroffset ; }} else {{ am(8) = 1 ; }} if (iqaxis==1) { p1x=.5*v1(1); p1y=.5*v2(1) # x-axis p2x=.5*v1(2); p2y=.5*v2(2) # y-axis p3x=.5*v1(3); p3y=.5*v2(3) # z-axis p4x=0.; p4y=0. # origin p1x=square(1) + (square(2)-square(1))*(p1x+s3)/(2*s3) + .999 p2x=square(1) + (square(2)-square(1))*(p2x+s3)/(2*s3) + .999 p3x=square(1) + (square(2)-square(1))*(p3x+s3)/(2*s3) + .999 p4x=square(1) + (square(2)-square(1))*(p4x+s3)/(2*s3) + .999 p1y=square(1) + (square(2)-square(1))*(p1y+s3)/(2*s3) + .999 p2y=square(1) + (square(2)-square(1))*(p2y+s3)/(2*s3) + .999 p3y=square(1) + (square(2)-square(1))*(p3y+s3)/(2*s3) + .999 p4y=square(1) + (square(2)-square(1))*(p4y+s3)/(2*s3) + .999 call segmtz(p1x,p1y,p4x,p4y,1) call segmtz(p2x,p2y,p4x,p4y,1) call segmtz(p3x,p3y,p4x,p4y,1) call textz(p1x,p1y,"x"//char(0)) call textz(p2x,p2y,"y"//char(0)) call textz(p3x,p3y,"z"//char(0)) } # draw the cube !!! if (icube==1) { p1x=0.; do j=1,3{p1x=p1x+v1(j)*cb(1,j)} p1y=0.; do j=1,3{p1y=p1y+v2(j)*cb(1,j)} do i = 2,16{ p2x=0.; do j=1,3{p2x=p2x+v1(j)*cb(i,j)} p2y=0.; do j=1,3{p2y=p2y+v2(j)*cb(i,j)} p1xx=square(1) + (square(2)-square(1))*(p1x+s3)/(2*s3) + .999 p2xx=square(1) + (square(2)-square(1))*(p2x+s3)/(2*s3) + .999 p1yy=square(1) + (square(2)-square(1))*(p1y+s3)/(2*s3) + .999 p2yy=square(1) + (square(2)-square(1))*(p2y+s3)/(2*s3) + .999 call segmtz(p1xx,p1yy,p2xx,p2yy,1) p1x=p2x; p1y=p2y } } inew = 0 # can use current ash for many parameter changes go to 5555 9005 return end subroutine qproject(c,ix,iy,iz,f1,f2,ns,txyz,tdel,iptr,v1,v2,v3,px,py,pz,xyz) dimension v1(3),v2(3),v3(3),txyz(3),tdel(3),px(12),py(12),pz(12),xyz(3),ixyz(3) # computes intersection of (f1,f2) at level c # txyz contains point of intersection # iptr indicates subscript in middle of a side # ixyz indicates beginning of line segment (== side) # txyz contains the (xyz) coordinates of origin # ns contains number of side # xyz should be the "ns-th" column of xyz above if( amax1(f1,f2)<=c | amin1(f1,f2)>=c ) return ixyz(1)=ix; ixyz(2)=iy;ixyz(3)=iz denom = f2-f1; if(denom==0.) {call intpr("denom========0",14,0,0);return} fc = (c-f1)/denom if (fc>0. & fc<1.) { # found contour do i=1,3{xyz(i)=txyz(i)+ixyz(i)*tdel(i)} xyz(iptr) = txyz(iptr) + fc*tdel(iptr) px(ns)=v1(1)*xyz(1) + v1(2)*xyz(2) + v1(3)*xyz(3) py(ns)=v2(1)*xyz(1) + v2(2)*xyz(2) + v2(3)*xyz(3) pz(ns)=v3(1)*xyz(1) + v3(2)*xyz(2) + v3(3)*xyz(3) ###call intpr(" intersection found on side",31,ns,1) } return end subroutine qoutward(e,vout) # compute a normal vector dimension e(3,3),vout(3), det(3) det(1) = e(2,1)*e(3,2) - e(3,1)*e(2,2) det(2) = e(1,1)*e(3,2) - e(3,1)*e(1,2) det(3) = e(1,1)*e(2,2) - e(2,1)*e(1,2) id=1; detm=abs(det(1)) do i=2,3{if(abs(det(i))>detm){id=i; detm=abs(det(i))}} if( detm==0.) { call intpr("determinant is zero!!!",22,0,0);return } ## id = 1 ####---------------TEMP FIX----------##### switch(id){ case 1: vout(1)=1. vout(2) = ( -e(1,1)*e(3,2) + e(1,2)*e(3,1) ) / det(1) vout(3) = ( -e(2,1)*e(1,2) + e(2,2)*e(1,1) ) / det(1) case 2: vout(1) = ( -e(2,1)*e(3,2) + e(2,2)*e(3,1) ) / det(2) vout(2)=1. vout(3) = ( -e(1,1)*e(2,2) + e(1,2)*e(2,1) ) / det(2) case 3: vout(1) = ( -e(3,1)*e(2,2) + e(3,2)*e(2,1) ) / det(3) vout(2) = ( -e(1,1)*e(3,2) + e(1,2)*e(3,1) ) / det(3) vout(3)=1. } pp=0.; do i=1,3{pp=pp+vout(i)**2}; pp=sqrt(pp) do i=1,3{vout(i)=vout(i)/pp} return end subroutine qtriout(xyorig,square) real xyorig(3,2), xy(3,2), y(3) # y(3) rather than y(2) if an error integer square(2) # find vertical lines intersecting triangle real am(200) common/bgrp/am common/encbfc/encbf0 common/encbfl/iqqbfl,iqqfil character*256 encbf0 ; integer iqqbfl,iqqfil common/decbfc/decbf0 character*512 decbf0 common/decbfl/jqqbfl,jqqinf,jqqflp,jqqfle,jqqlsc,jqdchk integer jqqbfl,jqqinf,jqqflp,jqqlsc; logical jqqfle,jqdchk do j=1,2{ do i=1,3{ xy(i,j)=square(1) + (xyorig(i,j)+sqrt(3.))*(square(2)-square(1))/(2*sqrt(3.)) + .999}} call zpolyz( xy(1,1), xy(1,2), 3 ) # 7-11/89 shaded polygon?? return ################-------------------replaced last part-----------------############# 987 continue x1=amin1(xy(1,1),xy(2,1),xy(3,1)); x2=amax1(xy(1,1),xy(2,1),xy(3,1)) k1 = x1 + .999 k2 = x2 ##############work on this so rounds rather than truncates xy(*,*) ###### k1 and k2 as well ###print(tstring(x1 x2 y1 y2 k1 k2),R(x1),R(x2),R(y1),R(y2),I(k1),I(k2)) do k=k1,k2 { iptr=0 x=k #### could add case where denom=0 (ie x1=x2) but x=x1=x2 (a vertical line, get 2 points) denom=(xy(2,1)-xy(1,1)) if(denom!=0.) { a = (x-xy(1,1))/denom if(a>=0.&a<=1.) { yy = (xy(2,2)-xy(1,2))*a + xy(1,2) iptr=iptr+1;y(iptr)=yy ####;print(I(1),R(x),R(yy)) } } denom=(xy(3,1)-xy(2,1)) if(denom!=0.) { a = (x-xy(2,1))/denom if(a>=0.&a<=1.) { yy = (xy(3,2)-xy(2,2))*a + xy(2,2) iptr=iptr+1;y(iptr)=yy ###;print(I(2),R(x),R(yy)) } } denom=(xy(1,1)-xy(3,1)) if(denom!=0.) { a = (x-xy(3,1))/denom if(a>=0.&a<=1.) { yy = (xy(1,2)-xy(3,2))*a + xy(3,2) iptr=iptr+1;y(iptr)=yy ###;print(I(3),R(x),R(yy)) } } if(iptr==2) { call segmtz(x,y(1),x,y(2),1) } else {call intpr("------error iptr==============",31,0,0)} } return end real function xkern(k,mv,iexp1,iexp2) xkern = (1.0 - (k/float(mv))**iexp1)**iexp2 return end