##               COPYRIGHT NOTICE

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

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




## rutherford_repulse.R
##  translation to R of numerical parts of rutherford_repulse.mac

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


    
##   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
##       chi = pi - 2*phi_inf = 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 = acos(w/sqrt(4+w^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 =  1.107149  rad, or  63.43495  deg
##   theta0 =  2.034444  rad, or  116.5651  deg
##   chi =  0.9272952  rad, or  53.1301  deg
##   rmin =  1.618034 

## analytic scattering angle in degrees for given
##     values of Ebar and bbar 
   
angle_a = function(e,b) 2*asin(1/sqrt(1 + 4*e^2*b^2))*180/pi
    
##  > angle_a(1,1)
##   [1] 53.1301
 



##    numerical scattering angle in degrees for given
##      values of Ebar and bbar for repulsive rutherford case

	  
angle_n = 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)*180/pi}
     
##  > angle_n(1,1)
##     [1] 53.1301



##   Veff_plot for repulsive Rutherford case 

Veff_plot = function(e,b,r0,r1,xmin,xmax,ymin,ymax) {    
    w = 1/e/b
    rmin = b*(w + sqrt(4+w^2))/2
    cat(" rmin = ",rmin,"\n")
    curve( 1/r + 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)}
    
    
    	  
##  > Veff_plot(1,1,0.9,5,0.8,5,0,2)
##    rmin =  1.618034 
 






##  TRAJECTORY PLOT FUNCTIONS


##  init(e,b)  specific to repulsive 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   # repulsive case
    cat(" rmin = ",rmin,"\n")
    phi_inf = acos(w/sqrt(4+w^2))   # repulsive 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 =  1.618034 
##   phi_inf =  1.107149  rad or  63.43495  deg
##   theta0 =  2.034444  rad or  116.5651  deg
##   chi =  0.9272952  rad, or  53.1301  deg
##   xc =  -0.7236068  yc =  1.447214 
##   vcx =  0.5527864  vcy =  0.2763932 
##   xa =  -0.5 
##  > chi
##  [1] 0.9272952
##  > xc
##  [1] -0.7236068
##  > yc
##  [1] 1.447214
##  > vcx
##  [1] 0.5527864
##  > vcy
##  [1] 0.2763932
##  > xa
##  [1] -0.5
##  > rmin
##  Error: object 'rmin' not found

##   feb. 5 version:
##   orbit_plot1 for repulsive 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   # repulsive case
        dy = vy
        dvy = y/2/e/r^3   # repulsive 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)}
    
##  > library(deSolve)
##  > orbit_plot1(1,1,5.5,5,0.01,5,-4.54,1.87,0,4)
##   rmin =  1.618034 
##   phi_inf =  1.107149  rad or  63.43495  deg
##   theta0 =  2.034444  rad or  116.5651  deg
##   chi =  0.9272952  rad, or  53.1301  deg
##   xc =  -0.7236068  yc =  1.447214 
##   vcx =  0.5527864  vcy =  0.2763932 
##   xa =  -0.5 
##    backwards from xc, yc 
##   xfirst =  -4.982187  
##    forwards from xc, yc 
##   xlast =  1.872902  ylast =  4.265932 
##   vx_last =  0.5421801  vy_last =  0.7009998 
##   xf =  2.5  yf =  5 




    

##   PLOT OF DIFFERENTIAL CROSS SECTION VS SCATT ANGLE    

## scattering angle in radians via numerical methods
## for repulsive 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}

     
##  > chi = achi(1,1)
##  > chi
##  [1] 0.9272952
##  > deg(chi)
##  [1] 53.1301


##    f1d(nv,hh,gL) returns 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}



## Example of use of f1d: getting numerical derivatives
## of sin(x) at grid points.

## divide interval (0,pi/2) into 100 subintervals
## h  = step size = (pi/2 - 0)/100 = pi/200
## number of data values = 101

## > h = pi/200; h
## [1] 0.01570796
## > xL = seq(0, pi/2, by = h); head(xL)
## [1] 0.00000000 0.01570796 0.03141593 0.04712389 0.06283185 0.07853982
## > length(xL)
## [1] 101

## > fL = sin(xL); head(fL)
## [1] 0.00000000 0.01570732 0.03141076 0.04710645 0.06279052 0.07845910
## > fpvec = f1d(101,h,fL); head(fpvec)
## [1] 1.0002056 0.9998355 0.9994655 0.9988488 0.9979857 0.9968763
## > tail(fpvec)
## [1] 7.845587e-02 6.278794e-02 4.710451e-02 3.140947e-02 1.570667e-02
## [6] 3.875386e-06

##  maxima output
##  (%i16) head(fplist);
##  (%o16) [1.00021,0.9998,0.9995,0.9988,0.998,0.9969]

##  (%i17) tail(fplist);
##  (%o17) [0.07846,0.06279,0.0471,0.03141,0.01571,3.875386E-6]

## > mylist = list(c(1,2),c(3,4))
## > mylist[[1]]
## [1] 1 2
## > mylist[[2]]
## [1] 3 4



##    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))}
    

    
## > spts = sigma_points(1,0.1,50,0.5)
## > chival = spts[[1]]
## > sigval = spts[[2]]
## > fll(chival)
##   2.746802   0.02016061   100 
## > fll(sigval)
##   -2.446814   15.6167   100 
##  > plot(chival, sigval, xlim = c(0,1), ylim= c(0,20),
##  +         type="l", col = "blue", lwd = 3,
##  +          xlab="chi", ylab="log(sigma)")
##  > abline(h=0, v=0)
##  > curve(log(1/16/sin(x/2)^4),0,1,n=200,add=TRUE,col = "red",lwd = 3)