/* COPYRIGHT NOTICE Copyright (C) 2014, 2015 Edwin L. (Ted) Woollett http://www.csulb.edu/~woollett rutherford_attract.mac is a utility file associated with Project 1 (Classical Scattering in a Central Potential) associated with Ch. 1 (Numerical Differentiation, Quadrature, and Roots) in the series Computational Physics with Maxima or R, See project1.pdf for more info. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details at http://www.gnu.org/copyleft/gpl.html */ /* rutherford_attract.mac V(r) = - alpha/r with alpha>0, attractive potential */ deg(z) := block([numer],numer:true, z*180/%pi)$ /* make_pts, list arithmetic version, fill in intermediate points along a straight line N = number of intervals, N+1 = number of points returns list [xL,yL] */ make_pts(x1,y1,xf,yf,N):= block([h, slope, xL,yL, numer], numer:true, h : (xf - x1)/N, slope : (yf-y1)/(xf-x1), xL : makelist(x,x,x1,xf,h), yL : y1 + slope*(xL - x1), [ xL, yL ])$ /* (%i14) make_pts(0,0,4,4,8); (%o14) [[0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0],[0,0.5,1.0,1.5,2.0,2.5,3.0,3.5, 4.0]] */ /* angles(e,b): attractive rutherford case some (e,b) dependent angles and rmin. theta0 defines the angle of closest approach (counter-clockwise from positive x axis) phi_inf defines the rotation angle of r-vec as the incident particle comes from theta=%pi, r = inf to the point of closest approach, so theta0 = %pi - phi_inf. chi is the scattering angle (attractive case) chi = %pi - 2*phi_inf = -2*asin(w/sqrt(w^2+4)) < 0. rmin is the distance of closest approach to the scattering center = ( x = 0, y = 0 ) */ angles(e,b) := block([w,phi_inf,theta0,chi,oldfp,rmin,numer],numer:true, oldfp : fpprintprec, fpprintprec : 8, w : 1/(e*b), phi_inf : asin(w/sqrt(w^2+4)) + %pi/2, /* attractive case */ print(" phi_inf = ",phi_inf," rad, or ",phi_inf*180/%pi," deg"), theta0 : %pi - phi_inf, print(" theta0 = ",theta0," rad, or ",theta0*180/%pi," deg"), chi : -2*asin(w/sqrt(w^2+4)), /* attractive case */ print(" chi = ",chi," rad, or ",chi*180/%pi," deg"), rmin : b*(-w + sqrt(4+w^2))/2, /* attractive case */ print(" rmin = ",rmin), fpprintprec : oldfp, done)$ /* (%i18) angles(1,1)$ phi_inf = 2.0344439 rad, or 116.56505 deg theta0 = 1.1071487 rad, or 63.434949 deg chi = -0.927295 rad, or -53.130102 deg rmin = 0.618034 */ /* analytic scattering angle in degrees for given values of Ebar and bbar */ angle_a(e,b) := block([numer],numer:true, -2*asin(1/sqrt(1 + 4*e^2*b^2))*180/%pi)$ /* numerical scattering angle in degrees for given values of Ebar and bbar uses achi() and deg(). */ angle_n(e,b) := (deg(achi(e,b)))$ /* (%i25) angle_a(1,1); (%o25) -53.130102 (%i26) angle_n(1,1); (%o26) -53.130102 */ /*********** plot of Veff **********************/ /* Veff_plot for attractive Rutherford case */ Veff_plot(e,b,r0,r1,xmin,xmax,ymin,ymax) := block([w,rmin,energy_line,numer],numer:true, w : 1/e/b, rmin : b*(-w + sqrt(4+w^2))/2, /* attractive case */ print(" rmin = ",rmin), energy_line : [discrete,[[rmin,e],[r1,e]]], plot2d([-1/r + e*b^2/r^2, energy_line],[r,r0,r1], [x,xmin,xmax], [y,ymin,ymax],[style,[lines,3,1],[lines,3,2]], [legend,false],[xlabel,"r"],[ylabel,"Veff"]))$ /* (%i51) Veff_plot(1,1,0.4,8,0.4,8,-2,2)$ rmin = 0.618034 plot2d: some values were clipped. */ Veff_plot_eps(e,b,r0,r1,xmin,xmax,ymin,ymax,fname) := block([w,rmin,energy_line,numer],numer:true, w : 1/e/b, rmin : b*(-w + sqrt(4+w^2))/2, /* attractive case */ energy_line : [discrete,[[rmin,e],[r1,e]]], plot2d([-1/r + e*b^2/r^2, energy_line],[r,r0,r1], [x,xmin,xmax], [y,ymin,ymax],[style,[lines,7,1],[lines,7,2]], [legend,false],[xlabel,"r"],[ylabel,"Veff"], [psfile, fname]))$ /* Veff_plot_eps(1,1,0.9,3,0.8,3,0,2,"c:/k1/rr9.eps"); */ /********* TRAJECTORY PLOT FUNCTIONS *********/ /* init(e,b) specific to attractive rutherford case, produces global definitions of chi,xc,yc, vcx,vcy,and xa. Besides printing out these values, init(e,b) also prints out the values of rmin, phi_inf, and theta0. */ init(e,b) := block([w,rmin,phi_inf,theta0,numer],numer:true, w : 1/e/b, rmin : b*(-w + sqrt(4+w^2))/2, /* attractive case */ print(" rmin = ",rmin), phi_inf : asin(w/sqrt(w^2+4)) + %pi/2, /* attractive case */ print(" phi_inf = ", phi_inf," rad or ", phi_inf*180/%pi," deg"), theta0 : %pi - phi_inf, print(" theta0 = ",theta0," rad or ", theta0*180/%pi," deg"), chi : (%pi - 2*phi_inf), /* global parameter chi in radians */ print(" chi = ", chi, " rad, or ", chi*180/%pi," deg"), /* point of closest approach */ xc : rmin*cos(theta0), /* global */ yc : rmin*sin(theta0), /* global */ vcx : b*sin(theta0)/rmin, /* global */ vcy : -b*cos(theta0)/rmin, /* global */ print(" xc = ", xc," yc = ", yc), print(" vcx = ", vcx," vcy = ", vcy), /* x-intersection of rmin line with y=b line */ xa : xc*b/yc, /* global */ print(" xa = ",xa))$ /* (%i53) init(1,1); rmin = 0.618034 phi_inf = 2.0344439 rad or 116.56505 deg theta0 = 1.1071487 rad or 63.434949 deg chi = -0.927295 rad, or -53.130102 deg xc = 0.276393 yc = 0.552786 vcx = 1.4472136 vcy = -0.723607 xa = 0.5 (%o53) 0.5 */ /* orbit_plot1 for attractive rutherford potential. calls init(e,b) to define local values of xc,yc,xa,vcx,vcy,chi. Calls make_pts to produce rmin_line extension. The final two args are used to control the number of points and their size for the rmin_line extension. Calls utility function take() to extract xL and yL from rkpts. */ orbit_plot1(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax,nextend,psize) := block([xc,yc,vcx,vcy,xa,chi,rkpts,rmag,dvxdt,dvydt, xL, yL,pm,pp,vxf,vyf,bline,rmin_line, xf,yf,chi_line,rmin_extend,numer],numer:true, /* define local xc, yc, vcx, vcy, xa, chi */ init(e,b), /* symbolic expressions for acceleration components */ rmag : sqrt(x^2 + y^2), /* rmag in terms of x and y */ dvxdt : -x/2/e/rmag^3, /* dvx/dt attractive rutherford case */ dvydt : -y/2/e/rmag^3, /* dvy/dt */ /* integrate backwards from xc, yc */ rkpts : rk([vx, dvxdt, vy, dvydt ], [x,vx,y,vy],[xc,vcx,yc,vcy],[t,0,-tm,-dt]), xL : take(rkpts,2), yL : take(rkpts,4), pm : [discrete, xL, yL], print(" backwards from xc, yc "), print(" xfirst = ",last(xL)), /* integrate forwards from xc,yc */ rkpts : rk([vx, dvxdt, vy, dvydt ], [x,vx,y,vy],[xc,vcx,yc,vcy],[t,0,tp,dt]), xL : take(rkpts,2), yL : take(rkpts,4), pp : [discrete, xL, yL], print(" forwards from xc, yc "), print(" xlast = ",last(xL)," ylast = ",last(yL)), vxf : last(rkpts)[3], vyf : last(rkpts)[5], print(" vx_last = ", vxf," vy_last = ", vyf), /* make plot of orbit */ bline : [discrete,[ [xmin,b],[xmax,b] ]], rmin_line : [discrete,[ [0,0],[xc,yc] ]], yf : b + rchi*sin(chi), xf : xa + rchi*cos(chi), print(" xf = ", xf," yf = ", yf), chi_line : [discrete, [ [xa,b],[xf,yf]]], rmin_extend : make_pts(xc,yc,xa,b,nextend), rmin_extend : [discrete, rmin_extend[1],rmin_extend[2]], plot2d( [ bline, rmin_line, chi_line, pm, pp, rmin_extend ], [x, xmin, xmax],[y, ymin, ymax], [style,[lines,2,5],[lines,2,2],[lines,2,5], [lines,3,1],[lines,3,1], [points, psize, 2, 1]],[legend,false]))$ /* new version orbit_plot1, feb 5: (%i6) orbit_plot1(1,1,5,5,0.01,5,-3,3.4,-2,2,5,1); rmin = 0.618034 phi_inf = 2.0344439 rad or 116.56505 deg theta0 = 1.1071487 rad or 63.434949 deg chi = -0.927295 rad, or -53.130102 deg xc = 0.276393 yc = 0.552786 vcx = 1.4472136 vcy = -0.723607 xa = 0.5 backwards from xc, yc xfirst = -5.7276639 forwards from xc, yc xlast = 4.2046373 ylast = -4.0061018 vx_last = 0.655096 vy_last = -0.861996 xf = 3.5 yf = -3.0 (%i7) orbit_plot1(1,0.6,5,5,0.01,5,-3,3.4,-2,2,5,1); rmin = 0.281025 phi_inf = 2.2655346 rad or 129.80557 deg theta0 = 0.876058 rad or 50.194429 deg chi = -1.3894766 rad, or -79.611142 deg xc = 0.179908 yc = 0.215889 vcx = 1.6401844 vcy = -1.3668203 xa = 0.5 backwards from xc, yc xfirst = -5.9196184 forwards from xc, yc xlast = 1.6346676 ylast = -5.7185896 vx_last = 0.198759 vy_last = -1.0623693 xf = 1.4016393 yf = -4.3180328 (%i8) orbit_plot1(1,0.3,5,5,0.01,5,-3,3.4,-2,2,5,1); rmin = 0.0830952 phi_inf = 2.6011732 rad or 149.03624 deg theta0 = 0.54042 rad or 30.963757 deg chi = -2.0607537 rad, or -118.07249 deg xc = 0.0712535 yc = 0.0427521 vcx = 1.8574929 vcy = -3.0958215 xa = 0.5 backwards from xc, yc xfirst = -6.0729577 forwards from xc, yc xlast = -2.6061534 ylast = -5.4927368 vx_last = -0.506288 vy_last = -0.951945 xf = -1.8529412 yf = -4.1117647 */ /* new orbit_plot1_eps feb.5 for attractive rutherford */ orbit_plot1_eps(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax,nextend,psize,fname) := block([xc,yc,vcx,vcy,xa,chi,rkpts,rmag,dvxdt,dvydt, xL, yL,pm,pp,vxf,vyf,bline,rmin_line, xf,yf,chi_line,rmin_extend,numer],numer:true, /* define local xc, yc, vcx, vcy, xa, chi */ init(e,b), /* symbolic expressions for acceleration components */ rmag : sqrt(x^2 + y^2), /* rmag in terms of x and y */ dvxdt : -x/2/e/rmag^3, /* dvx/dt attractive rutherford case */ dvydt : -y/2/e/rmag^3, /* dvy/dt */ /* integrate backwards from xc, yc */ rkpts : rk([vx, dvxdt, vy, dvydt ], [x,vx,y,vy],[xc,vcx,yc,vcy],[t,0,-tm,-dt]), xL : take(rkpts,2), yL : take(rkpts,4), pm : [discrete, xL, yL], print(" backwards from xc, yc "), print(" xfirst = ",last(xL)), /* integrate forwards from xc,yc */ rkpts : rk([vx, dvxdt, vy, dvydt ], [x,vx,y,vy],[xc,vcx,yc,vcy],[t,0,tp,dt]), xL : take(rkpts,2), yL : take(rkpts,4), pp : [discrete, xL, yL], print(" forwards from xc, yc "), print(" xlast = ",last(xL)," ylast = ",last(yL)), vxf : last(rkpts)[3], vyf : last(rkpts)[5], print(" vx_last = ", vxf," vy_last = ", vyf), /* make plot of orbit */ bline : [discrete,[ [xmin,b],[xmax,b] ]], rmin_line : [discrete,[ [0,0],[xc,yc] ]], yf : b + rchi*sin(chi), xf : xa + rchi*cos(chi), print(" xf = ", xf," yf = ", yf), chi_line : [discrete, [ [xa,b],[xf,yf]]], rmin_extend : make_pts(xc,yc,xa,b,nextend), rmin_extend : [discrete, rmin_extend[1],rmin_extend[2]], plot2d( [ bline, rmin_line, chi_line, pm, pp, rmin_extend ], [x, xmin, xmax],[y, ymin, ymax], [style,[lines,3,5],[lines,5,2],[lines,3,5], [lines,5,1],[lines,5,1], [points, psize, 2, 1]],[legend,false], [psfile, fname]))$ /* end new orbit_plot1_eps feb.5 for attractive rutherford */ /****** DIFFERENTIAL CROSS SECTION PLOT FUNCTIONS ********/ /* f1d(nv,hh,gL) returns list of first derivatives at the nv grid points for a function whose nv values at the grid points separated by hh are in the list gL. Would be more accurate if we used quadratic interpolation to define the end of grid derivatives. */ f1d(nv,hh,gL) := block([j,fpL:[],fp0,fpl,numer],numer:true, for j:2 thru nv-1 do fpL : cons( (gL[j+1] - gL[j-1])/2/hh,fpL), fpL : reverse(fpL), /* use linear interpolation to define first and last elements of fpL */ fp0 : 2*fpL[1] - fpL[2], fpl : 2*fpL[nv-2] - fpL[nv-3], fpL : cons(fp0,fpL), fpL : append(fpL,[fpl]), fpL)$ /* achi(e,b), specific to attractive rutherford case, returns the signed scattering angle in radians using numerical methods. */ achi(e,b) := block([z,w,zmin,phi_inf,rexpr,iexpr,numer],numer:true, w : 1/(b*e), rexpr : z^2 + w*z -1, /* root -> zmin for attractive case */ iexpr : 1/(z*sqrt(rexpr)), /* integrate to get phi_inf */ zmin : find_root(rexpr,z,1e-4,1e6), phi_inf : quad_qagi(iexpr,z,zmin,inf)[1], (%pi - 2*phi_inf))$ /* (%i15) achi(1,1); (%o15) -0.927295 */ /* sigma_points(e,b0,bmax,db) produces a list of two lists: [chi-list, log(d_sig/d_omega)-list] using numerical methods. Typical list arithmetic Maxima methods replace conventional loop methods here. calls achi() and f1d() */ sigma_points(ee,b0,bmax,db):= block([nb,bL,chiL,dchi_dbL,sigL,numer],numer:true, bL : makelist(b,b,b0,bmax,db), nb : length(bL), chiL : map(lambda([x], achi(ee,x)), bL), dchi_dbL : f1d(nb,db,chiL), sigL : abs(bL/sin(chiL)/dchi_dbL), [chiL,log(sigL)])$ /* (%i2) [chival,sigval] : sigma_points(1,0.1,50,0.1)$ (%i3) time(%); (%o3) [16.25] (%i4) fll(chival); (%o4) [-2.7468015,-0.0199993,500] (%i5) fll(sigval); (%o5) [-2.712689,15.648312,500] (%i6) chi_min : lmin(chival); (%o6) -2.7468015 (%i7) chi_max : lmax(chival); (%o7) -0.0199993 */ /* homemade list utilities */ fll(x) := [first(x),last(x),length(x)]$ head(mylist) := block([numL,nleft:6], numL : length(mylist), rest (mylist, -(numL - nleft)))$ tail(mylist) := block([numL,nleft:6], numL : length(mylist), rest (mylist, numL - nleft))$ take(%aL,%nn) := (map(lambda([x],part(x,%nn)), %aL))$ fpprintprec:8$