##               COPYRIGHT NOTICE

##  Copyright (C) 2014 Edwin L. (Ted) Woollett
##  http://www.csulb.edu/~woollett
##    example2.R is a supplement to Ch.2,
##      Computational Physics with Maxima or R: Initial Value Problems
##     See example2.pdf for discussion

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

##  example2.R, ted woollett, may 2014

##  global potential energy V(x,y)

V = function(x,y) (x^2 + y^2)/2 + y*x^2 - y^3/3

##  > V(2,2)
##  [1] 9.333333


##  yminmax(E) assumes KE = 0 and x = 0
##    to find limits on initial y from energy conservation 

yminmax = function(E) {
    mroots = sort(Re(polyroot(c(E,0,-1/2,1/3))))
    mroots[1:2]}
    
##  > sort(Re(polyroot(c(0.1,0,-1/2,1/3))))
##  [1] -0.3976099  0.5670689  1.3305410    
##  > yminmax(0.1)
##  [1] -0.3976099  0.5670689
    
gety0 = function(E) {
    ylim = yminmax(E)
    ymin = ylim[1]
    ymax = ylim[2]
    
    ok = FALSE
    while (! ok) {
       cat(" need ",ymin," < y0 < ",ymax,"\n")
       y0 = as.numeric(readline(" input y0: "))
       if ((y0 > ymin) && (y0 < ymax)) ok = TRUE
       if (y0 < -10) ok = TRUE}
    y0}
    
##  > gety0(0.1)
##  need  -0.3976099  < y0 <  0.5670689 
##   input y0: -1
##   need  -0.3976099  < y0 <  0.5670689 
##   input y0: 1
##   need  -0.3976099  < y0 <  0.5670689 
##   input y0: 0
##  [1] 0    


get_py0 = function(pymax) {
    ok = FALSE
    while (! ok) {
        cat(" need ",0.0," < py0 < ",pymax,"\n")
        py = as.numeric(readline(" input py0 "))
        if ((py > 0) && (py < pymax))  ok = TRUE
        if (py < 0)  ok =TRUE}
    py}
    
##  > get_py0(0.6)
##   need  0  < py0 <  0.6 
##   input py0 1
##   need  0  < py0 <  0.6 
##   input py0 0.4
##  [1] 0.4

##  some tools

##  > v1 = c(1,2,3)
##  > v2 = c(4,5,6,7)
##  > myL = list(v1,v2)
##  > test = sapply(myL,last)
##  > is.vector(test)
##  [1] TRUE
##  > test[1]
##  [1] 3
##  > test[2]
##  [1] 7





##  global derivatives function in proper form for myrk4
##  y[1] = x, y[2] = y, y[3] = px, y[4] = py

henon_heiles = function(t,y) {
  with( as.list(y),  {
     dx = y[3]
     dy = y[4]
     dpx = -2*y[1]*y[2] - y[1]
     dpy = y[2]^2 - y[2] - y[1]^2
     c(dx,dy,dpx,dpy)})}
     
## trajectory(E) returns list(xL,yL,pxL,pyL)


trajectory = function(E) {       
      x0 = 0
      if ((E <= 0) || (E > 1/6)) {
        cat(" need 0 <= E <= ", 1/6,"\n")
        cat(" aborting program \n")
        return()}
            
##   get y0  
      
     y0 = gety0(E)
     if (y0 < -10) {
         cat(" aborting program \n ")
         return()}
     
##   get py0. V(x,y) is global potential function

     py0max = sqrt(2*(E - V(x0, y0)))
     py0 = get_py0(py0max)
     if (py0 < 0) {
         cat(" aborting program \n")
         return()}
         
     px0 = sqrt(2*(E - V(x0,y0)) - py0^2)
     cat(" E = ",E," x0 = ",x0," y0 = ",y0,"\n")
     cat(" px0 = ", px0," py0 = ", py0, "\n")
     
## defined globally: tL and henon_heiles 
 
     initL = c(x0,y0,px0,py0)
     myrk4(initL,tL,henon_heiles)}

##  > tL = seq(0,3,0.1)
##  > rkvecs = trajectory(0.1)
##   need  -0.3976099  < y0 <  0.5670689 
##   input y0: 0.095
##   need  0  < py0 <  0.4376604 
##   input py0 0.096
##   E =  0.1  x0 =  0  y0 =  0.095 
##   px0 =  0.4270019  py0 =  0.096 
##  > is.list(rkpts)
##  [1] TRUE
##  > xL = rkvecs[[1]]
##  > yL = rkvecs[[2]]
##  > plot(xL,yL,type="l",lwd=2,xlab="x",ylab="y")

##  > get_last(rkvecs)
##  [1]  0.003305691 -0.243223862 -0.348455075 -0.099110096
##  > last(xL)
##  [1] 0.003305691
##  > last(yL)
##  [1] -0.2432239
##  > xLs = sign(xL)
##  > match(-1,xLs)
##  [1] NA
##  > if ( is.na(match(-1,xLs)) ) {cat(" yes \n")} else {cat(" no \n")}
##   yes

##  > t.x.y = data.frame(t=tL, x = xL, y = yL)
##  > head(t.x.y)
##      t          x         y
##  1 0.0 0.00000000 0.0950000
##  2 0.1 0.04261483 0.1041560
##  3 0.2 0.08471409 0.1123586
##  4 0.3 0.12577599 0.1194901
##  5 0.4 0.16528054 0.1254097
##  6 0.5 0.20271958 0.1299584

##  > ypy = getsection(rkvecs,0.01)
##  > ypy[[1]]
##  [1] 0.095
##  > ypy[[2]]
##  [1] 0.096

##  rk41 is almost the same as myrk4, except that
##  rk41 allows the step size to vary arbitrarily
##  as defined by grid: the vector of values  of
##  the independent variable.
##  If one dep var, rk41 returns vector yL 
##  yL is soln vector.
##  If 2+ dep var, rk41 returns  list(xL,yL,...)   


rk41 = function(init, grid ,func)  {
   num.var = length(init)
   solnList = list()   
   n.grid = length(grid)
   for (k in 1:num.var) {
        solnList[[k]] = vector(length = n.grid)
        solnList[[k]][1] = init[k] }   
   yL = vector(length = num.var)  # solution at beginning of each step
   for (j in 1:(n.grid-1)) {
      for (k in 1:num.var) yL[k] = solnList[[k]][j]
      h = grid[j+1] - grid[j]   # step size 
      k1 = func(grid[j], yL)  #  vector of derivatives
      k2 = func(grid[j] + h/2, yL + h*k1/2)  #  vector of derivatives
      k3 = func(grid[j] + h/2, yL + h*k2/2)  #  vector of derivatives
      k4 = func(grid[j] + h, yL + h*k3)  #  vector of derivatives
      for (k in 1:num.var) solnList[[k]][j+1] = solnList[[k]][j] +
                              h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 }
   if (num.var==1) solnList[[1]] else solnList}
   
   
##  mgrid adjusts the number of grid points and the
##  last interval size so that  the last point is the final point tf
   
mgrid = function(t0, tf, dt) {
      small = 1e-12
      aL = seq(t0, tf, dt)
      if (abs(tf - last(aL)) < small) return(aL)
      append(aL, tf)}
      
##  > last = function(xL) { xL[ length(xL)] }
##  > mgrid(0,1,0.1)
##   [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
##  > mgrid(0,1,0.12)
##   [1] 0.00 0.12 0.24 0.36 0.48 0.60 0.72 0.84 0.96 1.00
##  > mgrid(0,1,0.094)
##   [1] 0.000 0.094 0.188 0.282 0.376 0.470 0.564 0.658 0.752 0.846 0.940 1.000
##  > mgrid(0,-1,-0.1)
##   [1]  0.0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1.0
##  > mgrid(0,-1,-0.12)
##   [1]  0.00 -0.12 -0.24 -0.36 -0.48 -0.60 -0.72 -0.84 -0.96 -1.00
##  > mgrid(0,-1,-0.094)
##   [1]  0.000 -0.094 -0.188 -0.282 -0.376 -0.470 -0.564 -0.658 -0.752 -0.846
##  [11] -0.940 -1.000






     

##  deriv_gb is used by goback, and gives the
##  derivatives of y, px, and py with respect to x.        
##  y[1] = y, y[2] = px,  y[3] = py
        
deriv_gb = function(x,y) {
    with( as.list(y), {
      dy = y[3]/y[2]
      dpx = -(2*x*y[1] + x)/y[2]
      dpy = (y[1]^2 - y[1] - x^2)/y[2]
      c( dy, dpx, dpy )})}
      
##  goback modified may 30, 2014 to call mgrid and rk41      
##  goback: calls rk41 to integrate backwards with respect to
##    the independent variable x the three dependent
##     variables: (y, px, py). The goback arg
##      rkpoint is a vector c(xi,yi,pxi,pyi).
##      goback returns the vector c(y0,py0) corresponding to x=0.

goback = function(rkpoint,xprec) {
    init_gb = rkpoint[2:4]
    ##  cat(" init_gb = ",init_gb,"\n")
    xi = rkpoint[1]
    xprec = min (xprec, abs(xi)/10)
    ##  cat(" xprec = ",xprec,"\n")
    if (xi > 0)  {dx_gb = -xprec} else {dx_gb = xprec}
    ##  cat(" xi = ",xi," dx_gb = ",dx_gb,"\n")
    xL_gb = mgrid(xi,0,dx_gb)
    ##  cat(" fll(xL_gb) = ", first(xL_gb),last(xL_gb),length(xL_gb),"\n")
    out = rk41(init_gb,xL_gb,deriv_gb)
    yL = out[[1]]
    pyL = out[[3]]
    ylast = last(yL)
    pylast = last(pyL)
    c(ylast, pylast)}

## following test cases have debug printouts turned on.
##  case xi > 0 : pnt generated using trajectory
##  with E = 0.1, y0=0.095, py0=0.096, and integrating
##  from t=0 to t=3 with dt = 0.1
##  > pnt
##  [1]  0.0033057 -0.2432200 -0.3484600 -0.0991100    
##  > goback(pnt,0.01)
##   init_gb =  -0.24322 -0.34846 -0.09911 
##   xprec =  0.00033057 
##   xi =  0.0033057  dx_gb =  -0.00033057 
##  fll(xL_gb) =  0.0033057 0 11 
##  [1] -0.24414658 -0.09623499

##  case xi < 0 : pnt generated using trajectory
##  with E = 0.1, y0=0.095, py0=0.096, and integrating
##  from t=0 to t=3.1 with dt = 0.1
##  > get_last(rkvecs)
##  [1] -0.03151891 -0.25160089 -0.34774772 -0.06823856
##  > goback(get_last(rkvecs),0.01)
##  [1] -0.24415043 -0.09623503
##  > ypy = getsection(rkvecs,0.01)
##   nzero =  0 
##   nrle =  2 ns =  2
##  > ypy[[1]]
##  [1]  0.0950000 -0.2441504
##  > ypy[[2]]
##  [1]  0.09600000 -0.09623503


##  tools for getsection
##  > myvec = c(4,3,2,1,0,-1,-2,-3)
##  > myvec_sign = sign(myvec); myvec_sign
##  [1]  1  1  1  1  0 -1 -1 -1
##  > match(-1,myvec_sign)
##  [1] 6
##  > length(which(myvec_sign == 0))
##  [1] 1
##  > which(myvec_sign == 0,arr.ind=TRUE)
##  [1] 5
##  > length(which(myvec_sign == 0,arr.ind=TRUE))
##  [1] 1
##  > rle(myvec_sign)
##  Run Length Encoding
##    lengths: int [1:3] 4 1 3
##    values : num [1:3] 1 0 -1


    
    
##  getsection;  returns list(yLs,pyLs)
##  arg rkvecs has form list(xL,yL,pxL,pyL) as returned by
##  myrk4.
##  with px0 > 0, x0 = 0, x starts out positive. 

getsection = function(rkvecs,xtol) {
    xL = rkvecs[[1]]    
    xLs = sign(xL) # elements are 1, 0, and -1
    
    ## catch case x is never negative 
    if (is.na(match(-1, xLs))) return( list(rkvecs[[2]][1],rkvecs[[4]][1]))
    
    ##  We know x always starts out equal to 0.
    ##  nzero only counts additional x = 0 instances.
    
    nzero = length(which(xLs == 0)) -1
    ##  cat(" nzero = ", nzero, "\n")
    
    ## likewise nrle ignores initial 0
    nrle = length(rle(xLs)$lengths) -1
    ##  cat(" nrle = ", nrle,"\n")
    ns = nrle - nzero
    cat(" ns = ",ns,"\n")
    
    ysL = vector( length = ns )
    pysL = vector( length = ns)
    ysL[1] = rkvecs[[2]][1]
    pysL[1] = rkvecs[[4]][1]     
    
    nx = length(xL)
    nm = -1
    rvecs = rkvecs
    
    ## we ignore presence of 0 here
    for (k in 2:ns) {
       ##  cat(" k = ",k,"\n")
       nc = match(nm,sign(xL))
       ##  cat(" nc = ",nc,"\n")       
       rgoback = goback(get_pnt(rvecs,nc), xtol)
       ysL[k] = rgoback[1]
       pysL[k] = rgoback[2]
       rvecs = part(rvecs,nc+1,nx)
       xL = rvecs[[1]]       
       nx = length(xL)
       nm = - nm}
    list(ysL,pysL)}
    
## three (y,py) values if tf = 6.5
## maxima result: 
##  (%i3) ypy : getsection(rkpts,0.01);
##     working...  
##  (%o3) [[0.095,0.096],[-0.24415,-0.096235],[0.01341,-0.020683]]
##  > tL = seq(0,6.5,0.1)
##  > rkvecs = trajectory(0.1)
##   need  -0.3976099  < y0 <  0.5670689 
##   input y0: 0.095
##   need  0  < py0 <  0.4376604 
##   input py0 0.096
##   E =  0.1  x0 =  0  y0 =  0.095 
##   px0 =  0.4270019  py0 =  0.096 
##  > last.pnt = get_last(rkvecs)
##  > last.pnt
##  [1]  0.03284497  0.01185212  0.44529538 -0.02162857
##  > ypy = getsection(rkvecs,0.01)

##    ns =  3
  
##  > ypy[[1]]
##  [1]  0.0950000 -0.2441504  0.0134099
##  > ypy[[2]]
##  [1]  0.09600000 -0.09623503 -0.02068349

##  > goback(get_last(rkvecs),0.01)
##  [1]  0.01340990 -0.02068349

    
      
      
      
         
              
     



           




##  myrk4 returns a vector if there is only one dependent variable.

##  myrk4 returns a list of vectors if there are two or more
##     dependent variables. When called by trajectory, myrk4
##  returns list(xL,yL,pxL,pyL). When called by goback, myrk4
##  returns list(yL,pxL,pyL).     In the code below,
##  each element of solnList is a vector which contains the grid values
##  of one of the dependent variables. func returns a vector of derivatives.


myrk4 = function(init, grid ,func)  {
   num.var = length(init)
   solnList = list()      
   n.grid = length(grid)
   for (k in 1:num.var) {
        solnList[[k]] = vector(length = n.grid)
        solnList[[k]][1] = init[k] }
   h = grid[2] - grid[1]   # step size 
   yL = vector(length = num.var)  # solution at beginning of each step
   for (j in 1:(n.grid-1)) {
      for (k in 1:num.var) yL[k] = solnList[[k]][j]
      k1 = func(grid[j], yL)  #  vector of derivatives
      k2 = func(grid[j] + h/2, yL + h*k1/2)  #  vector of derivatives
      k3 = func(grid[j] + h/2, yL + h*k2/2)  #  vector of derivatives
      k4 = func(grid[j] + h, yL + h*k3)  #  vector of derivatives
      for (k in 1:num.var) solnList[[k]][j+1] = solnList[[k]][j] +
                              h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 }
   if (num.var==1) solnList[[1]] else solnList}

##  soln to dy/dt = -y with y(0) = 1 over [t,0,1]   
##  > derivs = function(t,y) {-y}
##  > tL = seq(0,1,0.1)
##  > head(tL)
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5
##  > yL = myrk4(init=1,grid=tL,func=derivs)
##  > head(yL)
##  [1] 1.0000000 0.9048375 0.8187309 0.7408184 0.6703203 0.6065309
##  > t.y = data.frame(t=tL,y=yL)
##  > head(t.y)
##      t         y
##  1 0.0 1.0000000
##  2 0.1 0.9048375
##  3 0.2 0.8187309
##  4 0.3 0.7408184
##  5 0.4 0.6703203
##  6 0.5 0.6065309



##  soln to pair of first order ode's, simple harmonic
##  oscillator with unit period, solving for x(t) and vx(t)
##  dx/dt = vx, dvx/dt = -4*pi^2*x, x(0) = 1, vx(0) = 0
##  > sho = function(t,y) with(as.list(y), c(y[2], -4*pi^2*y[1]) )
##  > out = myrk4(c(1,0),tL,sho)
##  > xL = out[[1]]
##  > fll(xL)
##    1   0.9959199   11 
##  > vxL = out[[2]]
##  > fll(vxL)
##    0   0.04406592   11 
##  > t.x.v = data.frame(t=tL,x=xL,vx=vxL)
##  > head(t.x.v)
##      t          x          vx
##  1 0.0  1.0000000  0.00000000
##  2 0.1  0.8091019 -3.68808418
##  3 0.2  0.3101040 -5.96807148
##  4 0.3 -0.3066331 -5.97246738
##  5 0.4 -0.8060469 -3.70144578
##  6 0.5 -0.9979641 -0.02207791


##  integrating forward and then backwards
##  > seq(0,-4,-1)
##  [1]  0 -1 -2 -3 -4
##  > deriv = function(t,y) {-y}
##  > tL = seq(0,1,0.01)
##  > rkpts = myrk4(1,tL,deriv)
##  > fll(rkpts)
##    1   0.3678794   101 
##  > tL = seq(1,0,-0.01)
##  > fll(tL)
##    1   0   101 
##  > rkpts = myrk4(exp(-1),tL,deriv)
##  > fll(rkpts)
##    0.3678794   1   101 
##  > seq(1,0,length=5)
##  [1] 1.00 0.75 0.50 0.25 0.00



mygrid = function() grid(lty = "solid", col = "darkgray")

##  Vector utilities first, last, fll work with R vectors, but
##   not with R lists

first = function(xL) { xL[1] }
last = function(xL) { xL[ length(xL)] }
   
fll = function(xL) {    
    cat(" ",first(xL), " ", last(xL), " ", length(xL), "\n") }
    
##  get_last(list.of.vecs) returns a vector whose elements
##  are the last components of the vectors    
    
get_last = function(rkvecs) sapply(rkvecs, last)

##  > myL = list(c(1,2,3,4),c(5,6,7,8),c(9,10,11,12))
##  > get_last(myL)
##  [1]  4  8 12

    
##  If myList = list(v1, v2, v3, ...), then get_pnt(myList,2)
##  returns the vector c(v1[2], v2[2], v3[2], ...)
##   in other words, using sapply,
##  get_pnt(myList,num) returns a vector whose elements
##  are myList[[k]][num]
    
get_pnt = function(myList, num) {

      ## check that the length of each vector in mylist
      ## is the same value.
      
      if ( length(rle(sapply(myList,length))$values) > 1 ) {
          cat(" vectors do not have the same length \n")
          return()}

      ##  num cannot be greater than the 
      ##  common length of the vectors.
      l1 = length(myList[[1]])
      if (num > l1) {
         cat(" num > ",l1,"\n")
         return()}
         
      ##  myList cannot be empty
      l2 = length(myList)
      if (l2 == 0) return()      
      
      sapply(myList, function(xL) xL[num])}

 
##  > myL = list(c(1,2),c(3,4),c(5,6))
##  > get_pnt(myL,2)
##  [1] 2 4 6

##  if myList = list(aL,bL,cL,,,), part(myList,n1,n2)
##  returns:   list(aL[n1:n2], bL[n1:n2], cL[n1:n2],...)
##  by using lapply.

part = function(myList, n1, n2) {

      if (n2 <= n1) return()

      ## check that the length of each vector in mylist
      ## is the same value.
      
      if ( length(rle(sapply(myList,length))$values) > 1 ) {
          cat(" vectors do not have the same length \n")
          return()}

      ##  n2 cannot be greater than the 
      ##  common length of the vectors.
      l1 = length(myList[[1]])
      if (n2 > l1) {
         cat(" n2 > ",l1,"\n")
         return()}
         
      ##  myList cannot be empty
      l2 = length(myList)  ## number of vectors in myList
      if (l2 == 0) return()      
      
      lapply(myList, function(xL) xL[n1:n2])}            
      
##  > myL = list(c(1,2,3,4),c(5,6,7,8),c(9,10,11,12))
##  > myL1 = part(myL,1,2)
##  > myL1
##  [[1]]
##  [1] 1 2

##  [[2]]
##  [1] 5 6

##  [[3]]
##  [1]  9 10

##   part code without checks

mypart = function(rkvecs,n1,n2) {   
         lapply(rkvecs, function(xL) xL[n1:n2])}
         
##  simple trajectory plot needs external tL defined
   
do_plot = function(E) {  
  rkvecs = trajectory(E)
  xL = rkvecs[[1]]
  yL = rkvecs[[2]]
  plot(xL,yL,type="l",lwd=2,xlab="x",ylab="y")
  mygrid()}
  
##  xsection_plot combines trajectory code with getsection to
##  do a plot and return surface of section list.
##  needs externally defined tL vector.
  
xsection_plot = function(E,y0,py0) {

    if ((E <= 0) || (E > 1/6)) {
        cat(" need 0 < E <= ", 1/6,"\n")
        cat(" aborting program \n")
        return()}
       
    ylim = yminmax(E)
    ymin = ylim[1]
    ymax = ylim[2]
    if ((y0 <= ymin) || (y0 >= ymax)) {
       cat(" need ", ymin, " <= y0 <= ",ymax,"\n")
       return()}
       
    x0 = 0
       
    py0max = sqrt(2*(E - V(x0,y0)))
    if ( (py0 < 0)  || (py0 >= py0max)) {
       cat("  need  0 < py0 < ", py0max,"\n")
       return()}
       
    px0 = sqrt(2*(E - V(0,y0)) - py0^2)
    
    ## defined globally: tL and henon_heiles 
 
     initL = c(x0,y0,px0,py0)
     rkvecs = myrk4(initL,tL,henon_heiles)
     
     ypyL = getsection(rkvecs, 0.01)
     
     plot(ypyL[[1]],ypyL[[2]],pch=19,cex=0.1,xlab="y",
         ylab="py",xlim=c(-0.5,0.5),ylim=c(-0.5,0.5))
    ##   mygrid()
     ypyL}
     

##  > tL = seq(0,1500,0.1)
##  > ypy1 = xsection_plot(0.1,0.095,0.096)
##   ns =  470 
##  > ypy2 = xsection_plot(0.1,0.095,0.03)
##   ns =  480 


## xsection does not check input, does not
##  make a plot, just returns list(yL,pyL)

xsection = function(E,y0,py0) {
       
    x0 = 0                  
    px0 = sqrt(2*(E - V(0,y0)) - py0^2)
    
    ## defined globally: tL and henon_heiles 
 
     initL = c(x0,y0,px0,py0)
     rkvecs = myrk4(initL,tL,henon_heiles)
     
     getsection(rkvecs, 0.01)}
     
     
##  > tL = seq(0,500,0.1)
##  > ypy = xsection(0.1,0.095,0.096)
##   ns =  157 

##  > plot(ypy[[1]],ypy[[2]],pch=19,cex=0.05,xlab="y",
##  +          ylab="py",xlim=c(-0.6,1),ylim=c(-0.6,0.6))

     
plot_sos = function(E,y0,py0, cex_val) {
       
    x0 = 0                  
    px0 = sqrt(2*(E - V(0,y0)) - py0^2)
    
    ## defined globally: tL and henon_heiles 
 
     initL = c(x0,y0,px0,py0)
     rkvecs = myrk4(initL,tL,henon_heiles)
     
     section_plot(rkvecs, 0.01,cex_val)}
     

section_plot = function(rkvecs,xtol, cex_val) {
    xL = rkvecs[[1]]    
    xLs = sign(xL) # elements are 1, 0, and -1
    
    ## catch case x is never negative 
    if (is.na(match(-1, xLs))) return( list(rkvecs[[2]][1],rkvecs[[4]][1]))
    
    ##  We know x always starts out equal to 0.
    ##  nzero only counts additional x = 0 instances.
    
    nzero = length(which(xLs == 0)) -1
    ##  cat(" nzero = ", nzero, "\n")
    
    ## likewise nrle ignores initial 0
    nrle = length(rle(xLs)$lengths) -1
    ##  cat(" nrle = ", nrle,"\n")
    ns = nrle - nzero
    cat(" ns = ",ns,"\n")
    
    ysL = vector( length = ns )
    pysL = vector( length = ns)
    ysL[1] = rkvecs[[2]][1]
    pysL[1] = rkvecs[[4]][1]     
    
    ## plot first SOS point
    plot(ysL[1], pysL[1], pch=19, cex = cex_val, xlab="y",
         ylab="py", xlim=c(-0.6,1), ylim=c(-0.6,0.6))
    
    nx = length(xL)
    nm = -1
    rvecs = rkvecs
    
    ## we ignore presence of 0 here
    for (k in 2:ns) {
       ##  cat(" k = ",k,"\n")
       nc = match(nm,sign(xL))
       ##  cat(" nc = ",nc,"\n")       
       rgoback = goback(get_pnt(rvecs,nc), xtol)
       ysL[k] = rgoback[1]
       pysL[k] = rgoback[2]
       points(ysL[k], pysL[k],pch=19, cex = cex_val)
       Sys.sleep(0.01)
       rvecs = part(rvecs,nc+1,nx)
       xL = rvecs[[1]]       
       nx = length(xL)
       nm = - nm}
    list(ysL,pysL)}