ODEint.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#include "ODEint.h"

namespace limpet {

void ODEint_RK(GlobalData_t *y, GlobalData_t *ydot, void f(GlobalData_t*,GlobalData_t*,void*), void *f_data, int len, GlobalData_t dt)
{
  int     i;
  GlobalData_t  k1[ODEint_MaxVecLen];
  GlobalData_t *ynpk1 = k1;
  GlobalData_t *k2    = k1;

  for (i=0;i<len;i++)  {
    k1[i]    = dt*ydot[i];
    ynpk1[i] = y[i]+0.5*k1[i];
  }

  f(ynpk1,ydot,f_data);

  for (i=0;i<len;i++) {
    k2[i] = dt*ydot[i];
    y[i] += k2[i];
  }
}

void ODEint_FE(GlobalData_t *y, GlobalData_t *ydot, void *f, void *f_data, int len, GlobalData_t dt)
{
  for (int i=0;i<len;i++)
    y[i] += ydot[i]*dt;
}


void d_ODEint_RK(d_OdeVec *y, d_OdeVec *ydot, void f(GlobalData_t*,d_OdeVec*,void*),
                 void *f_data, GlobalData_t dt)
{
  int     i;
  GlobalData_t  k1[ODEint_MaxVecLen];
  GlobalData_t *ynpk1 = k1;
  GlobalData_t *k2    = k1;

  for (i=0;i<y->len;i++)  {
    k1[i]    = dt*ydot->vec[i];
    ynpk1[i] = y->vec[i]+0.5*k1[i];
  }

  f(ynpk1,ydot,f_data);

  for (i=0;i<y->len;i++) {
    k2[i]      = dt*ydot->vec[i];
    y->vec[i] += k2[i];
  }
}

void d_ODEint_FE(d_OdeVec*y, d_OdeVec *ydot, void *f, void *f_data, GlobalData_t dt)
{
  for (int i=0;i<y->len;i++)
    y->vec[i] += ydot->vec[i]*dt;
}

}  // namespace limpet