I find an obscure error when calling rk() inside a
homemade function to integrate dr/dt > 0 part of
orbit for rutherford scattering (natural units).
The homemade function is points_plus().
using windows xp binary, xmaxima,
maxima version 5.28.0
--------------------------------------------
Maxima 5.28.0-2 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
load k1util.mac
fll(mylist), head(mylist), tail(mylist)
(%i1) points_plus(e,b,tmax,dt) :=
block([w,phi_inf,theta0,rmin,
th,r,pts,pair,xp,yp,xyL:[],j,numer],numer:true,
w : 1/(e*b),
phi_inf : acos(w/sqrt(4+w^2)),
theta0 : %pi - phi_inf,
rmin : b*(w + sqrt(4+w^2))/2,
print(" theta0 = ",theta0," rmin = ",rmin),
pts : rk([-b/r^2, sqrt(1 - 1/(e*r) -b^2/r^2)],[th,r],
[theta0,rmin],[t,0,tmax,dt]),
print(" pts = ",pts),
for j thru length(pts) do (
pair : rest(pts[j]),
xp : pair[2]*cos(pair[1]),
yp : pair[2]*sin(pair[1]),
xyL : cons([xp,yp],xyL)),
reverse(xyL))$
(%i2) pts_plus : points_plus(1,0.5,1,0.05)$
theta0 = 2.356194490192345 rmin = 1.207106781186548
Expecting a number when the initial state is replaced
in the equations, but instead found:
1.0536712127723509E-8*(-1)^0.5
#0: rk(odes=[-0.5/r^2,(-1/r-0.25/r^2+1)^0.5],state=[th,r],
initial=[2.356194490192345,1.207106781186548],
domain=[t,0,1,0.05])(dynamics.mac line 392)
#1: rk(?_l=[[-0.5/r^2,(-1/r-0.25/r^2+1)^0.5],[th,r],
[2.356194490192345,1.207106781186548],[t,0,1,0.05]])
#2: points_plus(e=1,b=0.5,tmax=1,dt=0.05)
-- an error. To debug this try: debugmode(true);
(%i3) pts: rk([-0.5/r^2,(-1/r-0.25/r^2+1)^0.5],[th,r],
[2.356194490192345,1.207106781186548],[t,0,1,0.05])$
(%i4) head(pts);
(%o4) [[0,2.356194490192345,1.207106781186548],
[0.05,2.339037376923019,1.207168112731729],
[0.1,2.321892342072049,1.208117390658052],
[0.15,2.304791481677257,1.210282757528192],
[0.2,2.287768700622341,1.213655110645663],
[0.25,2.27085695447911,1.218225863979847]]
------------------------------------------------------------------
Ted Woollett