src/sys_exc.c

   1 /* copyright 2016 Apache 2 sddekit authors */
   2
   3 #include "sddekit.h"
   4
   5 struct pars {
   6         double a, tau, D, k;
   7 };
   8
   9 struct pars pars_default = {.a=1.1, .tau=3.0, .D=1e-3, .k=0.01};
  10
  11 struct data {
  12         sd_sys sys_if;
  13         sd_sys_exc exc_if;
  14         struct pars pars;
  15 };
  16
  17 static void exc_free(sd_sys *s) { sd_free(s->ptr); }
  18
  19 /* get set */
  20 #define GETSET(field) \
  21 static double get_ ## field(sd_sys_exc *s) { \
  22         struct data *d = (struct data *) (*s->sys)(s)->ptr; \
  23         return d->pars.field; \
  24 }\
  25 static void set_ ## field(sd_sys_exc *s, double val) { \
  26         struct data *d = (struct data *) (*s->sys)(s)->ptr; \
  27         d->pars.field = val; \
  28 }
  29
  30 GETSET(a)
  31 GETSET(tau)
  32 GETSET(D)
  33 GETSET(k)
  34
  35 #undef GETSET
  36
  37 static sd_stat apply(sd_sys *s, sd_sys_in *in, sd_sys_out *out)
  38 {
  39         struct pars *d = &(((struct data*) s->ptr)->pars);
  40         double *x = in->x, *c = in->i, *f = out->f, *g = out->g, *o = out->o;
  41         f[0] = (x[0] - x[0]*x[0]*x[0]/3.0 + x[1]) * d->tau;
  42         f[1] = (d->a - x[0] + d->k*c[0]) / d->tau;
  43         g[0] = d->D;
  44         g[1] = d->D;
  45         o[0] = x[0];
  46         return SD_OK;
  47 }
  48
  49 uint32_t  ndim(sd_sys *s) { (void) s;  return 2; }
  50 uint32_t   ndc(sd_sys *s) { (void) s;  return 1; }
  51 uint32_t  nobs(sd_sys *s) { (void) s;  return 1; }
  52 uint32_t nrpar(sd_sys *s) { (void) s;  return 4; }
  53 uint32_t nipar(sd_sys *s) { (void) s;  return 0; }
  54
  55 sd_sys sys_default = {
  56         .ptr = NULL,
  57         .ndim = &ndim,
  58         .ndc = &ndc,
  59         .nobs = &nobs,
  60         .nrpar = &nrpar,
  61         .nipar = &nipar,
  62         .apply = &apply,
  63         .free = &exc_free
  64 };
  65
  66 sd_sys * get_sys(sd_sys_exc *e) { return &(((struct data *) e->ptr)->sys_if); }
  67
  68 sd_sys_exc exc_default = {
  69         .ptr = NULL,
  70         .sys = &get_sys,
  71         .get_a = &get_a,
  72         .set_a = &set_a,
  73         .get_tau = &get_tau,
  74         .set_tau = &set_tau,
  75         .get_D = &get_D,
  76         .set_D = &set_D,
  77         .get_k = &get_k,
  78         .set_k = &set_k
  79 };
  80
  81 sd_sys_exc *sd_sys_exc_new() {
  82         struct data *d, zero = {0};
  83         if ((d = sd_malloc (sizeof(struct data))) == NULL)
  84         {
  85                 sd_err("alloc exc sys failed.");
  86                 return NULL;
  87         }
  88         *d = zero;
  89         d->pars = pars_default;
  90         d->sys_if = sys_default;
  91         d->exc_if = exc_default;
  92         d->sys_if.ptr = d->exc_if.ptr = d;
  93         return &(d->exc_if);
  94 }