Add mathjax for dgesv

[maxima.git] / archive / src / plot3d.lisp

;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MAXIMA")
;;; PLOT3D Plotting Package (3D surfaces, Contours, etc.)

;ref: cacm algorithm 420 vol 15. (1972)

(declare-top(special maxdim $perspective $reverse scale-x scale-y max-xf min-xf $viewpt
                  $underside $howclose $crosshatch $xfun $yfun $zmax1 $zmin1
                  $zmax $zmin $labelcontours $plotnumprec $zigzag
                  x-3d y-3d xg-3d g-3d xh-3d h-3d x-arrv y-arrv z-arrv))

(setq  x-3d (make-array 1 :element-type 'long-float :adjustable t) 
       y-3d (make-array 1 :element-type 'long-float :adjustable t)
       xg-3d (make-array 2 :element-type 'long-float :adjustable t)
       g-3d (make-array 2 :element-type 'long-float :adjustable t :fill-pointer 2)
       xh-3d (make-array 1 :element-type 'long-float :adjustable t)
       h-3d (make-array 1 :element-type 'long-float :adjustable t))

(defun alter-array-size (aarray new-size)
  (setq aarray (if (< new-size (length aarray))
                  aarray
                  (apply 'adjust-array aarray (f+ new-size 2000.)
                         (if (array-has-fill-pointer-p aarray)
                             (list :fill-pointer (fill-pointer aarray))))))
  aarray)

(defun hide-init nil
  (let ((infin 1.e15))
    (setq xg-3d
          #-cl
          (make-array 2 :element-type 'long-float )
          #+cl (make-array 2 :element-type 'long-float :adjustable t)
           g-3d
           #-cl
           (zl-make-array 2 ':type '#+CADR art-float #-CADR 'art-q ':leader-list '(2))
           #+cl (make-array 2 :fill-pointer 2 :adjustable t )
           )
    (aset (-$ infin) xg-3d 0.) (aset infin xg-3d 1.)
    (aset (-$ infin) g-3d 0.) (aset (-$ infin) g-3d 1.)))

(defun howclose-line (mark)
  (let ((len (f- (car mark)
                 1 )) (xsta (cadr mark)) (ysta (caddr mark)))
    (setq mark (cddddr mark))
    (*$ 0.5 (+$ (funcall $howclose (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
              (funcall $howclose (aref x-arrv (f+ xsta (f* len (car mark))))
                       (aref y-arrv (f+ ysta (f* len (cadr mark)))) 0.0)))))

#+cl
(defun and-subr (&rest l)
  (sloop for v in l when (null v)
                     do (return v)
        finally (return v)))

;;; assumes all marks have the same x's and y's
(setq $zigzag nil)     ; Draws plots in zigzag pattern when CROSSHATCH:T$
(defun hide-drive (marks typel)
       (hide-init)
       (setq marks (copy-top-level marks ))
       (cond ((null (caar marks))
              (let ((marks1 nil)
                    (marks2 nil)
                    (border nil))
                (cond ((and $crosshatch $zigzag)
                       (setq marks1 (mapcar (function surf-expand2) marks))
                       (do ((sign 1.0 -1.0)) (nil)
                         (hide-init)
                         (do ((marks1 marks1 (mapcar (function cdr) marks1))
                              (i 0 (f1+ i)))
                             ((apply (function #+cl and-subr #-cl and)
                                     (mapcar (function null) marks1)))
                           (do ((marks1 marks1 (cdr marks1))
                                (typel1 typel (cdr typel1)))
                               ((null marks1))
                             (or typel1 (setq typel1 typel))
                             (let ((mark (caar marks1)) (type (car typel1)))
                               (and mark
                                    (cond ((> sign 0.0)
                                           (hide3d mark sign t type))
                                          (t (hide3d mark sign (> i 1) 0)))))))
                         (cond ((or (< sign 0.0) (not $underside))
                                (return nil)))))
                      (t
                (setq marks1 (mapcar (function (lambda (mark)
                                                 (surf-expand mark t)))
                                     marks)
                      marks2 (mapcar (function (lambda (mark)
                                                 (surf-expand mark nil)))
                                     marks))
                (cond ((> (howclose-line (caar marks1))
                          (howclose-line (car (last (car marks1)))))
                       (setq marks1 (mapcar (function nreverse) marks1))))
                (cond ((> (howclose-line (caar marks2))
                          (howclose-line (car (last (car marks2)))))
                       (setq marks2 (mapcar (function nreverse) marks2))))
                (setq border (list (caar marks1) (caar marks2)))
                (do ((xdir t nil)) (nil)
                    (do ((sign 1.0 -1.0) (ifplot t nil)) (nil)
                        (hide-init)
                        (cond (xdir (hide3d (cadr border) sign nil 0.))
                              (t (hide3d (car border) sign nil 0.)))
                        (do ((marks1 marks1 (mapcar (function cdr) marks1)))
                            ((apply (function #+cl and-subr #-cl and)
                                    (mapcar (function null) marks1)))
                            (do ((marks1 marks1 (cdr marks1))
                                 (typel1 typel (cdr typel1)))
                                ((null marks1))
                                (or typel1 (setq typel1 typel))
                                (and (caar marks1)
                                     (hide3d (caar marks1) sign ifplot
                                             (car typel1)))
                                (setq ifplot t)))
                        (cond ((or (< sign 0.0) (not $underside))
                               (return nil))))
                    (cond ((or (null xdir) (not $crosshatch)) (return nil)))
                    (setq marks1 marks2))))))
             (t (cond ((> (howclose-line (car marks))
                          (howclose-line (car (last marks))))
                       (setq marks (nreverse marks) typel (reverse typel))))
                (do ((sign 1.0 -1.0) (ifplot t nil)) (nil)
                    (do ((marks marks (cdr marks)) (typel1 typel (cdr typel)))
                        ((null marks))
                        (or typel1 (setq typel1 typel))
                        (hide3d (car marks) sign ifplot (car typel1))
                        (setq ifplot t))
                    (cond ((or (< sign 0.0) (not $underside)) (return nil))))))
;       (mapcar (function (lambda (arr) (alter-array-size (SYMBOL-VALUE arr) 1)))
;              '(x-3d y-3d xh-3d h-3d xg-3d g-3d))
       )

; Returns list of line descriptors (num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2)
; where num is the number of points, xs etc. give the starting point in
; the x-arr etc. arrays, dx1 and dx2 etc. give alternating steps to use through
; the arrays.
(declare-top(fixnum xlen ylen xstart ystart zstart xminc yinc zinc1 zinc2
                 i nummac num xs dx1 dx2 ys dy1 dy2 zs dz1 dz2
                 xl yl xj yj zi1 zi2)
         (flonum x0 x1 y0 y1 c0 c1 c2 c3 cmin))
(defun surf-expand2 (mark)
  (setq mark (reorientate (cdr mark)))
  (let ((xlen (f1- (car mark)))
        (ylen (f1- (cadr mark)))
        (xstart (cadddr mark))
        (ystart 0) (zstart 0) (xinc 0) (yinc 0) (zinc1 0) (zinc2 0))
    (setq mark (cddddr mark) ystart (car mark) zstart (cadr mark)
          xinc (caddr mark) yinc (cadddr mark) mark (cddddr mark)
          zinc1 (car mark) zinc2 (cadr mark))
    (cons (list (f1+ ylen) xstart ystart zstart 0 yinc zinc2)
     (cons (list (f1+ xlen) xstart ystart zstart xinc 0 zinc1)
      (append
       (do ((l) (i 0 (f1+ i)) (nummax (f1+ (f* 2 xlen)))
            (num 3 (min nummax (f+ 2 num)))
            (xs xstart) (dx1 xinc) (dx2 0)
            (ys (f+ ystart yinc) (f+ ys yinc)) (dy1 0) (dy2 (f- yinc))
            (zs (f+ zstart zinc2) (f+ zs zinc2)) (dz1 zinc1) (dz2 (f- zinc2)))
           ((= i ylen) (nreverse l))
         (setq l (cons (list nil num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2) l)))
       (do ((l) (i 1 (f1+ i)) (nummax (f1+ (f* 2 ylen)))
          (num 3 (min nummax (f+ 2 num)))
          (xs (f+ xstart (f* (f1- xlen) xinc)) (f- xs xinc)) (dx2 0) (dx1 xinc)
          (ys (f+ ystart (f* ylen yinc))) (dy2 (f- yinc)) (dy1 0)
          (zs (f+ zstart (f* (f1- xlen) zinc1) (f* ylen zinc2)) (f- zs zinc1))
          (dz2 (f- zinc2)) (dz1 zinc1))
         ((= i xlen) l)
       (setq l (cons (list nil num xs ys zs dx1 dy1 dz1 dx2 dy2 dz2) l))))))))

(defun reorientate (mark)
  ; Find closest corner and orientate
  (let ((xl (f1- (car mark)))
        (yl (f1- (cadr mark)))
        (xs (cadddr mark))
        (ys 0) (zs 0) (xj 0) (yj 0) (zi1 0) (zi2 0))
    (setq mark (cddddr mark) ys (car mark) zs (cadr mark)
          xj (caddr mark) yj (cadddr mark) mark (cddddr mark)
          zi1 (car mark) zi2 (cadr mark))
    (let ((x0 (aref x-arrv xs)) (x1 (aref x-arrv (f+ xs (f* xl xj))))
          (y0 (aref y-arrv ys)) (y1 (aref y-arrv (f+ ys (f* yl yj)))))
      (let ((c0 (funcall $howclose x0 y0 0.0))
            (c1 (funcall $howclose x1 y0 0.0))
            (c2 (funcall $howclose x1 y1 0.0))
            (c3 (funcall $howclose x0 y1 0.0))
            (cmin 0.0))
        (setq cmin (min c0 c1 c2 c3))
        (append (list (f1+ xl) (f1+ yl) 0)
                (cond ((= c0 cmin)
(list xs ys zs xj yj zi1 zi2))
                      ((= c1 cmin)
(list (f+ xs (f* xl xj)) ys (f+ zs (f* xl zi1)) (f- xj) yj (f- zi1) zi2))
                      ((= c2 cmin)
(list (f+ xs (f* xl xj)) (f+ ys (f* yl yj)) (f+ zs (f* xl zi1) (f* yl zi2))
      (f- xj) (f- yj) (f- zi1) (f- zi2)))
                      (t
(list xs (f+ ys (f* yl yj)) (f+ zs (f* yl zi2)) xj (f- yj) zi1 (f- zi2)))))))))

(defun lookupx (xx jj) 
  (do ((i (f1+ jj) (f1+ i)))
      ((not (< (aref x-3d i) xx))
       (cond ((= (aref x-3d i) xx) i)
             ((f1- i))))))

(defun lookupxg (xx jj) 
       (do ((i (f1+ jj) (f1+ i)))
           ((not (< (aref xg-3d i) xx)) (cond ((= (aref xg-3d i) xx) i) ((f1- i)))))) 

(defun enlarge-array nil 
  (setq maxdim (f+ 50. maxdim))
  (setq xh-3d (alter-array-size xh-3d (f1+ maxdim)))
  (setq h-3d (alter-array-size h-3d (f1+ maxdim)))
  T)

(defun f-intercept (xx xi yi xip1 yip1 &aux tem)
  (cond ((zerop (setq tem (-$ xip1 xi)))
         (+$ yi (*$ (-$ xx xi))))
        (t
         (+$ yi (*$ (-$ xx xi) (// (-$ yip1 yi)tem))))))

;; check over this stupid function!!

(defun hide3d (mark sign ifplot type)
 (let ((alt (null (car mark))))
  (and alt (setq mark (cdr mark)))
  (and (> (car mark) 1)
  (let ((eps (*$ 1.0e-5 (-$ max-xf min-xf)))
        (n1 (car mark))
        (ng (f1- #+lispm (array-leader g-3d 0)
                 #-lispm (fill-pointer g-3d)
                 ))
        (jj 0) (ig 0) (it 0) (igg 0) (itt 0) (indexg 0) (indext 0) (f1 0.0) (f2 0.0)
        (x1 0.0) (x2 0.0) (z1 0.0) (z2 0.0) (last) (maxdim -1))
    ($changedash (fixnum-remainder type 10.))
    (setq type (// type 10.))
    ;; Make them bigger than need be, this is expensive
    (setq x-3d (alter-array-size x-3d n1))
    (setq y-3d (alter-array-size y-3d n1))
    (cond (alt
                       (let ((xsta (cadr mark))
                             (ysta (caddr mark))
                             (zsta (cadddr mark))
                             (xinc1 0) (yinc1 0) (zinc1 0)
                             (xinc2 0) (yinc2 0) (zinc2 0))
                          (setq mark (cddddr mark) xinc1 (car mark)
                                yinc1 (cadr mark) zinc1 (caddr mark)
                                mark (cdddr mark) xinc2 (car mark)
                                yinc2 (cadr mark) zinc2 (caddr mark))
                          (let ((back
                                (> (call-x (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
                                   (call-x (aref x-arrv (f+ xsta
                                                     (f* (// n1 2) xinc1)
                                                     (f* (// (f1- n1) 2) xinc2)))
                                           (aref y-arrv (f+ ysta
                                                     (f* (// n1 2) yinc1)
                                                     (f* (// (f1- n1) 2) yinc2)))
                                           0.0)))
                                (dk 1) (ksta 0) (kend n1))
                            (and back (setq dk -1 ksta (f1- n1) kend -1))
                          (do ((kk ksta (f+ dk kk))
                               (i xsta (f+ xinc i)) (xinc xinc1 (f- qx xinc))
                               (qx (f+ xinc1 xinc2))
                               (j ysta (f+ yinc j)) (yinc yinc1 (f- qy yinc))
                               (qy (f+ yinc1 yinc2))
                               (k zsta (f+ zinc k)) (zinc zinc1 (f- qz zinc))
                               (qz (f+ zinc1 zinc2)))
                              ((= kk kend))
                            (aset (call-x (aref x-arrv i) (aref y-arrv j) (aref z-arrv k))
                                  x-3d kk)
                            (aset (*$ sign (call-y (aref x-arrv i) (aref y-arrv j)
                                                   (aref z-arrv k)))
                                  y-3d kk)))))
                      (t                                      
    (let ((xsta (cadr mark)) (ysta (caddr mark)) (zsta (cadddr mark))
          (xinc 0) (yinc 0) (zinc 0))
       (setq mark (cddddr mark) xinc (car mark) yinc (cadr mark)
             zinc (caddr mark))
       (cond ((> (call-x (aref x-arrv xsta) (aref y-arrv ysta) 0.0)
                 (call-x (aref x-arrv (f+ xsta (f* (f1- n1) xinc)))
                         (aref y-arrv (f+ ysta (f* (f1- n1) yinc))) 0.0))
              (setq xsta (f+ xsta (f* (f1- n1) xinc))
                    ysta (f+ ysta (f* (f1- n1) yinc))
                    zsta (f+ zsta (f* (f1- n1) zinc))
                    xinc (f- xinc) yinc (f- yinc) zinc (f- zinc))))
       (do ((kk 0. (f1+ kk)) (i xsta (f+ xinc i)) (j ysta (f+ yinc j))
            (k zsta (f+ zinc k)))
           ((= kk n1))
         (aset (call-x (aref x-arrv i) (aref y-arrv j) (aref z-arrv k)) x-3d kk)
         (aset (*$ sign (call-y (aref x-arrv i) (aref y-arrv j)
                               (aref z-arrv k))) y-3d kk)))))
    (setq n1 (f1- n1))
    (setq jj (lookupxg (aref x-3d 0.) 0.))
    (do nil ((< jj maxdim)) (enlarge-array))
    (do ((j 0. (f1+ j)))
        ((> j jj))
      (aset (aref xg-3d j) xh-3d j)
      (aset (aref g-3d j) h-3d j))
    (setq ig (f1+ jj))
    (aset (aref x-3d 0.) xh-3d ig)
    (aset (f-intercept (aref x-3d 0.) (aref xg-3d jj) (aref g-3d jj) (aref xg-3d ig)
                       (aref g-3d ig))
           h-3d ig)
    (setq indexg jj 
          indext 0. 
          z1 (aref x-3d 0.) 
          f1 (-$ (aref h-3d ig) (aref y-3d 0.)) 
          it 1. 
          jj ig)
    (cond ((< (aref h-3d ig) (aref y-3d 0.))
           (do nil ((< jj maxdim)) (enlarge-array))
           (setq jj (f1+ ig))
           (aset (aref y-3d 0.) h-3d jj)
           (aset (+$ z1 eps) xh-3d jj)))
    (setq x1 z1)
    (do ((zz 0.0) (k1 0.) (k2 0.) (n2 0.) (ngraph 0.)
         (relinc (// scale-x scale-y)))
        (nil)
      (do ((iwhich 0.) (slope 0.0))
          (last (setq z2 (aref x-3d n1) igg (lookupxg z2 indexg) itt (f1- n1))
                nil)
        (cond ((< (aref xg-3d ig) (aref x-3d it))
               (setq x2 (aref xg-3d ig) 
                     iwhich 1. 
                     f2 (-$ (aref g-3d ig)
                           (f-intercept x2 (aref x-3d (f1- it))
                                        (aref y-3d (f1- it)) (aref x-3d it)
                                        (aref y-3d it))) 
                     ig (f1+ ig)))
              (t (setq iwhich 0. 
                       x2 (aref x-3d it) 
                       f2 (-$ (f-intercept x2 (aref xg-3d (f1- ig))
                                          (aref g-3d (f1- ig)) (aref xg-3d ig)
                                          (aref g-3d ig))
                             (aref y-3d it)) 
                       it (f1+ it))))
        (and (> it n1) (setq last t))
        (cond ((or (and (> f1 0) (> f2 0))
                   (and (< f1 0) (< f2 0))
                   (and (= f1 0) (= f2 0)))
               (setq x1 x2 f1 f2))
              (t (setq slope (// (-$ f2 f1) (-$ x2 x1)) 
                       igg (f- ig 1. iwhich) 
                       itt (f+ it -2. iwhich))
                 (cond ((and (> (abs (*$ slope relinc)) 1.0e-6)
                             (> (-$ x2 x1) eps))
                        (setq z2 (-$ x1 (// f1 slope)))
                        (and (< (-$ z2 x1) eps) (setq z2 (+$ eps x1))))
                       (t (setq z2 x2)))
                 (return nil))))
      (setq zz (+$ z1 (*$ 0.01 (-$ z2 z1))) 
            k1 (lookupx zz indext) 
            k2 (lookupxg zz indexg))
      (cond ((> (cond ((= k1 n1) (aref x-3d k1))
                      (t (f-intercept zz (aref x-3d k1) (aref y-3d k1)
                                      (aref x-3d (f1+ k1)) (aref y-3d (f1+ k1)))))
                (f-intercept zz (aref xg-3d k2) (aref g-3d k2)
                             (aref xg-3d (f1+ k2)) (aref g-3d (f1+ k2))))
             (setq ngraph (f- itt indext -2.))
             (do nil
                 ((not (> (f+ jj ngraph -1.) maxdim)))
               (enlarge-array))
             (setq n2 jj)
             (do ((i (f1+ indext) (f1+ i)))
                 ((> i itt))
               (setq jj (f1+ jj))
               (aset (aref x-3d i) xh-3d jj)
               (aset (aref y-3d i) h-3d jj))
             (setq jj (f1+ jj))
             (aset z2 xh-3d jj)
             (aset (f-intercept z2
                                (aref x-3d itt) (aref y-3d itt)
                                (aref x-3d (f1+ itt)) (aref y-3d (f1+ itt)))
                    h-3d jj)
             (and ifplot (graph-hide n2 ngraph sign type)))
            (t (do nil
                   ((< (f+ jj igg (f- indexg)) maxdim))
                 (enlarge-array))
               (cond ((not (= indexg igg))
                      (do ((i (f1+ indexg) (f1+ i)))
                          ((> i igg))
                        (setq jj (f1+ jj))
                        (aset (aref xg-3d i) xh-3d jj)
                        (aset (aref g-3d i) h-3d jj))))
               (setq jj (f1+ jj))
               (aset z2 xh-3d jj)
               (aset (f-intercept z2 (aref x-3d itt) (aref y-3d itt)
                                  (aref x-3d (f1+ itt)) (aref y-3d (f1+ itt)))
                      h-3d jj)))
      (setq indext itt indexg igg)
      (and last (return nil))
      (setq x1 x2 f1 f2 z1 z2)
      (and (> it n1) (setq last t)))
    (do nil ((not (> (f+ jj 3. ng (f- igg)) maxdim))) (enlarge-array))
    (aset (+$ (aref xh-3d jj) eps) xh-3d (f1+ jj))
    (setq jj (f1+ jj))
    (aset (f-intercept (aref x-3d n1) (aref xg-3d igg) (aref g-3d igg)
                       (aref xg-3d (f1+ igg)) (aref g-3d (f1+ igg)))
           h-3d jj)
    (do ((j (f1+ igg) (f1+ j)))
        ((> j ng))
      (setq jj (f1+ jj))
      (aset (aref xg-3d j) xh-3d jj)
      (aset (aref g-3d j) h-3d jj))
    (setq xg-3d (alter-array-size xg-3d (f1+ jj)))
    (setq g-3d (alter-array-size g-3d (f1+ jj)))
    (do ((i 0. (f1+ i)) (j 0. (f1+ j)) (flg) (ox (*$ 2.0 (aref xh-3d 0.))))
        ((> i jj)
         (setq xg-3d (alter-array-size xg-3d j))
         (setq g-3d (alter-array-size g-3d j))
         nil)
      (cond ((not (> (aref xh-3d i) (+$ ox eps eps)))
             (setq ox (+$ ox eps))
             (cond (flg (setq j (f1- j))) (t (setq flg t))))
            (t (setq flg nil ox (aref xh-3d i))))
      (aset ox xg-3d j)
      (aset (aref h-3d i) g-3d j))))))

(defun graph-hide (n2 ngraph sign symtype) 
       (setq ngraph (f+ n2 ngraph) symtype (fixnum-remainder symtype 10.))
       ($setpoint (aref xh-3d n2) (*$ sign (aref h-3d n2)))
       (or (= symtype 0.)
           ($drawsymbol (aref xh-3d n2) (*$ sign (aref h-3d n2)) symtype))
       (do ((i (f1+ n2) (f1+ i)))
           ((= i ngraph))
           ($vector (aref xh-3d i) (*$ sign (aref h-3d i)))
           (or (= symtype 0.)
               ($drawsymbol (aref xh-3d i) (*$ sign (aref h-3d i)) symtype)))) 


;;ref: the computer journal vol 15 num 4 p 382 (1972)

(declare-top (special maxdim s zds ~n1 ~n2 ~xinc ~xstart ~yinc ~ystart ~zinc1 ~zinc2 ~zstart
                  ~symtype $diag xbd-cont cont-arr ybd-cont itg-cont))

(setq $diag t)

(setq xbd-cont (make-array 1 ':element-type '(signed-byte 16) :adjustable t)
      ybd-cont (make-array 1 ':element-type '(signed-byte 16) :adjustable t)
      cont-arr (make-array 1 :element-type 'long-float :adjustable t)
      itg-cont (make-array '(1 1) ':element-type '(signed-byte 16) :adjustable t))

(defun bdyp (i j) (or (= j 0.) (= j (f1- ~n2)) (= i 0) (= i (f1- ~n1))))

(defun phi-cont (i j) (aref z-arrv (f+ ~zstart (f* i ~zinc1) (f* j ~zinc2))))

(defun x-cont (i) (aref x-arrv (f+ ~xstart (f* i ~xinc))))

(defun y-cont (j) (aref y-arrv (f+ ~ystart (f* j ~yinc))))

(defun contour-set (contours cmin cmax)
 (cond (($listp contours)
        (adjust-array cont-arr (f1- (length contours)))
        (fillarray cont-arr (mapcar 'fmeval (cdr contours))))
       ((oldget contours 'array)
        (adjust-array cont-arr (cadr (arraydims contours)))
        (fillarray cont-arr contours))
       (t (let ((infin (^$ 8.0 42.)) (min  0.0) (max  0.0) (cnum  0.)
                (intflg  (eq contours '$integer)))
                   (or intflg (numberp contours) (setq contours 20.))
                   (setq min infin max (f- infin))
                   (cond ((not (and (numberp cmax) (numberp cmin)))
                          (do ((i 0. (f1+ i)) (zlen (cadr (arraydims z-arrv))) (pt))
                              ((= i zlen))
                              (setq pt (aref z-arrv i))
                              (and (> pt max) (setq max pt))
                              (and (< pt min) (setq min pt)))))
                   (and (numberp cmax) (setq max (float cmax)))
                   (and (numberp cmin) (setq min (float cmin)))
                   (cond (intflg
                          (setq max (float (fix max))
                                min (float (f- (fix (f- min)))))
                          (cond ((not (> max min)) (setq max (f+ 1.0 min))))
                          (setq cnum (fix (f+ 0.5 (f- max min)))))
                         (t (setq contours (fix contours))
                            (cond ((< max min)
                                   (setq max (prog2 nil min
                                                    (setq min max)))))
                            (cond ((not (> contours 0))
                                   (setq cnum (max 2. (f- contours))))
                                  (t (let ((ll (cdr (scale1 contours min max))))
                                              (setq min (cadr ll)
                                                    max (caddr ll)
                                                    cnum (fix (f+ (// (f- max min)
                                                                       (car ll))
                                                                  1.5))))))))
                   (adjust-array cont-arr cnum)
                   (setq $zmax1 max $zmin1 min)
                   (do ((i 0. (f1+ i)) (pt 0.0)) ((= i cnum))
                       (cond (intflg
                              (aset (f+ min (float i)) cont-arr i))
                             (t
                              (setq pt (// (float i)
                                            (float (f1- cnum))))
                              (aset (f+ (f* min (f- 1.0 pt))
                                         (f* max pt))
                                     cont-arr i))))))))

(defun contour-init (marks)
       (let ((mark (car marks)))
         (cond ((null (car mark))
                (setq ~n1 (cadr mark) ~n2 (caddr mark) mark (cddddr mark)
                      ~xstart (car mark) ~ystart (cadr mark) ~zstart (caddr mark)
                      ~xinc (cadddr mark) mark (cddddr mark) ~yinc (car mark)
                      ~zinc1 (cadr mark) ~zinc2 (caddr mark))
                (cdr marks))
               (t (setq ~n1 (car mark) ~n2 (length marks) ~xstart 0. ~ystart 0.
                        ~zstart 0. ~xinc 1. ~yinc 1. ~zinc1 1. ~zinc2 ~n1)
                  nil))))

(defun contour-drive (marks typel)
  (let ((~n1 0) (~n2 0) (~xstart 0) (~ystart 0) (~zstart 0)
        (~xinc 0) (~yinc 0) (~zinc1 0) (~zinc2 0)
        (n5 0) (ncn 0) (maxdim -1))
    (do ((typel1 typel (cdr typel1)) (type)) ((null marks))
      (setq marks (contour-init marks)
            n5 (f+ (f* 2. ~n1) (f* 2. ~n2) -3.))
      (cond ((null typel1) (setq typel1 typel)))
      (setq type (car typel1))
      (adjust-array xbd-cont n5)
      (adjust-array ybd-cont n5)
      (setq itg-cont (make-array `(,~n1 ,~n2)
                                    ':element-type '(signed-byte 16) :adjustable t))
      (do ((i 0 (f1+ i))) ((= i ~n2))
        (aset i ybd-cont i) (aset 0 xbd-cont i))
      (do ((i 1. (f1+ i))) ((= i ~n1))
        (aset (f1- ~n2) ybd-cont (f+ ~n2 i -1.))
        (aset i xbd-cont (f+ ~n2 i -1.)))
      (do ((i (f- ~n2 2) (f1- i))) ((< i 0.))
        (aset i ybd-cont (f- (f+ (f* 2 ~n2) ~n1 -3) i))
        (aset (f1- ~n1) xbd-cont (f- (f+ (f* 2 ~n2) ~n1 -3) i)))
      (do ((i (f- ~n1 2) (f1- i))) ((< i 0.))
        (aset 0. ybd-cont (f- (f+ (f* 2 ~n2) (f* 2 ~n1) -4.) i))
        (aset i xbd-cont (f- (f+ (f* 2 ~n2) (f* 2 ~n1) -4.) i)))
      (setq ncn (cadr (arraydims cont-arr)))
      ($changedash (fixnum-remainder type 10.))
      (setq type (fixnum-remainder (// type 10.) 10.))
      (enlarge-array)
      (contor ncn n5 type)
      (adjust-array xbd-cont 1)
      (adjust-array ybd-cont 1)
      (adjust-array xh-3d 1)
      (adjust-array h-3d 1)
      (adjust-array cont-arr 1)
      (setq itg-cont (make-array '(1 1)
                                    ':element-type '(signed-byte 16) :adjustable t)))))

(defun contor (ncn n5 ~symtype)
  (do ((cn 0. (f1+ cn)) (const) (i) (j) (ib) (jb))
      ((= cn ncn))
    (setq const (aref cont-arr cn))
    (do ((i 0. (f1+ i))) ((= i ~n1))
      (do ((j 0. (f1+ j))) ((= j ~n2))
        (aset 0. itg-cont i j)))
    (do ((k 1. (f1+ k))) ((= k n5))
      (setq i (aref xbd-cont k) j (aref ybd-cont k)
            ib (aref xbd-cont (f1- k)) jb (aref ybd-cont (f1- k)))
      (cond ((not (eq (< (f- (phi-cont i j) const) 0.0)
                      (< (f- (phi-cont ib jb) const) 0.0)))
             (look i j ib jb 1. const))))
    (do ((k 1. (f1+ k))) ((= k (f1- ~n1)))
      (do (( l 1. (f1+ l))) ((= l (f1- ~n2)))
        (setq i k ib k j l jb (f1- l))
        (cond ((and (not (and (bdyp i j) (bdyp ib jb)))
                    (not (eq (< (f- (phi-cont i j) const) 0.0)
                             (< (f- (phi-cont ib jb) const) 0.0))))
               (look i j ib jb 2. const)))))))

(defun look (i j ib jb qq const)
 (prog
  (jp ip jm im zds s ent) 
  (setq jp (f1+ j) ip (f1+ i) jm (f1- j) im (f1- i) zds 0.)
  (cond
   ((= jb jm) 
    (and (> (aref itg-cont i jm) 1.) (return nil))
    (setq ent 1.)
    (aset (x-cont i) xh-3d 0.)
    (aset (f-intercept const (phi-cont i jm) (y-cont jm) (phi-cont i j) (y-cont j))
           h-3d 0.))
   ((= ib im) 
    (and (or (= (aref itg-cont im j) 1.) (= (aref itg-cont im j) 3.)) (return nil))
    (setq ent 2.)
    (aset (y-cont j) h-3d 0.)
    (aset (f-intercept const (phi-cont im j) (x-cont im) (phi-cont i j) (x-cont i))
           xh-3d 0.))
   ((= jb jp) 
    (and (> (aref itg-cont i j) 1.) (return nil))
    (setq ent 3.)
    (aset (x-cont i) xh-3d 0.)
    (aset (f-intercept const (phi-cont i j) (y-cont j) (phi-cont i jp) (y-cont jp))
           h-3d 0.))
   (t (and (or (= (aref itg-cont i j) 1.) (= (aref itg-cont i j) 3.)) (return nil))
      (setq ent 4.)
      (aset (y-cont j) h-3d 0.)
      (aset (f-intercept const (phi-cont i j) (x-cont i)
                          (phi-cont ip j) (x-cont ip))
             xh-3d 0.)))
  (setq s 1.)
  (do
   nil (nil)
   (setq ip (f1+ i) jp (f1+ j) im (f1- i) jm (f1- j))
   (cond ((= ent 1.) (aset (f+ (aref itg-cont i jm) 2.) itg-cont i jm)
          (setq ent (ffnd i ip ip i j j jm jm ent qq const))
          (cond ((= ent 1.) (setq i ip)) ((= ent 2.) (setq i ip j jm))))
         ((= ent 2.) (aset (f1+ (aref itg-cont im j)) itg-cont im j)
          (setq ent (ffnd i i im im j jm jm j ent qq const))
          (cond ((= ent 2.) (setq j jm)) ((= ent 3.) (setq i im j jm))))
         ((= ent 3.) (aset (f+ (aref itg-cont i j) 2.) itg-cont i j)
          (setq ent (ffnd i im im i j j jp jp ent qq const))
          (cond ((= ent 3.) (setq i im)) ((= ent 4.) (setq i im j jp))))
         (t (aset (f1+ (aref itg-cont i j)) itg-cont i j)
            (setq ent (ffnd i i ip ip j jp jp j ent qq const))
            (cond ((= ent 4.) (setq j jp)) ((= ent 1.) (setq i ip j jp)))))
   (cond ((= zds 1.)
          (cond ((= ent 1.) (aset (f+ 2. (aref itg-cont i (f1- j))) itg-cont i (f1- j)))
                ((= ent 2.) (aset (f1+ (aref itg-cont (f1- i) j)) itg-cont (f1- i) j))
                ((= ent 3.) (aset (f+ 2 (aref itg-cont i j)) itg-cont i j))
                (t (aset (f1+ (aref itg-cont i j)) itg-cont i j)))
          (return nil)))
   (cond ((= qq 2.)
          (cond ((= ent 1.) (and (> (aref itg-cont i (f1- j)) 1.) (return nil)))
                ((= ent 2.) (and (oddp (aref itg-cont (f1- i) j)) (return nil)))
                ((= ent 3.) (and (> (aref itg-cont i j) 1.) (return nil)))
                (t (and (oddp (aref itg-cont i j)) (return nil)))))))
  (graph-contour s ~symtype const)))

;;pretend 0.0 phi's are +ve

(defun same-sign (phi1 phi2) (cond ((< phi1 0.0) (< phi2 0.0)) (t (not (< phi2 0.0)))))

(defun ffnd (i1 i2 i3 i4 j1 j2 j3 j4 ent qq const)
  (cond ((> (f+ s 4) maxdim) (enlarge-array)))
  (let ((phi1 (f- (phi-cont i1 j1) const)) (phi2 (f- (phi-cont i2 j2) const))
        (phi3 (f- (phi-cont i3 j3) const)) (phi4 (f- (phi-cont i4 j4) const))
        (phiav 0.0) (revflag nil) (xav (// (f+ (x-cont i1) (x-cont i3)) 2.0))
        (yav (// (f+ (y-cont j1) (y-cont j3)) 2.0)))
    (setq phiav (// (f+ phi1 phi2 phi3 phi4) 4.0))
    (cond ((not (same-sign phiav phi4))
           (setq revflag t
                 i1 (prog2 nil i4 (setq i4 i1))
                 j1 (prog2 nil j4 (setq j4 j1))
                 phi1 (prog2 nil phi4 (setq phi4 phi1))
                 i2 (prog2 nil i3 (setq i3 i2))
                 j2 (prog2 nil j3 (setq j3 j2))
                 phi2 (prog2 nil phi3 (setq phi3 phi2)))))
    (cond ($diag 
            (aset (f-intercept 0.0 phi1 (x-cont i1) phiav xav) xh-3d s)
            (aset (f-intercept 0.0 phi1 (y-cont j1) phiav yav) h-3d s)
            (setq s (f1+ s))))
    (do ((i 0. (f1+ i))) 
        ((not (same-sign phi1 phi2))
         (aset (f-intercept 0.0 phi1 (x-cont i1) phi2
                            (x-cont i2)) xh-3d s)
         (aset (f-intercept 0.0 phi1 (y-cont j1) phi2
                            (y-cont j2)) h-3d s)
         (setq s (f1+ s))
         (and (= qq 1.) (bdyp i1 j1) (bdyp i2 j2) (setq zds 1.))
         (and revflag (not (oddp i)) (setq i (f+ 2. i)))
         (f1+ (fixnum-remainder (f+ i ent 2.) 4.)))
      (cond ((and $diag (not (= phiav 0.0)))
             (aset (f-intercept 0.0 phi2 (x-cont i2) phiav
                                xav) xh-3d s)
             (aset (f-intercept 0.0 phi2 (y-cont j2) phiav yav) h-3d s)
             (setq s (f1+ s))))
      (setq i4 (prog2 nil i1 (setq i1 i2 i2 i3 i3 i4))
            j4 (prog2 nil j1 (setq j1 j2 j2 j3 j3 j4))
            phi4 (prog2 nil phi1 (setq phi1 phi2 phi2 phi3 phi3 phi4))))))

(defun graph-contour (ngraph symtype z)
       ($setpoint3 (aref xh-3d 0.) (aref h-3d 0.) z)
       (or (= symtype 0.)
           ($drawsymbol3 (aref xh-3d 0.) (aref h-3d 0.) z symtype))
       (do ((i 1. (f1+ i)))
           ((= i ngraph))
           ($vector3 (aref xh-3d i) (aref h-3d i) z)
           (or (= symtype 0.)
               ($drawsymbol3 (aref xh-3d i) (aref h-3d i) z symtype)))
       (and $labelcontours
            ($ghprint (format nil  `((f ,(fix $plotnumprec))) z)
                      (plot-x (call-x (aref xh-3d (// ngraph 2.))
                                     (aref h-3d (// ngraph 2.)) z))
                      (plot-y (call-y (aref xh-3d (// ngraph 2.))
                                     (aref h-3d (// ngraph 2.)) z))
                      1.)))

;;(declare (unspecial maxdim s zds ~n1 ~n2 scale-x scale-y max-xf min-xf
;;                  x-3d y-3d xg-3d g-3d xh-3d h-3d x-arrv y-arrv z-arrv))