##               COPYRIGHT NOTICE

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

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




## scatt.R generic scatt angle plot

##   scatt_plot1 for generic scattering plot
##    repulsive rutherford potential.
##   calls init(e,b) to define local values of
##     xc,yc,xa,vcx,vcy,chi. 


scatt_plot1 = function(e,b,tm,tp,dt,rchi,xmin,xmax,ymin,ymax) {

     n1 = 150  # vector element for incident rvec arrow
       
     # 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)
    ##  make xLb and yLb global so they can be looked at outside program
    xLb <<- out[,"x"]
    yLb <<- out[,"y"]
    
    x1 = xLb[n1]
    y1 = yLb[n1]
    cat(" x1 = ",x1," y1 = ",y1,"\n")
    xfirst = tail(xLb, n=1)
    cat(" xfirst = ", xfirst," \n")
    plot(xLb, yLb, 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 point at end of rvec arrow
    points(x1,y1,pch=19,col="blue",cex=2)
    
    #  add rvec arrow
    arrows(0,0,x1,y1, code=2, length=0.1,  lwd = 2 )        
        
    # add rmin vector
    arrows(0,0,xc,yc,code=2, length=0.1, col="red",lwd = 2)
    text(0,0.8,expression(r[min]),cex=1.5)
    
    #  show angle theta to rvec line, trial and error:
    plotcircle(mid=c(0,0), from=0, to= 2.1, r=0.5, lwd = 2, arrow=TRUE)
    text(0.8, 0.3, expression(theta),cex=2)   
    text(-1.2,0.6,"r",cex=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)
    
    #  show angle chi
    plotcircle(r=1,mid=c(xa,b), from=0,to=chi, arrow=TRUE)
    text(0.8, 1.4, expression(chi),cex=2)
    
    # show impact parameter b
    arrows(-4.2,0,-4.2,1,code=3,length=0.1,lwd=2)
    text(-4,0.5,"b",cex=1.7)
    }
    

##   production of Fig.3, Scattering Angle Illustrated:
##
##  > source("c:/k1/rutherford_repulse.R")
##  > library(deSolve)
##  > library(shape)
##  > source("c:/k1/scatt.R")
##  > scatt_plot1(1,1,5.5,5,0.01,5,-4.54,1.87,0,4)