home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / plot.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  43.0 KB  |  1,560 lines
  1. ;;Copyright William F. Schelter 1990, All Rights Reserved
  2. (in-package "MAXIMA")
  3. ;; see bottom of file for examples
  4.  
  5.  
  6. ;; I (RLT) don't think we should compile this with safety = 0.
  7.  
  8. ;;(eval-when (compile) (proclaim '(optimize (safety 0))))
  9.  
  10.  
  11. (eval-when (compile eval load)
  12.  (defmacro coerce-float (x)
  13.    `(lisp::float (meval* ,x) 1.d0))
  14.   )
  15.  
  16. (defvar *z-range* nil)
  17. (defvar *original-points* nil)
  18. (defvar $axes_length 4.0)
  19. (defvar *rot* (make-array 9 :element-type 'long-float))
  20. (defvar $rot nil)
  21.  
  22. (defvar $plot_options '((mlist)
  23.             ((mlist) |$x| -3 3)
  24.             ;; Make the default range on Y large.  Don't
  25.             ;; use most-positive-double-float because this
  26.             ;; causes overflow in the draw2d routine.
  27.             ((mlist) |$y| #.(- (/ most-positive-double-float 1024))
  28.                           #.(/ most-positive-double-float 1024))
  29.             ((mlist) |$t| -3 3)
  30.             ((mlist) $grid 30 30)
  31.             ((mlist) $view_direction 1 1 1)
  32.             ((mlist) $colour_z nil)
  33.             ((mlist) $transform_xy nil)
  34.             ((mlist) $run_viewer t)
  35.             ((mlist) $plot_format $openmath)
  36.             ((mlist) $nticks 100)
  37.             ))
  38.  
  39. (defun $get_plot_option (name &optional n)
  40.   (sloop for v in (cdr $plot_options)
  41.      when (eq (nth 1 v) name) do
  42.      (return (if n (nth n  v) v))))
  43.  
  44. (defun check-list-items (name lis type length)
  45.   (or (eql (length lis) length)
  46.       (merror "~M items were expected in the ~M list" length name))
  47.   `((mlist) , name ,@
  48.     (sloop for v in lis
  49.        do (setq v (meval* v))
  50.        when (not (typep v type))
  51.        do
  52.        (merror "~M is not of the right type in ~M list" v name)
  53.        collect v)
  54.     ))
  55.  
  56. (defun $set_plot_option ( value)
  57.   (setq $plot_options ($copylist $plot_options))
  58.   (unless (and  ($listp value)
  59.         (symbolp (setq name (nth 1 value))))
  60.       (merror "~M is not a plot option.  Must be [symbol,..data]" value))
  61.   (setq value
  62.     (case name
  63.       ((|$x| |$y|) (check-list-items name (cddr value) 'number 2)
  64.        (check-range value)
  65.        )
  66.       ($view_direction (check-list-items name (cddr value) 'number 3))
  67.       ($grid  (check-list-items name (cddr value) 'fixnum 2))
  68.       ($nticks  (check-list-items name (cddr value) 'fixnum 1))
  69.       (($colour_z $run_viewer $transform_xy)
  70.        (check-list-items name (cddr value) 't 1))
  71.       ($plot_format (or (member (nth 2 value) '($zic $geomview $ps
  72.                              $gnuplot
  73.                              $zplot
  74.                              $openmath
  75.                              ))
  76.                 (merror "Only [zic,geomview,ps,openmath,gnuplot] are available"))
  77.             value)
  78.       (t
  79.        (merror "Unknown plot option specified:  ~M" name))))
  80.   (sloop for v on (cdr $plot_options)
  81.      when (eq (nth 1 (car v)) name)
  82.      do (setf (car v) value))
  83.   $plot_options
  84.   )
  85.   
  86. (defvar $pstream nil)
  87.  
  88.  
  89.  
  90. (defun print-pt1 (f str)
  91.   (format str "~,3f " f))
  92.  
  93. (defmacro defbinop (name op type)
  94.   `(progn
  95.      (defun ,name (x y) (the ,type (,op  (the ,type x) (the ,type y))))
  96.      (eval-when (compile eval)     
  97.      (#+kcl si::DEFINE-COMPILER-MACRO #-kcl
  98.         defmacro ,name (x y)
  99.          `(the ,',type (,',op  (the ,',type ,x) (the ,',type ,y)))))))
  100.  
  101. (defbinop f+ + fixnum)
  102. (defbinop f- - fixnum)
  103.  
  104. (defbinop $+ + long-float)
  105. (defbinop $- - long-float)
  106. (defbinop $* * long-float)
  107. (defbinop $/ / long-float)
  108.  
  109.  
  110. (defstruct (polygon (:type list)
  111.             (:constructor %make-polygon (pts edges))
  112.             )
  113.   (dummy '($polygon simp))
  114.   pts edges)
  115.  
  116. (eval-when (compile eval)
  117.  
  118. (defmacro f* (a b &optional c)
  119.   (if c `(f* (f* ,a ,b) ,c)
  120.     `(the fixnum (* (the fixnum ,a) (the fixnum ,b)))))
  121.  
  122. (defmacro float-< (a b) `(< (the long-float ,a) (the long-float ,b)))
  123.   
  124.  
  125.  
  126.  
  127. (defmacro z-pt (ar i) `(aref ,ar (the fixnum (+ 2 (* ,i 3)))))
  128. (defmacro y-pt (ar i) `(aref ,ar (the fixnum (+ 1 (* ,i 3)))))
  129. (defmacro x-pt (ar i) `(aref ,ar (the fixnum (+ 0 (* ,i 3)))))
  130. (defmacro rot (m i j) `(aref ,m (the fixnum (+ ,i (the fixnum (* 3 ,j))))))
  131.  
  132.  
  133.  
  134. (defmacro print-pt (f)
  135.   `(print-pt1 ,f $pstream ))
  136.  
  137. (defmacro make-polygon (a b) `(list '($polygon) ,a ,b))
  138. )
  139.  
  140.  
  141.  
  142. (defun draw3d (f minx maxx miny maxy  nxint nyint)
  143.   (let* ((epsx (/ (- maxx minx) nxint))
  144.      (x 0.0)  ( y 0.0)
  145.      (epsy (/ (- maxy miny) nyint))
  146.      (nx (+ nxint 1))
  147.      (l 0)
  148.      (ny (+ nyint 1))
  149.      (ar (make-array  (+ 12  ; 12  for axes
  150.                  (f* (f* 3 nx) ny))  :fill-pointer (f* (f* 3 nx) ny)
  151.              :element-type 'long-float
  152.              :adjustable t
  153.              )))
  154.     (declare (long-float x y epsy epsx)
  155.          (fixnum nx  ny l)
  156.          (type (array long-float) ar))
  157.     (sloop for j below ny
  158.        initially (setq y miny)
  159.        do (setq x minx)
  160.        (sloop for i below nx
  161.           do
  162.           (setf (x-pt ar l) x)
  163.           (setf (y-pt ar l) y)
  164.           (setf (z-pt ar l) (funcall f x y))
  165.           (incf l)
  166.           (setq x (+ x epsx))
  167.           )
  168.        (setq y (+ y epsy)))
  169.     (make-polygon  ar  (make-grid-vertices nxint nyint))))
  170.  
  171. ;; The following is 3x2 = 6 rectangles
  172. ;; call (make-vertices 3 2)
  173. ;; there are 4x3 = 12 points.
  174. ;; ordering is x0,y0,z0,x1,y1,z1,....,x11,y11,z11
  175. ;; ----
  176. ;; ||||
  177. ;; ----
  178. ;; ||||
  179. ;; ----
  180.  
  181. (defun make-grid-vertices (nx ny)
  182.   (declare (fixnum nx ny))
  183.   (let* ((tem (make-array (+ 15 (f* 5 nx ny)) :fill-pointer (f* 5 nx ny)
  184.               :adjustable t
  185.               :element-type '(mod  65000)))
  186.      (m  nx )
  187.      (nxpt (+ nx 1))
  188.      (i 0)
  189.      )
  190.     (declare (fixnum i nxpt m)
  191.          (type (array (mod 65000)) tem))
  192.     (sloop for k below (length tem)
  193.        do
  194.        (setf (aref tem k) i)
  195.        (setf (aref tem (incf k))
  196.          (+ nxpt i))
  197.        (setf (aref tem (incf k))
  198.          (+ nxpt (incf i )))
  199.        (setf (aref tem (incf k)) i)
  200.        (setf (aref tem (incf k)) 0)    ;place for max
  201.        (setq m (- m 1))
  202.        (cond ((eql  m 0)
  203.           (setq m nx)
  204.           (setq i (+ i 1))))
  205.        )
  206.     tem))
  207.  
  208. (defun add-axes (pts vert)
  209.   (let ((origin (/ (length pts) 3)))
  210.     (sloop for i from -3 below 9
  211.        do
  212.        (vector-push-extend  (if (eql 0 (mod i 4)) $axes_length 0.0) pts))
  213.     (sloop for i below 15
  214.        do (vector-push-extend
  215.            (if (eql 1 (mod i 5)) (+ origin (ceiling i 5))  origin)
  216.            vert)
  217.        )))
  218.  
  219. (defun $rotation1 (phi th)
  220.   (let ((sinph (sin phi))
  221.     (cosph (cos phi))
  222.     (sinth (sin th))
  223.     (costh (cos th)))
  224.     `(($MATRIX SIMP)
  225.       ((MLIST SIMP) ,(* COSPH COSTH)
  226.        ,(* -1.0 COSPH SINTH)
  227.        ,SINPH)
  228.       ((MLIST SIMP) ,SINTH ,COSTH 0.0)
  229.       ((MLIST SIMP) ,(- (*  SINPH COSTH))
  230.        ,(* SINPH SINTH)
  231.        ,COSPH))))
  232.    
  233. ; pts is a vector of bts [x0,y0,z0,x1,y1,z1,...] and each tuple xi,yi,zi is rotated
  234. ; also the *z-range* is computed.
  235. (defun $rotate_pts(pts rotation-matrix)
  236.   (or ($matrixp rotation-matrix) (error "second arg not matrix"))
  237.   (let* ((rot *rot*)
  238.      (l (length pts))
  239.      (x 0.0) (y 0.0) (z 0.0)
  240.      )
  241.     (declare (long-float  x y z))
  242.     (declare (type (array long-float) rot))
  243.     ($copy_pts rotation-matrix *rot* 0)
  244.     
  245. ;    (setf (rot rot  0 0) (* cosphi costh))
  246. ;    (setf (rot rot  1 0) (- sinth))
  247. ;    (setf (rot rot 2 0) (* costh sinphi))
  248. ;
  249. ;    (setf (rot rot  0 1) (* cosphi sinth))
  250. ;    (setf (rot rot  1 1) costh)
  251. ;    (setf (rot rot 2 1) (* sinth sinphi))
  252. ;
  253. ;    (setf (rot rot  0 2) (- sinphi))
  254. ;    (setf (rot rot  1 2) 0.0)
  255. ;    (setf (rot rot 2 2) cosphi)
  256.  
  257.     (sloop with j = 0
  258.        while (< j l)
  259.        do
  260.        (setq x (aref pts j))
  261.        (setq y (aref pts (+ j 1)))
  262.        (setq z (aref pts (+ j 2)))
  263.        (sloop for i below 3 with a = 0.0
  264.           declare (long-float a)
  265.           do
  266.           (setq a (* x (aref rot (+ (* 3 i) 0))))
  267.           (setq a (+ a (* y (aref rot (+ (* 3 i) 1)))))
  268.           (setq a (+ a (* z (aref rot (+ (* 3 i) 2)))))
  269.           (setf (aref pts (+ j i )) a))
  270.        (setf j (+ j 3)))))
  271.  
  272. (defun $rotate_list (x)
  273.   (cond ((and ($listp x) (not (mbagp (nth 1 x))))
  274.          ($list_matrix_entries (ncmul2  $rot x)))
  275.         ((mbagp x) (cons (car x) (mapcar '$rotate_list (cdr x))))))
  276.  
  277. (defun $get_range (pts k &aux (z 0.0) (max most-negative-long-float) (min most-positive-long-float))
  278.   (declare (long-float z max min))
  279.   (declare (type (vector long-float) pts))
  280.   (sloop for i from k below (length pts) by 3
  281.      do (setq z (aref pts i))
  282.      (cond ((< z min) (setq min z)))
  283.      (cond ((> z max) (setq max z))))
  284.   (list min max (- max min)))
  285.  
  286.  
  287. ;; if this is '$polar_to_xy then the x y coordinates are interpreted as r theta
  288.  
  289. (defun add-ps-finish (opts)
  290.   (p (if opts
  291.      "/xr .30 def
  292. /xb 0.60 def
  293. /xg .60  def
  294. /myset { .005 mul dup xr add exch
  295.  dup xb add exch
  296.  xg add
  297. setrgbcolor} def
  298.  
  299. /myfinish { myset  gsave fill grestore 0 setgray stroke  } def"
  300.  
  301.        "/myfinish {.9 setgray gsave fill grestore .1 setgray stroke  } def")))
  302.  
  303.  
  304. (defun $draw_ngons(pts ngons number_edges &aux (i 0)(j 0) (s 0)
  305.                (opts *original-points*)
  306.                (maxz  most-negative-long-float))
  307.   (declare (type (array long-float) pts opts)
  308.        (type (array (mod 64000)) ngons)
  309.        (fixnum number_edges i s j number_edges)
  310.        (long-float maxz))
  311.   (setq j (length ngons))
  312.   (add-ps-finish opts)
  313.   (sloop while (< i j) 
  314.      do 
  315.      (sloop initially (setq s number_edges)
  316.         do
  317.         ;(print-pt (aref pts (f* 3 (aref ngons i))))
  318.         (print-pt (x-pt pts  (aref ngons i)))
  319.         ;(print-pt (aref pts (f+ 1 (f* 3 (aref ngons i)))))
  320.         (print-pt (y-pt pts  (aref ngons i)))
  321.         
  322.         (cond (opts (if (> (z-pt opts (aref ngons i)) maxz)
  323.                   (setq maxz (z-pt opts (aref ngons i))))))
  324.         (cond ((eql number_edges s) (p " moveto %"
  325.                            ;(aref pts (f+ 2 (f* 3 (aref ngons i))))
  326.                            
  327.                            ))
  328.               (t (p "lineto %" ;(aref pts (f+ 2 (f* 3 (aref ngons i))))
  329.                 )))
  330.         (setq i (f+ i 1))
  331.         while (> (setq s (f- s 1)) 0))
  332.      (setq i (f+ i 1))
  333.      (cond (opts
  334.         (p (f+ 1 (round ($* 100.0 ($/ ($- maxz (car *z-range*))
  335.                       (or (third *z-range*)
  336.                           ($- (second *z-range*) (car *z-range*))))))))
  337.         (setq maxz most-negative-long-float)
  338.         ))
  339.      (p " myfinish")
  340.      ))
  341.  
  342. ;; figure out the rotation to make pt the direction from which we view,
  343. ;; and to rotate z axis to vertical.
  344. ;; First get v, so p.v=0 then do u= p X v to give image of y axis
  345. ;; ans = transpose(matrix( v,u,p))
  346.  
  347. (defun $norm (pt) (sloop for v in (cdr pt) sum (* v v)))
  348. (defun $length_one (pt)
  349.   (let ((len (sqrt ($norm pt))))
  350.     (cons '(mlist) (sloop for v in (cdr pt) collect (/  (float v) len)))))
  351.  
  352. (defun $cross_product (u v)
  353.   (flet ((cp (i j)
  354.          (- (* (nth i u) (nth j v))
  355.         (* (nth i v) (nth j u)))))
  356.     `((mlist) ,(cp 2 3) ,(cp 3 1) ,(cp 1 2))))
  357.     
  358. (defun $GET_ROTATION (pt)
  359.   (setq pt ($length_one pt))
  360.   (let (v tem u)
  361.     (cond((setq tem (position 0.0 pt))
  362.        (setq v (cons '(mlist) (list 0.0 0.0 0.0)))
  363.        (setf (nth tem v) 1.0))
  364.       (t (setq v ($length_one `((mlist) ,(- (nth 2 pt))      , (nth 1 pt) 0.0)))))
  365.     (setq u ($cross_product pt v))
  366.     (let* (($rot   `(($matrix) ,v,u,pt))
  367.        (th (get-theta-for-vertical-z
  368.         (nth 3 (nth 1 $rot))
  369.         (nth 3 (nth 2 $rot)))))
  370.       (or (zerop th)
  371.       (setq $rot (ncmul2 ($rotation1 0.0 th)     $rot)))
  372.       $rot)))
  373.  
  374. (defun get-theta-for-vertical-z (z1 z2)
  375.   (cond ((eql z1 0.0) (if (> z2 0.0) 0.0  pi))
  376.     (t (lisp::atan  z2 z1 ))))
  377.  
  378. (defun $ps_axes ( rot )
  379.   (let ((tem (make-array 9 :element-type 'long-float)))
  380.     (setf (aref tem 0) 4.0)
  381.     (setf (aref tem 4) 4.0)
  382.     (setf (aref tem 8) 4.0)
  383.     ($rotate_pts tem rot)
  384.     (p 0 0 "moveto")
  385.     (p (aref tem 0) (aref tem 1) "lineto stroke")
  386.     (p  (aref tem 0) (aref tem 1) "moveto (x) show")
  387.     (p 0 0 "moveto")
  388.     (p (aref tem 3) (aref tem 4) "lineto stroke")
  389.     (p  (aref tem 3) (aref tem 4) "moveto (y) show")    
  390.     (p 0 0 "moveto")
  391.     (p (aref tem 6) (aref tem 7) "lineto stroke")
  392.     (p  (aref tem 6) (aref tem 7) "moveto (z) show")    
  393.     ))
  394.  
  395. (defun $polar_to_xy (pts &aux (r 0.0) (th 0.0))
  396.   (declare (long-float r th))
  397.   (declare (type (array long-float) pts))
  398.   (assert (typep pts '(vector long-float)))
  399.   (sloop for i below (length pts) by 3
  400.      do (setq r (aref pts i))
  401.      (setq th (aref pts (f+ i 1)))
  402.      (setf (aref pts i) (* r (cos th)))
  403.      (setf (aref pts (f+ i 1)) (* r (sin th)))))
  404.  
  405. (defun coerce-function-body (f lvars)
  406.   (setq f (coerce-float-fun f lvars))
  407.   (if (symbolp f) (symbol-function f) f))
  408.  
  409. ;; return a function suitable for the transform function in plot3d.
  410. (defun $make_transform (lvars fx fy fz  &aux ( $numer t))
  411.   (setq fx (coerce-function-body fx lvars))
  412.   (setq fy (coerce-function-body fy lvars))
  413.   (setq fz (coerce-function-body fz lvars))
  414.   (let ((sym (gensym "transform")))
  415.     (setf (symbol-function sym)
  416.   #'(lambda (pts &aux  (x1 0.0)(x2 0.0)(x3 0.0))
  417.       (declare (long-float  x1 x2 x3))
  418.       (declare (type (array long-float) pts))
  419.       (sloop for i below (length pts) by 3
  420.          do 
  421.      (setq x1 (aref pts i))
  422.      (setq x2 (aref pts (f+ i 1)))
  423.      (setq x3 (aref pts (f+ i 2)))
  424.      (setf (aref pts i) (funcall fx x1 x2 x3))
  425.      (setf (aref pts (f+ 1 i)) (funcall fy x1 x2 x3))
  426.      (setf (aref pts (f+ 2 i)) (funcall fz x1 x2 x3)))))
  427.   ))
  428.  
  429. (defun coerce-float-fun (expr &optional lvars)
  430.   (cond ((and (consp expr) (functionp expr))
  431.      expr)
  432.     ((and (symbolp expr) (not (member expr lvars)))
  433.      (cond ((fboundp expr) expr)
  434.            (t
  435.         (let* ((mexpr (mget expr 'mexpr))
  436.                (args (nth 1 mexpr)))
  437.           (or mexpr (merror "Undefined function ~a" expr))
  438.         (coerce `(lambda ,(cdr args)
  439.                  (declare (special ,@(cdr args)))
  440.                  (float ($realpart(meval* ',(nth 2 mexpr))) 1d0))
  441.               'function)))))
  442.     (t
  443.      (let ((vars (or lvars ($sort ($listofvars expr))))
  444.            ;(na (gensym "TMPF"))
  445.         )
  446.        (coerce `(lambda ,(cdr vars) (declare (special ,@(cdr vars)))
  447.             (float ($realpart (meval* ',expr)) 1d0))
  448.            'function)))))
  449.  
  450. (defmacro zval (points verts i) `(aref ,points (f+ 2 (f* 3 (aref ,verts ,i)))))
  451.  
  452. ;;sort the edges array so that drawing the edges will happen from the back towards
  453. ;; the front.   The if n==4 the edges array coming in looks like
  454. ;; v1 v2 v3 v4 0 w1 w2 w3 w4 0 ...
  455. ;; where vi,wi are indices pointint into the points array specifiying a point
  456. ;; in 3 space.   After the sorting is done, the 0 is filled in with the vertex
  457. ;; which is closer to us (ie highest z component after rotating towards the user)
  458. ;; and this is then they are sorted in groups of 5.   
  459. (defun sort-ngons (points edges n &aux lis )
  460.   (declare (type (array (long-float))  points)
  461.        (type (array (mod 65000)) edges)
  462.        (fixnum n))
  463.   (let ((new (make-array (length edges) :element-type  (array-element-type edges)))
  464.     (i 0)
  465.     (z 0.0)
  466.     (z1 0.0)
  467.     (n1 (- n 1))
  468.     (at 0)
  469.     (leng (length edges))
  470.     )
  471.     (declare (type (array (mod 65000)) new)
  472.          (fixnum i leng n1 at )
  473.          )
  474.     (declare (long-float z z1))
  475.     
  476.     (setq lis
  477.       (sloop  for i0 below leng by (+ n 1)
  478.           do 
  479.           (setq i i0)
  480.           (setq at 0)
  481.           (setq z (zval points edges i))
  482.           (setq i (+ i 1))
  483.           (sloop for j below n1
  484.              do (if (> (setq z1 (zval points edges i))  z)
  485.                 (setq z z1 at (aref edges i) ))
  486.              (setq i (+ i 1))
  487.              )
  488.           (setf (aref edges i) at)
  489.           collect (cons z i0)))
  490.     (setq lis (sortcar lis))
  491.     (setq i 0)
  492.     (sloop for v in lis
  493.        do
  494.        (sloop for j from (cdr v) 
  495.           for k to n
  496.           do (setf (aref new i) (aref edges j))
  497.           (incf i))
  498.        )
  499.     (copy-array-portion edges new  0 0 (length edges))
  500.     ))
  501.  
  502. (defun copy-array-portion (ar1 ar2 i1 i2 n1)
  503.  (declare (fixnum i1 i2 n1))
  504.  (sloop while (>= (setq n1 (- n1 1)) 0)
  505.         do (setf (aref ar1 i1) (aref ar2 i2))
  506.          (setq i1 (+ i1 1))
  507.         (setq i2 (+ i2 1))))
  508.  
  509.  
  510. (defun $concat_polygons (pl1 pl2 &aux tem new)
  511.   (setq new
  512.       (sloop for v in pl1 
  513.          for w in pl2
  514.          for l = (+ (length v) (length w))
  515.          do (setq tem (make-array l
  516.                       :element-type (array-element-type v)
  517.                       :fill-pointer  l
  518.                       )
  519.               )
  520.          collect tem))
  521.   (setq new (make-polygon (first new) (second new)) )
  522.  
  523.   (copy-array-portion (polygon-pts pl1) (polygon-pts new)
  524.               0 0 (length (polygon-pts pl1)))
  525.   (copy-array-portion (polygon-pts pl2) (polygon-pts new)
  526.               (length (polygon-pts pl1))
  527.               0 (length (polygon-pts pl2)))
  528.   (copy-array-portion (polygon-edges pl1) (polygon-edges new)
  529.               0 0 (length (polygon-edges pl1)))
  530.   (sloop for i from (length (polygon-edges pl1))
  531.      for j from 0 below (length (polygon-edges pl2))
  532.      with  lpts1  =  (length (polygon-pts pl1))
  533.      with ar2   =  (polygon-edges pl2)
  534.      with arnew =  (polygon-edges new)
  535.      do (setf (aref arnew i) (f+ lpts1 (aref ar2 j)))))
  536.  
  537. (defun $copy_pts(lis vec start)
  538.   (declare (fixnum start))
  539.   (let ((tem vec))
  540.     (declare (type (array long-float) tem))
  541.     (cond ((numberp lis)
  542.        (or (typep lis 'long-float) (setq lis (float lis 0.0)))
  543.        (setf (aref tem start) lis)
  544.        
  545.        (+ start 1))
  546.       ((typep lis 'cons)
  547.        ($copy_pts (cdr lis) vec  ($copy_pts (car lis) vec start)))
  548.  
  549.       ((symbolp lis) start)
  550.       (t (error "bad lis")))))
  551.   
  552.  
  553. ;; parametric ; [parametric,xfun,yfun,[t,tlow,thigh],[nticks ..]]
  554. ;; the rest of the parametric list after the list will be pushed plot_options
  555.  
  556. (defun draw2d-parametric (param range1 &aux range tem)
  557.   (cond ((and ($listp (setq tem (nth 4 param)))
  558.           (symbolp (cadr tem))
  559.           (eql ($length tem) 3)
  560.           (<= (length (symbol-name (cadr tem))) 2))
  561.           ;; sure looks like a range
  562.           (setq range tem)))
  563.   (let* (($plot_options ($append ($rest param 3)
  564.                  (if range1
  565.                      ($cons range1 $plot_options)
  566.                    $plot_options)))
  567.      (nticks (nth 2($get_plot_option '$nticks)))
  568.      (trange (or range ($get_plot_option '|$t|)))
  569.      (xrange ($get_plot_option '|$x|))
  570.      (yrange ($get_plot_option '|$y|))
  571.      ($numer t)
  572.      (tmin (coerce-float (nth 2 trange)))
  573.      (tmax (coerce-float (nth 3 trange)))
  574.      (xmin (coerce-float (nth 2 xrange)))
  575.      (xmax (coerce-float (nth 3 xrange)))
  576.      (ymin (coerce-float (nth 2 yrange)))
  577.      (ymax (coerce-float (nth 3 yrange)))
  578.      (x 0.0) ; have to initialize to some floating point..
  579.      (y 0.0)
  580.      (tt tmin)
  581.      (eps (/ (- tmax tmin) (- nticks 1)))
  582.      f1 f2 in-range-y in-range-x in-range last-ok 
  583.      )
  584.     (declare (long-float x y tt ymin ymax xmin xmax tmin tmax eps))
  585.     (setq f1 (coerce-float-fun (nth 2 param) `((mlist), (nth 1 trange))))
  586.     (setq f2 (coerce-float-fun (nth 3 param) `((mlist), (nth 1 trange))))
  587.    (cons '(mlist simp)    
  588.     (sloop 
  589.        do 
  590.        (setq x (funcall f1 tt))
  591.        (setq y (funcall f2 tt))
  592.        (setq in-range-y (and (<= y ymax) (>= y ymin)))
  593.        (setq in-range-x  (and  (<= x xmax) (>= x xmin)))
  594.        (setq in-range (and in-range-x in-range-y))
  595.        when (and (not in-range) (not last-ok))
  596.        collect  'moveto and collect 'moveto
  597.        do
  598.        (setq last-ok in-range)
  599.        collect (if in-range-x x (if (> x xmax) xmax xmin))
  600.        collect (if in-range-y y (if (> y ymax) ymax ymin))
  601.        when (>= tt tmax) do (sloop::loop-finish)
  602.        do (setq tt (+ tt eps))
  603.        (if (>= tt tmax) (setq tt tmax))
  604.        )))
  605. )
  606.  
  607. ;; arrange so that the list of points x0,y0,x1,y1,.. on the curve
  608. ;; never have abs(y1-y0 ) and (x1-x0) <= deltax
  609.  
  610. (defun draw2d (f range )
  611.   (if (and ($listp f) (equal '$parametric (cadr f)))
  612.       (return-from draw2d (draw2d-parametric f range)))
  613.   (let* ((nticks (nth 2($get_plot_option '$nticks)))
  614.      (yrange ($get_plot_option '|$y|))
  615.      ($numer t)
  616.      )
  617.  
  618.     (setq f (coerce-float-fun f `((mlist), (nth 1 range))))
  619.  
  620.     (let* ((x (coerce-float (nth 2 range)))
  621.        (xend (coerce-float (nth 3 range)))
  622.        (ymin (coerce-float (nth 2 yrange)))
  623.        (ymax (coerce-float (nth 3 yrange)))
  624.        (eps ($/ (- xend x) (coerce-float nticks)))
  625.        (x1 0.0)
  626.        (y1 0.0)
  627.        (y (funcall f x))
  628.        (dy 0.0)
  629.        (epsy ($/ (- ymax ymin) 1.0))
  630.        (eps2 (* eps eps))
  631.        in-range last-ok
  632.        )
  633.       (declare (long-float x1 y1 x y dy eps2 eps ymin ymax ))
  634.                     ;(print (list 'ymin ymin 'ymax ymax epsy))
  635.       (setq x ($- x eps))  
  636.       (cons '(mlist)
  637.         (sloop   do
  638.              (setq x1 ($+ eps x))
  639.              (setq y1 (funcall f x1))
  640.              (setq in-range (and (<= y1 ymax) (>= y1 ymin)))
  641.              (cond (in-range
  642.                 (setq dy (- y1 y))
  643.                 (cond ((< dy 0) (setq dy (- dy))))
  644.                 (cond ((> dy eps)
  645.                    (setq x1 (+ x (/ eps2 dy)))
  646.                    (setq y1 (funcall f x1))
  647.                    (setq in-range (and (<= y1 ymax) (>= y1 ymin)))
  648.                    (or in-range (setq x1 (+ eps x)))
  649.                    )
  650.                   ))
  651.                )
  652.              (setq x x1)
  653.              (setq y y1)
  654.              when (or (and (not last-ok)
  655.                    (not in-range))
  656.                   (> dy epsy))
  657.  
  658.              collect 'moveto and collect 'moveto
  659.              do
  660.              (setq last-ok in-range)
  661.              collect x1 
  662.              collect (if in-range y1 (if (> y1 ymax) ymax ymin))
  663.              when (>= x xend)
  664.              collect xend and
  665.              collect (let ((tem (funcall f xend)))
  666.                    (if (>= tem ymax) ymax (if (<= tem ymin) ymin tem)))
  667.              and do (sloop::loop-finish))))))
  668.  
  669. (defun get-range (lis)
  670.   (let ((ymin most-positive-long-float)
  671.     (ymax most-negative-long-float))
  672.     (declare (long-float ymin ymax))
  673.     (do ((l lis (cddr l)))
  674.     ((null l))
  675.     (or (floatp (car l)) (setf (car l) (float (car l) #. (coerce 2 'long-float))))
  676.       (cond ((float-<   (car l)ymin)
  677.          (setq ymin (car l))))
  678.       (cond ((float-<  ymax  (car l))
  679.           (setq ymax (car l)))))
  680.     (list '(mlist) ymin ymax)))
  681.  
  682. (defvar $window_size '((mlist)
  683.                #.(* 8.5 72) #. (* 11 72) ))
  684.  
  685. (defun $getrange(x &optional xrange &aux yrange)
  686.   (setq yrange  (cdr (get-range (cddr x))))
  687.   (or xrange (setq xrange (cdr (get-range (cdr x)))))
  688.   (setup-for-ps-range xrange yrange nil))
  689.  
  690. (defun $paramplot (f g range &optional (delta .1 supplied) &aux pts ($numer t)
  691.              )
  692.   (setq f (coerce-float-fun f))
  693.   (setq g (coerce-float-fun g))
  694.   (setq range (meval* range))
  695.   (or supplied (setq delta (/ (- (nth   2 range) (nth 1 range)) (nth 2 ($get_plot_option '$nticks)))))
  696.   (setq pts(cons '(Mlist)
  697.          (sloop with tt = (coerce-float (nth 1 range))
  698.             with end = (coerce-float (nth 2 range))
  699.             while (float-< tt end)
  700.             collect (funcall f tt)
  701.             collect (funcall g tt)
  702.             do (setq tt ($+ tt delta)))))
  703.   ($closeps)
  704.   ($getrange pts)
  705.   ($psdraw_curve pts)
  706.   ($closeps)
  707.   ($viewps))
  708.   
  709. (defun $plot2d_ps(fun range &rest options &aux ($numer t)  $display2d
  710.               ($plot_options $plot_options))
  711.   (dolist (v options) ($set_plot_option v))
  712.   (setq range (check-range range))
  713.   (let ((tem (draw2d fun range )))
  714.     ($closeps)
  715.     ($getrange tem (cddr range))
  716.     ($psdraw_curve tem)
  717.     ($psaxes ($rest range))
  718.     (p "showpage")
  719.     ($viewps)))
  720.  
  721.  
  722. (defvar $gnuplot_command "mgnuplot")
  723. (defvar $geomview_command "geomview maxout.geomview")
  724.  
  725. (defvar $openmath_plot_command "omplotdata")
  726.  
  727. (defun $plot2d(fun &optional range &rest options &aux ($numer t) $display2d
  728.                           (i 0) plot-format file plot-name
  729.               ($plot_options $plot_options))
  730.   (dolist (v options) ($set_plot_option v))
  731.   (or ($listp fun ) (setf fun `((mlist) ,fun)))
  732.   (cond ((eq (cadr fun) '$parametric)
  733.      (or range (setq range (nth 4 fun)))
  734.      (setf fun `((mlist) ,fun))))
  735.   (cond ((eq ($get_plot_option '$PLOT_FORMAT 2) '$openmath)
  736.          (return-from $plot2d (apply '$plot2dOpen fun range options))))
  737.   (check-range range)
  738.   (setq plot-format  ($get_plot_option '$plot_format 2))
  739.   (setq file (format nil "maxout.~(~a~)" (stripdollar plot-format)))
  740.   
  741.   (with-open-file (st  file :direction :output)  
  742.     (dolist (v (cdr fun))
  743.         (incf i)
  744.         (setq plot-name                (let ((string (coerce (mstring v) 'string)))
  745.              (cond ((< (length string) 9) string)
  746.                    (t (format nil "Fun~a" i)))))
  747.     (case plot-format
  748.           ($xgraph
  749.            (format st "~%~% \"~a\"~%" plot-name))
  750.           ($gnuplot
  751.            (format st "~%~%# \"~a\"~%" plot-name))
  752.            )
  753.       (sloop for (v w) on (cdr (draw2d v range )) by 'cddr
  754.      do
  755.      (cond ((eq v 'moveto)
  756.         (unless (equal plot-format '$gnuplot)
  757.             (format st "move ")))
  758.            (t  (format st "~,3f ~,3f ~%" v w))))))
  759.   (case plot-format
  760.     ($gnuplot 
  761.      ($system (concatenate 'string *maxima-plotdir* "/" $gnuplot_command) " -plot2d maxout.gnuplot -title '" plot-name "'"))
  762.     ($xgraph
  763.      ($system "xgraph -t 'Maxima Plot' < maxout.xgraph &"))
  764.     ))
  765.  
  766. (defun maxima-bin-search (command)
  767.   (or ($file_search command
  768.             `((mlist) , (maxima-path "bin" "###")))
  769.          command))
  770.     
  771.  
  772.  
  773. (defun $plot2dOpen(fun range &rest options &aux ($numer t)  $display2d
  774.                           (i 0)
  775.               ($plot_options $plot_options)
  776.     )
  777.   (declare (special linel))
  778.   (dolist (v options) ($set_plot_option v))
  779.   (setq range (check-range range))
  780.   (or ($listp fun ) (setf fun `((mlist) ,fun)))
  781.   (show-open-plot
  782.    (with-output-to-string
  783.     (st )
  784.     (format st "~%{plot2d ~%")
  785.     (sloop for f in (cdr fun)
  786.        do
  787.        (incf i)
  788.        (format st " {label \"~a\"}~% "
  789.            (let ((string (coerce (mstring f) 'string)))
  790.              (cond ((< (length string) 9) string)
  791.                (t (format nil "Fun~a" i))))
  792.            )
  793.        (format st "{xversusy ~%")
  794.        (let ((lis (cdr (draw2d f range ))))
  795.  
  796.          (sloop while lis
  797.             do
  798.  
  799.             (sloop while (and lis (not (eq (car lis) 'moveto)))
  800.  
  801.                collecting (car lis) into xx
  802.                collecting (cadr lis) into yy
  803.                do (setq lis (cddr lis))
  804.                finally
  805.                ;; only output if at least two points for line
  806.                (cond ((cdr xx)
  807.                   (tcl-output-list st xx)
  808.                   (tcl-output-list st yy)
  809.                   ))
  810.                ; remove the moveto
  811.                )
  812.             (setq lis (cddr lis))
  813.             ))
  814.        (format st "~%}"))
  815.     (format st "~%} ~%"))))
  816.  
  817.  
  818. (defvar $In_Netmath nil)
  819. (eval-when (load)
  820.   (cond ($In_Netmath
  821.      (setf (symbol-function '$plot2d) (symbol-function '$plot2dopen))
  822.      (setf $show_openplot nil)
  823.     )))
  824.  
  825. (defun $sprint(&rest args )
  826.   (sloop for v in args do
  827.      (cond ((symbolp v)
  828.         (setq v (string-left-trim '(#\$ #\&) (symbol-name v))))
  829.            ((numberp v) v)
  830.            (t (setq v (implode (strgrind v)))))
  831.      (princ v)
  832.      (princ " ")
  833.   )
  834.   (car args)
  835.   )
  836.  
  837. (defun $show_file(file)
  838.   (princ (file-to-string ($file_search file)))
  839.   '$done)
  840.  
  841.  
  842. (defun $tcl_output  (lis i &optional (skip 2))
  843.     (or (typep i 'fixnum) (error "~a should be an integer" i ))
  844.     ($listp_check 'list lis )
  845.     (format *standard-output* "~% {")
  846.     (cond (($listp (nth 1 lis))
  847.        (sloop for v in lis
  848.           do
  849.           (format *standard-output* "~,10f " (nth i v)))
  850.        )
  851.       (t
  852.        (setq lis (nthcdr i lis))
  853.        (sloop  with v = lis  while v
  854.           do
  855.           (format *standard-output* "~,10f " (car v))
  856.           (setq v (nthcdr skip v))
  857.           )
  858.        ))
  859.     (format *standard-output* "~% }")
  860.     )
  861.  
  862.            
  863.           
  864.        
  865.                 
  866.                 
  867.        
  868.  
  869.            
  870.  
  871.  
  872.  
  873. (defun tcl-output-list ( st lis )
  874.   (cond ((null lis) )
  875.     ((atom (car lis))
  876.      (princ " {  " st)
  877.      (sloop for v in lis
  878.         count t into n
  879.         when (eql 0 (mod n 5))
  880.         do (terpri st)
  881.         do
  882.         (format st "~,10f " v))
  883.      (format st  " }~%"))
  884.     (t (tcl-output-list st (car lis))
  885.        (tcl-output-list st (cdr lis)))))
  886.  
  887.  
  888.  
  889.  
  890. (defun $openplot_curves (lis &aux (linel 100000))
  891.   (declare (special linel))
  892.   (show-open-plot
  893.    (with-output-to-string
  894.     (st )
  895.     (format st "{plot2d ~%")
  896.     (or (and ($listp lis) ($listp (nth 1 lis)))
  897.     (merror "Need a list of curves, [[x1,y1,x2,y2,...],[u1,v1,u2,v2,...]] or [[[x1,y1],[x2,y2],...]"))
  898.     
  899.     (sloop for v in (cdr lis)
  900.        do
  901.        (or ($listp v) (merror "should be a list"))
  902.        (setq v (cdr v))
  903.        (format st "~%~%")
  904.        (sloop while (and (car v) (symbolp (car v)))
  905.           do   (mformat st "~M~%" ($concat (car v)))
  906.           (setq v (cdr v))
  907.           )
  908.            when v
  909.        do    
  910.        (format st "~%{ xversusy  ")
  911.        (sloop while v
  912.           with this with xvals with yvals
  913.           do
  914.           (cond   ((numberp (car v))
  915.                (setq this v) (setq v (cddr v))
  916.                (push (car this) xvals)
  917.                (push (second this) yvals)               
  918.                )
  919.               (($listp (car v))
  920.                (setq this (cdar v))
  921.                (push (car this) xvals)
  922.                (push (second this) yvals)
  923.                (and (third this) (push (third this ) yvals))
  924.                (setq v (cdr v)))
  925.               (t (princ (list 'bad (car v))) (setq v (cdr v))))
  926.  
  927.           finally
  928.           (tcl-output-list st (nreverse xvals))
  929.           (tcl-output-list st (nreverse yvals)))
  930.               (format st " }")
  931.        )
  932.     (format st "}~%")
  933.     )))
  934.  
  935.  
  936.  
  937. (defun $xgraph_curves (lis &rest options &aux w)
  938.   options
  939.   (with-open-file (st  "xgraph-out" :direction :output)
  940.     (format st "=600x600~%")
  941.     (sloop for v in (cdr lis)
  942.        do
  943.        (setq v (cdr v))
  944.        (format st "~%~%")
  945.        (sloop while v
  946.       do
  947.       (cond
  948.        ((symbolp (car v))
  949.         (mformat st "~M~%" ($concat (car v)))
  950.         (setq v (cdr v)))
  951.        (t (cond       ((numberp (car v))
  952.         (setq w v) (setq v (cddr v)))
  953.        (($listp (car v))
  954.         (setq w (cdar v))
  955.         (setq v (cdr v))))
  956.       (format st "~,3f ~,3f ~%" (car w) (second w)))))))
  957.   ($system "xgraph -t 'Maxima Plot' < xgraph-out &"))
  958.  
  959.  
  960.      
  961.  
  962. (defun average-slope (m1 m2)
  963.   (tan ($/ ($+ (lisp::atan m1) (lisp::atan m2)) 2.0)))
  964.  
  965. (defun slope (x1 y1 x2 y2 &aux (del ($- x2 x1)))
  966.   (declare (long-float x1 y1 x2 y2 del))
  967.   (cond ((eql del 0.0)
  968.      #. (expt 10 30))
  969.     (t ($/ ($- y2 y1) del))))
  970.        
  971.  
  972. ;;When we initialize we move the origin to the middle of $window_size
  973. ;;Then to offset from that use translate.
  974. (defvar $ps_translate '((mlist) 0 0))
  975.  
  976. ;; initially 1/72 of inch is the scale
  977. (defvar $ps_scale '((mlist) 72  72))
  978.  
  979.  
  980. (defun $pscom (&rest l)
  981.   (apply 'p l))
  982.  
  983. ;-10 to 10
  984. (defun psx (x) x)
  985. (defun psy (y) y)
  986.  
  987. (defun p (&rest l)
  988.   (assureps)
  989.   (sloop for v in l do
  990.      (if (symbolp v) (setq v (maxima-string v)))
  991.     ; (if (numberp v) (setq v (* 70 v)))
  992.      (cond ((consp v)
  993.         (sloop for w in (cdr v) do (p w)))
  994.            ((floatp v) (format $pstream "~,4f" v))
  995.            (t(princ v $pstream)))
  996.       (princ " " $pstream))
  997.   (terpri $pstream))
  998.  
  999. (defun psapply (f lis)
  1000.   (if ($listp lis) (setq lis (cdr lis)))
  1001.   (apply 'p lis)
  1002.   (p f))
  1003.  
  1004. (defun $moveto (x y)
  1005.   (p (psx x) (psy y) "moveto ")) 
  1006.  
  1007. (defun $join (x y)
  1008.   (cons '(mlist)(sloop for w in (cdr x) for u in (cdr y)
  1009.                collect w collect u)))
  1010.  
  1011. (defun $psline (a b c d)
  1012.   (p (psx a) (psy b) "moveto ")
  1013.   (p  (psx c) (psy d) "lineto"))
  1014.  
  1015. (defun setup-for-ps-range (xrange yrange do-prolog)
  1016.     (let* ((cy (/ (+ (nth 1 yrange) (nth 0 yrange)) 2.0))
  1017.        (CX (/ (+ (nth 1 xrange) (nth 0 xrange)) 2.0))
  1018.        (scaley (/ (nth 2 $window_size)
  1019.               (* 1.2  (- (nth 1 yrange) (nth 0 yrange)))))
  1020.        (scalex (/ (nth 1 $window_size)
  1021.               (* 1.2 (- (nth 1 xrange) (nth 0 xrange)))))
  1022.        ($ps_scale
  1023.          (progn
  1024.            (cond ((< scalex scaley)
  1025.                (setq scaley scalex)))
  1026.             `((mlist) , scaley ,scaley)))
  1027.        ($ps_translate `((mlist) , CX ,CY)))
  1028.       (ASSUREPS do-prolog)))
  1029.  
  1030. (defun assureps (&optional do-prolog)
  1031.   (cond ((streamp $pstream))
  1032.     (t (setq do-prolog t)))
  1033.   (or $pstream (setq $pstream (open "maxout.ps" :direction :output)))
  1034.   (cond (do-prolog
  1035.       (p "%!")
  1036.       (p  (* .5 (nth 1 $window_size))
  1037.          (* .5 (nth 2 $window_size))
  1038.          "translate"
  1039.          )
  1040.       (psapply "scale" $ps_scale)
  1041.       (p  (- (nth 1 $ps_translate))
  1042.          (- (nth 2 $ps_translate))
  1043.          "translate"
  1044.          )
  1045.       (p " 1.5 " (second $ps_scale) "div setlinewidth
  1046. /Helvetica findfont 14 " (second $ps_scale) " div scalefont setfont
  1047. /dotradius .05 def
  1048. /drawdot {
  1049.  /yy exch def
  1050.  /xx exch def
  1051.   gsave
  1052.   xx yy dotradius 0 360 arc
  1053.   fill
  1054.   grestore
  1055. }def
  1056.  
  1057. /ticklength  .03 def
  1058. /axiswidth  .01 def
  1059. /drawtick {
  1060.  /yy exch def
  1061.  /xx exch def
  1062.   gsave
  1063.   xx ticklength sub yy moveto
  1064.   xx ticklength add yy lineto
  1065.   stroke    
  1066.   xx yy  ticklength sub  moveto
  1067.   xx yy  ticklength add  lineto
  1068.   stroke
  1069.   grestore
  1070. } def
  1071. "))))
  1072.  
  1073.  
  1074. (defun $closeps ()
  1075.   (prog1
  1076.       (when (and (streamp $pstream)
  1077.          )
  1078.             (p "showpage")
  1079.         (close  $pstream))
  1080.     (setq $pstream nil)))
  1081.     
  1082. (defun ps-fixup-points(lis)
  1083.   (assert ($listp lis))
  1084.   (setq lis (cdr lis))
  1085.   (cond ((numberp (car lis)))
  1086.     ((and ($listp (car lis)) (numberp (nth 1 (car lis))))
  1087.      (setq lis (sloop for w in lis collect
  1088.               (nth 1 w)
  1089.               collect (nth 2 w))))
  1090.     (t (error
  1091.         "pscurve:Neither [x0,y0,x1,y1,...] nor [[x0,y0],[x1,y1],...]")))
  1092.   lis)
  1093.  
  1094. (defun $psdraw_curve (lis &aux (n 0))
  1095.   (declare (fixnum n))
  1096.   (setq lis (ps-fixup-points lis))
  1097.   (p "newpath" (nth 0 lis) (nth 1 lis) "moveto")
  1098.   (sloop while lis with second
  1099.      do
  1100.      (cond ((eq (car lis) 'moveto)
  1101.         (or (eql n 0) (p "stroke"))
  1102.         (setq n 0) (setq lis (cddr lis))))
  1103.      (or (setq second (cadr lis)) (error "odd length list of points"))
  1104.      (cond ((eql n 0)
  1105.         (p (car lis) second "moveto"))
  1106.        (t  (p (car lis) second "lineto")
  1107.            ))
  1108.      (setq n (+ n 1))
  1109.      (cond ((eql 0 (mod n 20))
  1110.         (p "stroke")
  1111.         (setq n 0))
  1112.        (t (setq lis (cddr lis)))))
  1113.   (or (eql n 0) (p "stroke")))
  1114.  
  1115. (defvar $viewps_command  "(ghostview  ~a)")
  1116.  
  1117. (defvar $viewps_command  "(gs -I. -Q  ~a)")
  1118.  
  1119. (defvar  $viewps_command "echo /def /show {pop} def |  cat - ~a | x11ps")
  1120.  
  1121. ;; If your gs has custom features for understanding mouse clicks
  1122.  
  1123.  
  1124. ;;Your gs will loop for ever if you don't have showpage at the end of it!!
  1125. (defvar $viewps_command   "echo '/showpage { .copypage readmouseclick /ke exch def ke 1 eq { erasepage initgraphics} {ke 5 ne {quit} if} ifelse} def  {(~a) run } loop' | gs  -title 'Maxima  (click left to exit,middle to redraw)' > /dev/null 2>/dev/null &")
  1126.  
  1127.     ;; allow this to be set in a system init file (sys-init.lsp)
  1128.  
  1129.  
  1130. (defun $viewps ( &optional file)
  1131.   (cond  ((and (streamp $pstream))
  1132.       ($pscom "showpage")
  1133.       (force-output $pstream)))
  1134.   (cond (file (setq file (maxima-string file)))
  1135.     (t(setq file "maxout.ps")
  1136.      
  1137.      (if (and (streamp $pstream))
  1138.          (force-output $pstream))))
  1139.   (if (equal $viewps_command "(gs -I. -Q  ~a)")
  1140.         (format t "~%type `quit' to exit back to affine or maxima
  1141.   To reprint a page do
  1142.   GS>showpage
  1143.   GS>(maxout.ps)run
  1144.   GS> -150 -150 translate 1.2 1.2 scale (maxout.ps)run
  1145.   would print it moved 150/72 inches to left, and down, and scaled by 1.2 times
  1146.   showpage clears scaling."))
  1147.  
  1148.   ($system    (format nil $viewps_command file)))
  1149.  
  1150. (defun $chkpt (a)
  1151.   (or (and ($listp a)
  1152.        (numberp (nth 1 a))
  1153.        (numberp (nth 2 a)))
  1154.       (error "illegal pt ~a" a))
  1155.   t)
  1156.        
  1157. (defvar $pslineno nil)
  1158. (defun $psdrawline (a &optional b c d)
  1159.   (cond ((null b)
  1160.      (assert (and ($listp a)
  1161.               ($chkpt (nth 1 a))))
  1162.      (setq b (nth 2 a))
  1163.      (setq a (nth 1 a))))
  1164.   (cond ((null c)
  1165.      ($chkpt b)
  1166.      (setq c (nth 1 b))
  1167.      (setq d (nth 2 b))))
  1168.   (cond (($listp a)
  1169.      (setq b (nth 2 a) a (nth 1 a))))
  1170.   (cond ((null c)
  1171.      (setq c (nth 1 b))
  1172.      (setq c (nth 2 b))))
  1173.  
  1174.   (when $pslineno
  1175.       (incf $pslineno)
  1176.       ($pslabelline a b c d $pslineno))
  1177.       
  1178.   (p a b "moveto")
  1179.   (p c d "lineto")
  1180.   (p "stroke"))
  1181.  
  1182. (defun $pslabelline (a b c d $pslineno)
  1183.    (p (/ (+ a c) 2)
  1184.       (/ (+ b d) 2)
  1185.       "moveto"
  1186.       (format nil "(<--L~a)show" $pslineno)))
  1187.  
  1188. (defun $sort_polys (lis)
  1189.   (let ((tem
  1190.      (sloop for v in (cdr lis)
  1191.         collect (cons
  1192.              (sloop for w in (cdr v)
  1193.                 maximize (nth 3 w))
  1194.              v))))
  1195.     (print 'next)
  1196.     (cons '(mlist) (mapcar 'cdr (sortcar tem '<)))))
  1197.  
  1198. (defun $drawpoly (x)
  1199.   (p "gsave")
  1200.   (p (cdadr x) "moveto")
  1201.   (setq x (cddr x))
  1202.   (sloop for edge in x
  1203.      do (p (cdr edge) "lineto")
  1204.      finally (p "1 setgray fill"))
  1205.   ($psdrawline x)
  1206.   (p "grestore"))
  1207.  
  1208. (defun $psdrawpolys (polys)
  1209.   (dolist (v (cdr polys))
  1210.       ($drawpoly v)))
  1211.       
  1212. ;; draw the axes through $psdef
  1213. (defun $psaxes (leng &optional (lengy leng))
  1214.   (p "gsave axiswidth setlinewidth")
  1215.   (let (begx begy endx endy)
  1216.     (cond ((numberp leng)
  1217.        (setq begx (- leng))
  1218.        (setq endx leng))
  1219.       (t (setq begx (nth 1 leng))
  1220.          (setq endx (nth 2 leng))))
  1221.     (cond ((numberp lengy)
  1222.        (setq begy (- lengy))
  1223.        (setq endy lengy))
  1224.       (t (setq begy (nth 1 lengy))
  1225.          (setq endy (nth 2 lengy))))
  1226.  
  1227.     (sloop for i from (floor begx) below (ceiling endx)
  1228.        do 
  1229.        ($psdrawline i 0 (+ i 1) 0)
  1230.        (p i 0 "drawtick")
  1231.        (p (+ i 1) 0 "drawtick"))
  1232.  
  1233.     (sloop for i from (floor begy) below (ceiling endy)
  1234.        do
  1235.        ($psdrawline 0 i 0 (+ i 1) )
  1236.        (p 0 i "drawtick")
  1237.        (p 0 (+ i 1)  "drawtick"))
  1238.  
  1239.   
  1240.     (p "grestore")))
  1241.  
  1242. (defun $psdraw_points(lis)
  1243.     (assert (and ($listp lis)
  1244.          ($listp (cadr lis))))
  1245.     (sloop for w in (cdr lis)
  1246.        do (p (nth 1 w)
  1247.          (nth 2 w)
  1248.          "drawdot")))
  1249.  
  1250.  
  1251. (defun $view_zic ()
  1252.   (let ((izdir (getenv "IZICDIR")))
  1253.     (or (probe-file
  1254.      (format nil "~a/tcl-files/maxima.tcl" izdir))
  1255.     (error
  1256.      "could not find file ~a :  Set environment variable IZICDIR" izdir))
  1257.     ($system "izic -interface ${IZICDIR}/tcl-files/maxima.tcl  1> /dev/null &")))
  1258.  
  1259. (defun $isend (x)
  1260.   ($system (format nil " izic -app izic -cmd '{~a}'" (string-trim '(#\&) (string x)))))
  1261.  
  1262.  
  1263.  
  1264. (defvar *some-colours*
  1265.   ;; from rgb.txt
  1266.   '(135 206 250        LightSkyBlue
  1267.     70 130 180        SteelBlue
  1268.     205  92  92        IndianRed
  1269.     178  34  34        firebrick
  1270.     176  48  96        maroon
  1271.     221 160 221        plum
  1272.     238 130 238        violet))
  1273.  
  1274.  
  1275. ;; one of $zic, $geomview, $ps
  1276.  
  1277. (defun plot-zic-colors (&aux (ncolors  (/ (length *some-colours*) 4))) 
  1278.   (format $pstream "couleurs ~% ~a ~% " ncolors  )
  1279.   (sloop for v in *some-colours*
  1280.      with do-ind = t with ind = -1 
  1281.      do
  1282.      (cond (do-ind
  1283.         (format $pstream "~a " (incf ind))
  1284.         (setq do-ind nil)
  1285.         ))
  1286.      (cond ((Numberp v)
  1287.         (print-pt (/ v 256.0)))
  1288.        (t 
  1289.         (setq do-ind t)
  1290.         (format $pstream "~%# ~(~a~)~%" v))))
  1291.   (format $pstream " 1
  1292. 0
  1293. ~a 0.801 0.359 0.359 
  1294. 128 .8 1 0
  1295. " ncolors))
  1296.  
  1297. (defun check-range (range &aux ($numer t) tem a b)
  1298.   (or (and ($listp range)
  1299.        (setq tem (cdr range))
  1300.        (symbolp (car tem))
  1301.        (numberp (setq a (meval* (second tem))))
  1302.        (numberp (setq b (meval* (third tem))))
  1303.        (< a b)
  1304.        )
  1305.       (merror "Bad Range ~%~M must be of the form [variable,min,max]" range))
  1306.   `((mlist) ,(car tem) ,(float a) ,(float b)))
  1307.  
  1308. (defun $zero_fun (x y) x y 0.0)
  1309.  
  1310. (defun output-points (pl &optional m)
  1311.   "If m is supplied print blank line every m lines"
  1312.   (let ((j -1))
  1313.     (declare (fixnum j))
  1314.   (sloop for i below (length (polygon-pts pl))
  1315.      with ar = (polygon-pts pl)
  1316.      do (print-pt (aref ar i))
  1317.      (setq i (+ i 1))
  1318.      (print-pt (aref ar i))
  1319.      (setq i (+ i 1))
  1320.      (print-pt (aref ar i))
  1321.      (terpri $pstream)
  1322.      (cond (m
  1323.         (setq j (+ j 1))
  1324.         (cond ((eql j (the fixnum m))
  1325.            (terpri $pstream)
  1326.            (setq j -1)))))
  1327.      )))
  1328.  
  1329. (defvar $show_openplot t)
  1330. (defun show-open-plot (ans)
  1331.   (cond ($show_openplot
  1332.      (with-open-file (st1 "maxout.openmath" :direction :output)
  1333.             (princ  ans st1))
  1334.      ($system (concatenate 'string *maxima-plotdir* "/" $openmath_plot_command) " maxout.openmath" ))
  1335.     (t (princ ans) "")))
  1336.  
  1337. (defun $plot3d ( fun &optional (xrange ($get_plot_option '$x))
  1338.              (yrange ($get_plot_option '|$y|) y-supplied)
  1339.              &rest options 
  1340.              &aux lvars trans *original-points*
  1341.              ($plot_options $plot_options)
  1342.              ($in_netmath $in_netmath)
  1343.              colour-z grid
  1344.              plot-format
  1345.              )
  1346.   (declare (special *original-points*))
  1347.   (cond (options
  1348.      (dolist (v options)
  1349.        ($set_plot_option v))))
  1350.   (setq plot-format  ($get_plot_option '$plot_format 2))
  1351.   (and $in_netmath (setq $in_netmath (eq plot-format '$openmath)))
  1352.   (setq xrange (check-range xrange))
  1353.   (setq yrange (check-range yrange))
  1354.   (cond ((not y-supplied)
  1355.      (let ((vars ($sort ($listofvars fun))))
  1356.        (or (eql ($length vars) 2)
  1357.            (merror "Please supply the range for variables eg [x,-3,3],[y,-3,4]"))
  1358.        (setq xrange ($cons (nth 1 vars) ($rest xrange)))
  1359.        (setq yrange ($cons (nth 2 vars) ($rest yrange))))))
  1360.   (setq grid ($get_plot_option '$grid))
  1361.   (setq colour-z  ($get_plot_option '$colour_z 2))
  1362.   (setq lvars `((mlist),(nth 1 xrange) ,(nth 1 yrange)))
  1363.   (cond (($listp fun)
  1364.      (or (eql 3 ($length fun)) (merror "List ~M is not of length 3" fun))
  1365.      (setq trans ($make_transform
  1366.               ($append lvars '((mlist) $z))
  1367.               (nth 1 fun)
  1368.               (nth 2 fun)
  1369.               (nth 3 fun)))
  1370.      (setq fun '$zero_fun))
  1371.     (t (setq fun (coerce-float-fun fun lvars))))
  1372.     (let* ((pl (draw3d fun
  1373.                (nth 2 xrange)
  1374.                (nth 3 xrange)
  1375.                (nth 2 yrange)
  1376.                (nth 3 yrange)
  1377.                (nth 2 grid)
  1378.                (nth 3 grid)))
  1379.        (ar (polygon-pts pl)) tem
  1380.        )
  1381.       (declare (type (array long-float) ar))
  1382.  
  1383.       (if trans  (mfuncall trans ar))
  1384.       (if (setq tem  ($get_plot_option '$transform_xy 2)) (mfuncall tem ar))
  1385.       ;; compute bounding box.
  1386.   (let (($pstream
  1387.      (cond ($in_netmath *standard-output*)
  1388.            (t (open (format nil "maxout.~(~a~)" (stripdollar plot-format))
  1389.             :direction :output)))))
  1390.     (unwind-protect
  1391.       (case  plot-format
  1392.     ($zic
  1393.      (let ((x-range ($get_range ar 0))
  1394.            (y-range ($get_range ar 1))
  1395.            (z-range ($get_range ar 2)))
  1396.        (plot-zic-colors)
  1397.        (format $pstream "domaine ~a ~a ~a ~a ~a ~a ~%"
  1398.            (nth 0 x-range)
  1399.            (nth 1 x-range)
  1400.            (nth 0 y-range)
  1401.            (nth 1 y-range)
  1402.            (nth 0 z-range)
  1403.            (nth 1 z-range))
  1404.        (format $pstream "surface ~a ~a ~%"
  1405.            (+ 1 (nth 3 grid))
  1406.            (+ 1 (nth 2 grid))
  1407.            )
  1408.        (output-points pl nil)))
  1409.     ($gnuplot
  1410.        (output-points pl (nth 2 grid)))
  1411.     ($openmath
  1412.      (progn
  1413.        (format $pstream "{plot3d {matrix_mesh ~%")
  1414.        ;; we do the x y z  separately:
  1415.        (sloop for off from 0 to 2
  1416.           with ar = (polygon-pts pl)
  1417.           with  i = 0
  1418.           declare (fixnum i)
  1419.           do (setq i off)
  1420.           (format $pstream "~%{")
  1421.           (sloop 
  1422.            while (< i (length ar))
  1423.            do (format $pstream "~% {")
  1424.            (sloop for j to (nth 2 grid)
  1425.              do (print-pt (aref ar i))
  1426.              (setq i (+ i 3)))
  1427.            (format $pstream "} ")
  1428.            )
  1429.           (format $pstream "} ")
  1430.           )
  1431.        (format $pstream "}}"))
  1432.      #+old
  1433.      (progn ; orig
  1434.        (print (list 'grid grid))
  1435.        (format $pstream "{plot3d {variable_grid ~%")
  1436.        (let* ((ar (polygon-pts pl))
  1437.           (x-coords
  1438.            (sloop for i to (nth 2 grid)
  1439.               collect (aref ar (* i 3))))
  1440.           (y-coords
  1441.            (sloop for i to (nth 3 grid)
  1442.               with m = (* 3 (+ 1 (nth 2 grid)))
  1443.               collect (aref ar (+ 1 (* i m)))))
  1444.           (z  (sloop for i to (nth 3 grid)
  1445.                  with k = 2
  1446.                  declare (fixnum k)
  1447.                  collect
  1448.                  (sloop for j to (nth 2 grid)
  1449.                     collect (aref ar k)
  1450.                     do(setq k (+ k 3))))))
  1451.          (tcl-output-list $pstream x-coords)
  1452.          (tcl-output-list $pstream y-coords)
  1453.          (format $pstream "~%{")
  1454.          (tcl-output-list $pstream z)
  1455.          (format $pstream "}}}"))
  1456.        )
  1457.      )
  1458.     ($geomview
  1459.      (format $pstream " MESH ~a ~a ~%"
  1460.          (+ 1 (nth 2 grid))
  1461.          (+ 1 (nth 3 grid))
  1462.          )
  1463.      (output-points pl nil))
  1464.     ($ps
  1465.      (let ((rot ($get_rotation ($rest ($get_plot_option '$view_direction)))))
  1466.        (cond (colour-z
  1467.           (setq *original-points* (copy-seq (polygon-pts pl)))
  1468.           (setq *z-range* ($get_range ar 2))))
  1469.        ($rotate_pts ar rot)
  1470.        (setup-for-ps-range ($get_range ar 0) ($get_range ar 1) t)
  1471.        (sort-ngons (polygon-pts pl) (polygon-edges pl) 4 )
  1472.            ;(print (polygon-edges pl))
  1473.        ($ps_axes rot)
  1474.        ($draw_ngons (polygon-pts pl) (polygon-edges pl) 4 )
  1475.        (p "[.1] 0 setdash")    ($ps_axes rot)
  1476.        (p "[] 0 setdash") (p " showpage "))))
  1477.       ;; close the stream and plot..
  1478.       (cond ($in_netmath (return-from $plot3d ""))
  1479.         (t (close $pstream)
  1480.            (setq $pstream nil)
  1481.            ))
  1482.       )
  1483.       (cond (($get_plot_option '$run_viewer 2)
  1484.          (case plot-format
  1485.            ($zic ($view_zic))
  1486.            ($ps ($viewps))
  1487.            ($openmath
  1488.         ($system (concatenate 'string *maxima-plotdir* "/" $openmath_plot_command) " maxout.openmath")
  1489.         )
  1490.            ($geomview ($system $geomview_command))
  1491.            ($gnuplot ($system (concatenate 'string *maxima-plotdir* "/" $gnuplot_command)
  1492.                   " -parametric3d maxout.gnuplot" ))
  1493.            )))
  1494.       )))
  1495.   
  1496.  
  1497.  
  1498. (setf (symbol-function '*$) (symbol-function '*))
  1499. (setf (symbol-function '+$) (symbol-function '+))
  1500. (setf (symbol-function '-$) (symbol-function '-))
  1501. (setf (symbol-function '/$) (symbol-function '/))
  1502.  
  1503. #|
  1504. Here is what the user inputs to draw the lattice picture.
  1505.  
  1506. /*Initially 1 unit = 1 pt = 1/72 inch
  1507. This makes 1 unit be 50/72 inch
  1508. */
  1509. ps_scale:[50,50];
  1510.  
  1511. /*This moves the origin to 400/72 inches up and over from bottom left corner
  1512. [ie roughly center of page]
  1513. */
  1514. ps_translate:[8,8];
  1515.  
  1516.  
  1517. f(x):=if (x = 0) then 100  else 1/x;
  1518.  
  1519. foo():=block([],
  1520. closeps(),
  1521. ps_translate:[6,6],
  1522. ps_scale:[50,50],
  1523. psdraw_curve(join(xcord,map(f,xcord))),
  1524. psdraw_curve(join(-xcord,map(f,xcord))),
  1525. psdraw_curve(join(xcord,-map(f,xcord))),
  1526. psdraw_curve(join(-xcord,-map(f,xcord))),
  1527. psdraw_points(lattice),
  1528. psaxes(8));
  1529.  
  1530.  
  1531. And here is the output .ps file which you should be
  1532. able to print on a laserwriter, or view on screen if you have
  1533. ghostscript (or another postscript screen previewer).
  1534.  
  1535. |#
  1536.  
  1537. #|
  1538. Examples
  1539.  
  1540. /* plot of z^(1/3)...*/
  1541. plot3d(r^.33*cos(th/3),[r,0,1],[th,0,6*%pi],['grid,12,80],['transform_xy,polar_to_xy],['view_direction,1,1,1.4],['colour_z,true],['plot_format,zic]);
  1542.  
  1543. /* plot of z^(1/2)...*/
  1544. plot3d(r^.5*cos(th/2),[r,0,1],[th,0,6*%pi],['grid,12,80],['transform_xy,polar_to_xy],['view_direction,1,1,1.4],['colour_z,true],['plot_format,zic]);
  1545.  
  1546. /* moebius */
  1547. plot3d([cos(x)*(3+y*cos(x/2)),sin(x)*(3+y*cos(x/2)),y*sin(x/2)],[x,-%pi,%pi],[y,-1,1],['grid,50,15]);
  1548.  
  1549. /* klein bottle */
  1550. plot3d([5*cos(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0) - 10.0,
  1551.           -5*sin(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0),
  1552.            5*(-sin(x/2)*cos(y)+cos(x/2)*sin(2*y))],[x,-%pi,%pi],[y,-%pi,%pi],
  1553.           ['grid,40,40]);
  1554. /* torus */
  1555. plot3d([cos(y)*(10.0+6*cos(x)),
  1556.            sin(y)*(10.0+6*cos(x)),
  1557.            -6*sin(x)], [x,0,2*%pi],[y,0,2*%pi],['grid,40,40]);
  1558. |#
  1559.  
  1560.