Richard Boegli's CnC_Generals_Zero_Hour Fork WIP
This is documentation of Richard Boegil's Zero Hour Fork
 
Loading...
Searching...
No Matches
ode.cpp
Go to the documentation of this file.
1/*
2** Command & Conquer Generals Zero Hour(tm)
3** Copyright 2025 Electronic Arts Inc.
4**
5** This program is free software: you can redistribute it and/or modify
6** it under the terms of the GNU General Public License as published by
7** the Free Software Foundation, either version 3 of the License, or
8** (at your option) any later version.
9**
10** This program is distributed in the hope that it will be useful,
11** but WITHOUT ANY WARRANTY; without even the implied warranty of
12** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13** GNU General Public License for more details.
14**
15** You should have received a copy of the GNU General Public License
16** along with this program. If not, see <http://www.gnu.org/licenses/>.
17*/
18
19/* $Header: /Commando/Code/wwmath/ODE.CPP 8 7/02/99 10:32a Greg_h $ */
20/***********************************************************************************************
21 *** Confidential - Westwood Studios ***
22 ***********************************************************************************************
23 * *
24 * Project Name : Commando *
25 * *
26 * $Archive:: /Commando/Code/wwmath/ODE.CPP $*
27 * *
28 * Author:: Greg_h *
29 * *
30 * $Modtime:: 6/25/99 6:23p $*
31 * *
32 * $Revision:: 8 $*
33 * *
34 *---------------------------------------------------------------------------------------------*
35 * Functions: *
36 * Euler_Integrate -- uses Eulers method to integrate a system of ODE's *
37 * Midpoint_Integrate -- midpoint method (Runge-Kutta 2) for integration *
38 * Runge_Kutta_Integrate -- Runge Kutta 4 method *
39 * Runge_Kutta5_Integrate -- 5th order Runge-Kutta *
40 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41
42#include "ode.h"
43#include <assert.h>
44
45static StateVectorClass Y0;
46static StateVectorClass Y1;
47static StateVectorClass _WorkVector0;
48static StateVectorClass _WorkVector1;
49static StateVectorClass _WorkVector2;
50static StateVectorClass _WorkVector3;
51static StateVectorClass _WorkVector4;
52static StateVectorClass _WorkVector5;
53static StateVectorClass _WorkVector6;
54static StateVectorClass _WorkVector7;
55
56/***********************************************************************************************
57 * Euler_Solve -- uses Eulers method to integrate a system of ODE's *
58 * *
59 * INPUT: *
60 * odesys - pointer to the ODE system to integrate *
61 * dt - size of the timestep *
62 * *
63 * OUTPUT: *
64 * state vector in odesys will be updated for the next timestep *
65 * *
66 * WARNINGS: *
67 * *
68 * HISTORY: *
69 * 08/11/1997 GH : Created. *
70 * 6/25/99 GTH : Updated to the new integrator system *
71 *=============================================================================================*/
72void IntegrationSystem::Euler_Integrate(ODESystemClass * sys, float dt)
73{
74 WWASSERT(sys != NULL);
75
76 /*
77 ** Get the current state
78 */
79 Y0.Reset();
80 sys->Get_State(Y0);
81 Y1.Resize(Y0.Count());
82
83 /*
84 ** make aliases to the work-vectors we need
85 */
86 StateVectorClass & dydt = _WorkVector0;
87 dydt.Resize(Y0.Count());
88
89 /*
90 ** Euler method, just evaluate the derivative, multiply
91 ** by the time-step and add to the current state vector.
92 */
93 sys->Compute_Derivatives(0,NULL,&dydt);
94
95 Y1.Resize(Y0.Count());
96 for (int i = 0; i < Y0.Count(); i++) {
97 Y1[i] = Y0[i] + dydt[i] * dt;
98 }
99
100 sys->Set_State(Y1);
101}
102
103/***********************************************************************************************
104 * Midpoint_Integrate -- midpoint method (Runge-Kutta 2) *
105 * *
106 * INPUT: *
107 * sys - pointer to the ODE system to integrate *
108 * dt - size of the timestep *
109 * *
110 * OUTPUT: *
111 * state vector in odesys will be updated for the next timestep *
112 * *
113 * WARNINGS: *
114 * *
115 * HISTORY: *
116 * 08/11/1997 GH : Created. *
117 * 6/25/99 GTH : Updated to the new integrator system *
118 *=============================================================================================*/
119void IntegrationSystem::Midpoint_Integrate(ODESystemClass * sys,float dt)
120{
121 int i;
122
123 /*
124 ** Get the current state
125 */
126 Y0.Reset();
127 sys->Get_State(Y0);
128 Y1.Resize(Y0.Count());
129
130 /*
131 ** make aliases to the work-vectors we need
132 */
133 StateVectorClass & dydt = _WorkVector0;
134 StateVectorClass & ymid = _WorkVector1;
135 dydt.Resize(Y0.Count());
136 ymid.Resize(Y0.Count());
137
138 /*
139 ** MidPoint method, first evaluate the derivitives of the
140 ** state vector just like the Euler method.
141 */
142 sys->Compute_Derivatives(0.0f,NULL,&dydt);
143
144 /*
145 ** Compute the midpoint between the Euler solution and
146 ** the input values.
147 */
148 for (i=0; i<Y0.Count(); i++) {
149 ymid[i] = Y0[i] + dt * dydt[i] / 2.0f;
150 }
151
152 /*
153 ** Re-compute derivatives at this point.
154 */
155 sys->Compute_Derivatives(dt/2.0f,&ymid,&dydt);
156
157 /*
158 ** Use these derivatives to compute the solution.
159 */
160 for (i=0; i<Y0.Count(); i++) {
161 Y1[i] = Y0[i] + dt * dydt[i];
162 }
163
164 sys->Set_State(Y1);
165}
166
167
168/***********************************************************************************************
169 * Runge_Kutta_Integrate -- Runge Kutta 4 method *
170 * *
171 * INPUT: *
172 * odesys - pointer to the ODE system to integrate *
173 * dt - size of the timestep *
174 * *
175 * OUTPUT: *
176 * state vector in odesys will be updated for the next timestep *
177 * *
178 * WARNINGS: *
179 * *
180 * HISTORY: *
181 * 08/11/1997 GH : Created. *
182 *=============================================================================================*/
183void IntegrationSystem::Runge_Kutta_Integrate(ODESystemClass * sys,float dt)
184{
185 int i;
186 float dt2 = dt / 2.0f;
187 float dt6 = dt / 6.0f;
188
189 /*
190 ** Get the current state
191 */
192 Y0.Reset();
193 sys->Get_State(Y0);
194 Y1.Resize(Y0.Count());
195
196 /*
197 ** make aliases to the work-vectors we need
198 */
199 StateVectorClass & dydt = _WorkVector0;
200 StateVectorClass & dym = _WorkVector1;
201 StateVectorClass & dyt = _WorkVector2;
202 StateVectorClass & yt = _WorkVector3;
203 dydt.Resize(Y0.Count());
204 dym.Resize(Y0.Count());
205 dyt.Resize(Y0.Count());
206 yt.Resize(Y0.Count());
207
208 /*
209 ** First Step
210 */
211 sys->Compute_Derivatives(0.0f,NULL,&dydt);
212 for (i=0; i<Y0.Count(); i++) {
213 yt[i] = Y0[i] + dt2 * dydt[i];
214 }
215
216 /*
217 ** Second Step
218 */
219 sys->Compute_Derivatives(dt2, &yt, &dyt);
220 for (i=0; i<Y0.Count(); i++) {
221 yt[i] = Y0[i] + dt2 * dyt[i];
222 }
223
224 /*
225 ** Third Step
226 */
227 sys->Compute_Derivatives(dt2, &yt, &dym);
228 for (i=0; i<Y0.Count(); i++) {
229 yt[i] = Y0[i] + dt*dym[i];
230 dym[i] += dyt[i];
231 }
232
233 /*
234 ** Fourth Step
235 */
236 sys->Compute_Derivatives(dt, &yt, &dyt);
237 for (i=0; i<Y0.Count(); i++) {
238 Y1[i] = Y0[i] + dt6 * (dydt[i] + dyt[i] + 2.0f*dym[i]);
239 }
240
241 sys->Set_State(Y1);
242}
243
244/***********************************************************************************************
245 * Runge_Kutta5_Integrate -- 5th order Runge-Kutta *
246 * *
247 * INPUT: *
248 * odesys - pointer to the ODE system to integrate *
249 * dt - size of the timestep *
250 * *
251 * OUTPUT: *
252 * state vector in odesys will be updated for the next timestep *
253 * *
254 * WARNINGS: *
255 * *
256 * HISTORY: *
257 * 08/11/1997 GH : Created. *
258 * 6/25/99 GTH : Converted to the new Integrator system *
259 *=============================================================================================*/
260void IntegrationSystem::Runge_Kutta5_Integrate(ODESystemClass * odesys,float dt)
261{
262 int i;
263 int veclen;
264 static const float a2 = 0.2f;
265 static const float a3 = 0.3f;
266 static const float a4 = 0.6f;
267 static const float a5 = 1.0f;
268 static const float a6 = 0.875f;
269 static const float b21 = 0.2f;
270 static const float b31 = 3.0f/40.0f;
271 static const float b32 = 9.0f/40.0f;
272 static const float b41 = 0.3f;
273 static const float b42 = -0.9f;
274 static const float b43 = 1.2f;
275 static const float b51 = -11.0f /54.0f;
276 static const float b52 = 2.5f;
277 static const float b53 = -70.0f/27.0f;
278 static const float b54 = 35.0f/27.0f;
279 static const float b61 = 1631.0f/55296.0f;
280 static const float b62 = 175.0f/512.0f;
281 static const float b63 = 575.0f/13824.0f;
282 static const float b64 = 44275.0f/110592.0f;
283 static const float b65 = 253.0f/4096.0f;
284 static const float c1 = 37.0f/378.0f;
285 static const float c3 = 250.0f/621.0f;
286 static const float c4 = 125.0f/594.0f;
287 static const float c6 = 512.0f/1771.0f;
288 static const float dc5 = -277.0f/14336.0f;
289 static const float dc1 = c1 - 2825.0f/27648.0f;
290 static const float dc3 = c3 - 18575.0f/48384.0f;
291 static const float dc4 = c4 - 13525.0f/55296.0f;
292 static const float dc6 = c6 - 0.25f;
293
294 /*
295 ** Get the current state
296 */
297 Y0.Reset();
298 odesys->Get_State(Y0);
299 veclen = Y0.Count();
300 Y1.Resize(veclen);
301
302 /*
303 ** make aliases to the work-vectors we need
304 */
305 StateVectorClass & dydt = _WorkVector0;
306 StateVectorClass & ak2 = _WorkVector1;
307 StateVectorClass & ak3 = _WorkVector2;
308 StateVectorClass & ak4 = _WorkVector3;
309 StateVectorClass & ak5 = _WorkVector4;
310 StateVectorClass & ak6 = _WorkVector5;
311 StateVectorClass & ytmp = _WorkVector6;
312 StateVectorClass & yerr = _WorkVector7;
313
314 dydt.Resize(veclen);
315 ak2.Resize(veclen);
316 ak3.Resize(veclen);
317 ak4.Resize(veclen);
318 ak5.Resize(veclen);
319 ak6.Resize(veclen);
320 ytmp.Resize(veclen);
321 yerr.Resize(veclen);
322
323 // First step
324 odesys->Compute_Derivatives(0.0f,NULL,&dydt);
325 for (i=0;i<veclen;i++) {
326 ytmp[i] = Y0[i] + b21*dt*dydt[i];
327 }
328
329 // Second step
330 odesys->Compute_Derivatives(a2*dt, &ytmp, &ak2);
331 for (i=0; i<veclen; i++) {
332 ytmp[i] = Y0[i] + dt*(b31*dydt[i] + b32*ak2[i]);
333 }
334
335 // Third step
336 odesys->Compute_Derivatives(a3*dt, &ytmp, &ak3);
337 for (i=0; i<veclen; i++) {
338 ytmp[i] = Y0[i] + dt*(b41*dydt[i] + b42*ak2[i] + b43*ak3[i]);
339 }
340
341 // Fourth step
342 odesys->Compute_Derivatives(a4*dt, &ytmp, &ak4);
343 for (i=0; i<veclen; i++) {
344 ytmp[i] = Y0[i] + dt*(b51*dydt[i] + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]);
345 }
346
347 // Fifth step
348 odesys->Compute_Derivatives(a5*dt, &ytmp, &ak5);
349 for (i=0; i<veclen; i++) {
350 ytmp[i] = Y0[i] + dt*(b61*dydt[i] + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]);
351 }
352
353 // Sixth step
354 odesys->Compute_Derivatives(a6*dt, &ytmp, &ak6);
355 for (i=0; i<veclen; i++) {
356 Y1[i] = Y0[i] + dt*(c1*dydt[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]);
357 }
358
359 // Error approximation!
360 // (maybe I should use this someday? nah not going to use this integrator anyway...)
361 for (i=0; i<veclen; i++) {
362 yerr[i] = dt*(dc1*dydt[i] + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i]);
363 }
364
365 odesys->Set_State(Y1);
366}
367
NULL
#define NULL
Definition BaseType.h:92
WWASSERT
#define WWASSERT
Definition aabtreebuilder.cpp:68
IntegrationSystem::Runge_Kutta5_Integrate
static void Runge_Kutta5_Integrate(ODESystemClass *odesys, float dt)
Definition ode.cpp:260
IntegrationSystem::Midpoint_Integrate
static void Midpoint_Integrate(ODESystemClass *sys, float dt)
Definition ode.cpp:119
IntegrationSystem::Euler_Integrate
static void Euler_Integrate(ODESystemClass *sys, float dt)
Definition ode.cpp:72
IntegrationSystem::Runge_Kutta_Integrate
static void Runge_Kutta_Integrate(ODESystemClass *sys, float dt)
Definition ode.cpp:183
ODESystemClass::Get_State
virtual void Get_State(StateVectorClass &set_state)=0
ODESystemClass::Set_State
virtual int Set_State(const StateVectorClass &new_state, int start_index=0)=0
ODESystemClass::Compute_Derivatives
virtual int Compute_Derivatives(float t, StateVectorClass *test_state, StateVectorClass *dydt, int start_index=0)=0
StateVectorClass::Resize
void Resize(int size)
Definition ode.h:63