001: /*
002:  *  rinaldi.c
003:  *  odemodels.cgi
004:  *
005:  *  Created by ashley on 09/03/2010.
006:  *  Copyright 2010 __MyCompanyName__. All rights reserved.
007:  *
008:  */
009: 
010: #include <math.h>
011: #include <stdio.h>
012: #include <stdlib.h>
013: #include <strings.h>
014: #include "control.h"
015: #include "rinaldi.h"
016: 
017: void fcnrinaldi(int n, double t, double *y, double *yprime, void *fdata);
018: 
019: control *makerinaldi()
020: {
021:         int nconstants = 13;
022:         int nvariables = 3;
023:         int nplots = 4;
024:         control *c = makecontrol(nvariables, nconstants, nplots);
025:         
026:         c->title = "Love Dynamics";
027:         c->description = "Rinaldi's model of the love between Petrarch & Laura";
028: 
029:         c->integrator = dverkintegrator;
030:         c->fcn = fcnrinaldi;
031: 
032:         int p = 0;
033:         c->parameters[p]->value = 3.0;
034:         c->parameters[p]->name = astrcpy("alpha1");
035:         c->parameters[p]->info = astrcpy("alpha1 parameter");
036:         p++;
037:         
038:         c->parameters[p]->value = 1.0;
039:         c->parameters[p]->name = astrcpy("alpha2");
040:         c->parameters[p]->info = astrcpy("alpha2 parameter");
041:         p++;
042:         
043:         c->parameters[p]->value = 0.1;
044:         c->parameters[p]->name = astrcpy("alpha3");
045:         c->parameters[p]->info = astrcpy("alpha3 parameter");
046:         p++;
047:         
048:         c->parameters[p]->value = 1.0;
049:         c->parameters[p]->name = astrcpy("beta1");
050:         c->parameters[p]->info = astrcpy("beta1 parameter");
051:         p++;
052:         
053:         c->parameters[p]->value = 5.0;
054:         c->parameters[p]->name = astrcpy("beta2");
055:         c->parameters[p]->info = astrcpy("beta2 parameter");
056:         p++;
057:         
058:         c->parameters[p]->value = 0.1e+02;
059:         c->parameters[p]->name = astrcpy("beta3");
060:         c->parameters[p]->info = astrcpy("beta3 parameter");
061:         p++;
062:         
063:         c->parameters[p]->value = 1.0;
064:         c->parameters[p]->name = astrcpy("gamma");
065:         c->parameters[p]->info = astrcpy("gamma parameter");
066:         p++;
067:         
068:         c->parameters[p]->value = 1.0;
069:         c->parameters[p]->name = astrcpy("delta");
070:         c->parameters[p]->info = astrcpy("delta parameter");
071:         p++;
072:         
073:         c->parameters[p]->value = 2.0;
074:         c->parameters[p]->name = astrcpy("A_L");
075:         c->parameters[p]->info = astrcpy("A_L parameter");
076:         p++;
077:         
078:         c->parameters[p]->value = -1.0;
079:         c->parameters[p]->name = astrcpy("A_P");
080:         c->parameters[p]->info = astrcpy("A_P parameter");
081:         p++;
082:         
083:         c->parameters[p]->value = 0.0;
084:         c->parameters[p]->name = astrcpy("L(0)");
085:         c->parameters[p]->info = astrcpy("Initial L");
086:         p++;
087:         
088:         c->parameters[p]->value = 0.0;
089:         c->parameters[p]->name = astrcpy("P(0)");
090:         c->parameters[p]->info = astrcpy("Initial P");
091:         p++;
092:         
093:         c->parameters[p]->value = 0.0;
094:         c->parameters[p]->name = astrcpy("Z(0)");
095:         c->parameters[p]->info = astrcpy("Initial Z");
096:         p++;
097:         
098:         c->variables[0]->name = astrcpy("L");
099:         c->variables[1]->name = astrcpy("P");
100:         c->variables[2]->name = astrcpy("Z");
101:         c->variables[0]->value = c->parameters[10]->value;
102:         c->variables[1]->value = c->parameters[11]->value;
103:         c->variables[2]->value = c->parameters[12]->value;
104:         
105:         int lines[] = {1, 2, 3};
106:         c->plots[0] = makeplot(lineplot, 3, 0, lines);
107:         c->plots[0]->title = astrcpy("Time series");
108:         c->plots[0]->label = astrcpy("x-y-z Time Series Plot");
109:         c->plots[0]->xlabel = astrcpy("y");
110:         c->plots[0]->ylabel = astrcpy("t");
111: 
112:         int lines2[] = {3};
113:         c->plots[1] = makeplot(lineplot, 1, 2, lines2);
114:         c->plots[1]->title = astrcpy("Phase Space");
115:         c->plots[1]->label = astrcpy("z-y Phase Space Plot");
116:         c->plots[1]->xlabel = astrcpy("w");
117:         c->plots[1]->ylabel = astrcpy("y");
118: 
119:         int lines3[] = {2};
120:         c->plots[2] = makeplot(lineplot, 1, 1, lines3);
121:         c->plots[2]->title = astrcpy("Phase Space");
122:         c->plots[2]->label = astrcpy("x-y Phase Space Plot");
123:         c->plots[2]->xlabel = astrcpy("x");
124:         c->plots[2]->ylabel = astrcpy("y");
125:         
126:         int lines4[] = {3};
127:         c->plots[3] = makeplot(lineplot, 1, 1, lines4);
128:         c->plots[3]->title = astrcpy("Phase Space");
129:         c->plots[3]->label = astrcpy("z-w Phase Space Plot");
130:         c->plots[3]->xlabel = astrcpy("w");
131:         c->plots[3]->ylabel = astrcpy("y");
132:         
133:         c->t0 = 0.0;
134:         c->t1 = 25.0;
135:         c->dt = 0.05;
136:         c->tol = 1.0e-07;
137: 
138:         
139:         return c;
140: }
141: 
142: void fcnrinaldi(int n, double t, double *y, double *yprime, void *fdata)
143: {
144:         control *c = (control *)fdata;
145:         constant **params = c->parameters;
146:         double alpha1 = params[0]->value;
147:         double alpha2 = params[1]->value;
148:         double alpha3 = params[2]->value;
149: 
150:         double beta1 = params[3]->value;
151:         double beta2 = params[4]->value;
152:         double beta3 = params[5]->value;
153: 
154:         double gamma = params[6]->value;
155:         double delta = params[7]->value;
156: 
157:         double A_L = params[8]->value;
158:         double A_P = params[9]->value;
159: 
160:         double L = y[0];
161:         double P = y[1];
162:         double Z = y[2];
163:         
164:         double pp = P/gamma;
165:         yprime[0] = -alpha1*L + beta1*(P*(1 - pp*pp) + A_P);
166:         yprime[1] = -alpha2*P + beta2*(L + A_L/(1 + delta*Z));
167:         yprime[2] = -alpha3*Z + beta3*P;
168: }
169: