1 /* modelu.c
2
3 Written by Don Maszle
4 7 October 1991
5
6 Copyright (c) 1991-2017 Free Software Foundation, Inc.
7
8 This file is part of GNU MCSim.
9
10 GNU MCSim is free software; you can redistribute it and/or
11 modify it under the terms of the GNU General Public License
12 as published by the Free Software Foundation; either version 3
13 of the License, or (at your option) any later version.
14
15 GNU MCSim is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
22
23 Model utilities.
24
25 This file contains a number of utility routines used by the model
26 file written by the code generator. These routines must know of
27 the structure of the model file, as they manipulate the model
28 variables. Since they do not vary with the model, however, it
29 makes sense to have them in their own file, instead of being
30 written each time.
31 */
32
33 #include <float.h> /* Floating point limits */
34 #include <math.h>
35 #include <stdlib.h>
36 #include <stdio.h>
37 #include <string.h>
38
39 #include "lexerr.h"
40 #include "modelu.h"
41 #include "strutil.h"
42 #include "hungtype.h"
43
44
45 /* -----------------------------------------------------------------------------
46 Model externs
47
48 The following will be defined inside the file created by the
49 code generator. These names must match with those that will be
50 written. These should NOT be in the include file created by the
51 generator as they are private variable to the model section of the
52 code.
53 */
54 extern int vnStates;
55 extern int vnOutputs;
56 extern int vnModelVars;
57 extern int vnInputs;
58 extern int vnParms;
59
60 extern double vrgModelVars[]; /* States and Outputs */
61 extern IFN vrgInputs[]; /* Input Function records */
62 extern VMMAPSTRCT vrgvmGlo[]; /* Global variable map */
63
64 /* WARNING! ===== GLOBAL VARIABLE ! */
65
66 /* This is used by UpdateInputs() which needs to know if the model
67 has been re-initialized in order to fixup dependent inputs, and
68 to set the start period heuristic. */
69
70 BOOL vbModelReinitd = FALSE; /* Model has been initialized */
71
72
73 /* ----------------------------------------------------------------------------
74 FixupDependentInputs
75
76 Fixes input definitions dependent on model variables.
77
78 This routine resolves input parameter that are dependent on model
79 parameters. If a parameter has a non-zero handle, the nominal
80 value of the variable associated with that handle is used for the
81 nominal value of the input parameter. As an example, an
82 exponential input might have a decay rate dependent on a model
83 variable called 'Decay', which changes with each experiment. The
84 handle to 'Decay' would be the hDecay parameter of that input
85 function.
86
87 This routine also adjust variant exposure times.
88 */
89
FixupDependentInputs(void)90 void FixupDependentInputs (void)
91 {
92 int i, j;
93
94 for (i = 0; i < vnInputs; i++) {
95
96 /* Resolve dependencies */
97
98 if (vrgInputs[i].hMag)
99 vrgInputs[i].dMag = GetVarValue (vrgInputs[i].hMag);
100 if (vrgInputs[i].hTper)
101 vrgInputs[i].dTper = GetVarValue (vrgInputs[i].hTper);
102 if (vrgInputs[i].hT0)
103 vrgInputs[i].dT0 = GetVarValue (vrgInputs[i].hT0);
104 if (vrgInputs[i].hTexp)
105 vrgInputs[i].dTexp = GetVarValue (vrgInputs[i].hTexp);
106 if (vrgInputs[i].hDecay)
107 vrgInputs[i].dDecay = GetVarValue (vrgInputs[i].hDecay);
108 if (vrgInputs[i].hNcpt)
109 vrgInputs[i].dNcpt = GetVarValue (vrgInputs[i].hNcpt);
110
111 /* NDoses or Spikes dependencies */
112 if ((vrgInputs[i].iType == IFN_NDOSES) || /* FB 12/10/96 */
113 (vrgInputs[i].iType == IFN_SPIKES) || /* FB 11/06/97 */
114 (vrgInputs[i].iType == IFN_EVENTS)) /* FB 15/10/16 */
115 for (j = 0; j < vrgInputs[i].nDoses; j++) {
116 if (vrgInputs[i].rghMags[j])
117 vrgInputs[i].rgMags[j] = GetVarValue (vrgInputs[i].rghMags[j]);
118 if (vrgInputs[i].rghT0s[j])
119 vrgInputs[i].rgT0s[j] = GetVarValue (vrgInputs[i].rghT0s[j]);
120
121 /* check that the times are in order */
122 if (j > 0) {
123 if (vrgInputs[i].rgT0s[j] <= vrgInputs[i].rgT0s[j-1]) {
124 printf ("\nError: unordered pair of times (%g, %g) in %s "
125 "statement - Exiting\n",
126 vrgInputs[i].rgT0s[j-1], vrgInputs[i].rgT0s[j],
127 (vrgInputs[i].iType == IFN_NDOSES ? "NDoses" : "Spikes"));
128 exit (0);
129 }
130 }
131 }
132
133 /* set exponential exposure times */
134 if ((vrgInputs[i].iType == IFN_PEREXP) ||
135 (vrgInputs[i].iType == IFN_PERTRANS)) {
136 vrgInputs[i].dTexp = vrgInputs[i].dTper;
137 }
138
139 /* For PerTransit the number of compartments must be > 0
140 and do not adjust exposure time, just set it at the period */
141 if (vrgInputs[i].iType == IFN_PERTRANS) {
142 if (vrgInputs[i].dNcpt <= 0) {
143 printf ("\nError: null or negative number of virtual compartment "
144 "in PerTransit input function - Exiting\n");
145 exit (0);
146 }
147 vrgInputs[i].dTexp = vrgInputs[i].dTper;
148 }
149
150 /* FB 9/9/97, modified 29/5/2006: make sure that the exposure does not
151 exceed the period */
152 if (vrgInputs[i].dTexp >= vrgInputs[i].dTper)
153 vrgInputs[i].dTexp = vrgInputs[i].dTper;
154
155 /* For doses and spikes allow scaling in Scale{} by use of dMag */
156 if (vrgInputs[i].iType == IFN_NDOSES || vrgInputs[i].iType == IFN_SPIKES)
157 vrgInputs[i].dMag = 1.0; /* Initially no scaling */
158
159 /* FB 29/5/2006 removed a condition leading to failure if the start time
160 was superior to the period */
161 /* If the input is ill-defined, turn it off */
162 if (vrgInputs[i].iType == IFN_PERDOSE ||
163 vrgInputs[i].iType == IFN_PERTRANS ||
164 vrgInputs[i].iType == IFN_PEREXP) {
165 if (vrgInputs[i].dTexp == 0.0 ||
166 vrgInputs[i].dT0 < 0.0 ||
167 vrgInputs[i].dTper < 0.0)
168 vrgInputs[i].dMag = 0.0;
169 } /* if */
170
171 } /* for */
172
173 } /* FixupDependentInputs */
174
175
176 /* ----------------------------------------------------------------------------
177 GetStartPeriods
178
179 Calculates the start of the current period for each input. This
180 is 0.0 if *pdTime is 0.0, otherwise we have
181 to figure out where we are. Note that this is only called
182 once at the begining of a simulation and for most simulations, the
183 start time will be 0.
184 FB - 22 Oct 94 replaced a modf(x,dTmp) call by x - (long) (x)
185 */
186
GetStartPeriods(PDOUBLE pdTime)187 void GetStartPeriods (PDOUBLE pdTime)
188 {
189 int i;
190 double dTmp, dDummy;
191
192 if (*pdTime == 0.0)
193 for (i = 0; i < vnInputs; i++) {
194 vrgInputs[i].dTStartPeriod = 0.0;
195 if (vrgInputs[i].iType == IFN_NDOSES ||
196 vrgInputs[i].iType == IFN_SPIKES)
197 vrgInputs[i].iDoseCur = 0;
198 } /* for */
199
200 else
201 for (i = 0; i < vnInputs; i++) {
202 if (vrgInputs[i].iType == IFN_NDOSES ||
203 vrgInputs[i].iType == IFN_SPIKES) {
204 for (vrgInputs[i].iDoseCur = 0;
205 vrgInputs[i].iDoseCur < vrgInputs[i].nDoses;
206 vrgInputs[i].iDoseCur++)
207 if (*pdTime < vrgInputs[i].rgT0s[vrgInputs[i].iDoseCur])
208 break;
209 vrgInputs[i].iDoseCur--; /* The one that fails is too far */
210
211 /* this is for possible negative times... Frederic 20/3/96 */
212 if (vrgInputs[i].iDoseCur < 0) vrgInputs[i].iDoseCur = 0;
213 } /* if */
214
215 else {
216 if (vrgInputs[i].dTper == 0.0) /* Not periodic - only one pulse */
217 vrgInputs[i].dTStartPeriod = 0.0;
218 else {
219 dTmp = (*pdTime / vrgInputs[i].dTper);
220 /* FB 6/4/99 introduced modf */
221 vrgInputs[i].dTStartPeriod = vrgInputs[i].dTper *
222 modf (dTmp, &dDummy);
223 } /* else */
224 } /* else */
225 } /* for */
226
227 } /* GetStartPeriods */
228
229
230 /* -----------------------------------------------------------------------------
231 UpdateNDoses
232
233 Update an NDoses type input.
234 */
235
UpdateNDoses(PIFN pifn,PDOUBLE pdTnext,PDOUBLE pdTime)236 void UpdateNDoses (PIFN pifn, PDOUBLE pdTnext, PDOUBLE pdTime)
237 {
238 int j;
239
240 j = pifn->iDoseCur; /* Current pulse */
241
242 if (j < pifn->nDoses) { /* More pulses left */
243 *pdTnext = pifn->rgT0s[j];
244 pifn->bOn = (*pdTime >= *pdTnext);
245
246 if (pifn->bOn) {
247 /* Look at end of pulse */
248 *pdTnext = pifn->rgT0s[j+1];
249 pifn->bOn = (*pdTime < *pdTnext);
250
251 if (!pifn->bOn)
252 if (++pifn->iDoseCur < pifn->nDoses) {
253 *pdTnext = pifn->rgT0s[pifn->iDoseCur + 1]; /* Use next time */
254 pifn->bOn = TRUE;
255 }
256
257 } /* if pulse is On */
258 } /* if more pulses */
259 else
260 /* No pulses left, make dT unusable */
261 *pdTnext = DBL_MAX;
262
263 /* If not exposure, clear input */
264 if (!pifn->bOn)
265 pifn->dVal = 0.0;
266
267 } /* UpdateNDoses */
268
269
270 /* -----------------------------------------------------------------------------
271 UpdateSpikes
272
273 Update a Spikes type input.
274 */
UpdateSpikes(PIFN pifn,PDOUBLE pdTnext,PDOUBLE pdTime)275 BOOL UpdateSpikes (PIFN pifn, PDOUBLE pdTnext, PDOUBLE pdTime)
276 {
277 register double *rgT0s = pifn->rgT0s; /* The schedule times */
278 register int j = pifn->iDoseCur; /* Index to current spike */
279
280 *pdTnext = DBL_MAX; /* If no spikes, make dTnext infinite */
281 pifn->bOn = FALSE; /* FB 9/9/97 */
282 if (j < pifn->nDoses) {
283 if (*pdTime < rgT0s[j]) {
284 *pdTnext = rgT0s[j]; /* Time of upcoming spike */
285 }
286 else {
287 if (*pdTime == rgT0s[j]) {
288 /* At a spike, return the time of the following spike */
289 pifn->bOn = TRUE;
290 if (j + 1 < pifn->nDoses)
291 *pdTnext = rgT0s[j + 1];
292 }
293 else {
294 /* try match with next IEEE representable pdTime value */
295 double nextRepTime = nextafter(*pdTime, rgT0s[j]);
296 if (nextRepTime == rgT0s[j]) {
297 printf ("\n UpdateSpikes: Discontinuity time failed match corrected "
298 "for floating point precision at\n"
299 "simulation time = %.17f \n"
300 "event time = %.17f \n", *pdTime, rgT0s[j]);
301
302 /* return the time of the following spike */
303 pifn->bOn = TRUE;
304 if (j + 1 < pifn->nDoses)
305 *pdTnext = rgT0s[j + 1];
306 }
307 else { /* Oops */
308 printf ("\n UpdateSpikes: Discontinuity was passed over at\n"
309 "simulation time = %.17f \n"
310 "event time = %.17f \n", *pdTime, rgT0s[j]);
311 }
312 }
313 } /* else */
314
315 } /* if */
316
317 return (pifn->bOn);
318
319 } /* UpdateSpikes */
320
321
322 /* -----------------------------------------------------------------------------
323 UpdateDefaultInput
324
325 Update default input types : PerDose, PerExp, PerTransit.
326 */
UpdateDefaultInput(PIFN pifn,PDOUBLE pdTnext,PDOUBLE pdTime)327 void UpdateDefaultInput (PIFN pifn, PDOUBLE pdTnext, PDOUBLE pdTime)
328 {
329 *pdTnext = pifn->dTStartPeriod + pifn->dT0;
330 pifn->bOn = (*pdTime >= *pdTnext); /* Time is after pulse goes on */
331
332 if (pifn->bOn) {
333 *pdTnext += pifn->dTexp; /* Check second transition of pulse */
334 /* if time is reasonably before pulse goes off then pulse is on */
335 pifn->bOn = (*pdTime < *pdTnext) &&
336 ((*pdTnext - *pdTime) > DBL_EPSILON * 2.0 * (*pdTnext));
337
338 if (!pifn->bOn) {
339 if (pifn->dTper != 0.0) /* Periodic */
340 *pdTnext = (pifn->dTStartPeriod += pifn->dTper); /* Next period */
341 else
342 *pdTnext = pifn->dTStartPeriod = DBL_MAX - pifn->dTper;
343 } /* if */
344 } /* if */
345
346 if (!pifn->bOn) /* If not exposure, clear input */
347 pifn->dVal = 0.0;
348 else
349 pifn->dVal = pifn->dMag; /* Set to exposure */
350
351 } /* UpdateDefaultInput */
352
353
354 /* -----------------------------------------------------------------------------
355 UpdateInputs
356
357 The controlling routine to integrator must know when inputs change
358 states abruptly to correct for integration across transition
359 boundaries.
360
361 Given the current time it returns (in *pdNextTransTime) the time of the
362 next closest input transition (which is a possible discontinuity).
363
364 It also updates information about whether inputs are on
365 or off, increments the start of periods for periodic inputs.
366
367 It helps optimize the calculation of inputs which is done inside
368 the model. It does not calculate the inputs. CalcInputs() does that.
369
370 UpdateInputs() assumes that time proceeds positively.
371 */
UpdateInputs(PDOUBLE pdTime,PDOUBLE pdNextTransTime)372 void UpdateInputs (PDOUBLE pdTime, PDOUBLE pdNextTransTime)
373 {
374 double dT; /* New time to try */
375 int i, j, stateID;
376
377 if (vbModelReinitd) { /* if model has been re-initialized */
378 ScaleModel(pdTime); /* scale the model parms */
379 FixupDependentInputs(); /* fixup input fields dependent on vars */
380 GetStartPeriods(pdTime); /* initialize the current period */
381 }
382
383 dT = *pdNextTransTime = DBL_MAX;
384 for (i = 0; i < vnInputs; i++) {
385 switch (vrgInputs[i].iType) {
386
387 case IFN_CONSTANT:
388 break;
389
390 case IFN_NDOSES:
391 UpdateNDoses (&vrgInputs[i], &dT, pdTime);
392 break;
393
394 case IFN_EVENTS:
395 if (vrgInputs[i].bOn)
396 vrgInputs[i].iDoseCur++;
397
398 UpdateSpikes (&vrgInputs[i], &dT, pdTime);
399
400 j = vrgInputs[i].iDoseCur;
401 stateID = HINDEX(vrgInputs[i].target_state);
402 if (vrgInputs[i].bOn && (j < vrgInputs[i].nDoses))
403 vrgModelVars[stateID] =
404 (vrgInputs[i].rgOper[j] == KM_REPLACE ?
405 vrgInputs[i].rgMags[j] :
406 vrgInputs[i].rgOper[j] == KM_ADD ?
407 vrgModelVars[stateID] + vrgInputs[i].rgMags[j] :
408 vrgModelVars[stateID] * vrgInputs[i].rgMags[j]);
409
410 break;
411
412
413 case IFN_SPIKES:
414 /* FB 9/9/97: make it advance only if this input spike is on.
415 At the start it is off. The on/off switch is only used here,
416 CalcInputs does not use it but checks that the current time is
417 equal to the spike time */
418 if (vrgInputs[i].bOn)
419 vrgInputs[i].iDoseCur++;
420
421 UpdateSpikes (&vrgInputs[i], &dT, pdTime);
422
423 break;
424
425 default:
426 case IFN_PERDOSE:
427 case IFN_PEREXP:
428 case IFN_PERTRANS:
429 if (vrgInputs[i].dMag != 0.0)
430 UpdateDefaultInput (&vrgInputs[i], &dT, pdTime);
431 break;
432
433 } /* switch */
434
435 if (dT < *pdNextTransTime)
436 /* Return minimum time */
437 *pdNextTransTime = dT;
438
439 } /* for */
440
441 if (vbModelReinitd) { /* If model has just been initialized */
442 vbModelReinitd = FALSE; /* don't do it again */
443 }
444
445 } /* UpdateInputs */
446
447
448 /* -----------------------------------------------------------------------------
449 GetVarPtr
450
451 Returns a PVMMAPSTRCT to the structure of the variable
452 specified by szName, or NULL if it does not exist.
453 */
454
GetVarPtr(PVMMAPSTRCT pvm,PSTR szName)455 PVMMAPSTRCT GetVarPtr (PVMMAPSTRCT pvm, PSTR szName)
456 {
457 while (*pvm->szName && MyStrcmp (szName, pvm->szName))
458 pvm++;
459
460 return (*pvm->szName ? pvm : NULL);
461
462 } /* GetVarPtr */
463
464
465 /* -----------------------------------------------------------------------------
466 GetVarHandle
467
468 Returns a handle to the variable szName given. If the string
469 does not reference a declared variable, returns 0.
470 */
471
GetVarHandle(PSTR szName)472 HVAR GetVarHandle (PSTR szName)
473 {
474 PVMMAPSTRCT pvm = GetVarPtr (vrgvmGlo, szName);
475
476 return (pvm ? pvm->hvar : ID_NULL);
477
478 } /* GetVarHandle */
479
480
481 /* -----------------------------------------------------------------------------
482 GetVarType
483
484 Returns the ID_ type of the HVAR given. If hvar is a valid
485 handle, returns the type of the variable, else returns 0.
486 */
487
GetVarType(HVAR hvar)488 int GetVarType (HVAR hvar)
489 {
490 BOOL bOK = FALSE;
491
492 switch (HTYPE(hvar)) {
493 case ID_INPUT:
494 bOK = HINDEX(hvar) < vnInputs;
495 break;
496
497 case ID_STATE:
498 bOK = HINDEX(hvar) < vnStates;
499 break;
500
501 case ID_OUTPUT:
502 bOK = HINDEX(hvar) >= vnStates
503 && HINDEX(hvar) < vnModelVars;
504 break;
505
506 case ID_PARM:
507 {
508 int nSOI = vnStates + vnOutputs + vnInputs;
509
510 bOK = (HINDEX(hvar) >= nSOI) && (HINDEX(hvar) < nSOI + vnParms);
511 break;
512 } /* ID_PARM */
513
514 default:
515 break;
516 } /* switch */
517
518 return (bOK ? HTYPE(hvar) : ID_NULL);
519
520 } /* GetVarType */
521
522
523 /* -----------------------------------------------------------------------------
524 GetVarName
525
526 returns the name of HVAR given, or "InvalidVariable?" if hvar is
527 not a valid handle.
528
529 This is embarrassingly inefficient because of the screwy manner
530 in which I setup the symbol table.
531 */
532
GetVarName(HVAR hvar)533 char *GetVarName (HVAR hvar)
534 {
535 static char szInvalid[] = "InvalidVariable?";
536 PVMMAPSTRCT pvm = vrgvmGlo;
537
538 while (*pvm->szName && hvar != pvm->hvar)
539 pvm++;
540
541 return (*pvm->szName ? pvm->szName : szInvalid);
542
543 } /* GetVarType */
544
545
546 /* -----------------------------------------------------------------------------
547 GetVarValue
548
549 Returns the value of the HVAR given. If hvar is not a valid
550 handle, returns 0.0.
551 */
552
GetVarValue(HVAR hvar)553 double GetVarValue (HVAR hvar)
554 {
555 double dReturn = 0.0;
556
557 switch (GetVarType(hvar)) {
558 case ID_INPUT:
559 dReturn = vrgInputs[HINDEX(hvar)].dVal;
560 break;
561
562 case ID_OUTPUT:
563 case ID_STATE:
564 dReturn = vrgModelVars[HINDEX(hvar)];
565 break;
566
567 case ID_PARM:
568 {
569 PVMMAPSTRCT pvm;
570
571 pvm = &vrgvmGlo[HINDEX(hvar)];
572 dReturn = *(PDOUBLE) pvm->pVar;
573 break;
574 } /* ID_PARM: */
575
576 default:
577 break;
578 } /* switch */
579
580 return (dReturn);
581
582 } /* GetVarValue */
583
584
585 /* -----------------------------------------------------------------------------
586 IsInput
587
588 Returns TRUE if hvar is an input
589 */
590
IsInput(HVAR hvar)591 BOOL IsInput (HVAR hvar)
592 {
593 return (GetVarType(hvar) == ID_INPUT);
594
595 } /* IsInput */
596
597 /* -----------------------------------------------------------------------------
598 IsState
599
600 Returns TRUE if hvar is a state
601 */
602
IsState(HVAR hvar)603 BOOL IsState (HVAR hvar)
604 {
605 return (GetVarType(hvar) == ID_STATE);
606
607 } /* IsState */
608
609
610 /* -----------------------------------------------------------------------------
611 IsOutput
612
613 Returns TRUE if hvar is a output
614 */
615
IsOutput(HVAR hvar)616 BOOL IsOutput (HVAR hvar)
617 {
618 return (GetVarType(hvar) == ID_OUTPUT);
619
620 } /* IsOutput */
621
622
623 /* -----------------------------------------------------------------------------
624 IsModelVar
625
626 Returns TRUE if hvar is a state or output
627 */
628
IsModelVar(HVAR hvar)629 BOOL IsModelVar (HVAR hvar)
630 {
631 int iType = GetVarType(hvar);
632 return (iType == ID_STATE || iType == ID_OUTPUT);
633
634 } /* IsState */
635
636
637 /* -----------------------------------------------------------------------------
638 IsParm
639
640 Returns TRUE if hvar is a parameter
641 */
642
IsParm(HVAR hvar)643 BOOL IsParm (HVAR hvar)
644 {
645 return (GetVarType(hvar) == ID_PARM);
646 } /* IsParm */
647
648
649 /* -----------------------------------------------------------------------------
650 ModelIndex
651
652 Returns the index of the model variable given or 0 if hvar is not
653 a model varible. Use this with caution.
654 */
655
ModelIndex(HVAR hvar)656 int ModelIndex (HVAR hvar)
657 {
658 return (IsModelVar(hvar) ? HINDEX(hvar) : 0);
659
660 } /* ModelIndex */
661
662
663 /* -----------------------------------------------------------------------------
664 SetVar
665
666 Sets the variable referenced by handle to the value dVal. Returns
667 TRUE on success, FALSE if szName is not declared, or is input.
668 */
669
SetVar(HVAR hvar,double dVal)670 BOOL SetVar (HVAR hvar, double dVal)
671 {
672 BOOL bReturn = TRUE;
673
674 switch (GetVarType (hvar)) {
675 default:
676 case ID_INPUT:
677 bReturn = FALSE;
678 break;
679
680 case ID_OUTPUT:
681 case ID_STATE:
682 vrgModelVars[HINDEX(hvar)] = dVal;
683 break;
684
685 case ID_PARM:
686 {
687 PVMMAPSTRCT pvm = &vrgvmGlo[HINDEX(hvar)];
688 *(PDOUBLE) pvm->pVar = dVal;
689 break;
690 } /* ID_PARM: */
691
692 } /* switch */
693
694 return (bReturn);
695
696 } /* SetVar */
697
698 /* -----------------------------------------------------------------------------
699 SetInput
700
701 Sets the input referenced by hvar to the input function record
702 given. Returns TRUE on success, FALSE if hvar does not reference
703 an input.
704 */
705
SetInput(HVAR hvar,PIFN pifn)706 BOOL SetInput (HVAR hvar, PIFN pifn)
707 {
708 if (!pifn || !IsInput(hvar))
709 return (FALSE);
710
711 memcpy (&vrgInputs[HINDEX(hvar)], pifn, sizeof(IFN));
712 return (TRUE);
713
714 } /* SetVar */
715
716
717 /* -----------------------------------------------------------------------------
718 GetModelVector
719
720 Returns a pointer to the model vector.
721 */
722
GetModelVector(void)723 PDOUBLE GetModelVector (void)
724 {
725 return ((PDOUBLE) vrgModelVars);
726
727 } /* GetModelVector */
728
729 /* -----------------------------------------------------------------------------
730 GetNModelVars
731
732 Returns the number of model variables,
733 i.e. the size of the model vector.
734 */
735
GetNModelVars(void)736 int GetNModelVars (void)
737 {
738 return (vnModelVars);
739
740 } /* GetNModelVars */
741
742 /* -----------------------------------------------------------------------------
743 GetNStates
744
745 Returns the number of states, i.e. the first number of vars in
746 the model vector that are states
747 */
748
GetNStates(void)749 int GetNStates (void)
750 {
751 return (vnStates);
752
753 } /* GetNStates */
754
755
756 /* -----------------------------------------------------------------------------
757 CalcInputs
758
759 Calculates the inputs which change with time at time dTime. This
760 routine assumes that the bExpoOn flag and dTStartPeriod time are
761 correctly set by the input scheduler UpdateInputs().
762
763 The scheduler will have set the flag during exposure periods and
764 also will have set dVal current value to 0.0 for inputs that are
765 off, and to dMag for constant exposure inputs.
766
767 This routine will only calculate "ON" inputs whose magnitude
768 changes with time, for example exponentials, transits and pulses.
769 */
770
771 #define SQRT2PI 2.506628274631000241612
772
CalcInputs(PDOUBLE pdTime)773 void CalcInputs (PDOUBLE pdTime)
774 {
775 int i;
776 double dTmp;
777
778 for (i = 0; i < vnInputs; i++) {
779
780 /* Constants do not change, exposure must be on, or it's a spike
781 (always on) */
782 if ((vrgInputs[i].iType != IFN_CONSTANT) &&
783 (vrgInputs[i].bOn || (vrgInputs[i].iType != IFN_SPIKES))) {
784
785 switch (vrgInputs[i].iType) {
786 case IFN_CONSTANT:
787 break; /* dVal does not change */
788
789 case IFN_NDOSES:
790 if (vrgInputs[i].iDoseCur < vrgInputs[i].nDoses)
791 vrgInputs[i].dVal = vrgInputs[i].rgMags[vrgInputs[i].iDoseCur] *
792 vrgInputs[i].dMag; /* Scaled */
793 break;
794
795 case IFN_PERDOSE: /* Does not change step-wise */
796 break;
797
798 case IFN_PEREXP: /* Calc current exp() */
799 vrgInputs[i].dVal = vrgInputs[i].dMag * ((int) vrgInputs[i].bOn) *
800 exp ((vrgInputs[i].dTStartPeriod +
801 vrgInputs[i].dT0 - *pdTime) *
802 vrgInputs[i].dDecay);
803 break;
804
805 case IFN_PERTRANS: /* Calc current transit input value */
806 /* compute k * t */
807 dTmp = vrgInputs[i].dDecay *
808 (*pdTime - vrgInputs[i].dTStartPeriod - vrgInputs[i].dT0);
809 vrgInputs[i].dVal = vrgInputs[i].dMag * ((int) vrgInputs[i].bOn) *
810 pow(dTmp, vrgInputs[i].dNcpt) * exp(-dTmp) /
811 (SQRT2PI *
812 pow(vrgInputs[i].dNcpt, vrgInputs[i].dNcpt+0.5) *
813 exp(-vrgInputs[i].dNcpt));
814
815 break;
816
817 case IFN_SPIKES:
818 if ((*pdTime == vrgInputs[i].rgT0s[vrgInputs[i].iDoseCur]) &&
819 (vrgInputs[i].iDoseCur < vrgInputs[i].nDoses))
820 vrgInputs[i].dVal = vrgInputs[i].rgMags[vrgInputs[i].iDoseCur] *
821 vrgInputs[i].dMag; /* Scaled */
822
823 else
824 vrgInputs[i].dVal = 0;
825
826 break;
827
828 } /* switch */
829 } /* if */
830 } /* for all inputs */
831
832 } /* CalcInputs */
833
834
835 /* -----------------------------------------------------------------------------
836 DumpSymbolTable
837
838 prints the entire table and values.
839 */
840
DumpSymbolTable(char * szFilename)841 void DumpSymbolTable (char *szFilename)
842 {
843 static char szStderr[] = "<stdout>";
844 FILE *pfile;
845 PVMMAPSTRCT pvm = &vrgvmGlo[0];
846
847 if (szFilename)
848 pfile = fopen(szFilename, "a");
849 else {
850 pfile = stdout;
851 szFilename = szStderr;
852 } /* else */
853
854 if (!pfile) {
855 printf ("Cannot dump symbol table to %s\n", szFilename);
856 return;
857 } /* if */
858
859 fprintf (pfile, "\nSymbol Table:\n");
860 if (!pvm) {
861 fprintf (pfile, "<null>\n");
862 return;
863 } /* if */
864
865 while (*pvm->szName) {
866 fprintf (pfile, "%s \t= ", pvm->szName);
867 if (IsInput(pvm->hvar)) {
868 PIFN pifn = (PIFN) pvm->pVar;
869
870 fprintf (pfile, "Mag=%g [Val=%g]\n", pifn->dMag, pifn->dVal);
871 } /* if */
872 else
873 fprintf (pfile, "%g\n", *(double *) pvm->pVar);
874 pvm++;
875 } /* while */
876
877 if (szFilename != szStderr) fclose (pfile);
878
879 } /* DumpSymbolTable */
880
881
882 /* ----------------------------------------------------------------------------
883 GetStateHandles
884
885 Sets up an array of state variable handles
886 */
887
GetStateHandles(HVAR * phvar)888 void GetStateHandles (HVAR *phvar)
889 {
890 int i = 0;
891
892 VMMAPSTRCT *pvm = vrgvmGlo;
893
894 while (pvm->pVar) {
895 if (IsState(pvm->hvar))
896 phvar[i++] = pvm->hvar;
897 ++pvm;
898 };
899
900 } /* GetStateHandles */
901
902