1 /* -- translated by f2c (version 20100827).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "libtinyf2c.h"
14
15 /* ----------------------------------------------------------------------- */
16 /* 4th order Runge-Kutta algorithm routine. */
17 /* Oct. 5, 1990 S.Sakai */
18 /* ----------------------------------------------------------------------- */
19 /* Copyright (C) 2000-2004 GFD Dennou Club. All rights reserved. */
20 /* ----------------------------------------------------------------------- */
odrk4_(integer * n,S_fp fcn,real * t,real * dt,real * x,real * dx,real * xout,real * work)21 /* Subroutine */ int odrk4_(integer *n, S_fp fcn, real *t, real *dt, real *x,
22 real *dx, real *xout, real *work)
23 {
24 /* System generated locals */
25 integer work_dim1, work_offset, i__1;
26
27 /* Local variables */
28 static integer i__;
29 static real tt, dtt;
30
31 /* Parameter adjustments */
32 work_dim1 = *n;
33 work_offset = 1 + work_dim1;
34 work -= work_offset;
35 --xout;
36 --dx;
37 --x;
38
39 /* Function Body */
40 dtt = *dt / 2.f;
41 tt = *t + dtt;
42 i__1 = *n;
43 for (i__ = 1; i__ <= i__1; ++i__) {
44 work[i__ + work_dim1] = x[i__] + dtt * dx[i__];
45 /* L10: */
46 }
47 (*fcn)(n, &tt, &work[work_dim1 + 1], &work[(work_dim1 << 1) + 1]);
48 i__1 = *n;
49 for (i__ = 1; i__ <= i__1; ++i__) {
50 work[i__ + work_dim1] = x[i__] + dtt * work[i__ + (work_dim1 << 1)];
51 /* L20: */
52 }
53 (*fcn)(n, &tt, &work[work_dim1 + 1], &work[work_dim1 * 3 + 1]);
54 tt = *t + *dt;
55 i__1 = *n;
56 for (i__ = 1; i__ <= i__1; ++i__) {
57 work[i__ + work_dim1] = x[i__] + *dt * work[i__ + work_dim1 * 3];
58 work[i__ + work_dim1 * 3] = work[i__ + (work_dim1 << 1)] + work[i__ +
59 work_dim1 * 3];
60 /* L30: */
61 }
62 (*fcn)(n, &tt, &work[work_dim1 + 1], &work[(work_dim1 << 1) + 1]);
63 dtt = *dt / 6.f;
64 i__1 = *n;
65 for (i__ = 1; i__ <= i__1; ++i__) {
66 xout[i__] = x[i__] + dtt * (dx[i__] + work[i__ + work_dim1 * 3] * 2.f
67 + work[i__ + (work_dim1 << 1)]);
68 /* L40: */
69 }
70 return 0;
71 } /* odrk4_ */
72
73