001: /*
002:  *  dde-exam.c
003:  *  dde
004:  *
005:  *  Created by ashley on 29/11/2009.
006:  *  Copyright 2009 __MyCompanyName__. All rights reserved.
007:  *
008:  */
009: 
010: /* Template file for coding d.d.e. models for solv95 - Example is a dde
011:    model - Nicholson's blowfly model of Gurney And Nisbet 1981
012: 
013: */
014: 
015: #include <math.h>
016: #include <stdio.h>
017: #include <stdlib.h>
018: #include <string.h>
019: #include "solv95.h"
020: 
021: /***************************************************************************/
022: /*  Put global variables  here. These should never be written to from      */
023: /* grad() or switchfunctions(), directly or indirectly.                         */
024: /***************************************************************************/
025: 
026: 
027: /***************************************************************************/
028: /*             Problem specific routines                                                                        */
029: /***************************************************************************/
030: 
031: void initcons(int *no_vars,int *no_cons)
032: /* this routine specifies the number of constants and number of variables in
033:         the model.
034: */
035: { 
036: //  printf("initcons\n");
037:   *no_cons=5;
038:   *no_vars=1;
039: }
040: 
041: void switchfunctions(sw,s,c,t)
042: double *sw,*s,*c,t;
043: 
044: 
045: /* This routine sets the values of the switch functions. When the switch
046:         functions pass through zero from positive to negative the state variables
047:         may be reset in function map(). The switch functions should pass smoothly
048:         through 0 and should never have both value and first derivative zero. The
049:         same switch must not pass through zero from positive to negative more than
050:         once in an integration timestep. An example of a switch function is:
051:                                                 sw[0]=sin(pi*t/30.0)
052:         which passes through zero every 60 time units. Switches may include state
053:         variables provided the above conditions are met. Note that to get 'Solver'
054:         style switches define twice as many switches and let e.g. sw[1]=-sw[0] */
055: 
056: 
057: {
058:   printf("switchfunctions\n");
059: }
060: 
061: 
062: 
063: void map(s,c,t,swno)
064: double *s,*c,t;int swno;
065: 
066: /* This routine is called whenever one of the switch functions passes through
067:         zero. 'swno' is the number of the switch function. The state variables
068:         can be changed discontinuously within this routine. eg:
069:    if (swno==1)
070:           { s[0]=coeff[1]*(s[0]);}
071:         time and the coefficients should not be changed.
072: */
073: 
074: {
075:   printf("map\n");
076: }
077: 
078: 
079: 
080: 
081: // Note that there is an extra argument to grad()!
082: void grad(g,s,c,t,pastvalue)
083: double *g,*s,*c,t;
084: double (*pastvalue)();
085: 
086: /* This routine must provide the gradients g for the state variables s.
087:         So ds[i]/dt=g[i]=fi(s,c,t) where c is the coefficient vector. lagged
088:    variables may be accessed here using pastvalue(i,x,j) which returns the
089:    ith (starting at zero) lagged variable at time x, using lag pointer k
090: 
091:    (lag pointers are used by pastvalue to store the history buffer location
092:     corresponding to a lag in order to save exectution time. For example if
093:     your code requires lagged varaible 0 at lags T and 2T for each gradient
094:     calculation then it is efficient to obtain these values using:
095:     pastvalue(0,t-T,0) and pastvalue(0,t-2*T,1) rather than
096:     pastvalue(0,t-T,0) and pastvalue(0,t-2*T,0). The latter works, it's just
097:     slower because more time is spent searching for lagged values)
098: */
099: 
100: {       double Alag;
101: 
102:         if (t>c[0]) {
103:                 Alag=pastvalue(0,t-c[0],0);
104:         }
105:         else {
106:                 Alag=0.0;
107:         }
108:         
109:         g[0]=c[2]*Alag*exp(-Alag/c[3])-c[1]*s[0];
110: }
111: 
112: void storehistory(his,ghis,g,s,c,t)
113: double *his,*ghis,*g,*s,*c,t;
114: 
115: /* This is the routine in which the values of the history variables at time
116:         t are calculated and put in the array his, along with gradients in ghis,
117:         using state variables s, gradients of s, g, and coefficients c
118:    e.g. if the state variable 2 is history variable 0, you would need the line:
119:    his[0]=s[2];ghis[0]=g[2];
120: */
121: 
122: { 
123: //  printf("storehistory\n");
124:         his[0]=s[0];
125:         ghis[0]=g[0];
126: }
127: 
128: 
129: 
130: void statescale(double *scale)
131: 
132: /* In this routine you can set scale factors for error control. For each
133:    state variable the maximum permisable error will be bounded below by the
134:    tolerance multiplied by scale[i]. If you don't supply values then zero will
135:    be used.
136:         Non-zero scale values are useful for variables that start at zero and
137:    leave zero without 3rd order continuity. */
138: 
139: { 
140: //  printf("statescale\n");
141:   scale[0]=0.0;
142: }
143: 
144: 
145: void initialstate(double *s, double *c, double t)
146: 
147: /* initialise state variables and any global constants here, you can use c */
148: 
149: { 
150: //   printf("initstate\n");
151:    s[0]=c[4];
152: 
153: }
154: 
155: void initout(usercontrol *out)
156: 
157: /* This routine is where the output windows are first set up. You have to
158:    fill in the global structure "out", which tells solv95 which
159:    variables to plot and in which windows */
160: 
161: { // int i;
162: 
163: //  printf("initout\n");
164:         g_printf("%08x %08x\n", (int)out, (int)(out->index));
165:   
166:  /* how many windows */
167:         out->no_windows=1;
168:         
169: /* how many lines in each window e.g. for two lines in window zero set
170:      out->lines[0]=2 ..... */
171:         out->lines[0]=1;
172:         
173: /* which state variables go in which windows and what are they called?
174:     e.g. to put state variable 4 in window 3 as curve 0, and label it `trash':
175:     out->index[4].win=3;out->index[4].cur=0;out->label[4]="trash";
176:     If you don't provide this information for a variable then it isn't plotted
177:  */
178:         out->label[0]="Adults";
179:         out->index[0].win=0;
180:         out->index[0].cur=0;
181:         
182:         /* labels for windows */
183:         out->wname[0]="Blowfly population";
184:         
185:         /* Now set up initial range of y axis variables for each window */
186:         out->range[0].y0=0.0;
187:         out->range[0].y1=-5000.0;
188:         
189:         out->xlabel="Time";
190: /* initialise file output. Set out->fileno to the number of columns to output
191:     (in addition to time). Set out->fout[i] to the number of the state variable
192:     that is to go in column i. Set out->fileno=0 for no output. */
193:         out->fileno=1;
194:         out->fout[0]=0;
195:         
196: /* integration details*/
197:         out->t0=0.0;        /* default start time */
198:         out->t1=300.0;      /* default stop time */
199:         out->dt=1.0;        /* initial timestep */
200:         out->tol=0.000005;  /* integration tolerance */
201:         out->dout = 1.0;    /* approximate average output timestep */
202:         out->hbsize=2000L;  /* how many past values to store for each history variable */
203:         out->nhv=1;         /* Number of history (lagged) variables */
204:         out->nlag=1;        /* Number of lag markers per history variable (set to 1 if unsure)*/
205:         out->nsw=0;         /* number of switch varaibles */
206:         
207: /* Initial Values and names for constants to be prompted for at run time */
208:         
209:         out->c[0]=12.0; out->cname[0]="tau";
210:         out->c[1]=0.25; out->cname[1]="delta";
211:         out->c[2]=10.0; out->cname[2]="P";
212:         out->c[3]=300.0;out->cname[3]="A_0";
213:         out->c[4]=100.0;out->cname[4]="N(0)";
214: /* Optional information text on parameters - split text between lines with `\'*/
215:         
216:         out->cinfo[0]="Development time.";
217:         out->cinfo[1]="Per capita death rate";
218:         out->cinfo[2]="Product of max. fecundity and juvenile survival.";
219:         out->cinfo[3]="Fecundity decay constant.";
220:         out->cinfo[4]="Initial population.";
221:         
222: /* Text and title for initial welcome box (optional) */
223:         out->initialtitle="DDE Model - Blowflies";
224:         out->initialtext=" This model is an example that came with solv95 ";
225:         
226: }
227: 
228: 
229: