# -*- coding: utf-8 -*-
# {{ model_name }}_cmdl.jl - DAE Command-line version

using DifferentialEquations
using Sundials  # For IDA solver

# Parameters
{% for name, value in parameters.items() %}
const {{ name }} = {{ value }}
{% endfor %}

# Time parameters
const t0 = {{ t0 }}
const t1 = {{ t1 }}
const dt = {{ dt }}

function dae!(out, du, u, p, t)
    # Extract state variables
    {% for name in variable_names %}
    {{ name }} = u[{{ loop.index }}]
    {% endfor %}

    # Extract derivatives
    {% for name in variable_names %}
    d{{ name }}_dt = du[{{ loop.index }}]
    {% endfor %}

    # Auxiliary equations
    {% for name, expr in auxiliary_equations.items() %}
    {{ name }} = {{ expr }}
    {% endfor %}

    # Compute f_<var> expressions
    {% for name, expr in derivative_computations %}
    {{ name }} = {{ expr }}
    {% endfor %}

    {% for name in variable_names %}
    out[{{ loop.index }}] = d{{ name }}_dt - f_{{ name }}
    {% endfor %}
end

# Initial conditions for state variables
u0 = [
    {% for name, value in initial_conditions.items() %}
    {{ value }}{% if not loop.last %},{% endif %}
    {% endfor %}
]

# Initial guess for derivatives (can be zeros)
du0 = zeros({{ variable_count }})

# Problem setup
tspan = (t0, t1)
# The IDA solver requires the `differential_vars` argument to specify which
# variables are differential (true) and which are algebraic (false).
# This assumes all variables are differential.
prob = DAEProblem(dae!, du0, u0, tspan, differential_vars = [{{ differential_vars_list | join(", ") }}])

# Output file
outfile = open("models/{{ model_name }}.csv", "w")
write(outfile, "t,{{ variable_names | join(",") }}\n")

# Callback for writing results
step_callback = function (integrator)
    t = integrator.t
    y = integrator.u
    write(outfile, string(t))
    {% for name in variable_names %}
    write(outfile, "," * string(y[{{ loop.index }}]))
    {% endfor %}
    write(outfile, "\n")
    flush(outfile)
    return false
end

# Solve the DAE
cb = DiscreteCallback((f,t,integrator)->true, step_callback)
sol = solve(prob, IDA(), dt=dt, adaptive=false, callback=cb, abstol=1e-8, reltol=1e-6)

# Cleanup
close(outfile)
println("Simulation completed successfully")