Re: How do I make a log-log plot?



--=-=-=
Content-Disposition: inline

>>>>> "Raymond" == Raymond Toy  writes:

    Raymond> I have a little hack now that seems to work reasonably well for
    Raymond> adaptive 2D plotting.

Here is the little hack in case anyone wants to play with it.  If you
specify the plot-option logx, it will assume you want log sampling of
the arg, and do the appropriate stuff.  It also assumes that you'll do
a log plot when plotting the data, but that's not necessary.  It seems
to be ok even if you don't.

Ray



Index: src/plot.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/plot.lisp,v
retrieving revision 1.49
diff -u -b -r1.49 plot.lisp
--- src/plot.lisp	12 Apr 2005 04:04:43 -0000	1.49
+++ src/plot.lisp	7 Jun 2005 18:33:08 -0000
@@ -68,6 +68,7 @@
 			 "set term dumb 79 22")
 			((mlist) $gnuplot_ps_term_command
 			 "set size 1.5, 1.5;set term postscript eps enhanced color solid 24")
+			((mlist) $logx nil)
 			))
 
 (defvar $viewps_command  "(ghostview  ~a)")
@@ -148,6 +149,7 @@
 	  ($gnuplot_dumb_term_command value)
 	  ($gnuplot_ps_term_command value)
 	  ($adapt_depth (check-list-items name (cddr value) 'fixnum 1))
+	  ($logx value)
 	  (t
 	   (merror "Unknown plot option specified:  ~M" name))))
   (loop for v on (cdr $plot_options)
@@ -868,7 +870,7 @@
 		 (right (adaptive-plot f b b1 c f-b f-b1 f-c (1- depth) (* 2 eps))))
 	     (append left (cddr right)))))))
 
-(defun draw2d (f range)
+(defun draw2d (f range &optional (log-x-p t))
   (if (and ($listp f) (equal '$parametric (cadr f)))
       (return-from draw2d (draw2d-parametric f range)))
   (if (and ($listp f) (equal '$discrete (cadr f)))
@@ -896,10 +898,20 @@
       ;; NTICKS too big.  Since adaptive plotting splits the sections
       ;; in half, it's also probably not a good idea to have NTICKS be
       ;; a power of two.
+      (when log-x-p
+	(setf x-start (log x-start))
+	(setf xend (log xend))
+	(setf x-step (/ (- xend x-start) (coerce-float nticks) 2)))
+
+      (flet ((fun (x)
+	       (if log-x-p
+		   (funcall f (exp x))
+		   (funcall f x))))
+	
       (dotimes (k (1+ (* 2 nticks)))
 	(let ((x (+ x-start (* k x-step))))
 	  (push x x-samples)
-	  (push (funcall f x) y-samples)))
+	    (push (fun x) y-samples)))
       (setf x-samples (nreverse x-samples))
       (setf y-samples (nreverse y-samples))
 
@@ -923,10 +935,10 @@
 	      (if result
 		  (append result
 			  (cddr
-			   (adaptive-plot f (car x-start) (car x-mid) (car x-end)
+			     (adaptive-plot #'fun (car x-start) (car x-mid) (car x-end)
 					  (car y-start) (car y-mid) (car y-end)
 					  depth 1d-5)))
-		  (adaptive-plot f (car x-start) (car x-mid) (car x-end)
+		    (adaptive-plot #'fun (car x-start) (car x-mid) (car x-end)
 				 (car y-start) (car y-mid) (car y-end)
 				 depth 1d-5))))
 	  
@@ -938,11 +950,13 @@
       (do ((x result (cddr x))
 	   (y (cdr result) (cddr y)))
 	  ((null y))
+	  (when log-x-p
+	    (setf (car x) (exp (car x))))
 	(unless (and (numberp (car y))
 		     (<= ymin (car y) ymax))
 	  (setf (car x) 'moveto)
 	  (setf (car y) 'moveto)))
-      (cons '(mlist) result))))
+	(cons '(mlist) result)))))
 
 (defun get-range (lis)
   (let ((ymin most-positive-double-float)
@@ -1085,7 +1099,8 @@
 (defun $plot2d (fun &optional range &rest options &aux ($numer t) $display2d
 		(i 0) plot-format gnuplot-term gnuplot-out-file
 		file plot-name
-		($plot_options $plot_options))
+		($plot_options $plot_options)
+		log-x)
   (dolist (v options) ($set_plot_option v))
   
   (cond ((and (consp fun) (eq (cadr fun) '$parametric))
@@ -1098,6 +1113,9 @@
   (cond ((eq ($get_plot_option '$plot_format 2) '$openmath)
          (return-from $plot2d (apply '$plot2dopen fun range options))))
 
+  (when ($get_plot_option '$logx 2)
+    (setf log-x t))
+  
   ;; this has to come after the checks for ps and openmath 
   ;; (see bug report #834729)
   (or ($listp fun ) (setf fun `((mlist) ,fun)))	
@@ -1171,7 +1189,7 @@
 	($mgnuplot
 	 (format st "~%~%# \"~a\"~%" plot-name))
 	)
-      (loop for (v w) on (cdr (draw2d v range )) by #'cddr
+      (loop for (v w) on (cdr (draw2d v range log-x)) by #'cddr
 	    do
 	    (cond ((eq v 'moveto)
 		   (cond 

--=-=-=--