##               COPYRIGHT NOTICE

##  Copyright (C) 2014, 2015 Edwin L. (Ted) Woollett
##  http://www.csulb.edu/~woollett

##  lennard_jones.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



##  list utility: 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



## make a plot of the arg of the radical -
##  largest positive zero is rmin            

fplot = function(e,b,xmin,xmax,ymin,ymax) {
    curve(e*x^12 - e*b^2*x^10 + 4*x^6 - 4, xmin, xmax, 
        n=200, col="blue", lwd=3, ylim = c(ymin,ymax),
        xlab = "r", ylab = "f")
    abline(h = 0)}

##  > fplot(1,1,0.01,2,-10,10)

##   rmin = largest positive root of arg of radical
          
##  > e = 1;  b = 1
##  > uniroot (function(x) e*x^12 - e*b^2*x^10 + 4*x^6 - 4, c(0.1,2),
##  +           tol = 1e-10 )$root
##  [1] 1

#########################################
acomment = "
  Maxima rearrangement of function whose root is sought:
  
(%i3) root_expr : 1 - 4*(1/r^12 - 1/r^6)/e - b^2/r^2$
(%i5) solve(root_expr,r);
(%o5) [0 = e*r^12-b^2*e*r^10+4*r^6-4]
(%i6) root_expr : rhs(%[1]);
(%o6) e*r^12-b^2*e*r^10+4*r^6-4 "
#########################################


##    angle_n(e,b,x1,x2) 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))  
##    For given e and b, one needs to use fplot (above)
##    to make sure the limits for rmin search (x1,x2)
##    bracket the largest real root. x1 should not be zero
##    since the first step used by uniroot is to evaluate
##    the given function at the end points to see if the
##    function sign is different.    

angle_n = function(e,b,x1,x2) {
     fun = function(r) 1-4*(1/r^12 -1/r^6)/e -b^2/r^2
     rmin = uniroot (function(r) e*r^12 - b^2*e*r^10 + 4*r^6 - 4,
               c(x1,x2), tol = 1e-16 )$root
     phi_inf = b*integrate(function(r) 1/r^2/sqrt(fun(r)), rmin, Inf)$val
     (pi - 2*phi_inf)*180/pi}
     
##  > angle_n(1,1,1e-6,3)
##   [1] 57.11997
##  > sapply(c(1,0.8), function(x) angle_n(1,x,1e-6,3))
##  [1] 57.11997 85.63091



##  Veff_plot for lennard-jones case 
##  last two args: x1 and x2: see angle_n notes above.

Veff_plot = function(e,b,r0,r1,xmin,xmax,ymin,ymax,x1,x2) {
    rmin = uniroot (function(r) e*r^12 - b^2*e*r^10 + 4*r^6 - 4,
               c(x1,x2), tol = 1e-16 )$root
    cat(" rmin = ",rmin,"\n")
    curve( 4*(r^(-12) - r^(-6)) + e*b^2/r^2, r0, r1,
       xname="r", n=200, col="blue", ylim = c(ymin,ymax),
       xlim = c(xmin,xmax), lwd=3, xlab = "r", ylab = "Veff")
    lines(c(rmin,xmax), c(e,e), col = "red", lwd = 3)
    abline(h = 0) }
    
    
##  > Veff_plot(1,1,0.9,3,0.8,3,-0.5,2,1e-6,3)
##   rmin =  1 

    

##  init(e,b,x1,x2)  specific to lennard-jones case,
##    uses numerical methods to produce global definitions
##    of chi,xc,yc, vcx,vcy,and xa. Besides printing out 
##    these values, init(e,b,x1,x2) also prints 
##    out the values of rmin, phi_inf, and theta0.
##    see angle_n notes above for x1 and x2
    

init = function(e,b,x1,x2) {
    rmin = uniroot (function(r) e*r^12 - b^2*e*r^10 + 4*r^6 - 4,
               c(x1,x2), tol = 1e-16 )$root  # lennard-jones case
    cat(" rmin = ",rmin,"\n")
    fun = function(r) 1-4*(1/r^12 -1/r^6)/e -b^2/r^2   # lennard-jones case
    phi_inf = b*integrate(function(r) 1/r^2/sqrt(fun(r)), rmin, Inf)$val    
    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") }
    
##########################################
acomment = "    
    
> init(1,1,1e-6,3)
 rmin =  1 
 phi_inf =  1.072331  rad or  61.44001  deg
 theta0 =  2.069262  rad or  118.56  deg
 chi =  0.9969316  rad, or  57.11997  deg
 xc =  -0.4780786  yc =  0.8783171 
 vcx =  0.8783171  vcy =  0.4780786 
 xa =  -0.5443121 
 > chi
[1] 0.9969316  "

####################################### 
      


##   orbit_plot1 for lennard-jones potential.
##   calls init(e,b,x1,x2) to define local values of
##     xc,yc,xa,vcx,vcy,chi. 
##   values of x1 and x2 are hardwired into this function.
##  If results are not reasonable, fiddle with x1 and x2.

orbit_plot1 = function(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax) {
       
     # define local xc, yc, vcx, vcy, xa, chi 
     x1 = 1e-8
     x2 = 10
    init(e,b,x1,x2)
    
    #  symbolic expressions for acceleration components 
    
    trajec = function(t, y, parms) {
     with( as.list(y), {
        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         
        dx = vx        
        dvx = -12*x*(1/r8 - 2/r14)/e   #   dvx/dt 
        dy = vy        
        dvy = -12*y*(1/r8 - 2/r14)/e   #   dvy/dt         
        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)}    



#############################################

acomment = "

Example of use for e = 1, b = 1

> orbit_plot1(1,1,5,5,0.01,5,-2,1.2,0,2)
 rmin =  1 
 phi_inf =  1.072331  rad or  61.44001  deg
 theta0 =  2.069262  rad or  118.56  deg
 chi =  0.9969316  rad, or  57.11997  deg
 xc =  -0.4780786  yc =  0.8783171 
 vcx =  0.8783171  vcy =  0.4780786 
 xa =  -0.5443121 
  backwards from xc, yc 
 xfirst =  -5.662456  
  forwards from xc, yc 
 xlast =  2.234246  ylast =  5.298241 
 vx_last =  0.542906  vy_last =  0.8398568 
 xf =  2.170097  yf =  5.199046   "
 
###############################################

##   PLOT OF DIFFERENTIAL CROSS SECTION VS SCATT ANGLE    

## scattering angle in radians via numerical methods
## for lennard-jones scattering
	  
achi = function(e,b,x1,x2) {
    rmin = uniroot (function(r) e*r^12 - b^2*e*r^10 + 4*r^6 - 4,
               c(x1,x2), tol = 1e-16 )$root  # lennard-jones case    
    fun = function(r) 1-4*(1/r^12 -1/r^6)/e -b^2/r^2   # lennard-jones case
    phi_inf = b*integrate(function(r) 1/r^2/sqrt(fun(r)), rmin, Inf)$val    
    pi - 2*phi_inf}

    
##  > achi(1,1,1e-6,3)
##  [1] 0.9969316

##   find_b( b_vec, chi_vec )
##    given two equal length vectors, b_vec and a vector
##   of the corresponding chi values, chi_vec, find the
##   value of b corresponding to the minimum value of chi. 

find_b = function(bvec,chi_vec) bvec[which.min(chi_vec)]

#################################################
acomment = "
> bL = seq(0.1,0.6,0.1); bL
[1] 0.1 0.2 0.3 0.4 0.5 0.6
> chiL = sin(bL); chiL
[1] 0.09983342 0.19866933 0.29552021 0.38941834 0.47942554 0.56464247
> min(chiL)
[1] 0.09983342
> find_b(bL,chiL)
[1] 0.1   "
######################################################

##   chi_sigma1(e,b,db)
##    calculate the  vector c(chival,log(sig)) for
##   one impact parameter, given the small
##   interval db around b. x1 and x2 hard-wired.
   
chi_sigma1 = function(e,b,db) {
   x1 = 1e-8
   x2 = 10
   chival = achi(e,b, x1, x2)
   dchi_db = (achi(e,b+db, x1, x2) - achi(e,b-db, x1, x2))/2/db
   sig = abs(b/sin(chival)/dchi_db)
   c(chival, log(sig)) }
   
############################################################
acomment = "
> chi_sigma1(1,0.1,0.01)
[1]  2.947505 -1.321683
> bL = c(0.1, 0.2)
> sapply (bL, function(x) chi_sigma1(1,x,0.01))
          [,1]      [,2]
[1,]  2.947505  2.752433
[2,] -1.321683 -1.312742  "
##########################################################



    
##    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(). x1 and x2 hard-wired into this code.

    
sigma_points = function(ee,b0,bmax,db) {
    x1 = 1e-8
    x2 = 10
    bL = seq(from=b0, to=bmax, by=db)
    nb = length(bL)
    chiL = sapply(bL, function(x) achi(ee,x,x1,x2))
    dchi_dbL = f1d(nb, db, chiL)
    sigL = abs(bL/sin(chiL)/dchi_dbL)
    list( chiL,log(sigL))}
    
##################################################
acomment = "    
> spts = sigma_points(1,0.1,1.3,0.025)
> chival = spts[[1]]
> fll(chival)
  2.947505   0.04908547   49 
> sigval = spts[[2]]
> fll(sigval)
  -1.321405   1.924964   49 
> chi_min = min(chival); chi_min
[1] 0.04908547
> chi_max = max(chival); chi_max
[1] 2.947505
> plot(chival, sigval, xlim = c(chi_min,chi_max),
+         type="l", col = "blue", lwd = 3,
+           xlab="chi", ylab="log(sigma)")
> abline(h=0, v=0) "
#####################################################