{% if (ODEmethod == "cvode") %}

{{modelName}}::{{modelName}}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}

{{modelName}}::{{modelName}}(const {{modelName}}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}


{{modelName}}::~{{modelName}}()
{

  if ( B_.sys_ )
  {
    /* Free y vector */
    N_VDestroy_Serial(B_.y);

    /* Free integrator memory */
    if (B_.sys_ != NULL)
    {
      CVodeFree(&B_.sys_);
      B_.sys_ = NULL;
    }
  }
}


void {{modelName}}::init_node_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}


void {{modelName}}::init_state_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    S_ = State_(pr.P_);
}


void {{modelName}}::init_buffers_()
{

   {% for synapticEvent in synapticEventDefs %}
   B_.spike_{{synapticEvent.pscId}}.clear();
   {% endfor %}

   B_.currents_.clear();           
   Archiving_Node::clear_history();

   B_.logger_.reset();

   B_.step_ = Time::get_resolution().get_ms();
   B_.IntegrationStep_ = B_.step_;

   B_.I_stim_ = 0.0;

   int status, N, rootdir;

   N = {{stateSize}};
   // only positive direction (rising edge) of spike events will be detected
   rootdir = 1;

   /* Creates serial vector of length N */
   B_.y = N_VNew_Serial(N);
   if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw CVodeSolverFailure (get_name(), 0);

   for (int i = 0; i < N; i++)
   {
      Ith(B_.y,i) = S_.y_[i];
   }
 
   /* Calls CVodeCreate to create the solver memory and specify the 
    * Backward Differentiation Formula and the use of a Newton iteration */
   B_.sys_ = CVodeCreate(CV_BDF, CV_NEWTON);
   if (check_flag((void *)B_.sys_, "CVodeCreate", 0)) throw CVodeSolverFailure (get_name(), 0);

   /* Calls CVodeInit to initialize the integrator memory and specifies the
    * right hand side function in y'=f(t,y), the initial time, and
    * the initial values. */
   status = CVodeInit (B_.sys_, {{modelName}}_dynamics, 0.0, B_.y);
   if (check_flag(&status, "CVodeInit", 1)) throw CVodeSolverFailure (get_name(), status);

{% if ("V_t" in defaultDefs) %}

   /* Spike event handler (detects zero-crossing of V-V_t) */
   status = CVodeRootInit(B_.sys_, 1, (CVRootFn){{modelName}}_event);
   if (check_flag(&status, "CVodeRootInit", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Detect only the rising edge of spikes */
   status = CVodeSetRootDirection(B_.sys_, &rootdir);
   if (check_flag(&status, "CVodeSetRootDirection", 1)) throw CVodeSolverFailure (get_name(), status);

{% endif %}

   /* Sets the relative and absolute error tolerances of CVode  */
   status = CVodeSStolerances (B_.sys_, 
                               {% if abstol %}{{abstol}}{% else %}1e-7{% endif %}, 
                               {% if reltol %}{{reltol}}{% else %}1e-7{% endif %});
   if (check_flag(&status, "CVodeSStolerances", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Turns on CVode stability limit detection (only applicable for order 3 and greater) */
   status = CVodeSetStabLimDet (B_.sys_,true);
   if (check_flag(&status, "CVodeSetStabLimDet", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Sets the maximum order of CVode  */
   status = CVodeSetMaxOrd (B_.sys_,5);
   if (check_flag(&status, "CVodeSetMaxOrd", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Sets maximum step size. */
   status = CVodeSetMaxStep (B_.sys_,{% if maxstep %}{{maxstep}}{% else %}B_.step_{% endif %});
   if (check_flag(&status, "CVodeSetMaxStep", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Configures the integrator to pass the params structure to the right-hand function */
   status = CVodeSetUserData(B_.sys_, reinterpret_cast<void*>(this));
   if (check_flag(&status, "CVodeSetUserData", 1)) throw CVodeSolverFailure (get_name(), status);

   /* Initializes diagonal linear solver. */
   status = CVDiag (B_.sys_);
   if (check_flag(&status, "CVDiag", 1)) throw CVodeSolverFailure (get_name(), status);
  }


void {{modelName}}::calibrate()
{
    B_.logger_.init();  
}


{% elif (ODEmethod == "ida") %}

{{modelName}}::{{modelName}}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}


{{modelName}}::{{modelName}}(const {{modelName}}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}


{{modelName}}::~{{modelName}}()
{

  if ( B_.sys_ )
  {
    /* Free y vector */
    N_VDestroy_Serial(B_.y);
    N_VDestroy_Serial(B_.yp);

    /* Free integrator memory */
    if (B_.sys_ != NULL)
    {
      IDAFree(&B_.sys_);
      B_.sys_ = NULL;
    }

  }
}


void {{modelName}}::init_node_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}


void {{modelName}}::init_state_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    S_ = State_(pr.P_);
}


void {{modelName}}::init_buffers_()
{

   {% for synapticEvent in synapticEventDefs %}
   B_.spike_{{synapticEvent.pscId}}.clear();
   {% endfor %}

   B_.currents_.clear();           
   Archiving_Node::clear_history();

   B_.logger_.reset();

   B_.step_ = Time::get_resolution().get_ms();
   B_.IntegrationStep_ = B_.step_;

   B_.I_stim_ = 0.0;


   int status, N, rootdir;

   N = {{stateSize}};

   // only positive direction (rising edge) of spike events will be detected
   rootdir = 1;

   /* Creates serial vectors of length N */
   B_.y = N_VNew_Serial(N);
   B_.y1 = N_VNew_Serial(N);
   B_.yp = N_VNew_Serial(N);
   if (check_flag((void *)B_.y, "N_VNew_Serial", 0)) throw IDASolverFailure (get_name(), 0);

   for (int i = 0; i < N; i++)
   {
      Ith(B_.y,i) = S_.y_[i];
   }
 
   {{modelName}}_dynamics (0.0, B_.y, B_.yp, reinterpret_cast<void*>(this));

   /* Calls IDACreate to create the solver memory */ 
   B_.sys_ = IDACreate();
   if (check_flag((void *)B_.sys_, "IDACreate", 0)) throw IDASolverFailure (get_name(), 0);

  /* Calls IDAInit to initialize the integrator memory and specify the
   * resdual function, the initial time, and the initial values. */
   status = IDAInit (B_.sys_, {{modelName}}_residual, 0.0, B_.y, B_.yp);

   if (check_flag(&status, "IDAInit", 1)) throw IDASolverFailure (get_name(), status);

{% if ("V_t" in defaultDefs) %}

   /* Spike event handler (detects zero-crossing of V-V_t) */
   status = IDARootInit(B_.sys_, 1, (IDARootFn){{modelName}}_event);
   if (check_flag(&status, "IDARootInit", 1)) throw IDASolverFailure (get_name(), status);

   /* Detect only the rising edge of spikes */
   status = IDASetRootDirection(B_.sys_, &rootdir);
   if (check_flag(&status, "IDASetRootDirection", 1)) throw IDASolverFailure (get_name(), status);

{% endif %}

   /* Sets the relative and absolute error tolerances of IDA  */
   status = IDASStolerances (B_.sys_, 
                             {% if abstol %}{{abstol}}{% else %}1e-7{% endif %}, 
                             {% if reltol %}{{reltol}}{% else %}1e-7{% endif %});
   if (check_flag(&status, "IDASStolerances", 1)) throw IDASolverFailure (get_name(), status);

   /* Sets the maximum order of IDA  */
   status = IDASetMaxOrd (B_.sys_,5);
   if (check_flag(&status, "IDASetMaxOrd", 1)) throw IDASolverFailure (get_name(), status);

   /* Sets maximum step size. */
   status = IDASetMaxStep (B_.sys_,{% if maxstep %}{{maxstep}}{% else %}B_.step_{% endif %});
   if (check_flag(&status, "IDASetMaxStep", 1)) throw IDASolverFailure (get_name(), status);

   /* Calls IDASetUserData to configure the integrator to pass the 
    * params structure to the right-hand function */
   status = IDASetUserData(B_.sys_, reinterpret_cast<void*>(this));
   if (check_flag(&status, "IDASetUserData", 1)) throw IDASolverFailure (get_name(), status);

   /* Initializes dense linear solver. */
   status = IDADense (B_.sys_, N);
   if (check_flag(&status, "IDADense", 1)) throw IDASolverFailure (get_name(), status);

   status = IDACalcIC(B_.sys_, IDA_Y_INIT, 0.0);
   if (check_flag(&status, "IDACalcIC", 1)) throw IDASolverFailure (get_name(), status);

}

void {{modelName}}::calibrate()
{
   B_.logger_.init();  
}

{% elif (ODEmethod == "gsl") %}

{{modelName}}::{{modelName}}()
    : Archiving_Node(), 
      P_(), 
      S_(P_),
      B_(*this)
{
    recordablesMap_.create();
}


{{modelName}}::{{modelName}}(const {{modelName}}& n)
    : Archiving_Node(n), 
      P_(n.P_), 
      S_(n.S_),
      B_(n.B_, *this)
{
}


{{modelName}}::~{{modelName}} ()
{
    // GSL structs only allocated by init_nodes_(), so we need to protect destruction
    if ( B_.s_ != NULL) gsl_odeiv2_step_free (B_.s_);
    if ( B_.c_ != NULL) gsl_odeiv2_control_free (B_.c_);
    if ( B_.e_ != NULL) gsl_odeiv2_evolve_free (B_.e_);
    if ( B_.u != NULL) free (B_.u);
    if ( B_.jac != NULL) free (B_.jac);
}


void {{modelName}}::init_node_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    P_ = pr.P_;
    S_ = State_(P_);
}


void {{modelName}}::init_state_(const Node& proto)
{
    const {{modelName}}& pr = downcast<{{modelName}}>(proto);
    S_ = State_(pr.P_);
}


void {{modelName}}::init_buffers_()
{
   {% for synapticEvent in synapticEventDefs %}
   B_.spike_{{synapticEvent.pscId}}.clear();
   {% endfor %}

   B_.currents_.clear();           
   Archiving_Node::clear_history();

   B_.logger_.reset();

   B_.step_ = Time::get_resolution().get_ms();
   B_.IntegrationStep_ = B_.step_;

   B_.I_stim_ = 0.0;


   static const gsl_odeiv2_step_type* T1 = gsl_odeiv2_step_rk2;
   B_.N = {{stateSize}};
  
   if ( B_.s_ == 0 )
     B_.s_ = gsl_odeiv2_step_alloc (T1, B_.N);
   else 
     gsl_odeiv2_step_reset(B_.s_);
    
   if ( B_.c_ == 0 )  
     B_.c_ = gsl_odeiv2_control_standard_new ({% if abstol %}{{abstol}}{% else %}1e-7{% endif %}, 
                                              {% if reltol %}{{reltol}}{% else %}1e-7{% endif %}, 1.0, 0.0);
   else
     gsl_odeiv2_control_init(B_.c_, {% if abstol %}{{abstol}}{% else %}1e-7{% endif %}, 
                             {% if reltol %}{{reltol}}{% else %}1e-7{% endif %}, 1.0, 0.0);
    
   if ( B_.e_ == 0 )  
     B_.e_ = gsl_odeiv2_evolve_alloc(B_.N);
   else 
     gsl_odeiv2_evolve_reset(B_.e_);
   
   B_.sys_.function  = {{modelName}}_dynamics; 
   B_.sys_.jacobian  = {{modelName}}_jacobian;
   B_.sys_.dimension = B_.N;
   B_.sys_.params    = reinterpret_cast<void*>(this);

   B_.u = (double *)malloc(sizeof(double) * B_.N);
   assert (B_.u);
   B_.jac = (double *)malloc(sizeof(double) * B_.N);
   assert (B_.jac);

}


void {{modelName}}::calibrate()
{
    B_.logger_.init();  
    {% if ("v" in stateIndexMap) %}V_.U_old_ = S_.y_[{{stateIndexMap | attr("v")}}];{% endif %}
{% if ("t_ref" in defaultDefs) %}
    V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref)).get_steps();
{% else %}
    V_.RefractoryCounts_ = 0;
{% endif %}
}


{% endif %}