1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
|
/* file lorenz.c */
#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
#include <R_ext/Rdynload.h>
#include <Rmath.h>
#include <stdio.h>
#include <math.h>
static double parms[10];
#define K parms[0]
#define kapa parms[1]
#define t0 parms[2]
#define alpha parms[3]
#define sigma parms[4]
#define ts parms[5]
#define te parms[6]
#define gamma parms[7]
#define beta parms[8]
#define c parms[9]
/* initializer */
void initmod(void (* odeparms)(int *, double *))
{
int N = 10;
odeparms(&N, parms);
}
#define S y[0]
#define I y[1]
#define R y[2]
#define dSdt ydot[0]
#define dIdt ydot[1]
#define dRdt ydot[2]
/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
double *yout, int *ip)
{
if (ip[0] < 1) error("nout should be at least 1") ;
double gammabis ;
if(*t<=te) gammabis=0* *t;
if(*t>te) gammabis=gamma;
double betabis ;
if(*t<=te) betabis=0* *t;
if(*t>te) betabis=beta;
double f1;
if(*t<=te) f1=0* *t;
if(*t>te) f1=powl(- R ,c)+1;
double tau ;
if(*t<=ts) tau = alpha * exp(-kapa * (*t-t0));
if(*t>ts) tau=0* *t;
double delta ;
if(*t<=ts) delta=0* *t;
if(*t>ts) delta=sigma *(1- exp(-kapa * (*t-ts)));
dSdt = -betabis*S*I *f1 + tau * S * (1 - S / K )- delta*S ;
dIdt = betabis * S * I *f1 - (gammabis + delta )*I;
dRdt= delta*S +(gammabis + delta)*I;
yout[0] = S + I + R ;
}
/* END file mymod.c */ |
Partager