archive/src/plot3d.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8
   9 (in-package "MAXIMA")
  10 ;;; PLOT3D Plotting Package (3D surfaces, Contours, etc.)
  11
  12 ;ref: cacm algorithm 420 vol 15. (1972)
  13
  14 (declare-top(special maxdim $perspective $reverse scale-x scale-y max-xf min-xf $viewpt
  15                   $underside $howclose $crosshatch $xfun $yfun $zmax1 $zmin1
  16                   $zmax $zmin $labelcontours $plotnumprec $zigzag
  17                   x-3d y-3d xg-3d g-3d xh-3d h-3d x-arrv y-arrv z-arrv))
  18
  19 (setq  x-3d (make-array 1 :element-type 'long-float :adjustable t)
  20        y-3d (make-array 1 :element-type 'long-float :adjustable t)
  21        xg-3d (make-array 2 :element-type 'long-float :adjustable t)
  22        g-3d (make-array 2 :element-type 'long-float :adjustable t :fill-pointer 2)
  23        xh-3d (make-array 1 :element-type 'long-float :adjustable t)
  24        h-3d (make-array 1 :element-type 'long-float :adjustable t))
  25
  26 (defun alter-array-size (aarray new-size)
  27   (setq aarray (if (< new-size (length aarray))
  28                   aarray
  29                   (apply 'adjust-array aarray (f+ new-size 2000.)
  30                          (if (array-has-fill-pointer-p aarray)
  31                              (list :fill-pointer (fill-pointer aarray))))))
  32   aarray)
  33
  34 (defun hide-init nil
  35   (let ((infin 1.e15))
  36     (setq xg-3d
  37           #-cl
  38           (make-array 2 :element-type 'long-float )
  39           #+cl (make-array 2 :element-type 'long-float :adjustable t)
  40            g-3d
  41            #-cl
  42            (zl-make-array 2 ':type '#+CADR art-float #-CADR 'art-q ':leader-list '(2))
  43            #+cl (make-array 2 :fill-pointer 2 :adjustable t )
  44            )
  45     (aset (-$ infin) xg-3d 0.) (aset infin xg-3d 1.)
  46     (aset (-$ infin) g-3d 0.) (aset (-$ infin) g-3d 1.)))
  47
  48 (defun howclose-line (mark)
  49   (let ((len (f- (car mark)
  50                  1 )) (xsta (cadr mark)) (ysta (caddr mark)))
  51     (setq mark (cddddr mark))
  52     (*$ 0.5 (+$ (funcall $howclose (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
  53               (funcall $howclose (aref x-arrv (f+ xsta (f* len (car mark))))
  54                        (aref y-arrv (f+ ysta (f* len (cadr mark)))) 0.0)))))
  55
  56 #+cl
  57 (defun and-subr (&rest l)
  58   (sloop for v in l when (null v)
  59                      do (return v)
  60         finally (return v)))
  61
  62 ;;; assumes all marks have the same x's and y's
  63 (setq $zigzag nil)     ; Draws plots in zigzag pattern when CROSSHATCH:T$
  64 (defun hide-drive (marks typel)
  65        (hide-init)
  66        (setq marks (copy-top-level marks ))
  67        (cond ((null (caar marks))
  68               (let ((marks1 nil)
  69                     (marks2 nil)
  70                     (border nil))
  71                 (cond ((and $crosshatch $zigzag)
  72                        (setq marks1 (mapcar (function surf-expand2) marks))
  73                        (do ((sign 1.0 -1.0)) (nil)
  74                          (hide-init)
  75                          (do ((marks1 marks1 (mapcar (function cdr) marks1))
  76                               (i 0 (f1+ i)))
  77                              ((apply (function #+cl and-subr #-cl and)
  78                                      (mapcar (function null) marks1)))
  79                            (do ((marks1 marks1 (cdr marks1))
  80                                 (typel1 typel (cdr typel1)))
  81                                ((null marks1))
  82                              (or typel1 (setq typel1 typel))
  83                              (let ((mark (caar marks1)) (type (car typel1)))
  84                                (and mark
  85                                     (cond ((> sign 0.0)
  86                                            (hide3d mark sign t type))
  87                                           (t (hide3d mark sign (> i 1) 0)))))))
  88                          (cond ((or (< sign 0.0) (not $underside))
  89                                 (return nil)))))
  90                       (t
  91                 (setq marks1 (mapcar (function (lambda (mark)
  92                                                  (surf-expand mark t)))
  93                                      marks)
  94                       marks2 (mapcar (function (lambda (mark)
  95                                                  (surf-expand mark nil)))
  96                                      marks))
  97                 (cond ((> (howclose-line (caar marks1))
  98                           (howclose-line (car (last (car marks1)))))
  99                        (setq marks1 (mapcar (function nreverse) marks1))))
 100                 (cond ((> (howclose-line (caar marks2))
 101                           (howclose-line (car (last (car marks2)))))
 102                        (setq marks2 (mapcar (function nreverse) marks2))))
 103                 (setq border (list (caar marks1) (caar marks2)))
 104                 (do ((xdir t nil)) (nil)
 105                     (do ((sign 1.0 -1.0) (ifplot t nil)) (nil)
 106                         (hide-init)
 107                         (cond (xdir (hide3d (cadr border) sign nil 0.))
 108                               (t (hide3d (car border) sign nil 0.)))
 109                         (do ((marks1 marks1 (mapcar (function cdr) marks1)))
 110                             ((apply (function #+cl and-subr #-cl and)
 111                                     (mapcar (function null) marks1)))
 112                             (do ((marks1 marks1 (cdr marks1))
 113                                  (typel1 typel (cdr typel1)))
 114                                 ((null marks1))
 115                                 (or typel1 (setq typel1 typel))
 116                                 (and (caar marks1)
 117                                      (hide3d (caar marks1) sign ifplot
 118                                              (car typel1)))
 119                                 (setq ifplot t)))
 120                         (cond ((or (< sign 0.0) (not $underside))
 121                                (return nil))))
 122                     (cond ((or (null xdir) (not $crosshatch)) (return nil)))
 123                     (setq marks1 marks2))))))
 124              (t (cond ((> (howclose-line (car marks))
 125                           (howclose-line (car (last marks))))
 126                        (setq marks (nreverse marks) typel (reverse typel))))
 127                 (do ((sign 1.0 -1.0) (ifplot t nil)) (nil)
 128                     (do ((marks marks (cdr marks)) (typel1 typel (cdr typel)))
 129                         ((null marks))
 130                         (or typel1 (setq typel1 typel))
 131                         (hide3d (car marks) sign ifplot (car typel1))
 132                         (setq ifplot t))
 133                     (cond ((or (< sign 0.0) (not $underside)) (return nil))))))
 134 ;       (mapcar (function (lambda (arr) (alter-array-size (SYMBOL-VALUE arr) 1)))
 135 ;              '(x-3d y-3d xh-3d h-3d xg-3d g-3d))
 136        )
 137
 138 ; Returns list of line descriptors (num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2)
 139 ; where num is the number of points, xs etc. give the starting point in
 140 ; the x-arr etc. arrays, dx1 and dx2 etc. give alternating steps to use through
 141 ; the arrays.
 142 (declare-top(fixnum xlen ylen xstart ystart zstart xminc yinc zinc1 zinc2
 143                  i nummac num xs dx1 dx2 ys dy1 dy2 zs dz1 dz2
 144                  xl yl xj yj zi1 zi2)
 145          (flonum x0 x1 y0 y1 c0 c1 c2 c3 cmin))
 146 (defun surf-expand2 (mark)
 147   (setq mark (reorientate (cdr mark)))
 148   (let ((xlen (f1- (car mark)))
 149         (ylen (f1- (cadr mark)))
 150         (xstart (cadddr mark))
 151         (ystart 0) (zstart 0) (xinc 0) (yinc 0) (zinc1 0) (zinc2 0))
 152     (setq mark (cddddr mark) ystart (car mark) zstart (cadr mark)
 153           xinc (caddr mark) yinc (cadddr mark) mark (cddddr mark)
 154           zinc1 (car mark) zinc2 (cadr mark))
 155     (cons (list (f1+ ylen) xstart ystart zstart 0 yinc zinc2)
 156      (cons (list (f1+ xlen) xstart ystart zstart xinc 0 zinc1)
 157       (append
 158        (do ((l) (i 0 (f1+ i)) (nummax (f1+ (f* 2 xlen)))
 159             (num 3 (min nummax (f+ 2 num)))
 160             (xs xstart) (dx1 xinc) (dx2 0)
 161             (ys (f+ ystart yinc) (f+ ys yinc)) (dy1 0) (dy2 (f- yinc))
 162             (zs (f+ zstart zinc2) (f+ zs zinc2)) (dz1 zinc1) (dz2 (f- zinc2)))
 163            ((= i ylen) (nreverse l))
 164          (setq l (cons (list nil num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2) l)))
 165        (do ((l) (i 1 (f1+ i)) (nummax (f1+ (f* 2 ylen)))
 166           (num 3 (min nummax (f+ 2 num)))
 167           (xs (f+ xstart (f* (f1- xlen) xinc)) (f- xs xinc)) (dx2 0) (dx1 xinc)
 168           (ys (f+ ystart (f* ylen yinc))) (dy2 (f- yinc)) (dy1 0)
 169           (zs (f+ zstart (f* (f1- xlen) zinc1) (f* ylen zinc2)) (f- zs zinc1))
 170           (dz2 (f- zinc2)) (dz1 zinc1))
 171          ((= i xlen) l)
 172        (setq l (cons (list nil num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2) l))))))))
 173
 174 (defun reorientate (mark)
 175   ; Find closest corner and orientate
 176   (let ((xl (f1- (car mark)))
 177         (yl (f1- (cadr mark)))
 178         (xs (cadddr mark))
 179         (ys 0) (zs 0) (xj 0) (yj 0) (zi1 0) (zi2 0))
 180     (setq mark (cddddr mark) ys (car mark) zs (cadr mark)
 181           xj (caddr mark) yj (cadddr mark) mark (cddddr mark)
 182           zi1 (car mark) zi2 (cadr mark))
 183     (let ((x0 (aref x-arrv xs)) (x1 (aref x-arrv (f+ xs (f* xl xj))))
 184           (y0 (aref y-arrv ys)) (y1 (aref y-arrv (f+ ys (f* yl yj)))))
 185       (let ((c0 (funcall $howclose x0 y0 0.0))
 186             (c1 (funcall $howclose x1 y0 0.0))
 187             (c2 (funcall $howclose x1 y1 0.0))
 188             (c3 (funcall $howclose x0 y1 0.0))
 189             (cmin 0.0))
 190         (setq cmin (min c0 c1 c2 c3))
 191         (append (list (f1+ xl) (f1+ yl) 0)
 192                 (cond ((= c0 cmin)
 193 (list xs ys zs xj yj zi1 zi2))
 194                       ((= c1 cmin)
 195 (list (f+ xs (f* xl xj)) ys (f+ zs (f* xl zi1)) (f- xj) yj (f- zi1) zi2))
 196                       ((= c2 cmin)
 197 (list (f+ xs (f* xl xj)) (f+ ys (f* yl yj)) (f+ zs (f* xl zi1) (f* yl zi2))
 198       (f- xj) (f- yj) (f- zi1) (f- zi2)))
 199                       (t
 200 (list xs (f+ ys (f* yl yj)) (f+ zs (f* yl zi2)) xj (f- yj) zi1 (f- zi2)))))))))
 201
 202 (defun lookupx (xx jj)
 203   (do ((i (f1+ jj) (f1+ i)))
 204       ((not (< (aref x-3d i) xx))
 205        (cond ((= (aref x-3d i) xx) i)
 206              ((f1- i))))))
 207
 208 (defun lookupxg (xx jj)
 209        (do ((i (f1+ jj) (f1+ i)))
 210            ((not (< (aref xg-3d i) xx)) (cond ((= (aref xg-3d i) xx) i) ((f1- i))))))
 211
 212 (defun enlarge-array nil
 213   (setq maxdim (f+ 50. maxdim))
 214   (setq xh-3d (alter-array-size xh-3d (f1+ maxdim)))
 215   (setq h-3d (alter-array-size h-3d (f1+ maxdim)))
 216   T)
 217
 218 (defun f-intercept (xx xi yi xip1 yip1 &aux tem)
 219   (cond ((zerop (setq tem (-$ xip1 xi)))
 220          (+$ yi (*$ (-$ xx xi))))
 221         (t
 222          (+$ yi (*$ (-$ xx xi) (// (-$ yip1 yi)tem))))))
 223
 224 ;; check over this stupid function!!
 225
 226 (defun hide3d (mark sign ifplot type)
 227  (let ((alt (null (car mark))))
 228   (and alt (setq mark (cdr mark)))
 229   (and (> (car mark) 1)
 230   (let ((eps (*$ 1.0e-5 (-$ max-xf min-xf)))
 231         (n1 (car mark))
 232         (ng (f1- #+lispm (array-leader g-3d 0)
 233                  #-lispm (fill-pointer g-3d)
 234                  ))
 235         (jj 0) (ig 0) (it 0) (igg 0) (itt 0) (indexg 0) (indext 0) (f1 0.0) (f2 0.0)
 236         (x1 0.0) (x2 0.0) (z1 0.0) (z2 0.0) (last) (maxdim -1))
 237     ($changedash (fixnum-remainder type 10.))
 238     (setq type (// type 10.))
 239     ;; Make them bigger than need be, this is expensive
 240     (setq x-3d (alter-array-size x-3d n1))
 241     (setq y-3d (alter-array-size y-3d n1))
 242     (cond (alt
 243                        (let ((xsta (cadr mark))
 244                              (ysta (caddr mark))
 245                              (zsta (cadddr mark))
 246                              (xinc1 0) (yinc1 0) (zinc1 0)
 247                              (xinc2 0) (yinc2 0) (zinc2 0))
 248                           (setq mark (cddddr mark) xinc1 (car mark)
 249                                 yinc1 (cadr mark) zinc1 (caddr mark)
 250                                 mark (cdddr mark) xinc2 (car mark)
 251                                 yinc2 (cadr mark) zinc2 (caddr mark))
 252                           (let ((back
 253                                 (> (call-x (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
 254                                    (call-x (aref x-arrv (f+ xsta
 255                                                      (f* (// n1 2) xinc1)
 256                                                      (f* (// (f1- n1) 2) xinc2)))
 257                                            (aref y-arrv (f+ ysta
 258                                                      (f* (// n1 2) yinc1)
 259                                                      (f* (// (f1- n1) 2) yinc2)))
 260                                            0.0)))
 261                                 (dk 1) (ksta 0) (kend n1))
 262                             (and back (setq dk -1 ksta (f1- n1) kend -1))
 263                           (do ((kk ksta (f+ dk kk))
 264                                (i xsta (f+ xinc i)) (xinc xinc1 (f- qx xinc))
 265                                (qx (f+ xinc1 xinc2))
 266                                (j ysta (f+ yinc j)) (yinc yinc1 (f- qy yinc))
 267                                (qy (f+ yinc1 yinc2))
 268                                (k zsta (f+ zinc k)) (zinc zinc1 (f- qz zinc))
 269                                (qz (f+ zinc1 zinc2)))
 270                               ((= kk kend))
 271                             (aset (call-x (aref x-arrv i) (aref y-arrv j) (aref z-arrv k))
 272                                   x-3d kk)
 273                             (aset (*$ sign (call-y (aref x-arrv i) (aref y-arrv j)
 274                                                    (aref z-arrv k)))
 275                                   y-3d kk)))))
 276                       (t
 277     (let ((xsta (cadr mark)) (ysta (caddr mark)) (zsta (cadddr mark))
 278           (xinc 0) (yinc 0) (zinc 0))
 279        (setq mark (cddddr mark) xinc (car mark) yinc (cadr mark)
 280              zinc (caddr mark))
 281        (cond ((> (call-x (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
 282                  (call-x (aref x-arrv (f+ xsta (f* (f1- n1) xinc)))
 283                          (aref y-arrv (f+ ysta (f* (f1- n1) yinc))) 0.0))
 284               (setq xsta (f+ xsta (f* (f1- n1) xinc))
 285                     ysta (f+ ysta (f* (f1- n1) yinc))
 286                     zsta (f+ zsta (f* (f1- n1) zinc))
 287                     xinc (f- xinc) yinc (f- yinc) zinc (f- zinc))))
 288        (do ((kk 0. (f1+ kk)) (i xsta (f+ xinc i)) (j ysta (f+ yinc j))
 289             (k zsta (f+ zinc k)))
 290            ((= kk n1))
 291          (aset (call-x (aref x-arrv i) (aref y-arrv j) (aref z-arrv k)) x-3d kk)
 292          (aset (*$ sign (call-y (aref x-arrv i) (aref y-arrv j)
 293                                (aref z-arrv k))) y-3d kk)))))
 294     (setq n1 (f1- n1))
 295     (setq jj (lookupxg (aref x-3d 0.) 0.))
 296     (do nil ((< jj maxdim)) (enlarge-array))
 297     (do ((j 0. (f1+ j)))
 298         ((> j jj))
 299       (aset (aref xg-3d j) xh-3d j)
 300       (aset (aref g-3d j) h-3d j))
 301     (setq ig (f1+ jj))
 302     (aset (aref x-3d 0.) xh-3d ig)
 303     (aset (f-intercept (aref x-3d 0.) (aref xg-3d jj) (aref g-3d jj) (aref xg-3d ig)
 304                        (aref g-3d ig))
 305            h-3d ig)
 306     (setq indexg jj
 307           indext 0.
 308           z1 (aref x-3d 0.)
 309           f1 (-$ (aref h-3d ig) (aref y-3d 0.))
 310           it 1.
 311           jj ig)
 312     (cond ((< (aref h-3d ig) (aref y-3d 0.))
 313            (do nil ((< jj maxdim)) (enlarge-array))
 314            (setq jj (f1+ ig))
 315            (aset (aref y-3d 0.) h-3d jj)
 316            (aset (+$ z1 eps) xh-3d jj)))
 317     (setq x1 z1)
 318     (do ((zz 0.0) (k1 0.) (k2 0.) (n2 0.) (ngraph 0.)
 319          (relinc (// scale-x scale-y)))
 320         (nil)
 321       (do ((iwhich 0.) (slope 0.0))
 322           (last (setq z2 (aref x-3d n1) igg (lookupxg z2 indexg) itt (f1- n1))
 323                 nil)
 324         (cond ((< (aref xg-3d ig) (aref x-3d it))
 325                (setq x2 (aref xg-3d ig)
 326                      iwhich 1.
 327                      f2 (-$ (aref g-3d ig)
 328                            (f-intercept x2 (aref x-3d (f1- it))
 329                                         (aref y-3d (f1- it)) (aref x-3d it)
 330                                         (aref y-3d it)))
 331                      ig (f1+ ig)))
 332               (t (setq iwhich 0.
 333                        x2 (aref x-3d it)
 334                        f2 (-$ (f-intercept x2 (aref xg-3d (f1- ig))
 335                                           (aref g-3d (f1- ig)) (aref xg-3d ig)
 336                                           (aref g-3d ig))
 337                              (aref y-3d it))
 338                        it (f1+ it))))
 339         (and (> it n1) (setq last t))
 340         (cond ((or (and (> f1 0) (> f2 0))
 341                    (and (< f1 0) (< f2 0))
 342                    (and (= f1 0) (= f2 0)))
 343                (setq x1 x2 f1 f2))
 344               (t (setq slope (// (-$ f2 f1) (-$ x2 x1))
 345                        igg (f- ig 1. iwhich)
 346                        itt (f+ it -2. iwhich))
 347                  (cond ((and (> (abs (*$ slope relinc)) 1.0e-6)
 348                              (> (-$ x2 x1) eps))
 349                         (setq z2 (-$ x1 (// f1 slope)))
 350                         (and (< (-$ z2 x1) eps) (setq z2 (+$ eps x1))))
 351                        (t (setq z2 x2)))
 352                  (return nil))))
 353       (setq zz (+$ z1 (*$ 0.01 (-$ z2 z1)))
 354             k1 (lookupx zz indext)
 355             k2 (lookupxg zz indexg))
 356       (cond ((> (cond ((= k1 n1) (aref x-3d k1))
 357                       (t (f-intercept zz (aref x-3d k1) (aref y-3d k1)
 358                                       (aref x-3d (f1+ k1)) (aref y-3d (f1+ k1)))))
 359                 (f-intercept zz (aref xg-3d k2) (aref g-3d k2)
 360                              (aref xg-3d (f1+ k2)) (aref g-3d (f1+ k2))))
 361              (setq ngraph (f- itt indext -2.))
 362              (do nil
 363                  ((not (> (f+ jj ngraph -1.) maxdim)))
 364                (enlarge-array))
 365              (setq n2 jj)
 366              (do ((i (f1+ indext) (f1+ i)))
 367                  ((> i itt))
 368                (setq jj (f1+ jj))
 369                (aset (aref x-3d i) xh-3d jj)
 370                (aset (aref y-3d i) h-3d jj))
 371              (setq jj (f1+ jj))
 372              (aset z2 xh-3d jj)
 373              (aset (f-intercept z2
 374                                 (aref x-3d itt) (aref y-3d itt)
 375                                 (aref x-3d (f1+ itt)) (aref y-3d (f1+ itt)))
 376                     h-3d jj)
 377              (and ifplot (graph-hide n2 ngraph sign type)))
 378             (t (do nil
 379                    ((< (f+ jj igg (f- indexg)) maxdim))
 380                  (enlarge-array))
 381                (cond ((not (= indexg igg))
 382                       (do ((i (f1+ indexg) (f1+ i)))
 383                           ((> i igg))
 384                         (setq jj (f1+ jj))
 385                         (aset (aref xg-3d i) xh-3d jj)
 386                         (aset (aref g-3d i) h-3d jj))))
 387                (setq jj (f1+ jj))
 388                (aset z2 xh-3d jj)
 389                (aset (f-intercept z2 (aref x-3d itt) (aref y-3d itt)
 390                                   (aref x-3d (f1+ itt)) (aref y-3d (f1+ itt)))
 391                       h-3d jj)))
 392       (setq indext itt indexg igg)
 393       (and last (return nil))
 394       (setq x1 x2 f1 f2 z1 z2)
 395       (and (> it n1) (setq last t)))
 396     (do nil ((not (> (f+ jj 3. ng (f- igg)) maxdim))) (enlarge-array))
 397     (aset (+$ (aref xh-3d jj) eps) xh-3d (f1+ jj))
 398     (setq jj (f1+ jj))
 399     (aset (f-intercept (aref x-3d n1) (aref xg-3d igg) (aref g-3d igg)
 400                        (aref xg-3d (f1+ igg)) (aref g-3d (f1+ igg)))
 401            h-3d jj)
 402     (do ((j (f1+ igg) (f1+ j)))
 403         ((> j ng))
 404       (setq jj (f1+ jj))
 405       (aset (aref xg-3d j) xh-3d jj)
 406       (aset (aref g-3d j) h-3d jj))
 407     (setq xg-3d (alter-array-size xg-3d (f1+ jj)))
 408     (setq g-3d (alter-array-size g-3d (f1+ jj)))
 409     (do ((i 0. (f1+ i)) (j 0. (f1+ j)) (flg) (ox (*$ 2.0 (aref xh-3d 0.))))
 410         ((> i jj)
 411          (setq xg-3d (alter-array-size xg-3d j))
 412          (setq g-3d (alter-array-size g-3d j))
 413          nil)
 414       (cond ((not (> (aref xh-3d i) (+$ ox eps eps)))
 415              (setq ox (+$ ox eps))
 416              (cond (flg (setq j (f1- j))) (t (setq flg t))))
 417             (t (setq flg nil ox (aref xh-3d i))))
 418       (aset ox xg-3d j)
 419       (aset (aref h-3d i) g-3d j))))))
 420
 421 (defun graph-hide (n2 ngraph sign symtype)
 422        (setq ngraph (f+ n2 ngraph) symtype (fixnum-remainder symtype 10.))
 423        ($setpoint (aref xh-3d n2) (*$ sign (aref h-3d n2)))
 424        (or (= symtype 0.)
 425            ($drawsymbol (aref xh-3d n2) (*$ sign (aref h-3d n2)) symtype))
 426        (do ((i (f1+ n2) (f1+ i)))
 427            ((= i ngraph))
 428            ($vector (aref xh-3d i) (*$ sign (aref h-3d i)))
 429            (or (= symtype 0.)
 430                ($drawsymbol (aref xh-3d i) (*$ sign (aref h-3d i)) symtype))))
 431
 432
 433 ;;ref: the computer journal vol 15 num 4 p 382 (1972)
 434
 435 (declare-top (special maxdim s zds ~n1 ~n2 ~xinc ~xstart ~yinc ~ystart ~zinc1 ~zinc2 ~zstart
 436                   ~symtype $diag xbd-cont cont-arr ybd-cont itg-cont))
 437
 438 (setq $diag t)
 439
 440 (setq xbd-cont (make-array 1 ':element-type '(signed-byte 16) :adjustable t)
 441       ybd-cont (make-array 1 ':element-type '(signed-byte 16) :adjustable t)
 442       cont-arr (make-array 1 :element-type 'long-float :adjustable t)
 443       itg-cont (make-array '(1 1) ':element-type '(signed-byte 16) :adjustable t))
 444
 445 (defun bdyp (i j) (or (= j 0.) (= j (f1- ~n2)) (= i 0) (= i (f1- ~n1))))
 446
 447 (defun phi-cont (i j) (aref z-arrv (f+ ~zstart (f* i ~zinc1) (f* j ~zinc2))))
 448
 449 (defun x-cont (i) (aref x-arrv (f+ ~xstart (f* i ~xinc))))
 450
 451 (defun y-cont (j) (aref y-arrv (f+ ~ystart (f* j ~yinc))))
 452
 453 (defun contour-set (contours cmin cmax)
 454  (cond (($listp contours)
 455         (adjust-array cont-arr (f1- (length contours)))
 456         (fillarray cont-arr (mapcar 'fmeval (cdr contours))))
 457        ((oldget contours 'array)
 458         (adjust-array cont-arr (cadr (arraydims contours)))
 459         (fillarray cont-arr contours))
 460        (t (let ((infin (^$ 8.0 42.)) (min  0.0) (max  0.0) (cnum  0.)
 461                 (intflg  (eq contours '$integer)))
 462                    (or intflg (numberp contours) (setq contours 20.))
 463                    (setq min infin max (f- infin))
 464                    (cond ((not (and (numberp cmax) (numberp cmin)))
 465                           (do ((i 0. (f1+ i)) (zlen (cadr (arraydims z-arrv))) (pt))
 466                               ((= i zlen))
 467                               (setq pt (aref z-arrv i))
 468                               (and (> pt max) (setq max pt))
 469                               (and (< pt min) (setq min pt)))))
 470                    (and (numberp cmax) (setq max (float cmax)))
 471                    (and (numberp cmin) (setq min (float cmin)))
 472                    (cond (intflg
 473                           (setq max (float (fix max))
 474                                 min (float (f- (fix (f- min)))))
 475                           (cond ((not (> max min)) (setq max (f+ 1.0 min))))
 476                           (setq cnum (fix (f+ 0.5 (f- max min)))))
 477                          (t (setq contours (fix contours))
 478                             (cond ((< max min)
 479                                    (setq max (prog2 nil min
 480                                                     (setq min max)))))
 481                             (cond ((not (> contours 0))
 482                                    (setq cnum (max 2. (f- contours))))
 483                                   (t (let ((ll (cdr (scale1 contours min max))))
 484                                               (setq min (cadr ll)
 485                                                     max (caddr ll)
 486                                                     cnum (fix (f+ (// (f- max min)
 487                                                                        (car ll))
 488                                                                   1.5))))))))
 489                    (adjust-array cont-arr cnum)
 490                    (setq $zmax1 max $zmin1 min)
 491                    (do ((i 0. (f1+ i)) (pt 0.0)) ((= i cnum))
 492                        (cond (intflg
 493                               (aset (f+ min (float i)) cont-arr i))
 494                              (t
 495                               (setq pt (// (float i)
 496                                             (float (f1- cnum))))
 497                               (aset (f+ (f* min (f- 1.0 pt))
 498                                          (f* max pt))
 499                                      cont-arr i))))))))
 500
 501 (defun contour-init (marks)
 502        (let ((mark (car marks)))
 503          (cond ((null (car mark))
 504                 (setq ~n1 (cadr mark) ~n2 (caddr mark) mark (cddddr mark)
 505                       ~xstart (car mark) ~ystart (cadr mark) ~zstart (caddr mark)
 506                       ~xinc (cadddr mark) mark (cddddr mark) ~yinc (car mark)
 507                       ~zinc1 (cadr mark) ~zinc2 (caddr mark))
 508                 (cdr marks))
 509                (t (setq ~n1 (car mark) ~n2 (length marks) ~xstart 0. ~ystart 0.
 510                         ~zstart 0. ~xinc 1. ~yinc 1. ~zinc1 1. ~zinc2 ~n1)
 511                   nil))))
 512
 513 (defun contour-drive (marks typel)
 514   (let ((~n1 0) (~n2 0) (~xstart 0) (~ystart 0) (~zstart 0)
 515         (~xinc 0) (~yinc 0) (~zinc1 0) (~zinc2 0)
 516         (n5 0) (ncn 0) (maxdim -1))
 517     (do ((typel1 typel (cdr typel1)) (type)) ((null marks))
 518       (setq marks (contour-init marks)
 519             n5 (f+ (f* 2. ~n1) (f* 2. ~n2) -3.))
 520       (cond ((null typel1) (setq typel1 typel)))
 521       (setq type (car typel1))
 522       (adjust-array xbd-cont n5)
 523       (adjust-array ybd-cont n5)
 524       (setq itg-cont (make-array `(,~n1 ,~n2)
 525                                     ':element-type '(signed-byte 16) :adjustable t))
 526       (do ((i 0 (f1+ i))) ((= i ~n2))
 527         (aset i ybd-cont i) (aset 0 xbd-cont i))
 528       (do ((i 1. (f1+ i))) ((= i ~n1))
 529         (aset (f1- ~n2) ybd-cont (f+ ~n2 i -1.))
 530         (aset i xbd-cont (f+ ~n2 i -1.)))
 531       (do ((i (f- ~n2 2) (f1- i))) ((< i 0.))
 532         (aset i ybd-cont (f- (f+ (f* 2 ~n2) ~n1 -3) i))
 533         (aset (f1- ~n1) xbd-cont (f- (f+ (f* 2 ~n2) ~n1 -3) i)))
 534       (do ((i (f- ~n1 2) (f1- i))) ((< i 0.))
 535         (aset 0. ybd-cont (f- (f+ (f* 2 ~n2) (f* 2 ~n1) -4.) i))
 536         (aset i xbd-cont (f- (f+ (f* 2 ~n2) (f* 2 ~n1) -4.) i)))
 537       (setq ncn (cadr (arraydims cont-arr)))
 538       ($changedash (fixnum-remainder type 10.))
 539       (setq type (fixnum-remainder (// type 10.) 10.))
 540       (enlarge-array)
 541       (contor ncn n5 type)
 542       (adjust-array xbd-cont 1)
 543       (adjust-array ybd-cont 1)
 544       (adjust-array xh-3d 1)
 545       (adjust-array h-3d 1)
 546       (adjust-array cont-arr 1)
 547       (setq itg-cont (make-array '(1 1)
 548                                     ':element-type '(signed-byte 16) :adjustable t)))))
 549
 550 (defun contor (ncn n5 ~symtype)
 551   (do ((cn 0. (f1+ cn)) (const) (i) (j) (ib) (jb))
 552       ((= cn ncn))
 553     (setq const (aref cont-arr cn))
 554     (do ((i 0. (f1+ i))) ((= i ~n1))
 555       (do ((j 0. (f1+ j))) ((= j ~n2))
 556         (aset 0. itg-cont i j)))
 557     (do ((k 1. (f1+ k))) ((= k n5))
 558       (setq i (aref xbd-cont k) j (aref ybd-cont k)
 559             ib (aref xbd-cont (f1- k)) jb (aref ybd-cont (f1- k)))
 560       (cond ((not (eq (< (f- (phi-cont i j) const) 0.0)
 561                       (< (f- (phi-cont ib jb) const) 0.0)))
 562              (look i j ib jb 1. const))))
 563     (do ((k 1. (f1+ k))) ((= k (f1- ~n1)))
 564       (do (( l 1. (f1+ l))) ((= l (f1- ~n2)))
 565         (setq i k ib k j l jb (f1- l))
 566         (cond ((and (not (and (bdyp i j) (bdyp ib jb)))
 567                     (not (eq (< (f- (phi-cont i j) const) 0.0)
 568                              (< (f- (phi-cont ib jb) const) 0.0))))
 569                (look i j ib jb 2. const)))))))
 570
 571 (defun look (i j ib jb qq const)
 572  (prog
 573   (jp ip jm im zds s ent)
 574   (setq jp (f1+ j) ip (f1+ i) jm (f1- j) im (f1- i) zds 0.)
 575   (cond
 576    ((= jb jm)
 577     (and (> (aref itg-cont i jm) 1.) (return nil))
 578     (setq ent 1.)
 579     (aset (x-cont i) xh-3d 0.)
 580     (aset (f-intercept const (phi-cont i jm) (y-cont jm) (phi-cont i j) (y-cont j))
 581            h-3d 0.))
 582    ((= ib im)
 583     (and (or (= (aref itg-cont im j) 1.) (= (aref itg-cont im j) 3.)) (return nil))
 584     (setq ent 2.)
 585     (aset (y-cont j) h-3d 0.)
 586     (aset (f-intercept const (phi-cont im j) (x-cont im) (phi-cont i j) (x-cont i))
 587            xh-3d 0.))
 588    ((= jb jp)
 589     (and (> (aref itg-cont i j) 1.) (return nil))
 590     (setq ent 3.)
 591     (aset (x-cont i) xh-3d 0.)
 592     (aset (f-intercept const (phi-cont i j) (y-cont j) (phi-cont i jp) (y-cont jp))
 593            h-3d 0.))
 594    (t (and (or (= (aref itg-cont i j) 1.) (= (aref itg-cont i j) 3.)) (return nil))
 595       (setq ent 4.)
 596       (aset (y-cont j) h-3d 0.)
 597       (aset (f-intercept const (phi-cont i j) (x-cont i)
 598                           (phi-cont ip j) (x-cont ip))
 599              xh-3d 0.)))
 600   (setq s 1.)
 601   (do
 602    nil (nil)
 603    (setq ip (f1+ i) jp (f1+ j) im (f1- i) jm (f1- j))
 604    (cond ((= ent 1.) (aset (f+ (aref itg-cont i jm) 2.) itg-cont i jm)
 605           (setq ent (ffnd i ip ip i j j jm jm ent qq const))
 606           (cond ((= ent 1.) (setq i ip)) ((= ent 2.) (setq i ip j jm))))
 607          ((= ent 2.) (aset (f1+ (aref itg-cont im j)) itg-cont im j)
 608           (setq ent (ffnd i i im im j jm jm j ent qq const))
 609           (cond ((= ent 2.) (setq j jm)) ((= ent 3.) (setq i im j jm))))
 610          ((= ent 3.) (aset (f+ (aref itg-cont i j) 2.) itg-cont i j)
 611           (setq ent (ffnd i im im i j j jp jp ent qq const))
 612           (cond ((= ent 3.) (setq i im)) ((= ent 4.) (setq i im j jp))))
 613          (t (aset (f1+ (aref itg-cont i j)) itg-cont i j)
 614             (setq ent (ffnd i i ip ip j jp jp j ent qq const))
 615             (cond ((= ent 4.) (setq j jp)) ((= ent 1.) (setq i ip j jp)))))
 616    (cond ((= zds 1.)
 617           (cond ((= ent 1.) (aset (f+ 2. (aref itg-cont i (f1- j))) itg-cont i (f1- j)))
 618                 ((= ent 2.) (aset (f1+ (aref itg-cont (f1- i) j)) itg-cont (f1- i) j))
 619                 ((= ent 3.) (aset (f+ 2 (aref itg-cont i j)) itg-cont i j))
 620                 (t (aset (f1+ (aref itg-cont i j)) itg-cont i j)))
 621           (return nil)))
 622    (cond ((= qq 2.)
 623           (cond ((= ent 1.) (and (> (aref itg-cont i (f1- j)) 1.) (return nil)))
 624                 ((= ent 2.) (and (oddp (aref itg-cont (f1- i) j)) (return nil)))
 625                 ((= ent 3.) (and (> (aref itg-cont i j) 1.) (return nil)))
 626                 (t (and (oddp (aref itg-cont i j)) (return nil)))))))
 627   (graph-contour s ~symtype const)))
 628
 629 ;;pretend 0.0 phi's are +ve
 630
 631 (defun same-sign (phi1 phi2) (cond ((< phi1 0.0) (< phi2 0.0)) (t (not (< phi2 0.0)))))
 632
 633 (defun ffnd (i1 i2 i3 i4 j1 j2 j3 j4 ent qq const)
 634   (cond ((> (f+ s 4) maxdim) (enlarge-array)))
 635   (let ((phi1 (f- (phi-cont i1 j1) const)) (phi2 (f- (phi-cont i2 j2) const))
 636         (phi3 (f- (phi-cont i3 j3) const)) (phi4 (f- (phi-cont i4 j4) const))
 637         (phiav 0.0) (revflag nil) (xav (// (f+ (x-cont i1) (x-cont i3)) 2.0))
 638         (yav (// (f+ (y-cont j1) (y-cont j3)) 2.0)))
 639     (setq phiav (// (f+ phi1 phi2 phi3 phi4) 4.0))
 640     (cond ((not (same-sign phiav phi4))
 641            (setq revflag t
 642                  i1 (prog2 nil i4 (setq i4 i1))
 643                  j1 (prog2 nil j4 (setq j4 j1))
 644                  phi1 (prog2 nil phi4 (setq phi4 phi1))
 645                  i2 (prog2 nil i3 (setq i3 i2))
 646                  j2 (prog2 nil j3 (setq j3 j2))
 647                  phi2 (prog2 nil phi3 (setq phi3 phi2)))))
 648     (cond ($diag
 649             (aset (f-intercept 0.0 phi1 (x-cont i1) phiav xav) xh-3d s)
 650             (aset (f-intercept 0.0 phi1 (y-cont j1) phiav yav) h-3d s)
 651             (setq s (f1+ s))))
 652     (do ((i 0. (f1+ i)))
 653         ((not (same-sign phi1 phi2))
 654          (aset (f-intercept 0.0 phi1 (x-cont i1) phi2
 655                             (x-cont i2)) xh-3d s)
 656          (aset (f-intercept 0.0 phi1 (y-cont j1) phi2
 657                             (y-cont j2)) h-3d s)
 658          (setq s (f1+ s))
 659          (and (= qq 1.) (bdyp i1 j1) (bdyp i2 j2) (setq zds 1.))
 660          (and revflag (not (oddp i)) (setq i (f+ 2. i)))
 661          (f1+ (fixnum-remainder (f+ i ent 2.) 4.)))
 662       (cond ((and $diag (not (= phiav 0.0)))
 663              (aset (f-intercept 0.0 phi2 (x-cont i2) phiav
 664                                 xav) xh-3d s)
 665              (aset (f-intercept 0.0 phi2 (y-cont j2) phiav yav) h-3d s)
 666              (setq s (f1+ s))))
 667       (setq i4 (prog2 nil i1 (setq i1 i2 i2 i3 i3 i4))
 668             j4 (prog2 nil j1 (setq j1 j2 j2 j3 j3 j4))
 669             phi4 (prog2 nil phi1 (setq phi1 phi2 phi2 phi3 phi3 phi4))))))
 670
 671 (defun graph-contour (ngraph symtype z)
 672        ($setpoint3 (aref xh-3d 0.) (aref h-3d 0.) z)
 673        (or (= symtype 0.)
 674            ($drawsymbol3 (aref xh-3d 0.) (aref h-3d 0.) z symtype))
 675        (do ((i 1. (f1+ i)))
 676            ((= i ngraph))
 677            ($vector3 (aref xh-3d i) (aref h-3d i) z)
 678            (or (= symtype 0.)
 679                ($drawsymbol3 (aref xh-3d i) (aref h-3d i) z symtype)))
 680        (and $labelcontours
 681             ($ghprint (format nil  `((f ,(fix $plotnumprec))) z)
 682                       (plot-x (call-x (aref xh-3d (// ngraph 2.))
 683                                      (aref h-3d (// ngraph 2.)) z))
 684                       (plot-y (call-y (aref xh-3d (// ngraph 2.))
 685                                      (aref h-3d (// ngraph 2.)) z))
 686                       1.)))
 687
 688 ;;(declare (unspecial maxdim s zds ~n1 ~n2 scale-x scale-y max-xf min-xf
 689 ;;                  x-3d y-3d xg-3d g-3d xh-3d h-3d x-arrv y-arrv z-arrv))