## COPYRIGHT NOTICE ## Copyright (C) 2014, 2015 Edwin L. (Ted) Woollett ## http://www.csulb.edu/~woollett ## rutherford_attract.R 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 ## print out first, last and length of a vector fll = function(xL) { xlen = length(xL) cat(" ",xL[1]," ",xL[xlen]," ",xlen,"\n") } ## convert radians to degrees deg = function(z) z*180/pi ## > deg(1) ## [1] 57.29578 ## fill in intermediate points along a straight line ## N = number of intervals, N+1 = number of points make_pts = function(x1,y1,xf,yf,N) { h = (xf - x1)/N slope = (yf-y1)/(xf-x1) xL = seq(x1, xf, by = h) yL = y1 + slope*(xL - x1) list(xL, yL)} ## > pts = make_pts(0,0,4,4,8) ## > xL = pts[[1]]; xL ## [1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 ## > yL = pts[[2]]; yL ## [1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 ## > 1 + yL; ## [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 ## angles (e,b) ## attractive rutherford case ## some (e,b) dependent angles and rmin using analytic results. ## 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 = 2*phi_inf - %pi = 2*asin(w/sqrt(w^2+4)) ## rmin is the distance of closest approach ## to the scattering center = ( x = 0, y = 0 ) angles = function(e,b) { w = 1/e/b phi_inf = asin(w/sqrt(w^2+4)) + pi/2 cat(" phi_inf = ",phi_inf," rad, or ",phi_inf*180/pi," deg\n") theta0 = pi - phi_inf cat(" theta0 = ",theta0," rad, or ",theta0*180/pi," deg\n") chi = 2*asin(w/sqrt(w^2+4)) cat(" chi = ",chi," rad, or ",chi*180/pi," deg\n") rmin = b*(-w + sqrt(4+w^2))/2 cat(" rmin = ",rmin,"\n")} ## > angles(1,1) ## phi_inf = 2.034444 rad, or 116.5651 deg ## theta0 = 1.107149 rad, or 63.43495 deg ## chi = 0.9272952 rad, or 53.1301 deg ## rmin = 0.618034 ## angle_n(e,b) attractive rutherford case ## numerical scattering angle in degrees for given ## values of Ebar and bbar angle_n = function(e,b) deg(achi(e,b)) ## > angle_n(1,1) ## [1] -53.1301 ## TRAJECTORY PLOT FUNCTIONS ## init(e,b) specific to attractive rutherford case, ## uses analytic expressions to produce 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 = function(e,b) { w = 1/e/b rmin = b*(-w + sqrt(4+w^2))/2 # attractive case cat(" rmin = ",rmin,"\n") phi_inf = asin(w/sqrt(w^2+4)) + pi/2 # attractive case cat(" phi_inf = ", phi_inf," rad or ", phi_inf*180/pi," deg\n") theta0 = pi - phi_inf cat(" theta0 = ",theta0," rad or ", theta0*180/pi," deg\n") chi <<- pi - 2*phi_inf # global parameter chi in radians cat(" chi = ", chi, " rad, or ", chi*180/pi," deg\n") # 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 cat(" xc = ", xc," yc = ", yc,"\n") cat(" vcx = ", vcx," vcy = ", vcy,"\n") # x-intersection of rmin line with y=b line xa <<- xc*b/yc # global cat(" xa = ",xa,"\n") } ## > init(1,1) ## rmin = 0.618034 ## phi_inf = 2.034444 rad or 116.5651 deg ## theta0 = 1.107149 rad or 63.43495 deg ## chi = -0.9272952 rad, or -53.1301 deg ## xc = 0.2763932 yc = 0.5527864 ## vcx = 1.447214 vcy = -0.7236068 ## xa = 0.5 ## feb. 5 version: ## orbit_plot1 for attractive rutherford potential. ## calls init(e,b) to define local values of ## xc,yc,xa,vcx,vcy,chi. orbit_plot1 = function(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax) { # define local xc, yc, vcx, vcy, xa, chi init(e,b) # symbolic expressions for acceleration components trajec = function(t, y, parms) { with( as.list(y), { r = sqrt(x^2 + y^2) dx = vx dvx = -x/2/e/r^3 # attractive case dy = vy dvy = -y/2/e/r^3 # attractive case list( c(dx, dvx, dy, dvy) ) } ) } # integrate backwards from xc, yc cat(" backwards from xc, yc \n") yini = c(x = xc, vx = vcx, y = yc, vy = vcy) times = seq(0, -tm, -dt) out = ode(times = times, y = yini, func = trajec, parms = NULL) xL = out[,"x"] yL = out[,"y"] xfirst = tail(xL, n=1) cat(" xfirst = ", xfirst," \n") plot(xL, yL, xlim = c(xmin, xmax), ylim = c(ymin, ymax), type = "l", col = "blue", lwd = 3, xlab = "X", ylab = "Y") # integrate forwards from xc,yc cat(" forwards from xc, yc \n") times = seq(0, tp, dt) out = ode(times = times, y = yini, func = trajec, parms = NULL) xL = out[,"x"] yL = out[,"y"] cat(" xlast = ", tail(xL, n=1)," ylast = ",tail(yL, n=1),"\n") vxf = tail(out[,"vx"], n=1) vyf = tail(out[,"vy"], n=1) cat (" vx_last = ",vxf," vy_last = ",vyf,"\n") lines(xL, yL, lwd = 3, col = "blue") # add line y = b abline(h = b, lwd=2) abline(h = 0) # add tick mark at origin lines ( c(0,0), c(- 0.05,0.05), lwd=2) # add rmin line lines ( c(0,xc), c(0,yc), col = "red", lwd = 2) # add chi line yf = b + rchi*sin(chi) xf = xa + rchi*cos(chi) cat(" xf = ",xf," yf = ",yf,"\n") lines ( c(xa,xf), c(b, yf), lwd=2) # extend rmin line lines(c(xc,xa), c(yc,b),col = "red", lty = 3,lwd=2) points(xa, b, col="red", pch=19, cex = 0.5)} ## > source("rutherford_attract.R") ## > library(deSolve) ## > orbit_plot1(1,1,5,5,0.01,5,-3,3.4,-2,2) ## rmin = 0.618034 ## phi_inf = 2.034444 rad or 116.5651 deg ## theta0 = 1.107149 rad or 63.43495 deg ## chi = -0.9272952 rad, or -53.1301 deg ## xc = 0.2763932 yc = 0.5527864 ## vcx = 1.447214 vcy = -0.7236068 ## xa = 0.5 ## backwards from xc, yc ## xfirst = -5.727665 ## forwards from xc, yc ## xlast = 4.204644 ylast = -4.006098 ## vx_last = 0.6550977 vy_last = -0.861996 ## xf = 3.5 yf = -3 ## PLOT OF DIFFERENTIAL CROSS SECTION VS SCATT ANGLE ## scattering angle in radians via numerical methods ## for attractive Rutherford scattering achi = function(e,b) { w = 1/(b*e) fr = function(z) z^2 + w*z -1 # root -> zmin zmin = uniroot(fr, c(1e-4, 1e6),tol = 1e-10)$root phi_inf = integrate(function(z) 1/z/sqrt(fr(z)), zmin, Inf)$val pi - 2*phi_inf} ## > achi(1,1) ## [1] -0.9272952 ## f1d(nv,hh,gL) returns a vector of first derivatives ## at the nv grid points for a function whose ## nv values at the grid points separated by hh ## are in the vector gL. Would be more accurate if ## we used quadratic interpolation to define the ## end of grid derivatives. f1d = function(nv,hh,gL) { fpL = vector(mode = "numeric", length = nv) for (j in 2:(nv-1)) fpL[j] = (gL[j+1] - gL[j-1])/2/hh fpL[1] = 2*fpL[2] - fpL[3] fpL[nv] = 2*fpL[nv-1] - fpL[nv-2] fpL} ## sigma_points(e,b0,bmax,db) produces a list of two lists: ## [chi-list, log(d_sig/d_omega)-list] using numerical methods. ## calls achi() and f1d() sigma_points = function(ee,b0,bmax,db) { bL = seq(from=b0, to=bmax, by=db) nb = length(bL) chiL = sapply(bL, function(x) achi(ee,x)) dchi_dbL = f1d(nb, db, chiL) sigL = abs(bL/sin(chiL)/dchi_dbL) list( chiL,log(sigL))} #################################### acomment = " > source("rutherford_attract.R") > spts = sigma_points(1,0.1,50,0.1) > chival = spts[[1]] > fll(chival) -2.746802 -0.01999933 500 > sigval = spts[[2]] > fll(sigval) -2.712689 15.64831 500 > chi_min = min(chival); chi_min [1] -2.746802 > chi_max = max(chival); chi_max [1] -0.01999933 " #####################################