1 /* delays.c
2
3 Originally written by Frederic Bois
4
5 Copyright (c) 2015-2017 Free Software Foundation, Inc.
6
7 This file is part of GNU MCSim.
8
9 GNU MCSim is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License
11 as published by the Free Software Foundation; either version 3
12 of the License, or (at your option) any later version.
13
14 GNU MCSim is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
21
22 Contains routines for handling delay differential equations (essentially
23 initialization of memory arrays, storage and retrieval.
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <math.h>
30
31 #include "lexerr.h"
32 #include "simmonte.h"
33 #include "yourcode.h"
34
35
36 /* ----------------------------------------------------------------------------
37 Constants
38 */
39
40 /* Delay storage size */
41 #define MAX_DELAY 1000
42
43
44 /* -----------------------------------------------------------------------------
45 externs: delay management stuff
46 */
47 extern double vrgModelVars[]; /* States and Outputs */
48
49
50 /* WARNING! ===== GLOBAL VARIABLES ! */
51
52 /* -----------------------------------------------------------------------------
53 globals to this file
54 */
55 static double dInitialTime;
56 int iCurrentTime;
57 double *rgdTime = NULL;
58 long *rgiVars = NULL;
59 double **pdVar = NULL;
60
61
62 /* -----------------------------------------------------------------------------
63 CalcDelay
64
65 Returns the value of a state or output variable at a given (past) time.
66 Past values are stored in ... The step back in time is given "delay", the
67 current time is dTime. A loop is used.
68 */
69
CalcDelay(HVAR hvar,double dTime,double delay)70 double CalcDelay (HVAR hvar, double dTime, double delay)
71 {
72 int i;
73 int sentinel; /* stopping flag for circular array */
74 double dTmp, oldTime;
75
76 /* printf("\n\n in CalcDelay, time = %g", dTime); */
77
78 /* if the handle has not been seen, create the storage for that variable.
79 We should check that it's a state or output variable... */
80 if (!rgiVars[hvar]) {
81 /* printf("\n => creating past array for variable %ld...\n", hvar); */
82 pdVar[hvar] = InitdVector(MAX_DELAY);
83 pdVar[hvar][0] = vrgModelVars[hvar];
84 rgiVars[hvar] = 1; /* mark as initialized */
85 }
86
87 /* negative or null delays are forbidden */
88 /* that should be handled by mod, or stay dynamic? */
89 if (delay <= 0) {
90 printf ("\nError: negative or null delays aren't allowed - Exiting.\n");
91 exit(0);
92 /* for null delays the variable's current value could be returned... */
93 }
94 else { /* delay strictly positive, OK */
95
96 /* printf("\n back time: %e - %e = %e\n", dTime, delay, dTime - delay); */
97
98 oldTime = dTime - delay;
99 if (oldTime <= dInitialTime) {
100 /* printf("\n default var %ld at val %g\n", hvar, pdVar[hvar][0]); */
101 return (pdVar[hvar][0]);
102 }
103 else { /* past time sought is after initial time */
104
105 /* find the index of time (dTime - delay).
106 The current time is stored in rgdTime[iCurrentTime].
107 We could save on searching time for multiple calls at the same time */
108 i = iCurrentTime - 1;
109 if (i < 0)
110 i = MAX_DELAY - 1;
111 sentinel = 0;
112 while (rgdTime[i] > oldTime) {
113 i = i - 1;
114 if (i < 0)
115 i = MAX_DELAY - 1;
116 sentinel = sentinel + 1;
117 /* printf("\n sentinel: %d\n", sentinel); */
118 if (sentinel > MAX_DELAY - 1) {
119 printf ("Error: size MAX_DELAY of rgdTime array = "
120 "%ld too small.\n", (long) MAX_DELAY);
121 exit(0);
122 }
123 } /* end while */
124
125 /* compute delayed term */
126 if (i == (iCurrentTime - 1)) {
127 /* an integration step larger than delay is attempted; it's a special
128 case: interpolate between i-1 and i, rather than between i and
129 iCurrentTime (which points to an empty cell or outdated ones).
130 Note that deSolve forbids that... */
131 dTmp = pdVar[hvar][i-1] +
132 ((pdVar[hvar][i] - pdVar[hvar][i-1]) *
133 (oldTime - rgdTime[i-1]) / (rgdTime[i] - rgdTime[i-1]));
134 }
135 else {
136 /* time step < delay, OK, interpolate between i and i+1 */
137 dTmp = pdVar[hvar][i] +
138 ((pdVar[hvar][i+1] - pdVar[hvar][i]) *
139 (oldTime - rgdTime[i]) / (rgdTime[i+1] - rgdTime[i]));
140 }
141
142 /* printf("\n computed var %ld at val %g\n", hvar, dTmp); */
143 }
144
145 return dTmp;
146
147 } /* else */
148
149 } /* CalcDelay */
150
151
152 /* ---------------------------------------------------------------------------
153 InitDelays
154
155 routine called at start if delay differential equations are used.
156 Initializes storage and internal bookkeeping variables.
157 */
InitDelays(double dTime)158 void InitDelays (double dTime)
159 {
160 int i;
161
162 /* printf("\n\n in InitDelays"); */
163
164 /* initialize the time storage */
165 if (!rgdTime) {
166 rgdTime = InitdVector(MAX_DELAY);
167 iCurrentTime = -1;
168 dInitialTime = dTime;
169 /* printf("\n => past time array created."); */
170 }
171
172 /* initialize a table of the delayed variables
173 note that the storage vectors pdVar[x] will be init in CalcDelay
174 upon first call to that function for a given x (state or output)
175 That is messy and should be handled by the model generator */
176 if (!rgiVars) {
177 rgiVars = InitlVector(GetNModelVars());
178 pdVar = InitpdVector(GetNModelVars());
179 for (i = 0; i < GetNModelVars(); i++)
180 rgiVars[i] = 0;
181 }
182
183 } /* InitDelays */
184
185
186 /* ---------------------------------------------------------------------------
187 StoreDelayed
188
189 routine called after each successful step of the integrator.
190 It stores current time, states and outputs values for later retrieval.
191 */
StoreDelayed(double t)192 void StoreDelayed (double t)
193 {
194 int i;
195
196 /* printf("\n\n in StoreDelayed, time %g\n", t); */
197
198 iCurrentTime++;
199 if (iCurrentTime == MAX_DELAY) /* loop */
200 iCurrentTime = 0;
201
202 rgdTime[iCurrentTime] = t;
203
204 for (i = 0; i < GetNModelVars(); i++)
205 if (rgiVars[i]) {
206 pdVar[i][iCurrentTime] = vrgModelVars[i];
207 /* printf(" stored var %d, val %g\n", i, vrgModelVars[i]); */
208 }
209
210 } /* StoreDelayed */
211
212
213 /* End */
214