1 /*
2 * - - - - - - - - - - -
3 * g a l _ r k f s 7 8
4 * - - - - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * This routine takes a Runge-Kutte-Fehlberg 7(8) step
11 *
12 * Status:
13 *
14 * Math special function.
15 *
16 * Given:
17 *
18 * y d[] dependent variable vector
19 * dydx d[] derivative of dependent variable vector
20 * n i Number of equations to integrate
21 * x d Independent variable value
22 * h d Stepsize
23 * *derivs func User defined function for calculating the right hand side derivatives
24 * *derivsp i Pointer to parameters structure for derivs routine
25 *
26 * Returned:
27 *
28 * yout d[] Ending y values
29 * yerr d[] Errors
30 * gal_rkfs45 i Status
31 * 0 = Success
32 * 1 = Failed to allocate memory
33 *
34 * References:
35 *
36 * NASA Technical Report TR R-352
37 * Some Experimental Results Concerning The Error Propagation in
38 * Runge-Kutte type integration formulas
39 * by Erwin Fehlberg
40 * October 1970
41 *
42 * This revision:
43 *
44 * 2008 May 18
45 *
46 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
47 *
48 *-----------------------------------------------------------------------
49 */
50
51 #include <stdlib.h>
52 #include "gal_const.h"
53 #include "gal_rkfs78.h"
54
gal_rkfs78(double y[],double dydx[],int n,double x,double h,double yout[],double yerr[],void (* derivs)(double,double[],double[],int *),int * derivsp)55 int gal_rkfs78
56 (
57 double y[],
58 double dydx[],
59 int n,
60 double x,
61 double h,
62 double yout[],
63 double yerr[],
64 void ( *derivs ) ( double, double [], double [], int * ),
65 int *derivsp
66 )
67
68 {
69
70 struct {
71 double a[13] ;
72 double b[13][12] ;
73 double c[13] ;
74 double chat[13] ;
75 } param = {
76
77 {
78 0.0,
79 2.0 / 27.0,
80 1.0 / 9.0,
81 1.0 / 6.0,
82 5.0 / 12.0,
83 1.0 / 2.0,
84 5.0 / 6.0,
85 1.0 / 6.0,
86 2.0 / 3.0,
87 1.0 / 3.0,
88 1.0,
89 0.0,
90 1.0,
91 } ,
92
93 {
94 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
95 { 2.0 / 27.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
96 { 1.0 / 36.0, 1.0 / 12.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
97 { 1.0 / 24.0, 0.0, 1.0 / 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
98 { 5.0 / 12.0, 0.0, -25.0 / 16.0, 25.0 / 16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
99 { 1.0 / 20.0, 0.0, 0.0, 1.0 / 4.0, 1.0 / 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
100 { -25.0 / 108.0, 0.0, 0.0, 125.0 / 108.0, -65.0 / 27.0, 125.0 / 54.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
101 { 31.0 / 300.0, 0.0, 0.0, 0.0, 61.0 / 225.0, -2.0 / 9.0, 13.0 / 900.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ,
102 { 2.0, 0.0, 0.0, -53.0 / 6.0, 704.0 / 45.0, -107.0 / 9.0, 67.0 / 90.0, 3.0, 0.0, 0.0, 0.0, 0.0 } ,
103 { -91.0 / 108.0, 0.0, 0.0, 23.0 / 108.0, -976.0 / 135.0, 311.0 / 54.0, -19.0 / 60.0, 17.0 / 6.0, -1.0 / 12.0, 0.0, 0.0, 0.0 } ,
104 { 2383.0 / 4100.0, 0.0, 0.0, -341.0 / 164.0, 4496.0 / 1025.0, -301.0 / 82.0, 2133.0 / 4100, 45.0 / 82.0, 45.0 / 164.0, 18.0 / 41.0, 0.0, 0.0 } ,
105 { 3.0 / 205.0, 0.0, 0.0, 0.0, 0.0, -6.0 / 41.0, -3.0 / 205.0, -3.0 / 41.0, 3.0 / 41.0, 6.0 / 41.0, 0.0, 0.0 } ,
106 { -1777.0 / 4100.0, 0.0, 0.0, -341.0 / 164.0, 4496.0 / 1025.0, -289.0 / 82.0, 2193.0 / 4100.0, 51.0 / 82.0, 33.0 / 164.0, 12.0 / 41.0, 0.0, 1.0 } ,
107 } ,
108
109 {
110 41.0 / 840,
111 0.0,
112 0.0,
113 0.0,
114 0.0,
115 34.0 / 105.0,
116 9.0 / 35.0,
117 9.0 / 35.0,
118 9.0 / 280.0,
119 9.0 / 280.0,
120 41.0 / 840.0,
121 0.0,
122 0.0,
123 } ,
124
125 {
126 0.0,
127 0.0,
128 0.0,
129 0.0,
130 0.0,
131 34.0 / 105.0,
132 9.0 / 35.0,
133 9.0 / 35.0,
134 9.0 / 280.0,
135 9.0 / 280.0,
136 0.0,
137 41.0 / 840.0,
138 41.0 / 840.0,
139 }
140
141 } ;
142
143 #define A(i) param.a[i]
144 #define B(i,j) param.b[i][j]
145 #define C(i) param.c[i]
146 #define CHAT param.chat[i]
147 #define DC(i) (param.c[i] - param.chat[i])
148
149 int i ;
150 double *ak2, *ak3, *ak4, *ak5, *ak6, *ak7, *ak8, *ak9, *ak10, *ak11, *ak12, *ak13, *ytemp ;
151
152 /*
153 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154 */
155
156 /*
157 * Allocate Workspace
158 */
159
160 ak2 = malloc ( n * sizeof ( double ) ) ;
161 if ( ak2 == NULL ) {
162 return 1 ;
163 }
164 ak3 = malloc ( n * sizeof ( double ) ) ;
165 if ( ak3 == NULL ) {
166 free ( ak2 ) ;
167 return 1 ;
168 }
169 ak4 = malloc ( n * sizeof ( double ) ) ;
170 if ( ak4 == NULL ) {
171 free ( ak3 ) ;
172 free ( ak2 ) ;
173 return 1 ;
174 }
175 ak5 = malloc ( n * sizeof ( double ) ) ;
176 if ( ak5 == NULL ) {
177 free ( ak4 ) ;
178 free ( ak3 ) ;
179 free ( ak2 ) ;
180 return 1 ;
181 }
182 ak6 = malloc ( n * sizeof ( double ) ) ;
183 if ( ak6 == NULL ) {
184 free ( ak5 ) ;
185 free ( ak4 ) ;
186 free ( ak3 ) ;
187 free ( ak2 ) ;
188 return 1 ;
189 }
190 ak7 = malloc ( n * sizeof ( double ) ) ;
191 if ( ak7 == NULL ) {
192 free ( ak6 ) ;
193 free ( ak5 ) ;
194 free ( ak4 ) ;
195 free ( ak3 ) ;
196 free ( ak2 ) ;
197 return 1 ;
198 }
199 ak8 = malloc ( n * sizeof ( double ) ) ;
200 if ( ak8 == NULL ) {
201 free ( ak7 ) ;
202 free ( ak6 ) ;
203 free ( ak5 ) ;
204 free ( ak4 ) ;
205 free ( ak3 ) ;
206 free ( ak2 ) ;
207 return 1 ;
208 }
209 ak9 = malloc ( n * sizeof ( double ) ) ;
210 if ( ak9 == NULL ) {
211 free ( ak8 ) ;
212 free ( ak7 ) ;
213 free ( ak6 ) ;
214 free ( ak5 ) ;
215 free ( ak4 ) ;
216 free ( ak3 ) ;
217 free ( ak2 ) ;
218 return 1 ;
219 }
220 ak10 = malloc ( n * sizeof ( double ) ) ;
221 if ( ak10 == NULL ) {
222 free ( ak9 ) ;
223 free ( ak8 ) ;
224 free ( ak7 ) ;
225 free ( ak6 ) ;
226 free ( ak5 ) ;
227 free ( ak4 ) ;
228 free ( ak3 ) ;
229 free ( ak2 ) ;
230 return 1 ;
231 }
232 ak11 = malloc ( n * sizeof ( double ) ) ;
233 if ( ak11 == NULL ) {
234 free ( ak10 ) ;
235 free ( ak9 ) ;
236 free ( ak8 ) ;
237 free ( ak7 ) ;
238 free ( ak6 ) ;
239 free ( ak5 ) ;
240 free ( ak4 ) ;
241 free ( ak3 ) ;
242 free ( ak2 ) ;
243 return 1 ;
244 }
245 ak12 = malloc ( n * sizeof ( double ) ) ;
246 if ( ak12 == NULL ) {
247 free ( ak11 ) ;
248 free ( ak10 ) ;
249 free ( ak9 ) ;
250 free ( ak8 ) ;
251 free ( ak7 ) ;
252 free ( ak6 ) ;
253 free ( ak5 ) ;
254 free ( ak4 ) ;
255 free ( ak3 ) ;
256 free ( ak2 ) ;
257 return 1 ;
258 }
259 ak13 = malloc ( n * sizeof ( double ) ) ;
260 if ( ak13 == NULL ) {
261 free ( ak12 ) ;
262 free ( ak11 ) ;
263 free ( ak10 ) ;
264 free ( ak9 ) ;
265 free ( ak8 ) ;
266 free ( ak7 ) ;
267 free ( ak6 ) ;
268 free ( ak5 ) ;
269 free ( ak4 ) ;
270 free ( ak3 ) ;
271 free ( ak2 ) ;
272 return 1 ;
273 }
274
275 ytemp = malloc ( n * sizeof ( double ) ) ;
276 if ( ytemp == NULL ) {
277 free ( ak13 ) ;
278 free ( ak12 ) ;
279 free ( ak11 ) ;
280 free ( ak10 ) ;
281 free ( ak9 ) ;
282 free ( ak8 ) ;
283 free ( ak7 ) ;
284 free ( ak6 ) ;
285 free ( ak5 ) ;
286 free ( ak4 ) ;
287 free ( ak3 ) ;
288 free ( ak2 ) ;
289 return 1 ;
290 }
291
292 /*
293 * Do it!
294 */
295
296 for ( i = 0; i < n; i++ ) {
297 ytemp[i] = y[i] + B(1,0) * h * dydx[i] ;
298 }
299 ( *derivs ) ( x + A(1) * h, ytemp, ak2, derivsp ) ;
300
301 for ( i = 0; i < n; i++ ) {
302 ytemp[i] = y[i] + h * (B(2,0) * dydx[i] + B(2,1) * ak2[i]) ;
303 }
304 ( *derivs ) ( x + A(2) * h, ytemp, ak3, derivsp ) ;
305
306 for ( i = 0; i < n; i++ ) {
307 ytemp[i] = y[i] + h * (B(3,0) * dydx[i] + B(3,1) * ak2[i] + B(3,2) * ak3[i]) ;
308 }
309 ( *derivs ) ( x + A(3) * h, ytemp, ak4, derivsp ) ;
310
311 for ( i = 0; i < n; i++ ) {
312 ytemp[i] = y[i] + h * (B(4,0) * dydx[i] + B(4,1) * ak2[i] + B(4,2) * ak3[i] + B(4,3) * ak4[i]) ;
313 }
314 ( *derivs ) ( x + A(4) * h, ytemp, ak5, derivsp ) ;
315
316 for ( i = 0; i < n; i++ ) {
317 ytemp[i] = y[i] + h * (B(5,0) * dydx[i] + B(5,1) * ak2[i] + B(5,2) * ak3[i] + B(5,3) * ak4[i] + B(5,4) * ak5[i]) ;
318 }
319 ( *derivs ) ( x + A(5) * h, ytemp, ak6, derivsp ) ;
320
321 for ( i = 0; i < n; i++ ) {
322 ytemp[i] = y[i] + h * (B(6,0) * dydx[i] + B(6,1) * ak2[i] + B(6,2) * ak3[i] + B(6,3) * ak4[i] + B(6,4) * ak5[i] + B(6,5) * ak6[i]) ;
323 }
324 ( *derivs ) ( x + A(6) * h, ytemp, ak7, derivsp ) ;
325
326 for ( i = 0; i < n; i++ ) {
327 ytemp[i] = y[i] + h * (B(7,0) * dydx[i] + B(7,1) * ak2[i] + B(7,2) * ak3[i] + B(7,3) * ak4[i] + B(7,4) * ak5[i] + B(7,5) * ak6[i] + B(7,6) * ak7[i]) ;
328 }
329 ( *derivs ) ( x + A(7) * h, ytemp, ak8, derivsp ) ;
330
331 for ( i = 0; i < n; i++ ) {
332 ytemp[i] = y[i] + h * (B(8,0) * dydx[i] + B(8,1) * ak2[i] + B(8,2) * ak3[i] + B(8,3) * ak4[i] + B(8,4) * ak5[i] + B(8,5) * ak6[i] + B(8,6) * ak7[i] + B(8,7) * ak8[i]) ;
333 }
334 ( *derivs ) ( x + A(8) * h, ytemp, ak9, derivsp ) ;
335
336 for ( i = 0; i < n; i++ ) {
337 ytemp[i] = y[i] + h * (B(9,0) * dydx[i] + B(9,1) * ak2[i] + B(9,2) * ak3[i] + B(9,3) * ak4[i] + B(9,4) * ak5[i] + B(9,5) * ak6[i] + B(9,6) * ak7[i] + B(9,7) * ak8[i] + B(9,8) * ak9[i]) ;
338 }
339 ( *derivs ) ( x + A(9) * h, ytemp, ak10, derivsp ) ;
340
341 for ( i = 0; i < n; i++ ) {
342 ytemp[i] = y[i] + h * (B(10,0) * dydx[i] + B(10,1) * ak2[i] + B(10,2) * ak3[i] + B(10,3) * ak4[i] + B(10,4) * ak5[i] + B(10,5) * ak6[i] + B(10,6) * ak7[i] + B(10,7) * ak8[i] + B(10,8) * ak9[i] + B(10,9) * ak10[i]) ;
343 }
344 ( *derivs ) ( x + A(10) * h, ytemp, ak11, derivsp ) ;
345
346 for ( i = 0; i < n; i++ ) {
347 ytemp[i] = y[i] + h * (B(11,0) * dydx[i] + B(11,1) * ak2[i] + B(11,2) * ak3[i] + B(11,3) * ak4[i] + B(11,4) * ak5[i] + B(11,5) * ak6[i] + B(11,6) * ak7[i] + B(11,7) * ak8[i] + B(11,8) * ak9[i] + B(11,9) * ak10[i] + B(11,10) * ak11[i]) ;
348 }
349 ( *derivs ) ( x + A(11) * h, ytemp, ak12, derivsp ) ;
350
351 for ( i = 0; i < n; i++ ) {
352 ytemp[i] = y[i] + h * (B(12,0) * dydx[i] + B(12,1) * ak2[i] + B(12,2) * ak3[i] + B(12,3) * ak4[i] + B(12,4) * ak5[i] + B(12,5) * ak6[i] + B(12,6) * ak7[i] + B(12,7) * ak8[i] + B(12,8) * ak9[i] + B(12,9) * ak10[i] + B(12,10) * ak11[i] + B(12,11) * ak12[i]) ;
353 }
354 ( *derivs ) ( x + A(12) * h, ytemp, ak13, derivsp ) ;
355
356 for ( i = 0; i < n; i++ ) {
357 yout[i] = y[i] + h * (C(0) * dydx[i] + C(1) * ak2[i] + C(2) * ak3[i] + C(3) * ak4[i] + C(4) * ak5[i] + C(5) * ak6[i] + C(6) * ak7[i] + C(7) * ak8[i] + C(8) * ak9[i] + C(9) * ak10[i] + C(10) * ak11[i] + C(11) * ak12[i] + C(12) * ak13[i]) ;
358 }
359
360 for ( i = 0; i < n; i++ ) {
361 yerr[i] = h * (DC(0) * dydx[i] + DC(1) * ak2[i] + DC(2) * ak3[i] + DC(3) * ak4[i] + DC(4) * ak5[i] + DC(5) * ak6[i] + DC(6) * ak7[i] + DC(7) * ak8[i] + DC(8) * ak9[i] + DC(9) * ak10[i] + DC(10) * ak11[i] + DC(11) * ak12[i] + DC(12) * ak13[i]) ;
362 }
363
364 /*
365 * De-allocate Workspace Memory
366 */
367
368 free ( ytemp ) ;
369 free ( ak13 ) ;
370 free ( ak12 ) ;
371 free ( ak11 ) ;
372 free ( ak10 ) ;
373 free ( ak9 ) ;
374 free ( ak8 ) ;
375 free ( ak7 ) ;
376 free ( ak6 ) ;
377 free ( ak5 ) ;
378 free ( ak4 ) ;
379 free ( ak3 ) ;
380 free ( ak2 ) ;
381
382 return GAL_SUCCESS ;
383
384 /*
385 * Finished.
386 */
387
388 }
389
390 /*
391 * gal - General Astrodynamics Library
392 * Copyright (C) 2008 Paul C. L. Willmott
393 *
394 * This program is free software; you can redistribute it and/or modify
395 * it under the terms of the GNU General Public License as published by
396 * the Free Software Foundation; either version 2 of the License, or
397 * (at your option) any later version.
398 *
399 * This program is distributed in the hope that it will be useful,
400 * but WITHOUT ANY WARRANTY; without even the implied warranty of
401 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
402 * GNU General Public License for more details.
403 *
404 * You should have received a copy of the GNU General Public License along
405 * with this program; if not, write to the Free Software Foundation, Inc.,
406 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
407 *
408 * Contact:
409 *
410 * Paul Willmott
411 * vp9mu@amsat.org
412 */
413
414