module GR_Spherical

using DifferentialEquations
using BoundaryValueDiffEq
using OrdinaryDiffEq
using Fun1d
using DataFrames
using CSV
using Plots
using Roots

using BenchmarkTools
using InteractiveUtils
using RecursiveArrayTools
using StaticArrays
using LinearAlgebra

using Profile

struct Param{T}
    rtmin::T
    rtmax::T
    drt::T
    Mtot::T
    grid::Grid{T}
    r::Function
    drdrt::Function
    d2rdrt::Function
    rsamp::Vector{T}
    drdrtsamp::Vector{T}
    d2rdrtsamp::Vector{T}
    state::ArrayPartition{T, NTuple{12, Vector{T}}}
    drstate::ArrayPartition{T, NTuple{12, Vector{T}}}
    dr2state::ArrayPartition{T, NTuple{12, Vector{T}}}
    dtstate::ArrayPartition{T, NTuple{12, Vector{T}}}
    dissipation::ArrayPartition{T, NTuple{12, Vector{T}}}
    temp::ArrayPartition{T, NTuple{12, Vector{T}}}
end

numvar = 12

@inline function Base.similar(::Type{ArrayPartition},::Type{T},size::Int) where T
    return ArrayPartition(
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),
        similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size)),similar(Vector{T}(undef,size))
    )
end

function printlogo()
    # Just a fancy ASCII logo for the program
    println(
"\n",
"           /\\\\\\\\\\\\\\\\\\       /\\\\\\\\\\\\\\\\\\\\\\\\\n",
" ________/\\\\\\\\\\\\\\\\\\\\\\\\\\\\___/ /\\\\\\\\/ / /\\\\\\\\_____________________________\n",
"  ______//\\\\\\\\\\/ / / / /___\\/ /\\\\\\\\/_/ / /\\\\\\\\___________________________\n",
"   _____//\\\\\\\\/_/_/_/_/_____\\/ /\\\\\\\\ \\/ / /\\\\\\\\___________________________\n",
"    ____/ /\\\\\\\\____/\\\\\\\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\_______in___________________\n",
"     ___\\/ /\\\\\\\\__/ /  /\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\/________Spherical____________\n",
"      ___\\/ /\\\\\\\\_\\/__/ /\\\\\\___\\/ /\\\\\\\\ / /\\\\\\\\__________Symmetry____________\n",
"       ___\\/ /\\\\\\\\____\\/ /\\\\\\___\\/ /\\\\\\\\_/ //\\\\\\\\_____________________________\n",
"        ___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\___\\/ /\\\\\\\\\\// //\\\\\\\\____________________________\n",
"         ___\\/ / /\\\\\\\\\\\\\\\\\\  /____\\/ /\\\\\\\\_\\// //\\\\\\\\_____by Conner Dailey_____\n",
"          ___\\/_/ / / / / /_/______\\/ /  /___\\// /  /______________________________\n",
"           _____\\/_/_/_/_/__________\\/__/______\\/__/______________________________\n"
)

end

function sample!(f::Vector{T}, grid::Grid{S}, fun) where {S,T}

    f .= T[fun(location(grid, n)) for n in 1:(grid.ncells + 4)]

end

function init!(state::ArrayPartition, param) where T

    ############################################
    # Specifies the Initial Conditions
    ############################################

    αreg,A,βr,Br,χ,γtrrreg,γtθθreg,Arrreg,Kreg,Γr,𝜙,K𝜙 = state.x

    grid = param.grid
    drt = spacing(grid)
    r = param.r
    drdrt = param.drdrt
    d2rdrt = param.d2rdrt
    rtmin = param.rtmin
    rtmax = param.rtmax

    n = grid.ncells + 4
    m = 1.
    rtspan = (rtmin,rtmax)

    # Initial conditions for Schwarzschild metric (Kerr-Schild Coordinates)

    # Mass (no real reason not to use 1 here)
    #M = 1

    fα(M,rt) = real((1+2*M/(r(rt))+0im)^(-1/2))
    fA(rt) = 0.
    fβr(M,rt) = (2*M/r(rt))*fα(M,rt)^2
    fBr(rt) = 0.
    fχ(rt) = 1.
    fγtrr(M,rt) = 1+2*M/r(rt)
    fγtθθ(rt) = r(rt)^2
    fArr(M,∂M,rt) = (4/3)*(r(rt)*(M+r(rt))*∂M-M*(3*M+2*r(rt)))/real(((r(rt)^5)*(r(rt)+2*M)+0im)^(1/2))
    fK(M,∂M,rt) = (2*M*(3*M+r(rt))+2*r(rt)*∂M*(M+r(rt)))/real((r(rt)*(r(rt)+2*M)+0im)^(3/2))
    fΓr(M,∂M,rt) = (r(rt)*∂M-2*r(rt)-5*M)/(r(rt)+2*M)^2

    r0 = 10.
    σr = 0.5
    #Amp = 1.
    Amp = 0.

    f𝜙(rt) = Amp*(1/r(rt))*exp(-(1/2)*((r(rt)-r0)/σr)^2)
    f∂𝜙(rt) = Amp*exp(-(1/2)*((r(rt)-r0)/σr)^2)*(r(rt)*r0-r(rt)^2-σr^2)/(r(rt)^2*σr^2)
    fK𝜙(rt) = 0.

    fρ(M,rt) = (2*fK𝜙(rt)^2 + (1/2)*(fχ(rt)/fγtrr(M,rt))*f∂𝜙(rt)^2
        + (1/2)*m^2*f𝜙(rt)^2)

    function f∂rtM(M,rt)
         if r(rt) < 2*M
             return 0.
         else
             4*pi*(r(rt)^2)*fρ(M,rt)*drdrt(rt)
         end
    end

    function f∂M(M,rt)
         if r(rt) < 2*M
             return 0.
         else
             4*pi*(r(rt)^2)*fρ(M,rt)
         end
    end

    fαreg(M,rt) = 2*M
    fγtrrreg(M,rt) = 2*M
    fArrreg(M,∂M,rt) = sqrt(r(rt)^5)*fArr(M,∂M,rt)
    fKreg(M,∂M,rt) = sqrt(r(rt)^3)*fK(M,∂M,rt)

    # Constraint Equations

    function constraintSystem(M, param, rt)
        f∂rtM(M,rt)
    end

    atol = 1e-15

    BVP = ODEProblem(constraintSystem, 1., rtspan, param)
    M = solve(BVP, Tsit5(), abstol=atol, dt=drt, adaptive=false)

    ∂M(rt) = f∂M(M(rt),rt)

    sample!(αreg, grid, rt -> fαreg(M(rt),rt) )
    sample!(A, grid, fA)
    sample!(βr, grid, rt -> fβr(M(rt),rt) )
    sample!(Br, grid, fBr)
    sample!(χ, grid, fχ)
    sample!(γtrrreg, grid, rt -> fγtrrreg(M(rt),rt) )
    sample!(γtθθreg, grid, rt -> 0)
    sample!(Arrreg, grid, rt -> fArrreg(M(rt),∂M(rt),rt) )
    sample!(Kreg, grid, rt -> fKreg(M(rt),∂M(rt),rt) )
    sample!(Γr, grid, rt -> fΓr(M(rt),∂M(rt),rt))
    sample!(𝜙, grid, f𝜙)
    sample!(K𝜙, grid, fK𝜙)

end

@inline function deriv!(df::Vector{T}, f::Vector{T}, n::Int64, dx::T) where T

    # @inbounds @fastmath @simd

    df[1] = T((-25. *f[1] + 48. *f[2] - 36. *f[3] + 16. *f[4] - 3. *f[5])/(12. *dx))

    df[2] = T((-3. *f[1] - 10. *f[2] + 18. *f[3] - 6. *f[4] + f[5])/(12. *dx))

    for i in 3:(n - 2)
        df[i] = T(((f[i-2] - 8. *f[i-1] + 8. *f[i+1] - f[i+2])/(12. *dx)))
    end

    df[n-1] = T(-(-3. *f[n] - 10. *f[n-1] + 18. *f[n-2] - 6. *f[n-3] + f[n-4])/(12. *dx))

    df[n] = T(-(-25. *f[n] + 48. *f[n-1] - 36. *f[n-2] + 16. *f[n-3] - 3. *f[n-4])/(12. *dx))

end

@inline function dissipation!(df::Vector{T}, f::Vector{T},drdrt::Vector{T}, n::Int64) where T

    ############################################
    # Calculates the numerical dissipation terms
    #
    # Finite differencing fails to model high
    # frequency modes in the system, and these
    # modes can lead to instabilities. Adding
    # numerical dissipation terms damps these
    # high frequency modes by subtracting off
    # a high order derivative from each dynamical
    # variable in the system.
    ############################################

    # @simd @inbounds @fastmath

    df[1:2] .= 0.

    for i in 3:(n - 2)
        df[i] = T((f[i-2] - 4. *f[i-1] + 6. *f[i] - 4. *f[i+1] + f[i+2])/(drdrt[i]))
    end

    df[(n-1):n] .= 0.


end

function rhs_test(dtstate::ArrayPartition{T, NTuple{12, Vector{T}}},regstate::ArrayPartition{T, NTuple{12, Vector{T}}}, param::Param{T}, t) where T

    ############################################
    # Caculates the right hand ride of the
    # evolved variables
    #
    # This is the main meat of the program.
    # This function contains all of the boundary
    # conditions, coordinate conversions,
    # spatial derivative calculations,
    # evolution equations, and numerical
    # dissipation. Each time the Julia DiffEq
    # Solver moves one time step, it calls
    # this function to calculate the new
    # values of the evolved variables.
    ############################################

    # Unpack the parameters

    m = 1.
    Mtot = 1.

    grid = param.grid
    drt = param.drt
    r = param.rsamp
    drdrt = param.drdrtsamp
    d2rdrt = param.d2rdrtsamp
    rtmin = param.rtmin
    rtmax = param.rtmax

    state = param.state
    drstate = param.drstate
    dr2state = param.dr2state
    dissipation = param.dissipation
    dtstate2 = param.dtstate
    temp = param.temp

    n = grid.ncells + 4

    # Copy the state into the parameters
    # so that it can be changed

    #######################
    # Attention!
    #
    # Do not do the following:
    # state .= regstate
    #
    # This results in an intense slowdown
    # Do instead:
    for i in 1:12
        state.x[i] .= regstate.x[i]
    end

    # Give names to individual variables

    α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state.x
    ∂α,∂A,∂βr,∂Br,∂χ,∂γtrr,∂γtθθ,∂Arr,∂K,∂Γr,∂𝜙,∂K𝜙 = drstate.x
    ∂2α,∂2A,∂2βr,∂2Br,∂2χ,∂2γtrr,∂2γtθθ,∂2Arr,∂2K,∂2Γr,∂2𝜙,∂2K𝜙 = dr2state.x
    ∂tα,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrr,∂tγtθθ,∂tArr,∂tK,∂tΓr,∂t𝜙,∂tK𝜙 = dtstate.x
    ∂4α,∂4A,∂4βr,∂4Br,∂4χ,∂4γtrr,∂4γtθθ,∂4Arr,∂4K,∂4Γr,∂4𝜙,∂4K𝜙 = dissipation.x

    # Calculate first derivatives

    deriv!(∂α,α,n,drt)
    deriv!(∂A,A,n,drt)
    deriv!(∂βr,βr,n,drt)
    deriv!(∂Br,Br,n,drt)
    deriv!(∂χ,χ,n,drt)
    deriv!(∂γtrr,γtrr,n,drt)
    deriv!(∂γtθθ,γtθθ,n,drt)
    deriv!(∂Arr,Arr,n,drt)
    deriv!(∂K,K,n,drt)
    deriv!(∂Γr,Γr,n,drt)
    deriv!(∂𝜙,𝜙,n,drt)
    deriv!(∂K𝜙,K𝜙,n,drt)

    # Calculate second derivatives

    deriv!(∂2α,∂α,n,drt)
    deriv!(∂2βr,∂βr,n,drt)
    deriv!(∂2χ,∂χ,n,drt)
    deriv!(∂2γtrr,∂γtrr,n,drt)
    deriv!(∂2γtθθ,∂γtθθ,n,drt)
    deriv!(∂2𝜙,∂𝜙,n,drt)

    # Convert between computational rt coordinnate
    # and trasditional r coordinate

    ∂α ./= drdrt
    ∂A ./= drdrt
    ∂βr ./= drdrt
    ∂Br ./= drdrt
    ∂χ ./= drdrt
    ∂γtrr ./= drdrt
    ∂γtθθ ./= drdrt
    ∂Arr ./= drdrt
    ∂K ./= drdrt
    ∂Γr ./= drdrt
    ∂𝜙 ./= drdrt
    ∂K𝜙 ./= drdrt

    @. ∂2α = (∂2α - d2rdrt*∂α)/(drdrt^2)
    @. ∂2βr = (∂2βr - d2rdrt*∂βr)/(drdrt^2)
    @. ∂2χ = (∂2χ - d2rdrt*∂χ)/(drdrt^2)
    @. ∂2γtrr = (∂2γtrr - d2rdrt*∂γtrr)/(drdrt^2)
    @. ∂2γtθθ = (∂2γtθθ - d2rdrt*∂γtθθ)/(drdrt^2)
    @. ∂2𝜙 = (∂2𝜙 - d2rdrt*∂𝜙)/(drdrt^2)

    # Convert between regularized variables
    # and cannonical variables

    @. α = real((1 + α/r + 0im)^(-1/2))
    @. ∂α = α*(1 - (1 + ∂α)*α^2)/(2*r)
    @. ∂2α = (6*r*∂α^2 - 4*α*∂α - ∂2α*α^4)/(2*r*α)

    @. γtrr = γtrr/r + 1
    @. ∂γtrr = (1 - γtrr + ∂γtrr)/r
    @. ∂2γtrr = (∂2γtrr - 2*∂γtrr)/r

    @. γtθθ = (r^2)*(γtθθ + 1)
    @. ∂γtθθ = (2*γtθθ + ∂γtθθ*r^3)/r
    @. ∂2γtθθ = (4*∂γtθθ*r - 6*γtθθ + ∂2γtθθ*r^4)/(r^2)

    @. K = sqrt(r^(-3))*K
    @. ∂K = sqrt(r^(-3))*∂K - (3/2)*K/r

    @. Arr = sqrt(r^(-5))*Arr
    @. ∂Arr = sqrt(r^(-5))*∂Arr - (5/2)*Arr/r

    #########################################################
    # Evolution Equations
    #
    # This is the full suite of evolution equations
    # for GR in spherical symmetry in the BSSN framework.
    # I have tried to keep them looking as close to their
    # mathematically written counterpart as possible.
    #
    # They are written in the order they appear in the
    # reference (arXiv:0705.3845v2) except for the ∂tBr
    # equation since it contains a ∂tΓr term.
    #
    #########################################################

    # Lagrangian Gauge condition
    v = 1.

    @. ∂tχ = ((2/3)*K*α*χ - (1/3)*v*βr*χ*∂γtrr/γtrr - (2/3)*v*βr*χ*∂γtθθ/γtθθ
     - (2/3)*v*χ*∂βr + βr*∂χ)

    @. ∂tγtrr = (-2*Arr*α - (1/3)*v*βr*∂γtrr + βr*∂γtrr
     - (2/3)*v*γtrr*βr*∂γtθθ/γtθθ + 2*γtrr*∂βr - (2/3)*v*γtrr*∂βr)

    @. ∂tγtθθ = (Arr*γtθθ*α/γtrr - (1/3)*v*γtθθ*βr*∂γtrr/γtrr - (2/3)*v*βr*∂γtθθ
     + βr*∂γtθθ - (2/3)*v*γtθθ*∂βr)

    @. ∂tArr = (-2*α*(Arr^2)/γtrr + K*α*Arr - (1/3)*v*βr*Arr*∂γtrr/γtrr
     - (2/3)*v*βr*Arr*∂γtθθ/γtθθ - (2/3)*v*Arr*∂βr + 2*Arr*∂βr
     + (2/3)*α*χ*(∂γtrr/γtrr)^2 - (1/3)*α*χ*(∂γtθθ/γtθθ)^2
     - (1/6)*α*(∂χ^2)/χ - (2/3)*α*χ*γtrr/γtθθ + βr*∂Arr
     + (2/3)*α*χ*γtrr*∂Γr - (1/2)*α*χ*(∂γtrr/γtrr)*(∂γtθθ/γtθθ)
     + (1/3)*χ*∂γtrr*∂α/γtrr + (1/3)*χ*∂α*∂γtθθ/γtθθ - (1/6)*α*∂γtrr*∂χ/γtrr
     - (1/6)*α*∂γtθθ*∂χ/γtθθ - (2/3)*∂α*∂χ - (1/3)*α*χ*∂2γtrr/γtrr
     + (1/3)*α*χ*∂2γtθθ/γtθθ - (2/3)*χ*∂2α + (1/3)*α*∂2χ)

    @. ∂tK = ((3/2)*α*(Arr/γtrr)^2 + (1/3)*α*K^2 + βr*∂K
     + (1/2)*χ*∂γtrr*∂α/(γtrr^2) - χ*∂α*(∂γtθθ/γtθθ)/γtrr
     + (1/2)*∂α*∂χ/γtrr - χ*∂2α/γtrr)

    @. ∂tΓr = (-v*βr*((∂γtθθ/γtθθ)^2)/γtrr + α*Arr*(∂γtθθ/γtθθ)/(γtrr^2)
     - (1/3)*v*∂βr*(∂γtθθ/γtθθ)/γtrr + ∂βr*(∂γtθθ/γtθθ)/γtrr
     + βr*∂Γr + α*Arr*∂γtrr/(γtrr^3) - (4/3)*α*∂K/γtrr
     - 2*Arr*∂α/(γtrr^2) + (1/2)*v*∂βr*∂γtrr/(γtrr^2)
     - (1/2)*∂βr*∂γtrr/(γtrr^2) - 3*α*Arr*(∂χ/χ)/(γtrr^2)
     + (1/6)*v*βr*∂2γtrr/(γtrr^2) + (1/3)*v*βr*(∂2γtθθ/γtθθ)/γtrr
     + (1/3)*v*∂2βr/γtrr + ∂2βr/γtrr)

    #######################################################
    # Gauge Evolution

    @. ∂tα = -2*α*A
    @. ∂tA = ∂tK

    @. ∂tβr = (3/4)*Br
    @. ∂tBr = ∂tΓr

    # Gauge choices for the evolution of the
    # determinant of the conformal metric
    # (must have v = 1 to use this)

    # ∂tlnγt = temp.x[5]
    # ∂rt∂tlnγt = temp.x[6]
    #
    # ∂tlnγt .= 0
    #
    # #∂tlnγt = -8*pi*Sr.*real((γtθθ./γtrr .+ 0im).^(1/2))
    #
    # deriv!(∂rt∂tlnγt,∂tlnγt,n,drt)
    #
    # @. ∂r∂tlnγt = ∂rt∂tlnγt/drdrt
    #
    # # ∂tα = -(1/2)*α.*∂tlnγt
    # # ∂tβr = (χ./γtrr).*∂tlnγt
    #
    # @. ∂tχ += (1/3)*χ*∂tlnγt
    # @. ∂tγtrr += (1/3)*γtrr*∂tlnγt
    # @. ∂tγtθθ += (1/3)*γtθθ*∂tlnγt
    # @. ∂tArr += (1/3)*Arr*∂tlnγt
    # @. ∂tΓr += -(1/3)*Γr*∂tlnγt - (1/6)*(χ/γtrr)*∂r∂tlnγt

    #########################################################
    # Source Terms and Source Evolution
    #
    # This currently includes the addition of source terms
    # to GR that come from a Klein-Gordon scalar field
    #
    #########################################################

    # Klein-Gordon System

    @. ∂t𝜙 = βr*∂𝜙 - 2*α*K𝜙

    @. ∂tK𝜙 = (βr*∂K𝜙 + α*K*K𝜙 + (1/2)*m^2*α*𝜙 - (1/2)*χ.*∂α.*∂𝜙./γtrr
        - (1/2)*(α*χ/γtrr)*(∂2𝜙 + ∂𝜙*(∂γtθθ/γtθθ - (1/2)*∂γtrr/γtrr
        - (1/2)*(∂χ/χ))))

    ρ = temp.x[1]
    Sr = temp.x[2]
    S = temp.x[3]
    Srr = temp.x[4]

    @. ρ = 2*K𝜙^2 + (1/2)*(χ/γtrr)*∂𝜙^2 + (1/2)*(m^2)*𝜙^2
    #Lower Index
    @. Sr = 2*K𝜙*∂𝜙
    @. S = 6*K𝜙^2 - (1/2)*(χ/γtrr)*∂𝜙^2 - (3/2)*(m^2)*𝜙^2
    @. Srr = (γtrr/χ)*(2*K𝜙^2 + (1/2)*(χ/γtrr)*∂𝜙^2 - (1/2)*(m^2)*𝜙^2)

    @. ∂tArr += -8*pi*α*(χ*Srr - (1/3)*S*γtrr)
    @. ∂tK += 4*pi*α*(ρ + S)
    @. ∂tΓr += -16*pi*α*Sr/γtrr


    # fr = param.r
    #
    # fα(M,rt) = real((1+2*M/(fr(rt))+0im)^(-1/2))
    # fβr(M,rt) = (2*M/fr(rt))*fα(M,rt)^2
    # fγtrr(M,rt) = 1+2*M/fr(rt)
    # fγtθθ(rt) = fr(rt)^2
    # fArr(M,∂M,rt) = (4/3)*(fr(rt)*(M+fr(rt))*∂M-M*(3*M+2*fr(rt)))/real(((fr(rt)^5)*(fr(rt)+2*M)+0im)^(1/2))
    # fK(M,∂M,rt) = (2*M*(3*M+fr(rt))+2*fr(rt)*∂M*(M+fr(rt)))/real((fr(rt)*(fr(rt)+2*M)+0im)^(3/2))
    # fΓr(M,∂M,rt) = (fr(rt)*∂M-2*fr(rt)-5*M)/(fr(rt)+2*M)^2
    #
    # f∂α(M,rt) = M*real((fr(rt)*(fr(rt)+2*M+0im)^3)^(-1/2))
    # f∂βr(M,rt) = -2*M/(fr(rt)+2*M)^2
    # f∂γtrr(M,rt) = -2*M/(fr(rt)^2)
    # f∂γtθθ(rt) = 2*fr(rt)
    # f∂Arr(M,rt) = (4*M/3)*(15*M^2+15*M*fr(rt)+4*fr(rt)^2)/real(((fr(rt)^7)*((fr(rt)+2*M)^3)+0im)^(1/2))
    # f∂K(M,rt) = -2*M*(9*M^2+10*M*fr(rt)+2*fr(rt)^2)/real((fr(rt)*(fr(rt)+2*M)+0im)^(5/2))
    # f∂Γr(M,rt) = 2*(fr(rt)+3*M)/(fr(rt)+2*M)^3
    #
    # rt = sample(Float64, A.grid, rt->rt)
    #
    # # for i=1:2
    #
    # ∂tα[1:2] .= (α[1:2] .- fα.(1.,rt[1:2]))./r[1:2] + ∂α[1:2] - f∂α.(1.,rt[1:2])
    # ∂tA[1:2] .= (A[1:2] .- 0.)./r[1:2] + ∂A[1:2]
    # ∂tβr[1:2] .= (βr[1:2] .- fβr.(1.,rt[1:2]))./r[1:2] + ∂βr[1:2] - f∂βr.(1.,rt[1:2])
    # ∂tBr[1:2] .= (Br[1:2] .- 0.)./r[1:2] + ∂Br[1:2]
    # ∂tχ[1:2] .= (χ[1:2] .- 1.)./r[1:2] + ∂χ[1:2]
    # ∂tγtrr[1:2] .= (γtrr[1:2] .- fγtrr.(1.,rt[1:2]))./r[1:2] + ∂γtrr[1:2] - f∂γtrr.(1.,rt[1:2])
    # ∂tγtθθ[1:2] .= (γtθθ[1:2] .- fγtθθ.(rt[1:2]))./r[1:2] + ∂γtθθ[1:2] - f∂γtθθ.(rt[1:2])
    # ∂tArr[1:2] .= (Arr[1:2] .- fArr.(1.,0.,rt[1:2]))./r[1:2] + ∂Arr[1:2] - f∂Arr.(1.,rt[1:2])
    # ∂tK[1:2] .= (K[1:2] .- fK.(1.,0.,rt[1:2]))./r[1:2] + ∂K[1:2] - f∂K.(1.,rt[1:2])
    # ∂tΓr[1:2] .= (Γr[1:2] .- fΓr.(1.,0.,rt[1:2]))./r[1:2] + ∂Γr[1:2] - f∂Γr.(1.,rt[1:2])
    # ∂t𝜙[1:2] .= (𝜙[1:2] .- 0.)./r[1:2] + ∂𝜙[1:2]
    # ∂tK𝜙[1:2] .= (K𝜙[1:2] .- 0.)./r[1:2] + ∂K𝜙[1:2]

    ######################################################
    # Calculate the numerical dissipation

    dissipation!(∂4α,α,drdrt,n)
    dissipation!(∂4A,A,drdrt,n)
    dissipation!(∂4βr,βr,drdrt,n)
    dissipation!(∂4Br,Br,drdrt,n)
    dissipation!(∂4χ,χ,drdrt,n)
    dissipation!(∂4γtrr,γtrr,drdrt,n)
    dissipation!(∂4γtθθ,γtθθ,drdrt,n)
    dissipation!(∂4Arr,Arr,drdrt,n)
    dissipation!(∂4K,K,drdrt,n)
    dissipation!(∂4Γr,Γr,drdrt,n)
    dissipation!(∂4𝜙,𝜙,drdrt,n)
    dissipation!(∂4K𝜙,K𝜙,drdrt,n)

    # Magnitude of dissipation
    σ = 0.3

    @. ∂tα -= (1/16)*σ*∂4α/drt
    @. ∂tA -= (1/16)*σ*∂4A/drt
    @. ∂tβr -= (1/16)*σ*∂4βr/drt
    @. ∂tBr -= (1/16)*σ*∂4Br/drt
    @. ∂tχ -= (1/16)*σ*∂4χ/drt
    @. ∂tγtrr -= (1/16)*σ*∂4γtrr/drt
    @. ∂tγtθθ -= (1/16)*σ*∂4γtθθ/drt
    @. ∂tArr -= (1/16)*σ*∂4Arr/drt
    @. ∂tK -= (1/16)*σ*∂4K/drt
    @. ∂tΓr -= (1/16)*σ*∂4Γr/drt
    @. ∂t𝜙 -= (1/16)*σ*∂4𝜙/drt
    @. ∂tK𝜙 -= (1/16)*σ*∂4K𝜙/drt

    # Convert back to regularized variables
    # for the time derivatives

    @. ∂tα = -2*r*∂tα/(α^3)
    @. ∂tγtrr = r*∂tγtrr
    @. ∂tγtθθ = ∂tγtθθ/(r^2)
    @. ∂tArr = sqrt(r^5)*∂tArr
    @. ∂tK = sqrt(r^3)*∂tK

    # Specify the inner temporal boundary conditions

    for i in 1:12
        dtstate.x[i][1:2] .= 0.
    end

    # Specify the outer temporal boundary conditions

    for i in 1:12
        dtstate.x[i][(n-1):n] .= 0.
    end

    # Store the calculated state into the param
    # so that we can print it to the screen

    for i in 1:12
        dtstate2.x[i] .= dtstate.x[i]
    end

end

function rhs_all(regstate::ArrayPartition{Float64, NTuple{12, Vector{Float64}}}, param::Param{Float64}, t)

    n = param.grid.ncells + 4

    dtstate = similar(ArrayPartition,Float64,n)

    rhs_test(dtstate,regstate,param,t)

    return dtstate

end

function constraints(state::ArrayPartition,drstate::ArrayPartition,dr2state::ArrayPartition,param)

    ############################################
    # Caculates the constraints of the system
    #
    # These are outputed at every saved time step.
    # They should limit to zero as the spatial
    # resolution is increased. They are important to
    # monitor to make sure the physics of the system
    # is being properly modeled.
    ############################################

    # Unpack Variables

    α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state.x

    ∂α,∂A,∂βr,∂Br,∂χ,∂γtrr,∂γtθθ,∂Arr,∂K,∂Γr,∂𝜙,∂K𝜙 = drstate.x

    ∂2α,∂2A,∂2βr,∂2Br,∂2χ,∂2γtrr,∂2γtθθ,∂2Arr,∂2K,∂2Γr,∂2𝜙,∂2K𝜙 = dr2state.x

    m = 1.
    n = param.grid.ncells + 4
    drt = param.drt
    r = param.rsamp
    drdrt = param.drdrtsamp
    d2rdrt = param.d2rdrtsamp

    deriv!(∂χ,χ,n,drt)
    deriv!(∂γtrr,γtrr,n,drt)
    deriv!(∂γtθθ,γtθθ,n,drt)
    deriv!(∂Arr,Arr,n,drt)
    deriv!(∂K,K,n,drt)
    deriv!(∂Γr,Γr,n,drt)
    deriv!(∂𝜙,𝜙,n,drt)
    deriv!(∂K𝜙,K𝜙,n,drt)

    deriv!(∂2χ,∂χ,n,drt)
    deriv!(∂2γtrr,∂γtrr,n,drt)
    deriv!(∂2γtθθ,∂γtθθ,n,drt)
    deriv!(∂2𝜙,∂𝜙,n,drt)

    ∂χ ./= drdrt
    ∂γtrr ./= drdrt
    ∂γtθθ ./= drdrt
    ∂Arr ./= drdrt
    ∂K ./= drdrt
    ∂Γr ./= drdrt
    ∂𝜙 ./= drdrt
    ∂K𝜙 ./= drdrt

    @. ∂2χ = (∂2χ - d2rdrt*∂χ)/(drdrt^2)
    @. ∂2γtrr = (∂2γtrr - d2rdrt*∂γtrr)/(drdrt^2)
    @. ∂2γtθθ = (∂2γtθθ - d2rdrt*∂γtθθ)/(drdrt^2)
    @. ∂2𝜙 = (∂2𝜙 - d2rdrt*∂𝜙)/(drdrt^2)

    @. γtrr = γtrr/r + 1
    @. ∂γtrr = (1 - γtrr + ∂γtrr)/r
    @. ∂2γtrr = (∂2γtrr - 2*∂γtrr)/r

    @. γtθθ = (r^2)*(γtθθ + 1)
    @. ∂γtθθ = (2*γtθθ + ∂γtθθ*r^3)/r
    @. ∂2γtθθ = (4*∂γtθθ*r - 6*γtθθ + ∂2γtθθ*r^4)/(r^2)

    @. K = sqrt(r^(-3))*K
    @. ∂K = sqrt(r^(-3))*∂K - (3/2)*K/r

    @. Arr = sqrt(r^(-5))*Arr
    @. ∂Arr = sqrt(r^(-5))*∂Arr - (5/2)*Arr/r


    ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2
    Sr = 2*K𝜙.*∂𝜙

    # Constraint Equations

    𝓗 = (-(3/2)*(Arr./γtrr).^2 + (2/3)*K.^2 - (5/2)*((∂χ.^2)./χ)./γtrr
     + 2*∂2χ./γtrr + 2*χ./γtθθ - 2*χ.*(∂2γtθθ./γtθθ)./γtrr + 2*∂χ.*(∂γtθθ./γtθθ)./γtrr
     + χ.*(∂γtrr./(γtrr.^2)).*(∂γtθθ./γtθθ) - ∂χ.*∂γtrr./(γtrr.^2)
     + (1/2)*χ.*((∂γtθθ./γtθθ).^2)./γtrr - 16*pi*ρ)

    𝓜r = (∂Arr./γtrr - (2/3)*∂K - (3/2)*Arr.*(∂χ./χ)./γtrr
     + (3/2)*Arr.*(∂γtθθ./γtθθ)./γtrr - Arr.*∂γtrr./(γtrr.^2)
     - 8*pi*Sr)

    𝓖r = -(1/2)*∂γtrr./(γtrr.^2) + Γr + (∂γtθθ./γtθθ)./γtrr

    return (𝓗, 𝓜r, 𝓖r, ρ)

end

function horizon(T,state::ArrayPartition,param)

    ############################################
    # Caculates the apparent horizon
    #
    # Where the function crosses zero is the
    # apparent horizon of the black hole.
    ############################################

    #v = param[3]

    # Unpack Variables

    χ = state.χ
    γtrr = state.γtrr
    γtθθ = state.γtθθ
    Arr = state.Arr
    K = state.K

    # Gauge condition

    # if v == 1 # Lagrangian Condition
    #     γtθθreg[2] = ((-315*γtθθreg[3] + 210*γtθθreg[4] - 126*γtθθreg[5]
    #     + 45*γtθθreg[6] - 7*γtθθreg[7])/63)
    # end

    #r = sample(T, χ.grid, param[4])
    drdrt = param.drdrtsamp

    # Conversions from regularized variables to canonical variables

    # γtrr = γtrrreg./r .+ 1
    #
    # γtθθ = (r.^2).*(γtθθreg .+ 1)
    #
    # K = real((r .+ 0im).^(-3/2)).*Kreg
    #
    # Arr = real((r .+ 0im).^(-5/2)).*Arrreg

    # Intermediate calculations

    Kθθ = ((1/3)*γtθθ.*K - (1/2)*Arr.*γtθθ./γtrr)./χ

    grr =  γtrr./χ

    gθθ =  γtθθ./χ

    # Spatial Derivatives

    ∂rtgθθ = deriv(gθθ,4,1)

    # Coordinate transformations from computational rt coordinate
    # to physical r coordinate

    ∂gθθ = ∂rtgθθ./drdrt

    # Apparent horizon function

    Θ = (∂gθθ./gθθ)./real((grr .+ 0im).^(1/2)) - 2*Kθθ./gθθ

    # cross = GridFun(χ.grid, sign.(Θ))

    return Θ

end

# function crossings(Θ::GridFun)
#
#     GridFun(χ.grid, sign.(Θ))
#
#     for
#
# end


function custom_progress_message(dt,state::ArrayPartition{T, NTuple{numvar, Vector{T}}},param,t) where T

    ###############################################
    # Outputs status numbers while the program runs
    ###############################################

    dtstate = param.dtstate::ArrayPartition{T, NTuple{numvar, Vector{T}}}

    ∂tα,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrr,∂tγtθθ,∂tArr,∂tK,∂tΓr,∂t𝜙,∂tK𝜙 = dtstate.x

    println("  ",
    #rpad(string(param[1]),6," "),
    #rpad(string(round(dt,digits=3)),10," "),
    rpad(string(round(t,digits=3)),10," "),
    rpad(string(round(maximum(abs.(∂tα)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tχ)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tγtrr)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(∂tγtθθ)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(∂tArr)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.𝜙)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.K𝜙)),digits=3)),14," ")
    rpad(string(round(maximum(abs.(∂tK)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(∂tΓr)),digits=3)),14," ")
    )

    #PrettyTables.jl

    # param[1] += param[2]

    return

end


function solution_saver(T,grid,sol,param,folder)

    ###############################################
    # Saves all of the variables in nice CSV files
    # in the choosen data folder directory
    ###############################################

    vars = (["α","A","βr","Br","χ","γtrr","γtθθ","Arr","K","Γr","𝜙","K𝜙",
    "∂tα","∂tA","∂tβr","∂tBr","∂tχ","∂tγtrr","∂tγtθθ","∂tArr","∂tK","∂tΓr","∂t𝜙","∂tK𝜙",
    "H","Mr","Gr","ρ","appHorizon"])
    varlen = length(vars)
    #mkdir(string("data\\",folder))
    tlen = size(sol)[2]
    rlen = grid.ncells + 4
    r = param.rsamp
    #cons = Array{GridFun,2}(undef,tlen,4)
    #state = Array{ArrayPartition,1}(undef,tlen)
    #dstate = Array{ArrayPartition,1}(undef,tlen)
    #d2state = Array{ArrayPartition,1}(undef,tlen)
    #derivs = Array{ArrayPartition,1}(undef,tlen)
    #apphorizon = Array{GridFun,1}(undef,tlen)

    # for i in 1:tlen
    #
    #     αreg,A,βr,Br,χ,γtrrreg,γtθθreg,Arrreg,Kreg,Γr,𝜙,K𝜙 = sol[i].x
    #
    #     # Conversions from regularized variables to canonical variables
    #     α = real((1 .+ αreg./r .+ 0im).^(-1/2))
    #     γtrr = γtrrreg./r .+ 1
    #     γtθθ = (r.^2).*(γtθθreg .+ 1)
    #     K = sqrt.(r.^(-3)).*Kreg
    #     Arr = sqrt.(r.^(-5)).*Arrreg
    #
    #     state[i] = ArrayPartition(α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙)
    #     derivs[i] = param.dtstate
    #     rhs_test(derivs[i],sol[i],param,0)
    #     #cons[i,1:4] .= constraints!(T,state[i],dstate[i],d2state[i],param)
    #     #apphorizon[i] = horizon(T,state[i],param)
    # end

    drstate = param.drstate
    dr2state = param.dr2state

    dtstate = [rhs_all(sol[i],param,0.) for i = 1:tlen]

    cons = [constraints(sol[i],drstate,dr2state,param) for i = 1:tlen]

    array = Array{T,2}(undef,tlen+1,rlen+1)

    array[1,1] = 0
    array[1,2:end] .= r

    for j = 1:varlen
        if j < 13
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= sol[i-1].x[j]
                # if j==1
                #     println(array[i,2:10])
                # end
            end
        elseif j < 25

            for i = 2:tlen+1

                array[i,1] = sol.t[i-1]
                array[i,2:end] .= dtstate[i-1].x[j-12]
            end
        elseif j < 29
            for i = 2:tlen+1
                #println(size(cons[i-1][j-24]))
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= cons[i-1][j-24]
            end
        # elseif j == varlen
        #     for i = 2:tlen+1
        #         array[i,1] = sol.t[i-1]
        #         array[i,2:end] .= apphorizon[i-1]
        #     end
        end

        CSV.write(
        string("data/",folder,"/",vars[j],".csv"),
        DataFrame(array, :auto),
        header=false
        )

    end

    # for i = 2:tlen+1
    #     array[i,1] = sol.t[i-1]
    #     array[i,2:end] .= sol[i-1].x[1][:]
    #     println(sol[i-1].x[1][1:10])
    # end
    #
    # CSV.write(
    # string("data/",folder,"/",vars[1],".csv"),
    # DataFrame(array, :auto),
    # header=false
    # )

end

function main(points)

    ###############################################
    # Main Program
    #
    # Calls each of the above functions to run a
    # simulation. Sets up the numerical grid,
    # sets the gauge conditions, sets up a new
    # computational rt coordinate that makes the
    # physical r coordinate step larger near the
    # outer boundary, sets the initial conditions,
    # and finally runs the numerical DiffEq
    # package to run the time integration.
    #
    # All data is saved in the folder specified to
    # the solution_saver, each in their own CSV
    # file.
    ###############################################

    T = Float64
    rspan = T[1.,210.]
    rtspan = T[1.,21.]

    rtmin, rtmax = rtspan

    domain = Domain{S}(rtmin, rtmax)
    grid = Grid(domain, points)

    n = grid.ncells + 4

    drt = spacing(grid)
    dt = drt/4.

    tspan = T[0.,20. ]
    tmin, tmax = tspan

    printtimes = 1.

    v = 1.

    m = 0.

    Mtot = 1.

    f(b) = b*tan(rtspan[2]/b)-rspan[2]

    scale = find_zero(f, 0.64*rtspan[2])

    r(rt) = scale*tan(rt/scale)
    drdrt(rt) = sec(rt/scale)^2
    d2rdrt(rt) = (2/scale)*(sec(rt/scale)^2)*tan(rt/scale)

    # r(rt) = rt
    # drdrt(rt) = 1
    # d2rdrt(rt) = 0

    atol = eps(T)^(T(3) / 4)

    #alg = KuttaPRK2p5()
    alg = RK4()
    #alg = ImplicitMidpoint()

    #printlogo()

    custom_progress_step = round(Int, printtimes/dt)
    step_iterator = custom_progress_step

    regstate = similar(ArrayPartition,T,n)
    state = similar(ArrayPartition,T,n)
    drstate = similar(ArrayPartition,T,n)
    dr2state = similar(ArrayPartition,T,n)
    dtstate = similar(ArrayPartition,T,n)
    dissipation = similar(ArrayPartition,T,n)
    temp = similar(ArrayPartition,T,n)

    # regstate = similar(Array{T}(undef,numvar,n))
    # state = similar(Array{T}(undef,numvar,n))
    # drstate = similar(Array{T}(undef,numvar,n))
    # dr2state = similar(Array{T}(undef,numvar,n))
    # dtstate = similar(Array{T}(undef,numvar,n))
    # dissipation = similar(Array{T}(undef,numvar,n))
    # temp = similar(Array{T}(undef,numvar,n))

    #println("Defining Problem...")
    rsamp = similar(Vector{T}(undef,n))
    drdrtsamp = similar(Vector{T}(undef,n))
    d2rdrtsamp = similar(Vector{T}(undef,n))

    sample!(rsamp, grid, r)
    sample!(drdrtsamp, grid, drdrt)
    sample!(d2rdrtsamp, grid, d2rdrt)

    param = Param(rtmin,rtmax,drt,Mtot,grid,r,drdrt,d2rdrt,rsamp,drdrtsamp,d2rdrtsamp,state,drstate,dr2state,dtstate,dissipation,temp)
    #
    init!(regstate, param)

    # @btime rhs!($regstate,$regstate,$param, 0.0)
    # # @btime radial_derivatives!($regstate,$regstate,$regstate, $param, 0.0)
    # # @btime regular_to_canonical!($regstate,$regstate,$regstate, $param, 0.0)
    # # # #
    # return

    prob = ODEProblem(rhs_test, regstate, tspan, param)

    println("Starting Solution...")

    #return

    println("")
    println("| Time | max α'(t) | max χ'(t) | max γtrr'(t) | max γtθθ'(t) | max Arr'(t) | max K'(t) | max Γr'(t) |")
    println("|______|___________|___________|______________|______________|_____________|___________|____________|")
    println("")


    sol = solve(
        prob, alg,
        abstol = atol,
        dt = drt/4,
        adaptive = false,
        saveat = printtimes,
        alias_u0 = true,
        progress = true,
        progress_steps = custom_progress_step,
        progress_message = custom_progress_message
    )

    #@code_native deriv!([1.0,1.0,1.0,1.0,1.0,1.0], [1.0,1.0,1.0,1.0,1.0,1.0], 6, 1.)

    #@btime α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = $state

    #@btime deriv!($dtstate[1],$regstate[1], $n, 1.)

    # rhs_test(dtstate,regstate, param, 0)
    #
    # Profile.clear_malloc_data()
    #
    # rhs_test(dtstate,regstate, param, 0)

    #@btime rhs_test($dtstate,$regstate, $param, 1.)

    #@code_warntype rhs_test(dtstate,regstate, param, 1.)

    # @profile rhs_test(dtstate,regstate, param, 1.)
    #
    # Profile.print()

    #@btime α,A,βr,Br,χ,γtrr,γtθθ,Arr,K,Γr,𝜙,K𝜙 = state

    #@btime radial_derivatives!($drstate, $dr2state, $regstate, $param, 0)

    # df = similar(Vector{T}(undef,n))
    # f = similar(Vector{T}(undef,n))
    #
    # a = GridFun(grid,Vector{T}(zeros(n)))
    # b = GridFun(grid,Vector{T}(zeros(n)))
    #
    # test = similar(VectorOfArray,T,n)
    # )



    #return

    solution_saver(T,grid,sol,param,"test")


    return

end


end

# k1 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k2 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k3 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
# k4 = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
#
# temp = zero(GBSSN_Variables{T,T}, grid)::GBSSN_Variables
#
# saveat = dt
#
# trange = tmin:dt:tmax
#
# tlen = length(trange)
#
# save_range = tmin:saveat:tmax
#
# save_len = length(save_range)
#
# save = Array{GBSSN_Variables,1}(undef,save_len)
#
# save[1] = regstate
#
# for t = trange
#
#     rhs!(dtstate, regstate, param, t)
#
#     k1 .= dtstate
#
#     temp .= regstate.+dt.*k1./2
#
#     rhs!(dtstate, temp, param , t + dt/2)
#
#     k2 .= dtstate
#
#     temp .= regstate.+dt.*k2./2
#
#     rhs!(dtstate, temp, param , t + dt/2)
#
#     k3 .= dtstate
#
#     temp .= regstate.+dt.*k3
#
#     rhs!(dtstate, temp, param , t + dt)
#
#     k4 .= dtstate
#
#     @. regstate = regstate + dt*(k1+2*k2+2 *k3+k4)/6
#
#     #save
#
#     if t in save_range
#
#         i = Int(round((tmax-tmin)/saveat) + 1)
#
#         save[i] = regstate
#
#     end
#
# end
#
# function M_init(::Type{T}, grid::Grid, param) where {T}
#
#     ############################################
#     # Specifies the Initial Conditions
#     ############################################
#
#     n = grid.ncells + 4
#     domain = grid.domain
#     initgrid = grid
#     drt = spacing(grid)
#     r = param[4]
#     drdrt = param[5]
#     d2rdrt = param[6]
#     m = param[7]
#     rtspan = param[8]
#
#     num = 0
#
#     fρ(M,rt) = (2*fK𝜙(rt)^2 + (1/2)*(fχ(rt)/fγtrr(M,rt))*f∂𝜙(rt)^2
#         + (1/2)*m^2*f𝜙(rt)^2)
#
#     fSr(rt) = 2*fK𝜙(rt)*f∂𝜙(rt)
#
#     #f∂M(M,rt) = 4*pi*(r(rt)^2)*fρ(M,rt)
#
#     function f∂M(M,rt)
#          if rt < 2
#              return 0.
#          else
#              4*pi*(r(rt)^2)*fρ(M,rt)
#          end
#     end
#
#     function f𝓗(M,∂M,rt)
#          (-(3/2)*(fArr(M,∂M,rt)/fγtrr(M,rt))^2 + (2/3)*fK(M,∂M,rt)^2
#          - (5/2)*((f∂χ(rt)^2)/fχ(rt))/fγtrr(M,rt) + 2*f∂2χ(rt)/fγtrr(M,rt)
#          + 2*fχ(rt)/fγtθθ(rt) - 2*fχ(rt)*(f∂2γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
#          + 2*f∂χ(rt)*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
#          + fχ(rt)*(f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2))*(f∂γtθθ(rt)/fγtθθ(rt))
#          - f∂χ(rt)*f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2)
#          + (1/2)*fχ(rt)*((f∂γtθθ(rt)/fγtθθ(rt))^2)/fγtrr(M,rt) - 16*pi*fρ(M,rt))
#     end
#
#     fαreg(M,rt) = 2*M
#     fγtrrreg(M,rt) = 2*M
#     fArrreg(M,∂M,rt) = real((r(rt)+ 0im)^(5/2))*fArr(M,∂M,rt)
#     fKreg(M,∂M,rt) = real((r(rt)+ 0im)^(3/2))*fK(M,∂M,rt)
#
#     # Constraint Equations
#
#     rtspan = (rtspan[1], rtspan[2])
#     #rtspan = (rtspan[2], 0.5)
#
#     function constraintSystem(M, param, rt)
#         f∂M(M,rt)
#     end
#
#     # function boundaryCondition!(residual, M, param, rt)
#     #     residual = M[end] - 1. #inner boundary condition
#     # end
#
#     atol = 1e-15
#
#     BVP = ODEProblem(constraintSystem, 1., rtspan, param)
#     M = solve(BVP, Tsit5(), abstol=atol, dt=drt, adaptive=false)
#
#     ∂M(rt) = f∂M(M(rt),rt)
#
#     # M(rt) = 1.
#     # ∂M(rt) = 0
#
#
#     state = GridFun(grid, M)
#
#     return
#
# end
#
# function M_rhs(M::GridFun, param, t)
#
#     #########################################################
#     # Source Terms and Source Evolution
#     #
#     # This currently includes the addition of source terms
#     # to GR that come from a Klein-Gordon scalar field
#     #
#     #########################################################
#
#     # Klein-Gordon System
#
#     ∂t𝜙 = βr.*∂𝜙 - 2*α.*K𝜙
#     ∂tK𝜙 = (βr.*∂K𝜙 + α.*K.*K𝜙 - (1/2)*α.*χ.*∂2𝜙./γtrr
#         + (1/4)*α.*χ.*∂γtrr.*∂𝜙./γtrr.^2 - (1/4)*α.*∂χ.*∂𝜙./γtrr
#         - (1/2)*χ.*∂α.*∂𝜙./γtrr - (1/2)*χ.*∂γtθθ.*∂𝜙./(γtrr.*γtθθ)
#         + (1/2)*∂χ.*∂𝜙./(γtrr) + (1/2)*m^2*𝜙)
#
#     ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2
#     #Lower Index
#     Sr = 2*γtrr.*K𝜙.*∂𝜙./χ
#     # S = 6*K𝜙.^2 - (1/2)*(χ./γtrr).*∂𝜙.^2 - (3/2)*m^2*𝜙.^2
#     # Srr = (γtrr./χ).*(2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 - (1/2)*m^2*𝜙.^2)
#
#     # ∂tArr .+= -8*pi*α.*(χ.*Srr - (1/3)*S.*γtrr)
#     # ∂tK .+= 4*pi*α.*(ρ + S)
#     # ∂tΓr .+= -16*pi*α.*Sr./γtrr
#
#     # Inner temporal boundary Conditions
#
#     # ∂tα[1:2] .= 0.
#     # ∂tA[1:2] .= 0.
#     # ∂tβr[1:2] .= 0.
#     # ∂tBr[1:2] .= 0.
#     # ∂tχ[1:2] .= 0.
#     # ∂tγtrr[1:2] .= 0.
#     # ∂tγtθθ[1:2] .= 0.
#     # ∂tArr[1:2] .= 0.
#     # ∂tK[1:2] .= 0.
#     # ∂tΓr[1:2] .= 0.
#     # ∂t𝜙[1:2] .= 0.
#     # ∂tK𝜙[1:2] .= 0.
#
#     return GBSSN_Variables(∂tαreg,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrrreg,∂tγtθθreg,∂tArrreg,∂tKreg,∂tΓr,∂t𝜙,∂tK𝜙)
#
# end
