using JuMP
using Ipopt
using DifferentialEquations
using Plots

# ------Define the model parameters----------

# Define Constants
α = 1700.0               #Heat transfer coefficient between reactor and jacket                        [kJ/Khm^2]
r = 0.092              #Radius of cross section inner part                                            [m]
ρ = 1000.0               #Density of the reactor contents                                             [g/L]
Cp = 4.2*1e-3          #Specific heat capacity of the reactor contents                                [kJ/gK]
Cbin = 3.0               #Input concentration of reactant B                                           [mol/L]
Vj = 2.22              #Volume of the contents inside the cooling jacket                              [L]
Tin = 300              #Temperature of the flows entering the reactor                                 [K]

# Define initial conditions
Vr0 = 3.5                                                                                              #[L]
Ca0 = 2.0                                                                                              #[mol/L]
Cb0 = 0.0                                                                                              #[mol/L]
Cc0 = 0.0                                                                                              #[mol/L]
Tr0 = 325.0                                                                                            #[K]
Tj0 = 325.0  
Aw0 = (π*(r^2))+((0.002*Vr0)/(r))                                                                      #[m^2]                
x0 = [Vr0,Ca0,Cb0,Tr0,Tj0]

# ------Define the plant model-----------
function dotStates(inputStates,p,t)

    Vr,Ca,Cb,Tr,Tj = inputStates

    u = p[1:2]
    par = p[3:4] # par[1] = K par[2] = H
                                                                                                                    
    Aw = (π*(r^2))+((0.002*Vr)/(r))

    VrDot = u[1]
    CaDot = -((u[1]*Ca)/Vr) - (par[1]*Ca*Cb)
    CbDot = ((u[1]/Vr)*(Cbin-Cb)) - (par[1]*Ca*Cb)
    TrDot = ((u[1]/Vr)*(Tin-Tr)) - ((α*Aw*(Tr-Tj))/(ρ*Vr*Cp)) - ((par[1]*Ca*Cb*par[2])/(ρ*Cp))
    TjDot = (u[2] + α*Aw*(Tr-Tj))/(ρ*Vj*Cp)

    return [VrDot,CaDot,CbDot,TrDot,TjDot]
end

# ------Define the ODEsolver-----------
function ODEmodel(x,u,t)
    return dotStates(x,u,t)
end

#------Define the Plant-------------
function PlantODE(x0,u0,dt)
    tspan = (0.0,dt)
    function f(x,u,t)
        return ODEmodel(x,u,t)
    end

    prob = ODEProblem(f,x0,tspan,u0)

    sol = DifferentialEquations.solve(prob)
    x = sol.u
    xk = sol.u[end]

    return xk
end

# Since we can calculate Cc and Aw explicitly from the states we do not need to code a DAEsolver.