github.com/PDP-10.its
1
0
Fork 0
mirror of https://github.com/PDP-10/its.git synced 2026-01-16 00:14:18 +00:00
Code Issues Releases Wiki Activity
PDP-10.its/src/share/plot3d.66
Eric Swenson 85994ed770 Added files to support building and running Macsyma. 
Resolves #284.

Commented out uses of time-origin in maxtul; mcldmp (init) until we
can figure out why it gives arithmetic overflows under the emulators.

Updated the expect script statements in build_macsyma_portion to not
attempt to match expected strings, but simply sleep for some time
since in some cases the matching appears not to work.
2018-03-11 13:10:19 -07:00

808 lines
30 KiB
Common Lisp
Raw Blame History

;;;-*-Lisp-*-
; (c) Copyright 1976, 1983 Massachusetts Institute of Technology.
; All Rights Reserved. Enhancements (c) Copyright 1983 Symbolics Inc.
; All Rights Reserved.
;
; The data and information in the Enhancements is proprietary to, and a
; valuable trade secret of, SYMBOLICS, INC., a Delaware corporation. It is
; given in confidence by SYMBOLICS, and may not be used as the basis of
; manufacture, or be reproduced or copied, or distributed to any other
; party, in whole or in part, without the prior written consent of SYMBOLICS.
;ref: cacm algorithm 420 vol 15. (1972)
(declare (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)
(flonum eps z f1 f2 x1 x2 z1 z2 xx xi yi xip1 yip1 sign zz relinc
slope scale-x scale-y max-xf min-xf ox infin $zmax1 $zmin1
(f-intercept flonum flonum flonum flonum flonum)
(call-x flonum flonum flonum) (call-y flonum flonum flonum))
(fixnum maxdim ng jj ig it igg itt indexg indext i j k kk k1 k2 len
n1 n2 xinc yinc iwhich symtype type xsta ysta zinc zsta ksta
dk kend (lookupx flonum fixnum) (lookupxg flonum fixnum))
(notype (enlarge-array) (surf-expand2 notype))
(array* (flonum x-arr 1. y-arr 1. z-arr 1. x-3d 1. y-3d 1.
xg-3d 1. g-3d 1. xh-3d 1. h-3d 1.)))
(progn (array x-3d flonum 1.) (array y-3d flonum 1.)
(array xg-3d flonum 2.) (array g-3d flonum 2.)
(array xh-3d flonum 1.) (array h-3d flonum 1.))
(defun hide-init nil
((lambda (infin)
(*rearray 'xg-3d 'flonum 2.) (*rearray 'g-3d 'flonum 2.)
(store (xg-3d 0.) (-$ infin)) (store (xg-3d 1.) infin)
(store (g-3d 0.) (-$ infin)) (store (g-3d 1.) (-$ infin)))
1.e15))
(defun howclose-line (mark)
((lambda (len xsta ysta)
(setq mark (cddddr mark))
(*$ 0.5 (+$ (funcall $howclose (x-arr xsta) (y-arr ysta) 0.0)
(funcall $howclose (x-arr (+ xsta (* len (car mark))))
(y-arr (+ ysta (* len (cadr mark)))) 0.0))))
(1- (car mark)) (cadr mark) (caddr mark)))
;;; 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 (append marks nil))
(cond ((null (caar marks))
((lambda (marks1 marks2 border)
(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 (1+ i)))
((apply (function 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 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)))))
nil nil nil))
(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) (*rearray arr 'flonum 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 (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 (1- (car mark)))
(ylen (1- (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 (1+ ylen) xstart ystart zstart 0 yinc zinc2)
(cons (list (1+ xlen) xstart ystart zstart xinc 0 zinc1)
(append
(do ((l) (i 0 (1+ i)) (nummax (1+ (* 2 xlen)))
(num 3 (min nummax (+ 2 num)))
(xs xstart) (dx1 xinc) (dx2 0)
(ys (+ ystart yinc) (+ ys yinc)) (dy1 0) (dy2 (- yinc))
(zs (+ zstart zinc2) (+ zs zinc2)) (dz1 zinc1) (dz2 (- 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 (1+ i)) (nummax (1+ (* 2 ylen)))
(num 3 (min nummax (+ 2 num)))
(xs (+ xstart (* (1- xlen) xinc)) (- xs xinc)) (dx2 0) (dx1 xinc)
(ys (+ ystart (* ylen yinc))) (dy2 (- yinc)) (dy1 0)
(zs (+ zstart (* (1- xlen) zinc1) (* ylen zinc2)) (- zs zinc1))
(dz2 (- 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 (1- (car mark)))
(yl (1- (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 (x-arr xs)) (x1 (x-arr (+ xs (* xl xj))))
(y0 (y-arr ys)) (y1 (y-arr (+ ys (* 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 (1+ xl) (1+ yl) 0)
(cond ((= c0 cmin)
(list xs ys zs xj yj zi1 zi2))
((= c1 cmin)
(list (+ xs (* xl xj)) ys (+ zs (* xl zi1)) (- xj) yj (- zi1) zi2))
((= c2 cmin)
(list (+ xs (* xl xj)) (+ ys (* yl yj)) (+ zs (* xl zi1) (* yl zi2))
(- xj) (- yj) (- zi1) (- zi2)))
(t
(list xs (+ ys (* yl yj)) (+ zs (* yl zi2)) xj (- yj) zi1 (- zi2)))))))))
(defun lookupx (xx jj)
(do ((i (1+ jj) (1+ i)))
((not (< (x-3d i) xx)) (cond ((= (x-3d i) xx) i) ((1- i))))))
(defun lookupxg (xx jj)
(do ((i (1+ jj) (1+ i)))
((not (< (xg-3d i) xx)) (cond ((= (xg-3d i) xx) i) ((1- i))))))
(defun enlarge-array nil
(setq maxdim (+ 50. maxdim))
(*rearray 'xh-3d 'flonum (1+ maxdim))
(*rearray 'h-3d 'flonum (1+ maxdim))
nil)
(defun f-intercept (xx xi yi xip1 yip1)
(+$ yi (*$ (-$ xx xi) (//$ (-$ yip1 yi) (-$ xip1 xi)))))
;; 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)
((lambda (eps n1 ng jj ig it igg itt indexg
indext f1 f2 x1 x2 z1 z2 last maxdim)
($changedash (\ type 10.))
(setq type (// type 10.))
(*rearray 'x-3d 'flonum n1)
(*rearray 'y-3d 'flonum 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 (x-arr xsta) (y-arr ysta) 0.0)
(call-x (x-arr (+ xsta
(* (// n1 2) xinc1)
(* (// (1- n1) 2) xinc2)))
(y-arr (+ ysta
(* (// n1 2) yinc1)
(* (// (1- n1) 2) yinc2)))
0.0)))
(dk 1) (ksta 0) (kend n1))
(and back (setq dk -1 ksta (1- n1) kend -1))
(do ((kk ksta (+ dk kk))
(i xsta (+ xinc i)) (xinc xinc1 (- qx xinc))
(qx (+ xinc1 xinc2))
(j ysta (+ yinc j)) (yinc yinc1 (- qy yinc))
(qy (+ yinc1 yinc2))
(k zsta (+ zinc k)) (zinc zinc1 (- qz zinc))
(qz (+ zinc1 zinc2)))
((= kk kend))
(store (x-3d kk)
(call-x (x-arr i) (y-arr j) (z-arr k)))
(store (y-3d kk)
(*$ sign (call-y (x-arr i) (y-arr j)
(z-arr k))))))))
(t
((lambda (xsta ysta zsta xinc yinc zinc)
(setq mark (cddddr mark) xinc (car mark) yinc (cadr mark)
zinc (caddr mark))
(cond ((> (call-x (x-arr xsta) (y-arr ysta) 0.0)
(call-x (x-arr (+ xsta (* (1- n1) xinc)))
(y-arr (+ ysta (* (1- n1) yinc))) 0.0))
(setq xsta (+ xsta (* (1- n1) xinc))
ysta (+ ysta (* (1- n1) yinc))
zsta (+ zsta (* (1- n1) zinc))
xinc (- xinc) yinc (- yinc) zinc (- zinc))))
(do ((kk 0. (1+ kk)) (i xsta (+ xinc i)) (j ysta (+ yinc j))
(k zsta (+ zinc k)))
((= kk n1))
(store (x-3d kk) (call-x (x-arr i) (y-arr j) (z-arr k)))
(store (y-3d kk) (*$ sign
(call-y (x-arr i) (y-arr j)
(z-arr k))))))
(cadr mark) (caddr mark) (cadddr mark) 0. 0. 0.)))
(setq n1 (1- n1))
(setq jj (lookupxg (x-3d 0.) 0.))
(do nil ((< jj maxdim)) (enlarge-array))
(do ((j 0. (1+ j)))
((> j jj))
(store (xh-3d j) (xg-3d j))
(store (h-3d j) (g-3d j)))
(setq ig (1+ jj))
(store (xh-3d ig) (x-3d 0.))
(store (h-3d ig)
(f-intercept (x-3d 0.) (xg-3d jj) (g-3d jj) (xg-3d ig)
(g-3d ig)))
(setq indexg jj
indext 0.
z1 (x-3d 0.)
f1 (-$ (h-3d ig) (y-3d 0.))
it 1.
jj ig)
(cond ((< (h-3d ig) (y-3d 0.))
(do nil ((< jj maxdim)) (enlarge-array))
(setq jj (1+ ig))
(store (h-3d jj) (y-3d 0.))
(store (xh-3d jj) (+$ z1 eps))))
(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 (x-3d n1) igg (lookupxg z2 indexg) itt (1- n1))
nil)
(cond ((< (xg-3d ig) (x-3d it))
(setq x2 (xg-3d ig)
iwhich 1.
f2 (-$ (g-3d ig)
(f-intercept x2 (x-3d (1- it))
(y-3d (1- it)) (x-3d it)
(y-3d it)))
ig (1+ ig)))
(t (setq iwhich 0.
x2 (x-3d it)
f2 (-$ (f-intercept x2 (xg-3d (1- ig))
(g-3d (1- ig)) (xg-3d ig)
(g-3d ig))
(y-3d it))
it (1+ it))))
(and (> it n1) (setq last t))
(cond ((or (and (> f1 0.0) (> f2 0.0))
(and (< f1 0.0) (< f2 0.0))
(and (= f1 0.0) (= f2 0.0)))
(setq x1 x2 f1 f2))
(t (setq slope (//$ (-$ f2 f1) (-$ x2 x1))
igg (- ig 1. iwhich)
itt (+ 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) (x-3d k1))
(t (f-intercept zz (x-3d k1) (y-3d k1)
(x-3d (1+ k1)) (y-3d (1+ k1)))))
(f-intercept zz (xg-3d k2) (g-3d k2)
(xg-3d (1+ k2)) (g-3d (1+ k2))))
(setq ngraph (- itt indext -2.))
(do nil
((not (> (+ jj ngraph -1.) maxdim)))
(enlarge-array))
(setq n2 jj)
(do ((i (1+ indext) (1+ i)))
((> i itt))
(setq jj (1+ jj))
(store (xh-3d jj) (x-3d i))
(store (h-3d jj) (y-3d i)))
(setq jj (1+ jj))
(store (xh-3d jj) z2)
(store (h-3d jj)
(f-intercept z2
(x-3d itt) (y-3d itt)
(x-3d (1+ itt)) (y-3d (1+ itt))))
(and ifplot (graph-hide n2 ngraph sign type)))
(t (do nil
((< (+ jj igg (- indexg)) maxdim))
(enlarge-array))
(cond ((not (= indexg igg))
(do ((i (1+ indexg) (1+ i)))
((> i igg))
(setq jj (1+ jj))
(store (xh-3d jj) (xg-3d i))
(store (h-3d jj) (g-3d i)))))
(setq jj (1+ jj))
(store (xh-3d jj) z2)
(store (h-3d jj)
(f-intercept z2 (x-3d itt) (y-3d itt)
(x-3d (1+ itt)) (y-3d (1+ itt))))))
(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 (> (+ jj 3. ng (- igg)) maxdim))) (enlarge-array))
(store (xh-3d (1+ jj)) (+$ (xh-3d jj) eps))
(setq jj (1+ jj))
(store (h-3d jj)
(f-intercept (x-3d n1) (xg-3d igg) (g-3d igg)
(xg-3d (1+ igg)) (g-3d (1+ igg))))
(do ((j (1+ igg) (1+ j)))
((> j ng))
(setq jj (1+ jj))
(store (xh-3d jj) (xg-3d j))
(store (h-3d jj) (g-3d j)))
(*rearray 'xg-3d 'flonum (1+ jj))
(*rearray 'g-3d 'flonum (1+ jj))
(do ((i 0. (1+ i)) (j 0. (1+ j)) (flg) (ox (*$ 2.0 (xh-3d 0.))))
((> i jj) (*rearray 'xg-3d 'flonum j)
(*rearray 'g-3d 'flonum j) nil)
(cond ((not (> (xh-3d i) (+$ ox eps eps)))
(setq ox (+$ ox eps))
(cond (flg (setq j (1- j))) (t (setq flg t))))
(t (setq flg nil ox (xh-3d i))))
(store (xg-3d j) ox)
(store (g-3d j) (h-3d i))))
(*$ 1.0e-5 (-$ max-xf min-xf))
(car mark)
(1- (cadr (arraydims 'g-3d)))
0. 0. 0. 0. 0. 0. 0.
0.0 0.0 0.0 0.0 0.0 0.0 nil -1.))))
(defun graph-hide (n2 ngraph sign symtype)
(setq ngraph (+ n2 ngraph) symtype (\ symtype 10.))
($setpoint (xh-3d n2) (*$ sign (h-3d n2)))
(or (= symtype 0.)
($drawsymbol (xh-3d n2) (*$ sign (h-3d n2)) symtype))
(do ((i (1+ n2) (1+ i)))
((= i ngraph))
($vector (xh-3d i) (*$ sign (h-3d i)))
(or (= symtype 0.)
($drawsymbol (xh-3d i) (*$ sign (h-3d i)) symtype))))
;;ref: the computer journal vol 15 num 4 p 382 (1972)
(declare (special maxdim s zds @n1 @n2 @xinc @xstart @yinc @ystart @zinc1 @zinc2 @zstart
@symtype $diag)
(flonum xx xi yi xip1 yip1 const phi1 phi2 phi3 phi4 phiav xav yav infin max
min pt z (f-intercept flonum flonum flonum flonum flonum)
(phi-cont fixnum fixnum) (x-cont fixnum) (y-cont fixnum))
(fixnum i j n1 n2 n5 @n1 @n2 @symtype type nm cn ent i1 i2 i3 i4 ib im ip j1 j2
j3 j4 jb jm jp k l ncn qq s zds maxdim ngraph zlen cnum
(ffnd fixnum fixnum fixnum fixnum fixnum fixnum fixnum fixnum fixnum
fixnum flonum))
(notype (bdyp fixnum fixnum) (look fixnum fixnum fixnum fixnum fixnum flonum)
(same-sign flonum flonum))
(array* (flonum cont-arr 1.) (fixnum xbd-cont 1. ybd-cont 1. itg-cont 2.)))
(setq $diag t)
(progn (array xbd-cont fixnum 1.) (array ybd-cont fixnum 1.)
(array cont-arr flonum 1.) (array itg-cont fixnum 1. 1.))
(defun bdyp (i j) (or (= j 0.) (= j (1- @n2)) (= i 0) (= i (1- @n1))))
(defun phi-cont (i j) (z-arr (+ @zstart (* i @zinc1) (* j @zinc2))))
(defun x-cont (i) (x-arr (+ @xstart (* i @xinc))))
(defun y-cont (j) (y-arr (+ @ystart (* j @yinc))))
(defun contour-set (contours cmin cmax)
(cond (($listp contours)
(*rearray 'cont-arr 'flonum (1- (length contours)))
(fillarray 'cont-arr (mapcar 'fmeval (cdr contours))))
((and (get contours 'array) (eq (car (arraydims contours)) 'flonum))
(*rearray 'cont-arr 'flonum (cadr (arraydims contours)))
(fillarray 'cont-arr contours))
(t ((lambda (infin min max cnum intflg)
(or intflg (numberp contours) (setq contours 20.))
(setq min infin max (-$ infin))
(cond ((not (and (numberp cmax) (numberp cmin)))
(do ((i 0. (1+ i)) (zlen (cadr (arraydims 'z-arr))) (pt))
((= i zlen))
(setq pt (z-arr 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 (- (fix (-$ min)))))
(cond ((not (> max min)) (setq max (+$ 1.0 min))))
(setq cnum (fix (+$ 0.5 (-$ 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. (- contours))))
(t ((lambda (ll)
(setq min (cadr ll)
max (caddr ll)
cnum (fix (+$ (//$ (-$ max min)
(car ll))
1.5))))
(cdr (progn (cond ((< max min) (setq max (prog2 nil min (setq min max)))))
(cond ((= max min)
(cond ((= max 0.0) (setq max 1.0 min -1.0))
(t (setq max (+$ max (*$ (abs max) 0.1))
min (-$ min (*$ (abs min) 0.1)))))))
(scale1 contours min max))))))))
(*rearray 'cont-arr 'flonum cnum)
(setq $zmax1 max $zmin1 min)
(do ((i 0. (1+ i)) (pt 0.0)) ((= i cnum))
(cond (intflg
(store (cont-arr i) (+$ min (float i))))
(t
(setq pt (//$ (float i)
(float (1- cnum))))
(store (cont-arr i)
(+$ (*$ min (-$ 1.0 pt))
(*$ max pt)))))))
(^$ 8.0 42.) 0.0 0.0 0. (eq contours '$integer)))))
(defun contour-init (marks)
((lambda (mark)
(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)))
(car marks)))
(defun contour-drive (marks typel)
((lambda (@n1 @n2 @xstart @ystart @zstart @xinc @yinc @zinc1 @zinc2 n5 ncn
maxdim)
(do ((typel1 typel (cdr typel1)) (type)) ((null marks))
(setq marks (contour-init marks)
n5 (+ (* 2. @n1) (* 2. @n2) -3.))
(cond ((null typel1) (setq typel1 typel)))
(setq type (car typel1))
(*rearray 'xbd-cont 'fixnum n5)
(*rearray 'ybd-cont 'fixnum n5)
(*rearray 'itg-cont 'fixnum @n1 @n2)
(do ((i 0 (1+ i))) ((= i @n2))
(store (ybd-cont i) i) (store (xbd-cont i) 0))
(do ((i 1. (1+ i))) ((= i @n1))
(store (ybd-cont (+ @n2 i -1.)) (1- @n2))
(store (xbd-cont (+ @n2 i -1.)) i))
(do ((i (- @n2 2) (1- i))) ((< i 0.))
(store (ybd-cont (- (+ (* 2 @n2) @n1 -3) i)) i)
(store (xbd-cont (- (+ (* 2 @n2) @n1 -3) i)) (1- @n1)))
(do ((i (- @n1 2) (1- i))) ((< i 0.))
(store (ybd-cont (- (+ (* 2 @n2) (* 2 @n1) -4.) i)) 0.)
(store (xbd-cont (- (+ (* 2 @n2) (* 2 @n1) -4.) i)) i))
(setq ncn (cadr (arraydims 'cont-arr)))
($changedash (\ type 10.))
(setq type (\ (// type 10.) 10.))
(enlarge-array)
(contor ncn n5 type)
(*rearray 'xbd-cont 'fixnum 1.) (*rearray 'ybd-cont 'fixnum 1.)
(*rearray 'itg-cont 'fixnum 1. 1.) (*rearray 'xh-3d 'flonum 1.)
(*rearray 'h-3d 'flonum 1.))
(*rearray 'cont-arr 'flonum 1.))
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1.))
(defun contor (ncn n5 @symtype)
(do ((cn 0. (1+ cn)) (const) (i) (j) (ib) (jb))
((= cn ncn))
(setq const (cont-arr cn))
(do ((i 0. (1+ i))) ((= i @n1))
(do ((j 0. (1+ j))) ((= j @n2))
(store (itg-cont i j) 0.)))
(do ((k 1. (1+ k))) ((= k n5))
(setq i (xbd-cont k) j (ybd-cont k)
ib (xbd-cont (1- k)) jb (ybd-cont (1- k)))
(cond ((not (eq (< (-$ (phi-cont i j) const) 0.0)
(< (-$ (phi-cont ib jb) const) 0.0)))
(look i j ib jb 1. const))))
(do ((k 1. (1+ k))) ((= k (1- @n1)))
(do (( l 1. (1+ l))) ((= l (1- @n2)))
(setq i k ib k j l jb (1- l))
(cond ((and (not (and (bdyp i j) (bdyp ib jb)))
(not (eq (< (-$ (phi-cont i j) const) 0.0)
(< (-$ (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 (1+ j) ip (1+ i) jm (1- j) im (1- i) zds 0.)
(cond
((= jb jm)
(and (> (itg-cont i jm) 1.) (return nil))
(setq ent 1.)
(store (xh-3d 0.) (x-cont i))
(store (h-3d 0.)
(f-intercept const (phi-cont i jm) (y-cont jm) (phi-cont i j) (y-cont j))))
((= ib im)
(and (or (= (itg-cont im j) 1.) (= (itg-cont im j) 3.)) (return nil))
(setq ent 2.)
(store (h-3d 0.) (y-cont j))
(store (xh-3d 0.)
(f-intercept const (phi-cont im j) (x-cont im) (phi-cont i j) (x-cont i))))
((= jb jp)
(and (> (itg-cont i j) 1.) (return nil))
(setq ent 3.)
(store (xh-3d 0.) (x-cont i))
(store (h-3d 0.)
(f-intercept const (phi-cont i j) (y-cont j) (phi-cont i jp) (y-cont jp))))
(t (and (or (= (itg-cont i j) 1.) (= (itg-cont i j) 3.)) (return nil))
(setq ent 4.)
(store (h-3d 0.) (y-cont j))
(store (xh-3d 0.)
(f-intercept const (phi-cont i j) (x-cont i)
(phi-cont ip j) (x-cont ip)))))
(setq s 1.)
(do
nil (nil)
(setq ip (1+ i) jp (1+ j) im (1- i) jm (1- j))
(cond ((= ent 1.) (store (itg-cont i jm) (+ (itg-cont i jm) 2.))
(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.) (store (itg-cont im j) (1+ (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.) (store (itg-cont i j) (+ (itg-cont i j) 2.))
(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 (store (itg-cont i j) (1+ (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.) (store (itg-cont i (1- j)) (+ 2. (itg-cont i (1- j)))))
((= ent 2.) (store (itg-cont (1- i) j) (1+ (itg-cont (1- i) j))))
((= ent 3.) (store (itg-cont i j) (+ 2 (itg-cont i j))))
(t (store (itg-cont i j) (1+ (itg-cont i j)))))
(return nil)))
(cond ((= qq 2.)
(cond ((= ent 1.) (and (> (itg-cont i (1- j)) 1.) (return nil)))
((= ent 2.) (and (oddp (itg-cont (1- i) j)) (return nil)))
((= ent 3.) (and (> (itg-cont i j) 1.) (return nil)))
(t (and (oddp (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 ((> (+ s 4) maxdim) (enlarge-array)))
((lambda (phi1 phi2 phi3 phi4 phiav revflag xav yav)
(setq phiav (//$ (+$ 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
(store (xh-3d s) (f-intercept 0.0 phi1 (x-cont i1) phiav xav))
(store (h-3d s) (f-intercept 0.0 phi1 (y-cont j1) phiav yav))
(setq s (1+ s))))
(do ((i 0. (1+ i)))
((not (same-sign phi1 phi2))
(store (xh-3d s) (f-intercept 0.0 phi1 (x-cont i1) phi2
(x-cont i2)))
(store (h-3d s) (f-intercept 0.0 phi1 (y-cont j1) phi2
(y-cont j2)))
(setq s (1+ s))
(and (= qq 1.) (bdyp i1 j1) (bdyp i2 j2) (setq zds 1.))
(and revflag (not (oddp i)) (setq i (+ 2. i)))
(1+ (\ (+ i ent 2.) 4.)))
(cond ((and $diag (not (= phiav 0.0)))
(store (xh-3d s) (f-intercept 0.0 phi2 (x-cont i2) phiav
xav))
(store (h-3d s) (f-intercept 0.0 phi2 (y-cont j2) phiav yav))
(setq s (1+ 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)))))
(-$ (phi-cont i1 j1) const) (-$ (phi-cont i2 j2) const)
(-$ (phi-cont i3 j3) const) (-$ (phi-cont i4 j4) const)
0.0 nil (//$ (+$ (x-cont i1) (x-cont i3)) 2.0)
(//$ (+$ (y-cont j1) (y-cont j3)) 2.0)))
(defun graph-contour (ngraph symtype z)
($setpoint3 (xh-3d 0.) (h-3d 0.) z)
(or (= symtype 0.)
($drawsymbol3 (xh-3d 0.) (h-3d 0.) z symtype))
(do ((i 1. (1+ i)))
((= i ngraph))
($vector3 (xh-3d i) (h-3d i) z)
(or (= symtype 0.)
($drawsymbol3 (xh-3d i) (h-3d i) z symtype)))
(and $labelcontours
($ghprint (pfp1 z (fix $plotnumprec))
(tek-x (call-x (xh-3d (// ngraph 2.)) (h-3d (// ngraph 2.)) z))
(tek-y (call-y (xh-3d (// ngraph 2.)) (h-3d (// ngraph 2.)) z))
1.)))
(comment
(defun $test (ifplot)
((lambda (eps n1 ng jj ig it igg itt indexg indext f1 f2 x1 x2 z1 z2 last)
(setq maxdim -1. n1 (1- n1))
(setq jj (lookupxg (x-3d 0.) 0.))
(do nil ((< jj maxdim)) (enlarge-array))
(do ((j 0. (1+ j)))
((> j jj))
(store (xh-3d j) (xg-3d j))
(store (h-3d j) (g-3d j)))
(setq ig (1+ jj))
(store (xh-3d ig) (x-3d 0.))
(store (h-3d ig)
(f-intercept (x-3d 0.) (xg-3d jj) (g-3d jj) (xg-3d ig)
(g-3d ig)))
(setq indexg jj
indext 0.
z1 (x-3d 0.)
f1 (-$ (h-3d ig) (y-3d 0.))
it 1.
jj ig)
(cond ((< (h-3d ig) (y-3d 0.))
(do nil ((< jj maxdim)) (enlarge-array))
(setq jj (1+ ig))
(store (h-3d jj) (y-3d 0.))
(store (xh-3d jj) (+$ z1 eps))))
(setq x1 z1)
(do ((zz 0.0) (k1 0.) (k2 0.) (n2 0.) (ngraph 0.))
(nil)
(do ((iwhich 0.) (slope 0.0) (relinc (//$ scale-x scale-y)))
(last (setq z2 (x-3d n1) igg (lookupxg z2 indexg) itt (1- n1))
nil)
(cond ((< (xg-3d ig) (x-3d it))
(setq x2 (xg-3d ig)
iwhich 1.
f2 (-$ (g-3d ig)
(f-intercept x2 (x-3d (1- it))
(y-3d (1- it)) (x-3d it)
(y-3d it)))
ig (1+ ig)))
(t (setq iwhich 0.
x2 (x-3d it)
f2 (-$ (f-intercept x2 (xg-3d (1- ig))
(g-3d (1- ig)) (xg-3d ig)
(g-3d ig))
(y-3d it))
it (1+ it))))
(and (> it n1) (setq last t))
(cond ((and (= (lsh f1 -35.) (lsh f2 -35.))
(= (lsh (-$ f1) -35.) (lsh (-$ f2) -35.)))
(setq x1 x2 f1 f2))
(t (setq slope (//$ (-$ f2 f1) (-$ x2 x1))
igg (- ig 1. iwhich)
itt (+ 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) (x-3d k1))
(t (f-intercept zz (x-3d k1) (y-3d k1)
(x-3d (1+ k1)) (y-3d (1+ k1)))))
(f-intercept zz (xg-3d k2) (g-3d k2)
(xg-3d (1+ k2)) (g-3d (1+ k2))))
(setq ngraph (- itt indext -2.))
(do nil
((not (> (+ jj ngraph -1.) maxdim)))
(enlarge-array))
(setq n2 jj)
(do ((i (1+ indext) (1+ i)))
((> i itt))
(setq jj (1+ jj))
(store (xh-3d jj) (x-3d i))
(store (h-3d jj) (y-3d i)))
(setq jj (1+ jj))
(store (xh-3d jj) z2)
(store (h-3d jj)
(f-intercept z2
(x-3d itt) (y-3d itt)
(x-3d (1+ itt)) (y-3d (1+ itt))))
(and ifplot (graph-hide n2 ngraph sign type)))
(t (do nil
((< (+ jj igg (- indexg)) maxdim))
(enlarge-array))
(cond ((not (= indexg igg))
(do ((i (1+ indexg) (1+ i)))
((> i igg))
(setq jj (1+ jj))
(store (xh-3d jj) (xg-3d i))
(store (h-3d jj) (g-3d i)))))
(setq jj (1+ jj))
(store (xh-3d jj) z2)
(store (h-3d jj)
(f-intercept z2 (x-3d itt) (y-3d itt)
(x-3d (1+ itt)) (y-3d (1+ itt))))))
(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 (> (+ jj 3. ng (- igg)) maxdim))) (enlarge-array))
(store (xh-3d (1+ jj)) (+$ (xh-3d jj) eps))
(setq jj (1+ jj))
(store (h-3d jj)
(f-intercept (x-3d n1) (xg-3d igg) (g-3d igg)
(xg-3d (1+ igg)) (g-3d (1+ igg))))
(do ((j (1+ igg) (1+ j)))
((> j ng))
(setq jj (1+ jj))
(store (xh-3d jj) (xg-3d j))
(store (h-3d jj) (g-3d j)))
(*rearray 'xg-3d 'flonum (1+ jj))
(*rearray 'g-3d 'flonum (1+ jj))
(do ((i 0. (1+ i)) (j 0. (1+ j)) (flg) (ox (*$ 2.0 (xh-3d 0.))))
((> i jj) (*rearray 'xg-3d 'flonum j)
(*rearray 'g-3d 'flonum j) nil)
(cond ((not (> (xh-3d i) (+$ ox eps eps)))
(setq ox (+$ ox eps))
(cond (flg (setq j (1- j))) (t (setq flg t))))
(t (setq flg nil ox (xh-3d i))))
))
(*$ 1.0e-5 (-$ max-xf min-xf))
(cadr (arraydims 'x-3d))
(1- (cadr (arraydims 'g-3d)))
0. 0. 0. 0. 0. 0. 0.
0.0 0.0 0.0 0.0 0.0 0.0 nil)))
Reference in New Issue View Git Blame Copy Permalink