1 /*
2 * - - - - - - - -
3 * g a l _ r k f
4 * - - - - - - - -
5 *
6 * This routine is part of the General Astrodynamics Library
7 *
8 * Description:
9 *
10 * This routine integrates an ordinary deferential equation
11 * using the Runge-Kutte-Fehlberg method.
12 *
13 * Status:
14 *
15 * Math special function.
16 *
17 * Given:
18 *
19 * ystart d[] Starting y values
20 * nvar i Number of equations to integrate
21 * x1 d Starting x value
22 * x2 d Ending x value
23 * eps d Accuracy
24 * h1 d First guess step-size
25 * hmin d Minimum step-size
26 * *derivs func User defined function for calculating the right hand side derivatives
27 * *rkfs func Required Runge-Kutte-Fehlberg stepper routine
28 * *derivsp i Pointer to parameters structure for derivs routine
29 *
30 * Returned:
31 *
32 * ystart d[] Ending y values
33 * gal_rkf i Status
34 * 0 = Success
35 * 1 = Failed to allocate workspace memory
36 * 2 = Stepsize underflow
37 * 3 = Maximum steps exceeded
38 * 4 = Stepsize too small
39 *
40 * References:
41 *
42 * NASA Technical Report TR R-352
43 * Some Experimental Results Concerning The Error Propagation in
44 * Runge-Kutte type integration formulas
45 * by Erwin Fehlberg
46 * October 1970
47 *
48 * This revision:
49 *
50 * 2008 May 18
51 *
52 * Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
53 *
54 *-----------------------------------------------------------------------
55 */
56
57 #include <stdlib.h>
58 #include <math.h>
59 #include "gal_const.h"
60 #include "gal_rkf.h"
61 #include "gal_rkfqs.h"
62
gal_rkf(double ystart[],int nvar,double x1,double x2,double eps,double h1,double hmin,void (* derivs)(double,double[],double[],int *),int (* rkfs)(double[],double[],int,double,double,double[],double[],void (*)(double,double[],double[],int *),int *),int * derivsp)63 int gal_rkf
64 (
65 double ystart[],
66 int nvar,
67 double x1,
68 double x2,
69 double eps,
70 double h1,
71 double hmin,
72 void ( *derivs ) ( double, double [], double [], int * ),
73 int ( *rkfs ) ( double [], double [], int, double, double, double [], double [], void ( * ) ( double, double [], double [], int * ), int * ) ,
74 int *derivsp
75 )
76 {
77
78 #define RKF_MAXSTEP 1000000
79 #define RKF_TINY 1.0e-30
80
81 int nstp, i, jstat ;
82 double x, hnext, hdid, h, *yscal, *y, *dydx ;
83
84 /*
85 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 */
87
88 /*
89 * Allocate Workspace
90 */
91
92 yscal = malloc ( sizeof ( double ) * nvar ) ;
93 if ( yscal == NULL ) {
94 return 1 ;
95 }
96 y = malloc ( sizeof ( double ) * nvar ) ;
97 if ( y == NULL ) {
98 free ( yscal ) ;
99 return 1 ;
100 }
101 dydx = malloc ( sizeof ( double ) * nvar ) ;
102 if ( dydx == NULL ) {
103 free ( y ) ;
104 free ( yscal ) ;
105 return 1 ;
106 }
107
108 /*
109 * Do the stuff
110 */
111
112 x = x1 ;
113 h = GAL_SIGN ( h1, x2-x1 ) ;
114
115 for ( i = 0 ; i < nvar ; i++ ) {
116 y[i] = ystart[i] ;
117 }
118
119 for ( nstp = 0; nstp < RKF_MAXSTEP; nstp++ ) {
120
121 ( *derivs ) ( x, y, dydx, derivsp ) ;
122
123 for ( i = 0; i < nvar; i++ ) {
124 yscal[i] = fabs ( dydx[i] * h ) + RKF_TINY ;
125 }
126
127 if ( ( x + h - x2 ) * ( x + h - x1 ) > 0.0) {
128 h = x2 - x ;
129 }
130
131 if ( ( jstat = gal_rkfqs ( y, dydx, nvar, &x, h, eps, yscal, &hdid, &hnext, derivs, rkfs, derivsp ) ) != GAL_SUCCESS ) {
132 free ( dydx ) ;
133 free ( y ) ;
134 free ( yscal ) ;
135 return jstat ;
136 } ;
137
138 if ( ( x - x2 ) * ( x2 - x1 ) >= 0.0 ) {
139
140 for ( i = 0; i < nvar; i++ ) {
141 ystart[i] = y[i] ;
142 }
143
144 free ( dydx ) ;
145 free ( y ) ;
146 free ( yscal ) ;
147 return GAL_SUCCESS ;
148
149 }
150
151 if ( fabs ( hnext ) <= hmin ) {
152 free ( dydx ) ;
153 free ( y ) ;
154 free ( yscal ) ;
155 return 4 ;
156 }
157
158 h = hnext ;
159
160 }
161
162 /*
163 * De-allocate workspace memory
164 */
165
166 free ( dydx ) ;
167 free ( y ) ;
168 free ( yscal ) ;
169 return 3 ;
170
171 /*
172 * Finished.
173 */
174
175 }
176
177 /*
178 * gal - General Astrodynamics Library
179 * Copyright (C) 2008 Paul C. L. Willmott
180 *
181 * This program is free software; you can redistribute it and/or modify
182 * it under the terms of the GNU General Public License as published by
183 * the Free Software Foundation; either version 2 of the License, or
184 * (at your option) any later version.
185 *
186 * This program is distributed in the hope that it will be useful,
187 * but WITHOUT ANY WARRANTY; without even the implied warranty of
188 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
189 * GNU General Public License for more details.
190 *
191 * You should have received a copy of the GNU General Public License along
192 * with this program; if not, write to the Free Software Foundation, Inc.,
193 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
194 *
195 * Contact:
196 *
197 * Paul Willmott
198 * vp9mu@amsat.org
199 */
200
201