# The function accel() adds up all forces and returns the acceleration
# There is a small inperfection in how the frictional force ffri is treated:
# - initially when the mass is at rest, sign(v)=0 and there will be no friction
# - at the following step, the mass has velocity > 0, and ffri < 0
# - the small v > 0 of the first step will then rapidly decrease
# - in some cases, a small negative velocity may result from the discretization
#    - but then the frictional force will be > 0 and again act to impede the motion
# The result of the imperfection is that the mass initially, when it should
#   really stay at rest, has some very low, in some case oscillating, velocity
# The behavior is "self correcting" for small dt, diminishing when dt is decreased
# - because of the self-correction there is no need to treat v=0 as a special case

# The spring has native (relaxed) length lspr and can be elongated without bounds
#   and compressed to length 0. The final compression force is divergent to avoid the
#   spring length becoming negative (in principle a max length should also be imposed)
# The friction force if fr0 at rest, decays rapidly to fr1 as v increases

 function accel(x,v,t,mass,kspr,lspr,xwall,vwall,kwall,fr0,fr1,fr2)
    fspr=(xwall-x-lspr)*kspr                    # force from spring of native length lspr
    fspr=fspr-kwall/((xwall-x)/lspr+0.0001)^6   # nonlinear repulsive force when close to wall  
    ffri=fr1+(fr0-fr1)*exp(-fr2*v^2)            # friction force magnitude
    ffri=-ffri*sign(v)                          # sign of friction depends on velocity
    return (fspr+ffri)/mass                     # return the acceleration
 end

# The function integrate1() carries out nt steps of the standard leapfrog algorithm,
# while integrate12() uses the modified leapfrog algorithm with better trated friction.
# Results are written to the file "x.dat" every wt steps and also pushed to vectors.
# The step error is not very good with the standard leapfrog method in the presence of,
# but good results can still be obtained quickly with reasonable dt values

 function integrate1(dt,wt,nt,mass,kspr,lspr,vwall,kwall,fr0,fr1,fr2)
    v=0.
    x=0.
    fi=open("x.dat","w")
    vect=Vector{Float64}()
    vecx=Vector{Float64}()
    vecv=Vector{Float64}()
    vecl=Vector{Float64}()
    for i=0:nt
        t=i*dt        
        xwall=lspr+t*vwall 
        if mod(i,wt)==0
           println(fi,t,"  ",x,"  ",xwall-x,"  ",v)           
           push!(vect,t)
           push!(vecx,x)
           push!(vecv,v)
           push!(vecl,xwall-x)
        end	
        v=v+dt*accel(x,v,t,mass,kspr,lspr,xwall,vwall,kwall,fr0,fr1,fr2)
        x=x+dt*v        
    end
    close(fi)
    return vect,vecx,vecv,vecl
 end

 function integrate2(dt,wt,nt,mass,kspr,lspr,vspr,kwall,fr0,fr1,fr2)
    v=0.
    x=0.
    x1=0.
    fi=open("x.dat","w")
    vect=Vector{Float64}()
    vecx=Vector{Float64}()
    vecv=Vector{Float64}()
    vecl=Vector{Float64}()
    for i=0:nt
        t=i*dt        
        xwall=lspr+t*vspr 
        if mod(i,wt)==0
           println(fi,t,"  ",x,"  ",xwall-x,"  ",v)           
           push!(vect,t)
           push!(vecx,x)
           push!(vecv,v)
           push!(vecl,xwall-x)
        end
	u=v+dt*accel(x,v,t,mass,kspr,lspr,xwall,vwall,kwall,fr0,fr1,fr2)
        y=x+dt*u        
        u=0.5*(y-x1)/dt
        v=v+dt*accel(x,u,t,mass,kspr,lspr,xwall,vwall,kwall,fr0,fr1,fr2)
        x1=x
        x=x+dt*v                
    end
    close(fi)
    return vect,vecx,vecv,vecl
 end

# Setting parameters related to the leapfrog integrationloop
 dt=0.05         # time step
 nt=2000         # number of time steps
 wt=5           # write to file every wt step


# Setting parameters related to the system 
# - the mass is connected to a spring, the end of which is connected to a moving wall

 mass=1.          # block mass
 kspr=1.          # spring constant
 lspr=1.          # native spring length
 vwall=0.1        # constant velocity of spring end (at moving wall); a kind of driving
 kwall=0.0001     # constant in repulsive force of spring compressed very close to 0
 fr0=0.5          # friction force at rest
 fr1=0.1          # friction force while in motion
 fr2=20.0         # coefficient for cross-over between fr0 and fr1

# Doing the integration for time step dt and then the same for dt/10,
# plotting t,x,v, or l (the latter being the spring length) after each time.

 t,x,v,l=integrate1(dt,wt,nt,mass,kspr,lspr,vwall,kwall,fr0,fr1,fr2)

 using Plots

 display(plot(t,x))

 dt=dt/10
 nt=nt*10
 wt=wt*10
 
 t,x,v,l=integrate1(dt,wt,nt,mass,kspr,lspr,vwall,kwall,fr0,fr1,fr2)
 display(plot!(t,x))