archive/src/plot2d_ps.lisp

   1 #| Existing documentation from the reference manual:
   2
   3  -- Function: plot2d_ps (<expr>, <range>)
   4      Writes to pstream a sequence of PostScript commands which plot
   5      <expr> over <range>.
   6
   7      <expr> is an expression.  <range> is a list of the form `[<x>,
   8      <min>, <max>]' in which <x> is a variable which appears in <expr>.
   9
  10      See also `closeps'.
  11
  12  -- Function: closeps ()
  13      This should usually becalled at the end of a sequence of plotting
  14      commands.   It closes the current output stream <pstream>, and sets
  15      it to nil.   It also may be called at the start of a plot, to
  16      ensure pstream is closed if it was open.    All commands which
  17      write to pstream, open it if necessary.   `closeps' is separate
  18      from the other plotting commands, since we may want to plot 2
  19      ranges or superimpose several plots, and so must keep the stream
  20      open.
  21
  22  -- Function: psdraw_curve (<ptlist>)
  23      Draws a curve connecting the points in <ptlist>.   The latter may
  24      be of the form `[x0, y0, x1, y1, ...]' or `[[x0, y0], [x1, y1],
  25      ...]'
  26
  27      The function `join' is handy for taking a list of x's and a list
  28      of y's and splicing them together.
  29
  30      <psdraw_curve> simply invokes the more primitive function
  31      <pscurve>.   Here is the definition:
  32
  33           (defun $psdraw_curve (lis)
  34             (p "newpath")
  35             ($pscurve lis)
  36             (p "stroke"))
  37
  38  -- Function: pscom (<cmd>)
  39      <cmd> is inserted in the PostScript file.  Example:
  40           pscom ("4.5 72 mul 5.5 72 mul translate 14 14 scale");
  41
  42    End of stuff from the reference manual.
  43  |#
  44
  45                                         #|
  46 Here is what the user inputs to draw the lattice picture.
  47
  48 /*Initially 1 unit = 1 pt = 1/72 inch
  49 This makes 1 unit be 50/72 inch
  50 */
  51 ps_scale:[50,50]                        ;
  52
  53 /*This moves the origin to 400/72 inches up and over from bottom left corner
  54 [ie roughly center of page]
  55 */
  56 ps_translate:[8,8]                      ;
  57
  58
  59 f(x):=if (x = 0) then 100  else 1/x     ;
  60
  61 foo():=block([],
  62 closeps(),
  63 ps_translate:[6,6],
  64 ps_scale:[50,50],
  65 psdraw_curve(join(xcord,map(f,xcord))),
  66 psdraw_curve(join(-xcord,map(f,xcord))),
  67 psdraw_curve(join(xcord,-map(f,xcord))),
  68 psdraw_curve(join(-xcord,-map(f,xcord))),
  69 psdraw_points(lattice),
  70 psaxes(8))                              ;
  71
  72
  73 And here is the output .ps file which you should be
  74 able to print on a laserwriter, or view on screen if you have
  75 ghostscript (or another postscript screen previewer).
  76
  77 |#
  78
  79 ;;(defvar $viewps_command  "(gs -I. -Q  ~a)")
  80
  81 ;;(defvar  $viewps_command "echo /def /show {pop} def |  cat - ~a | x11ps")
  82
  83 ;; If your gs has custom features for understanding mouse clicks
  84
  85
  86 ;;Your gs will loop for ever if you don't have showpage at the end of it!!
  87 ;;(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 &")
  88
  89 (defvar $viewps_command  "(ghostview \"~a\")")
  90
  91 (defvar $window_size '((mlist)
  92                        #.(* 8.5 72) #. (* 11 72) ))
  93
  94 (defun $getrange (x &optional xrange &aux yrange)
  95   (setq yrange  (cdr (get-range (cddr x))))
  96   (or xrange (setq xrange (cdr (get-range (cdr x)))))
  97   (setup-for-ps-range xrange yrange nil))
  98
  99 (defun $paramplot (f g range &optional (delta .1 supplied) &aux pts ($numer t))
 100   (setq f (coerce-float-fun f))
 101   (setq g (coerce-float-fun g))
 102   (setq range (meval* range))
 103   (or supplied (setq delta (/ (- (nth   2 range) (nth 1 range)) (nth 2 ($get_plot_option '$nticks)))))
 104   (setq pts(cons '(mlist)
 105                  (loop with tt = (coerce-float (nth 1 range))
 106                         with end = (coerce-float (nth 2 range))
 107                         while (float-< tt end)
 108                         collect (funcall f tt)
 109                         collect (funcall g tt)
 110                         do (setq tt ($+ tt delta)))))
 111   ($closeps)
 112   ($getrange pts)
 113   ($psdraw_curve pts)
 114   ($closeps)
 115   ($viewps))
 116
 117 (defun $plot2d_ps(fun range &rest options &aux ($numer t)  $display2d
 118                   ($plot_options $plot_options))
 119   (dolist (v options) ($set_plot_option v))
 120   (setq range (check-range range))
 121   (let ((tem (draw2d fun range )))
 122     ($closeps)
 123     ($getrange tem (cddr range))
 124     ($psdraw_curve tem)
 125     ($psaxes ($rest range))
 126     (p "showpage") ;; IS THIS NEEDED ??? ($viewps WILL APPEND A showpage TOO.)
 127     ($viewps)))
 128
 129 ;; do-ps-created-date was cribbed from the Common Lisp Cookbook -- Dates and Times.
 130 ;; See: http://cl-cookbook.sourceforge.net/dates_and_times.html
 131
 132 (defun do-ps-created-date (my-stream)
 133   (let ((day-names #("Mon" "Tue" "Wed" "Thu" "Fri" "Sat" "Sun")))
 134     (multiple-value-bind
 135           (second minute hour date month year day-of-week dst-p tz)
 136         (get-decoded-time)
 137       (declare (ignore dst-p))
 138       (format my-stream "%%CreatedDate: ~2,'0d:~2,'0d:~2,'0d ~a, ~d/~2,'0d/~d (GMT~@d)~%"
 139               hour minute second (aref day-names day-of-week) month date year (- tz)))))
 140
 141
 142 (defun do-ps-trailer ()
 143   (p "showpage")
 144   (p "%%Trailer")
 145   (p "%%EOF"))
 146
 147 ;;When we initialize we move the origin to the middle of $window_size
 148 ;;Then to offset from that use translate.
 149 (defvar $ps_translate '((mlist) 0 0))
 150
 151 ;; initially 1/72 of inch is the scale
 152 (defvar $ps_scale '((mlist) 72  72))
 153
 154
 155 (defun $pscom (&rest l)
 156   (apply 'p l))
 157
 158 ;;-10 to 10
 159 (defun psx (x) x)
 160 (defun psy (y) y)
 161
 162 (defun p (&rest l)
 163   (assureps)
 164   (loop for v in l do
 165          (if (symbolp v) (setq v (maxima-string v)))
 166                                         ; (if (numberp v) (setq v (* 70 v)))
 167          (cond ((consp v)
 168                 (loop for w in (cdr v) do (p w)))
 169                ((floatp v) (format $pstream "~,4g" v))
 170                (t(princ v $pstream)))
 171          (princ " " $pstream))
 172   (terpri $pstream))
 173
 174 (defun psapply (f lis)
 175   (if ($listp lis) (setq lis (cdr lis)))
 176   (apply 'p lis)
 177   (p f))
 178
 179 (defun $moveto (x y)
 180   (p (psx x) (psy y) "moveto "))
 181
 182 (defun $psline (a b c d)
 183   (p (psx a) (psy b) "moveto ")
 184   (p  (psx c) (psy d) "lineto"))
 185
 186 (defun setup-for-ps-range (xrange yrange do-prolog)
 187   (let* ((cy (/ (+ (nth 1 yrange) (nth 0 yrange)) 2.0))
 188          (cx (/ (+ (nth 1 xrange) (nth 0 xrange)) 2.0))
 189          (scaley (/ (nth 2 $window_size)
 190                     (* 1.2  (- (nth 1 yrange) (nth 0 yrange)))))
 191          (scalex (/ (nth 1 $window_size)
 192                     (* 1.2 (- (nth 1 xrange) (nth 0 xrange)))))
 193          ($ps_scale
 194           (progn
 195             (cond ((< scalex scaley)
 196                    (setq scaley scalex)))
 197             `((mlist) , scaley ,scaley)))
 198          ($ps_translate `((mlist) , cx ,cy)))
 199     (assureps do-prolog)))
 200
 201 (defun assureps (&optional do-prolog)
 202   (cond ((streamp $pstream))
 203         (t (setq do-prolog t)))
 204   (or $pstream (setq $pstream (open (plot-temp-file "maxout.ps") :direction :output :if-exists :supersede)))
 205   (cond (do-prolog
 206             (p "%!PS-Adobe-2.0")
 207           (p "%%Title: Maxima 2d plot") ;; title could be filename and/or plot equation
 208           (p "%%Creator: Maxima version:" *autoconf-version* "on" (lisp-implementation-type))
 209           (do-ps-created-date $pstream)
 210           (p "%%DocumentFonts: Helvetica")
 211
 212           ;; Put the lower left corner of the bounding box at $ps_translate,
 213           ;; and put the upper right corner at ($ps_translate + $window_size).
 214           (p "%%BoundingBox:"
 215              (round (nth 1 $ps_translate))
 216              (round (nth 2 $ps_translate))
 217              (round (+ (nth 1 $ps_translate) (nth 1 $window_size)))
 218              (round (+ (nth 2 $ps_translate) (nth 2 $window_size))))
 219
 220           (p "%%Orientation: Portrait")
 221           (p "%%Pages: 1")
 222           (p "%%EndComments")
 223           (p "%%BeginPrologue")
 224           (p "%%EndPrologue")
 225           (p "%%BeginSetup")
 226           (p "%%PaperSize: Letter")
 227           (p "%%EndSetup")
 228           (p "%%Page: 1 1")
 229           (p  (* .5 (nth 1 $window_size))
 230               (* .5 (nth 2 $window_size))
 231               "translate"
 232               )
 233           (psapply "scale" $ps_scale)
 234           (p  (- (nth 1 $ps_translate))
 235               (- (nth 2 $ps_translate))
 236               "translate"
 237               )
 238           (p " 1.5 " (second $ps_scale) "div setlinewidth
 239 /Helvetica findfont 14 " (second $ps_scale) " div scalefont setfont
 240 /dotradius .05 def
 241 /drawdot {
 242  /yy exch def
 243  /xx exch def
 244   gsave
 245   xx yy dotradius 0 360 arc
 246   fill
 247   grestore
 248 }def
 249
 250 /ticklength  .03 def
 251 /axiswidth  .01 def
 252 /drawtick {
 253  /yy exch def
 254  /xx exch def
 255   gsave
 256   xx ticklength sub yy moveto
 257   xx ticklength add yy lineto
 258   stroke
 259   xx yy  ticklength sub  moveto
 260   xx yy  ticklength add  lineto
 261   stroke
 262   grestore
 263 } def
 264 "))))
 265
 266
 267 (defun $closeps ()
 268   (prog1
 269       (when (and (streamp $pstream)
 270                  )
 271         (p "showpage")
 272         (close  $pstream))
 273     (setq $pstream nil)))
 274
 275 (defun ps-fixup-points(lis)
 276   (assert ($listp lis))
 277   (setq lis (cdr lis))
 278   (cond ((numberp (car lis)))
 279         ((and ($listp (car lis)) (numberp (nth 1 (car lis))))
 280          (setq lis (loop for w in lis collect
 281                           (nth 1 w)
 282                           collect (nth 2 w))))
 283         (t (error
 284             "pscurve:Neither [x0,y0,x1,y1,...] nor [[x0,y0],[x1,y1],...]")))
 285   lis)
 286
 287 (defun $psdraw_curve (lis &aux (n 0))
 288   (declare (fixnum n))
 289   (setq lis (ps-fixup-points lis))
 290   (p "newpath" (nth 0 lis) (nth 1 lis) "moveto")
 291   (loop while lis with second
 292          do
 293          (cond ((eq (car lis) 'moveto)
 294                 (or (eql n 0) (p "stroke"))
 295                 (setq n 0) (setq lis (cddr lis))))
 296          (or (setq second (cadr lis)) (error "odd length list of points"))
 297          (cond ((eql n 0)
 298                 (p (car lis) second "moveto"))
 299                (t  (p (car lis) second "lineto")
 300                    ))
 301          (setq n (+ n 1))
 302          (cond ((eql 0 (mod n 20))
 303                 (p "stroke")
 304                 (setq n 0))
 305                (t (setq lis (cddr lis)))))
 306   (or (eql n 0) (p "stroke")))
 307
 308
 309
 310 (defun $viewps ( &optional file)
 311   (cond  ((and (streamp $pstream))
 312           (do-ps-trailer)
 313           (force-output $pstream)))
 314   (cond (file (setq file (maxima-string file)))
 315         (t(setq file (plot-temp-file "maxout.ps"))
 316
 317           (if (and (streamp $pstream))
 318               (force-output $pstream))))
 319   (if (equal $viewps_command "(gs -I. -Q  ~a)")
 320       (format t "~%type `quit' to exit back to affine or maxima
 321   To reprint a page do
 322   GS>showpage
 323   GS>(maxout.ps)run
 324   GS> -150 -150 translate 1.2 1.2 scale (maxout.ps)run
 325   would print it moved 150/72 inches to left, and down, and scaled by 1.2 times
 326   showpage clears scaling."))
 327
 328   ($system    (format nil $viewps_command file)))
 329
 330 (defun $chkpt (a)
 331   (or (and ($listp a)
 332            (numberp (nth 1 a))
 333            (numberp (nth 2 a)))
 334       (error "illegal pt ~a" a))
 335   t)
 336
 337 (defvar $pslineno nil)
 338 (defun $psdrawline (a &optional b c d)
 339   (cond ((null b)
 340          (assert (and ($listp a)
 341                       ($chkpt (nth 1 a))))
 342          (setq b (nth 2 a))
 343          (setq a (nth 1 a))))
 344   (cond ((null c)
 345          ($chkpt b)
 346          (setq c (nth 1 b))
 347          (setq d (nth 2 b))))
 348   (cond (($listp a)
 349          (setq b (nth 2 a) a (nth 1 a))))
 350   (cond ((null c)
 351          (setq c (nth 1 b))
 352          (setq c (nth 2 b))))
 353
 354   (when $pslineno
 355     (incf $pslineno)
 356     ($pslabelline a b c d $pslineno))
 357
 358   (p a b "moveto")
 359   (p c d "lineto")
 360   (p "stroke"))
 361
 362 (defun $pslabelline (a b c d $pslineno)
 363   (p (/ (+ a c) 2)
 364      (/ (+ b d) 2)
 365      "moveto"
 366      (format nil "(<--L~a)show" $pslineno)))
 367
 368 (defun $sort_polys (lis)
 369   (let ((tem
 370          (loop for v in (cdr lis)
 371                 collect (cons
 372                          (loop for w in (cdr v)
 373                                 maximize (nth 3 w))
 374                          v))))
 375     (print 'next)
 376     (cons '(mlist) (mapcar 'cdr (sortcar tem '<)))))
 377
 378 (defun $drawpoly (x)
 379   (p "gsave")
 380   (p (cdadr x) "moveto")
 381   (setq x (cddr x))
 382   (loop for edge in x
 383          do (p (cdr edge) "lineto")
 384          finally (p "1 setgray fill"))
 385   ($psdrawline x)
 386   (p "grestore"))
 387
 388 (defun $psdrawpolys (polys)
 389   (dolist (v (cdr polys))
 390     ($drawpoly v)))
 391
 392 ;; draw the axes through $psdef
 393 (defun $psaxes (leng &optional (lengy leng))
 394   (p "gsave axiswidth setlinewidth")
 395   (let (begx begy endx endy)
 396     (cond ((numberp leng)
 397            (setq begx (- leng))
 398            (setq endx leng))
 399           (t (setq begx (nth 1 leng))
 400              (setq endx (nth 2 leng))))
 401     (cond ((numberp lengy)
 402            (setq begy (- lengy))
 403            (setq endy lengy))
 404           (t (setq begy (nth 1 lengy))
 405              (setq endy (nth 2 lengy))))
 406
 407     (loop for i from (floor begx) below (ceiling endx)
 408            do
 409            ($psdrawline i 0 (+ i 1) 0)
 410            (p i 0 "drawtick")
 411            (p (+ i 1) 0 "drawtick"))
 412
 413     (loop for i from (floor begy) below (ceiling endy)
 414            do
 415            ($psdrawline 0 i 0 (+ i 1) )
 416            (p 0 i "drawtick")
 417            (p 0 (+ i 1)  "drawtick"))
 418
 419
 420     (p "grestore")))
 421
 422 (defun $psdraw_points(lis)
 423   (assert (and ($listp lis)
 424                ($listp (cadr lis))))
 425   (loop for w in (cdr lis)
 426          do (p (nth 1 w)
 427                (nth 2 w)
 428                "drawdot")))
 429
 430 (defun $ps_axes ( rot )
 431   (let ((tem (make-array 9 :element-type 'double-float :initial-element 0d0)))
 432     (setf (aref tem 0) 4.0)
 433     (setf (aref tem 4) 4.0)
 434     (setf (aref tem 8) 4.0)
 435     ($rotate_pts tem rot)
 436     (p 0 0 "moveto")
 437     (p (aref tem 0) (aref tem 1) "lineto stroke")
 438     (p  (aref tem 0) (aref tem 1) "moveto (x) show")
 439     (p 0 0 "moveto")
 440     (p (aref tem 3) (aref tem 4) "lineto stroke")
 441     (p  (aref tem 3) (aref tem 4) "moveto (y) show")
 442     (p 0 0 "moveto")
 443     (p (aref tem 6) (aref tem 7) "lineto stroke")
 444     (p  (aref tem 6) (aref tem 7) "moveto (z) show")
 445     ))
 446
 447 ;; if this is '$polar_to_xy then the x y coordinates are interpreted as r theta
 448
 449 (defun add-ps-finish (opts)
 450   (p (if opts
 451          "/xr .30 def
 452 /xb 0.60 def
 453 /xg .60  def
 454 /myset { .005 mul dup xr add exch
 455  dup xb add exch
 456  xg add
 457 setrgbcolor} def
 458
 459 /myfinish { myset  gsave fill grestore 0 setgray stroke  } def"
 460
 461          "/myfinish {.9 setgray gsave fill grestore .1 setgray stroke  } def")))
 462
 463
 464 (defun $draw_ngons (pts ngons number_edges &aux (i 0) (j 0) (s 0)
 465                     (opts *original-points*)
 466                     (maxz  most-negative-double-float))
 467   (declare (type (cl:array double-float) pts)
 468            #-(or cmu scl sbcl) (type (cl:array double-float) opts)
 469            (type (cl:array (mod 64000)) ngons)
 470            (fixnum number_edges i s j number_edges)
 471            (double-float maxz))
 472   (setq j (length ngons))
 473   (add-ps-finish opts)
 474   (loop while (< i j)
 475          do
 476          (loop initially (setq s number_edges)
 477                 do
 478                                         ;(print-pt (aref pts (f* 3 (aref ngons i))))
 479                 (print-pt (x-pt pts  (aref ngons i)))
 480                                         ;(print-pt (aref pts (f+ 1 (f* 3 (aref ngons i)))))
 481                 (print-pt (y-pt pts  (aref ngons i)))
 482
 483                 (cond (opts (if (> (z-pt opts (aref ngons i)) maxz)
 484                                 (setq maxz (z-pt opts (aref ngons i))))))
 485                 (cond ((eql number_edges s) (p " moveto %"
 486                                         ;(aref pts (f+ 2 (f* 3 (aref ngons i))))
 487
 488                                                ))
 489                       (t (p "lineto %" ;(aref pts (f+ 2 (f* 3 (aref ngons i))))
 490                             )))
 491                 (setq i (f+ i 1))
 492                 while (> (setq s (f- s 1)) 0))
 493          (setq i (f+ i 1))
 494          (cond (opts
 495                 (p (f+ 1 (round ($* 100.0 ($/ ($- maxz (car *z-range*))
 496                                               (or (third *z-range*)
 497                                                   ($- (second *z-range*) (car *z-range*))))))))
 498                 (setq maxz most-negative-double-float)
 499                 ))
 500          (p " myfinish")
 501          ))
 502