001: /*
002:  *  dynamo.c
003:  *  dynamo.cgi
004:  *
005:  *  Created by ashley on 18/05/2009.
006:  *  Copyright 2009 __MyCompanyName__. All rights reserved.
007:  *
008:  */
009: #include <stdio.h>
010: #include <stdlib.h>
011: #include <math.h>
012: #include <strings.h>
013: #include "control.h"
014: #include "dynamo.h"
015: 
016: //void initialise(int *nconstants, int *nvariables, int *nplots);
017: // void initcontrol(control *c);
018: void fcndynamo(int n, double t, double *y, double *yprime, void *fdata);
019: 
020: /*void initialise(int *nconstants, int *nvariables, int *nplots)
021: {
022:         *nconstants = 4;
023:         *nvariables = 3;
024:         *nplots = 4;
025: } */
026: 
027: control *makedynamo()
028: {
029:         int nconstants = 4;
030:         int nvariables = 3;
031:         int nplots = 4;
032:         control *c = makecontrol(nvariables, nconstants, nplots);
033:         
034:         c->title = "Double Disc Dynamo ODE model";
035:         c->description = "Adapted from the gnu graph 'ode' package.";
036:         
037:         c->integrator = dverkintegrator;
038: //      c->init = initcontrol;
039:         c->fcn = fcndynamo;
040:         
041:         c->parameters[0]->value = 1.0;
042:         c->parameters[0]->name = astrcpy("A");
043:         c->parameters[0]->info = astrcpy("A parameter");
044:         
045:         c->parameters[1]->value = 1.0;
046:         c->parameters[1]->name = astrcpy("V");
047:         c->parameters[1]->info = astrcpy("V parameter");
048:         
049:         c->parameters[2]->value = 14.625;
050:         c->parameters[2]->name = astrcpy("Q");
051:         c->parameters[2]->info = astrcpy("Q parameter");
052:         
053:         c->parameters[3]->value = 5.0;
054:         c->parameters[3]->name = astrcpy("S");
055:         c->parameters[3]->info = astrcpy("S parameter");
056: 
057:         c->variables[0]->name = astrcpy("w");
058:         c->variables[1]->name = astrcpy("y");
059:         c->variables[2]->name = astrcpy("z");
060:         c->variables[0]->value = 1.0;
061:         c->variables[1]->value = 1.0;
062:         c->variables[2]->value = 1.0;
063:         
064:         int lines[] = {1, 2, 3};
065:         c->plots[0] = makeplot(lineplot, 3, 0, lines);
066:         c->plots[0]->title = astrcpy("Time series");
067:         c->plots[0]->label = astrcpy("w-y-z Time Series Plot");
068:         c->plots[0]->xlabel = astrcpy("y");
069:         c->plots[0]->ylabel = astrcpy("t");
070: 
071:         int lines2[] = {3};
072:         c->plots[1] = makeplot(lineplot, 1, 2, lines2);
073:         c->plots[1]->title = astrcpy("Phase Space");
074:         c->plots[1]->label = astrcpy("z-y Phase Space Plot");
075:         c->plots[1]->xlabel = astrcpy("w");
076:         c->plots[1]->ylabel = astrcpy("y");
077: 
078:         int lines3[] = {2};
079:         c->plots[2] = makeplot(lineplot, 1, 1, lines3);
080:         c->plots[2]->title = astrcpy("Phase Space");
081:         c->plots[2]->label = astrcpy("w-y Phase Space Plot");
082:         c->plots[2]->xlabel = astrcpy("w");
083:         c->plots[2]->ylabel = astrcpy("y");
084:         
085:         int lines4[] = {3};
086:         c->plots[3] = makeplot(lineplot, 1, 1, lines4);
087:         c->plots[3]->title = astrcpy("Phase Space");
088:         c->plots[3]->label = astrcpy("z-w Phase Space Plot");
089:         c->plots[3]->xlabel = astrcpy("w");
090:         c->plots[3]->ylabel = astrcpy("y");
091:         
092:         c->t0 = 0.0;
093:         c->t1 = 25.0;
094:         c->dt = 0.05;
095:         c->tol = 1.0e-07;
096:         
097:         return c;
098: }
099: 
100: /*void initcontrol(control *c)
101: {
102:         c->title = "Double Disc Dynamo ODE model";
103:         c->description = "Adapted from the gnu graph 'ode' package.";
104:         
105:         c->integrator = dverkintegrator;
106: 
107:         c->parametervalues[0] = 1.0;
108:         c->parameternames[0] = astrcpy("A");
109:         c->parameterinfo[0] = astrcpy("A parameter");
110:         
111:         c->parametervalues[1] = 1.0;
112:         c->parameternames[1] = astrcpy("V");
113:         c->parameterinfo[1] = astrcpy("V parameter");
114:         
115:         c->parametervalues[2] = 14.625;
116:         c->parameternames[2] = astrcpy("Q");
117:         c->parameterinfo[2] = astrcpy("Q parameter");
118:         
119:         c->parametervalues[3] = 5.0;
120:         c->parameternames[3] = astrcpy("S");
121:         c->parameterinfo[3] = astrcpy("S parameter");
122: 
123:         c->variablenames[0] = astrcpy("w");
124:         c->variablenames[1] = astrcpy("y");
125:         c->variablenames[2] = astrcpy("z");
126:         c->variablevalues[0] = 1.0;
127:         c->variablevalues[1] = 1.0;
128:         c->variablevalues[2] = 1.0;
129:         
130:         int lines[] = {1, 2, 3};
131:         c->plots[0] = makeplot(lineplot, 3, 0, lines);
132:         c->plots[0]->title = astrcpy("Time series");
133:         c->plots[0]->label = astrcpy("w-y-z Time Series Plot");
134:         c->plots[0]->xlabel = astrcpy("y");
135:         c->plots[0]->ylabel = astrcpy("t");
136: 
137:         int lines2[] = {3};
138:         c->plots[1] = makeplot(lineplot, 1, 2, lines2);
139:         c->plots[1]->title = astrcpy("Phase Space");
140:         c->plots[1]->label = astrcpy("z-y Phase Space Plot");
141:         c->plots[1]->xlabel = astrcpy("w");
142:         c->plots[1]->ylabel = astrcpy("y");
143: 
144:         int lines3[] = {2};
145:         c->plots[2] = makeplot(lineplot, 1, 1, lines3);
146:         c->plots[2]->title = astrcpy("Phase Space");
147:         c->plots[2]->label = astrcpy("w-y Phase Space Plot");
148:         c->plots[2]->xlabel = astrcpy("w");
149:         c->plots[2]->ylabel = astrcpy("y");
150:         
151:         int lines4[] = {3};
152:         c->plots[3] = makeplot(lineplot, 1, 1, lines4);
153:         c->plots[3]->title = astrcpy("Phase Space");
154:         c->plots[3]->label = astrcpy("z-w Phase Space Plot");
155:         c->plots[3]->xlabel = astrcpy("w");
156:         c->plots[3]->ylabel = astrcpy("y");
157:         
158:         c->t0 = 0.0;
159:         c->t1 = 25.0;
160:         c->dt = 0.05;
161:         c->tol = 1.0e-07;
162: } */
163: 
164: void fcndynamo(int n, double t, double *y, double *yprime, void *fdata)
165: {
166: //      yprime[0] = -y[1] - 0.0*y[0];
167: //      yprime[1] = y[0];
168: //      yprime[0] = -y[0];
169:         control *c = (control *)fdata;
170:         constant **params = c->parameters;
171:         double A = params[0]->value;
172:         double V = params[1]->value;
173:         double Q = params[2]->value;
174:         double S = params[3]->value;
175: 
176:         yprime[0] = Q - y[2] * y[1] - V * y[0];                                                                 // w = Q - z * y - V * w
177:         yprime[1] = S * ( A * y[2] - y[1]);                                                                             // y = S * ( A * z - y)
178:         yprime[2] = y[0] * y[1] - y[2];                                                                                 // z = w * y - z
179: }
180: