#0----------------------------USEFUL INFORMATION--------------------------------
#To run only part of the code, select the desired lines and press (ctrl + enter)
#To run the whole code, press (ctrl + shift + enter)

#There is no equivalent to MATLAB's 'clc', 'clear all', 'close all', but
clearconsole()
#Works more or less like the 'clc'
#You can also clear the console by pressing (ctrl + L) in the command window
#0------------------------------------------------------------------------------

#System of ODES - Lorenz Equation

#Packages necessary for this file
using Printf
using Statistics
using LinearAlgebra
using DifferentialEquations
using Plots
#This is the 'backend' for the plots package, see more in http://docs.juliaplots.org/latest/install/
plotlyjs()

#One can see the equations here (Example 2): https://diffeq.sciml.ai/stable/tutorials/ode_example/

#0------------------Use one or the other, of course-----------------------------
#Defining the Integrator Function - Directly from equations
function Lorenz!(du,u,p,t)
    sigma, rho, beta = p
    du[1] = sigma*(u[2] - u[1])
    du[2] = u[1]*(rho - u[3]) - u[2]
    du[3] = u[1]*u[2] - beta*u[3]
end

#Alternative way in defining the Integrator Function - Matricial form
function Lorenz!(du,u,p,t)
    sigma, rho, beta = p
    u[1],u[2],u[3] = u
    du[1],du[2],du[3] = du
    Amatrix = [-sigma sigma 0; rho -1 -u[1]; 0 u[1] -beta]
    du .= Amatrix*u
end
#0------------------------------------------------------------------------------

#Initial conditions
u0 = [1.0,0.0,0.0]
#Defining parameters
p = (10,28,8/3)
#Time span
tspan = (0.0,100.0)

#Here we define a problem, using the initial conditions, tspan and parameters vector
prob = ODEProblem(Lorenz!,u0,tspan,p)

#Solving the problem------------------------------------------------------------
#Simplest way -> Default
sol = solve(prob)

#Defining max iterations and tolerances
#abstol -> accuracy when u is near zero (good rule ->
#actual tolerance one to two digits smaller than reltol)
#Reltol -> accuracy is 'most cases'
sol = solve(prob,maxiters = 1e8,abstol = 1e-8,reltol = 1e-8)

#I can also ask for specific time instants:
sol = solve(prob,saveat = 0.1)

#IT can even be an uneven 'grid':
sol = solve(prob,saveat = [0.0, 0.2, 0.7, 0.9, 1.0])

#We can tell the algorithm that the problem is going to be stiff:
sol = solve(prob, alg_hints = [:stiff])

#Only use when required, takes longer
#We can also specify one in particular:
sol = solve(prob, Tsit5())
sol = solve(prob, RK4())
sol = solve(prob, BS3())
sol = solve(prob, BS5())
#More information and full list of solvers at -> https://diffeq.sciml.ai/stable/solvers/ode_solve/

#Acessing and visualizing the solution------------------------------------------

#Time vector
println(sol.t)
#Solutions vector
println(sol.u)
#Only first DOF
println(sol[1,:])
#Only Second DOF and so on...
println(sol[2,:])
#In the : you can asign a specific time interval through the vectors indexes
println(sol[1,1:2])

#The solution is always an interpolating function, so i can ask for values
#outside of the discrete time intervals, e.g.:
println(sol(1.0))


#Plotting the solution----------------------------------------------------------
#Simple plot
fig = plot(sol)

#Plotting only the first DOF
fig = plot(sol.t,sol[1,:])
#Only the second DOF and so on...
fig = plot(sol.t,sol[2,:])

#'0' is the independent variable - Time
fig = plot(sol,vars = (0,1))

#We can do a phase plot also using the vars command
#We are saying plot x versus y versus z

fig = plot(sol,vars = (1,2,3))