/* COPYRIGHT NOTICE Copyright (C) 2014, 2015 Edwin L. (Ted) Woollett http://www.csulb.edu/~woollett lennard_jones.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 */ deg(z) := block([numer],numer:true, z*180/%pi)\$ /* make_pts(x1,y1,xf,yf,N), an example of Maxima list arithmetic methods and use of makelist, fills in intermediate points along a straight line between (x1,y1) and (xf,yf), and returns the list [xL,yL]. N = number of intervals, N+1 = number of points */ 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 ])\$ /* (%i7) make_pts(0,0,4,4,8); (%o7) [[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]] */ /* fplot(e,b,xmin,xmax,ymin,ymax) makes a simple plot of the argument of the radical whose zero is the parameter rmin */ fplot(e,b,xmin,xmax,ymin,ymax) := block([f,numer], numer : true, f : e*x^12 - e*b^2*x^10 + 4*x^6 - 4, plot2d(f,[x, xmin, xmax],[y,ymin,ymax], [style, [lines, 3]], [ylabel, "f"], [xlabel, "r"],[nticks,200]))\$ /* (%i2) fplot(1,1,0.01,2,-10,10)\$ plot2d: some values were clipped. */ fplot_eps(e,b,xmin,xmax,ymin,ymax,fname) := block([f,numer], numer : true, f : e*x^12 - e*b^2*x^10 + 4*x^6 - 4, plot2d(f,[x, xmin, xmax],[y,ymin,ymax], [style, [lines, 5]], [ylabel, "f"], [xlabel, "r"],[nticks,200],[psfile,fname]))\$ /* (%i4) fplot_eps(1,1,0.01,2,-10,10,"c:/k1/lj1.eps"); plot2d: some values were clipped. (%o4) "c:/k1/lj1.eps" */ /* angle_n(e,b) uses numerical methods to return the scattering angle in degrees for the lennard-jones potential. The scattering angle can be either positive or negative (see init(e,b)) */ angle_n(e,b) := block([root_expr,rmin,phi_inf,numer],numer:true, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : map('rhs, realroots(root_expr,1e-15)), rmin : lmax(rmin), phi_inf : b*quad_qagi(1/r^2/sqrt(root_expr),r,rmin,inf)[1], (%pi - 2*phi_inf)*180/%pi)\$ /* (%i24) angle_n(1,1); (%o24) 57.119973 (%i25) bval : [1,0.8]; (%o25) [1,0.8] (%i26) map('lambda([x],angle_n(1,x)),bval); (%o26) [57.119973,85.630912] */ /* anglef(e,b,xm,xp) uses find_root to find signed scattering angle in degrees */ anglef(e,b,xm,xp) := block([root_expr,rmin,phi_inf,numer],numer:true, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : find_root(root_expr,r,xm,xp), phi_inf : b*quad_qagi(1/r^2/sqrt(root_expr),r,rmin,inf)[1], (%pi - 2*phi_inf)*180/%pi)\$ /* (%i46) anglef(e,b,xm,xp) := block([root_expr,rmin,phi_inf,numer],numer:true, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : find_root(root_expr,r,xm,xp), phi_inf : b*quad_qagi(1/r^2/sqrt(root_expr),r,rmin,inf)[1], (%pi - 2*phi_inf)*180/%pi)\$ (%i47) anglef(1,1.6654,1.1,1.3); (%o47) -185.3681857689883 */ /* Veff_plot for lennard-jones case */ Veff_plot(e,b,r0,r1,xmin,xmax,ymin,ymax) := block([r,veff,root_expr,rmin,energy_line,numer],numer:true, veff : 4*(r^(-12) - r^(-6)) + e*b^2/r^2, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : map('rhs, realroots(root_expr,1e-15)), rmin : lmax(rmin), energy_line : [discrete,[[rmin,e],[r1,e]]], plot2d([veff, energy_line],[r,r0,r1], [x,xmin,xmax], [y,ymin,ymax],[style,[lines,3,1],[lines,3,2]], [legend,false],[xlabel,"r"],[ylabel,"Veff"]))\$ /* (%i41) Veff_plot(1,1,0.9,3,0.8,3,-0.5,2); */ Veff_plot_eps(e,b,r0,r1,xmin,xmax,ymin,ymax,fname) := block([r,veff,root_expr,rmin,energy_line,numer],numer:true, veff : 4*(r^(-12) - r^(-6)) + e*b^2/r^2, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : map('rhs, realroots(root_expr,1e-15)), rmin : lmax(rmin), energy_line : [discrete,[[rmin,e],[r1,e]]], plot2d([veff, energy_line],[r,r0,r1], [x,xmin,xmax], [y,ymin,ymax],[style,[lines,5,1],[lines,5,2]], [legend,false],[xlabel,"r"],[ylabel,"Veff"], [psfile, fname]))\$ /* (%i49) Veff_plot_eps(1,1.66597,1,5,1,5,-0.5,2,"c:/k1/lj8.eps"); */ /* init(e,b), with 'root_expr' specific to lennard-jones 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([root_expr,r,rmin,phi_inf,theta0,numer],numer:true, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : map('rhs, realroots(root_expr,1e-15)), rmin : lmax(rmin), print(" rmin = ",rmin), phi_inf : b*quad_qagi(1/r^2/sqrt(root_expr),r,rmin,inf)[1], 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))\$ /* restart maxima, load def: (%i3) init(1,1)\$ rmin = 1 phi_inf = 1.0723305 rad or 61.440014 deg theta0 = 2.0692621 rad or 118.55999 deg chi = 0.996932 rad, or 57.119973 deg xc = -0.478079 yc = 0.878317 vcx = 0.878317 vcy = 0.478079 xa = -0.544312 (%i4) rmin; (%o4) rmin (%i5) phi_inf; (%o5) phi_inf (%i6) theta0; (%o6) theta0 (%i7) chi; (%o7) 0.996932 (%i8) xc; (%o8) -0.478079 (%i9) yc; (%o9) 0.878317 (%i10) vcx; (%o10) 0.878317 (%i11) vcy; (%o11) 0.478079 (%i12) xa; (%o12) -0.544312 */ /* orbit_plot1 for lennard-jones 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,r8,r14,dvxdt,dvydt, xL, yL,pm,pp,vxf,vyf,bline,rmin_line, xf,yf,chi_line,rmin_extend,numer],numer:true, /* define xc, yc, vcx, vcy, xa, chi */ init(e,b), /* symbolic expressions for acceleration components */ r8 : (x^2 + y^2)^4, /* r^8 in terms of x and y */ r14 : (x^2 + y^2)^7, /* r^(14) in terms of x and y */ dvxdt : -12*x*(1/r8 - 2/r14)/e, /* dvx/dt */ dvydt : -12*y*(1/r8 - 2/r14)/e, /* 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]))\$ /* example e = 1, b = 1 (%i4) orbit_plot1(1,1,5,5,0.01,5,-2,1.2,0,2,4,0.4)\$ rmin = 1 phi_inf = 1.0723305 rad or 61.440014 deg theta0 = 2.0692621 rad or 118.55999 deg chi = 0.996932 rad, or 57.119973 deg xc = -0.478079 yc = 0.878317 vcx = 0.878317 vcy = 0.478079 xa = -0.544312 backwards from xc, yc xfirst = -5.6624496 forwards from xc, yc xlast = 2.2342414 ylast = 5.2982522 vx_last = 0.542905 vy_last = 0.839859 xf = 2.1700966 yf = 5.1990458 */ orbit_plot1_eps(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax,nextend,psize,fname) := block([xc,yc,vcx,vcy,xa,chi,rkpts,r8,r14,dvxdt,dvydt, xL, yL,pm,pp,vxf,vyf,bline,rmin_line, xf,yf,chi_line,rmin_extend,numer],numer:true, /* define xc, yc, vcx, vcy, xa, chi */ init(e,b), /* symbolic expressions for acceleration components */ r8 : (x^2 + y^2)^4, /* r^8 in terms of x and y */ r14 : (x^2 + y^2)^7, /* r^(14) in terms of x and y */ dvxdt : -12*x*(1/r8 - 2/r14)/e, /* dvx/dt */ dvydt : -12*y*(1/r8 - 2/r14)/e, /* 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]))\$ /****** DIFFERENTIAL CROSS SECTION PLOT FUNCTIONS ********/ /* f1d(nv,hh,gL) returns a 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. See rutherford_repulse.mac for notes on use. */ 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), here specific to the lennard-jones potential, returns the scattering angle in radians using numerical methods. */ achi(e,b) := block([root_expr,rmin,phi_inf,numer],numer:true, root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2, rmin : map('rhs, realroots(root_expr,1e-15)), rmin : lmax(rmin), phi_inf : b*quad_qagi(1/r^2/sqrt(root_expr),r,rmin,inf)[1], (%pi - 2*phi_inf))\$ /* (%i2) achi(1,1); (%o2) 0.996932 (%i12) find_root(lambda([x],achi(1,x)),x,1.3,1.4); (%o12) 1.3124992 */ /* find_b( b_list, chi_list ) given two equal length lists, a b_list and a list of the corresponding chi values, chi_list, find the value of b corresponding to the minimum value of chi. */ which_min(aL) := block([amin,j,jval:0,numer ],numer:true, amin:lmin(aL), for j thru length(aL) do if is(equal(aL[j], amin)) then ( jval : j, return()), jval)\$ /* (%i2) xL : makelist(x,x,1,10); (%o2) [1,2,3,4,5,6,7,8,9,10] (%i3) which_min(xL); (%o3) 1 */ find_b(xL,yL) := (xL[which_min(yL)])\$ /* (%i6) zL : sin(xL),numer; (%o6) [0.8414709848079,0.90929742682568,0.14112000805987,-0.75680249530793, -0.95892427466314,-0.27941549819893,0.65698659871879,0.98935824662338, 0.41211848524176,-0.54402111088937] (%i10) find_b(xL,zL); (%o10) 5 */ /* chi_sigma1(e,b,db) calculate the list [chival,log(sig)] for one impact parameter, given the small interval db around b. */ chi_sigma1(e,b,db):= block([chival,dchi_db,sig,numer],numer:true, chival : achi(e,b), dchi_db : (achi(e,b+db) - achi(e,b-db))/2/db, sig : abs(b/sin(chival)/dchi_db), [chival,log(sig)])\$ /* (%i12) chi_sigma1(1,0.1,0.01); (%o12) [2.9475046,-1.3216832] (%i13) bL : [0.1,0.2]\$ (%i14) map('lambda([x],chi_sigma1(1,x,0.01)), bL); (%o14) [[2.9475046,-1.3216832],[2.752433,-1.3127418]] */ sigma_plot(e,b0,bmax,db) := block([bL,pts,chiL,chi_min,chi_max, numer],numer:true, bL : makelist(b,b,b0,bmax,db), pts : map('lambda([x],chi_sigma1(1,x,db)), bL), chiL : take(pts,1), chi_min : lmin(chiL), chi_max : lmax(chiL), plot2d([discrete, pts],[x,chi_min,chi_max], [xlabel,"chi"], [ylabel,"ln(dsigma/do)"], [style,[lines,3]]))\$ sigma_plot_eps(e,b0,bmax,db,fname) := block([bL,pts,chiL,chi_min,chi_max, numer],numer:true, bL : makelist(b,b,b0,bmax,db), pts : map('lambda([x],chi_sigma1(1,x,db)), bL), chiL : take(pts,1), chi_min : lmin(chiL), chi_max : lmax(chiL), plot2d([discrete, pts],[x,chi_min,chi_max], [xlabel,"chi"], [ylabel,"ln(dsigma/do)"], [style,[lines,5]],[psfile,fname]))\$ /* 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)])\$ /* (%i20) [chival,sigval] : sigma_points(1,0.1,1.3,0.025)\$ fll(bL) = [0.1,1.3,49] nb = 49 fll(chiL) = [2.9475046,0.0490855,49] */ /* 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\$ ratprint : false\$ quad_control ('control,0)\$