using DifferentialEquations: include

#############################################################
# Running Simulations
#############################################################

#Include KG_Module to run simulation
include("../../src/KG_Module.jl")

#Helper function
function Kepler_Frequency( Rbh , R)
    return sqrt(Rbh/R^3)
end

function radius_parameters( Rbh)
    println("BH radius  : " , Rbh)
    println("LightRing  : " , 3*Rbh/2)
    println("ISCO radius: " , 3*Rbh)
    println("Safe Orbit : " , 3*Rbh)
end

#####################
# SIMULATION
#####################

#Setup masses
p0 = Param(

    #Temporal variables
    tmax           = 60000.0     , 
    out_every_t    = 2.5        ,  
    deltat         = 0.05       , 

    #Spacial variables
    xmin           =  -30.0     ,  
    xmax           =  30.0      , 
    xnodes         =  256       ,    

    #Spacial variables
    ymin           =  -30.0     ,  
    ymax           =  30.0      , 
    ynodes         =   256      , 
    
    #Initial Conditions
    A0              =   3.5     , 
    σ               =   1.5     , 
    r0              =   17.0    ,
    ω               =   0.5     , 
    m               =   2       , 
    
    #Klein Gordon Mass
    μ               =   0.0     ,

    #Dissipation terms
    dissipation     =  :true    ,   
    
    R_orbit         = 3.34716475041       ,
    R1              = 1.5       ,
    R2              = 0.0,#1.5       ,
       
    Ω               = Kepler_Frequency( 1.5 , 3.34716475041) , # Kepler's law
    alpha1          = 10.0       ,
    alpha2          = 0.0, #10.0       ,    

    #Periodic BC?
    Boundaries      = :radial    , # :square :periodic 
    
    #Data variables
    folder          = "/2020_07_01/"  ,       #must have a bar before the folder
    fname           = "default_name"

)

#Binary BH with a certain orbital radius
p1 = Param(     p0,

                R1              = 1.5       ,
                R2              = 1.5       ,
                R_orbit         = 4.5       ,
                
                alpha1          = 10.0          ,
                alpha2          = 10.0           ,   

                Ω               = Kepler_Frequency( 1.5 , 4.5) , # Kepler's law

                fname = "bbh_R1_1.5_Alpha_10.0_Rorbit_4.5_Rout_27.0_tmax_50000"
            )

p2 = Param(     p0,

            R1              = 1.5       ,
            R2              = 0.0       ,
            R_orbit         = 0.0       ,
            
            alpha1          = 10.0          ,
            alpha2          = 0.0           ,   

            Ω               = Kepler_Frequency( 1.5 , 4.5) , # Kepler's law

            fname = "sbh_R1_1.5_Alpha_10.0_Omega_0.13_Rout_27.0_tmax_50000"
        )

#s = [p1, p2]

#for x in s 
#    solve_wave_equation_2D(x)
#end

#include("../../src/KG_Plot_Module.jl")

#slices_observer(p2)
#slices_observer("./data/2020_07_01_bhb_sbh_comparison/default_name4.h5")