# install package:
#    julia -e 'using Pkg; Pkg.add("ForwardDiff")'
# best run from Julia console as
# include("Aziz-diff2.jl")

using ForwardDiff: derivative
using Printf


# 
# Definition of the function in file aziz.jl in current directory
#push!(LOAD_PATH,pwd()) # add current directory to module search path
#using Aziz
# 
# Another way, avoids mixing new and precompiled module 
include("Aziz.jl")


# define function, 1st and 2nd derivative
V   = x -> Aziz.aziz(float(x))
Vp  = x -> derivative(V, float(x))
Vpp = x -> derivative(Vp, float(x))

# compute values at points r
r  = [0.01:0.5:10.0;]
Vr  = V.(r)
Vpr = Vp.(r)
Vppr = Vpp.(r)

# output
@printf("%15s %15s %30s %30s\n","r","V(r)","V'(r)","V''(r)")
for (a,b,c,d) in zip(r,Vr,Vpr,Vppr)
    @printf("%15.5f %15.5f %30.20f %30.20f\n",a,b,c,d)   
end

# Another way:

#  project https://github.com/mth229, file https://github.com/mth229/MTH229.jl/blob/master/src/MTH229.jl
#  f'(x) will find the derivative of `f` using Automatic Differentation from the `ForwardDiff` package
#
Base.adjoint(f::Function) = x -> derivative(f, float(x)) #  overload ' (adjoint operator)
# also gradient can be overloaded,
# grad(f) = (x, xs...) -> ForwardDiff.gradient(f, vcat(x, xs...))

# V defined earlier
Vpr_new  = V'.(r)
Vppr_new = V''.(r)

# output
println("Compute derivatives using overloaded adjoint")
@printf("%15s %15s %30s %30s\n","r","V(r)","V'(r)","V''(r)")
for (a,b,c,d) in zip(r,Vr,Vpr_new,Vppr_new)
    @printf("%15.5f %15.5f %30.20f %30.20f\n",a,b,c,d)   
end


println("A small detail in V''(r)")
@printf("%15s %15s %30s\n","r","V'(r)","V''(r)")
for r in 4.268:0.0001:4.269    
    @printf("%15.5f %30.20f %30.20f\n",r,V'(r),V''(r))
end
println("... the potential has a poor 2nd derivative,")
println("caused by the if-clause. There's a discontinuity at r=1.438*2.9683 = 4.2684154.")


println("plotting...")
using Plots
pyplot() # pyplot backend
# more points for plotting
r  = [0.01:0.0001:10.0;]
Vr  = V.(r)
Vpr = V'.(r)
Vppr = V''.(r)
plot(r,Vr,label="V(r)",ylims=(-20,20))
plot!(r,Vpr,label="V'(r)") # update previous plot
plot!(r,Vppr,label="V''(r)") # update previous plot


#r  = [4.268:0.00001:4.269;]    
#plot(r,V''.(r))