#include
#include
#include
/*
compile by: tcc RK4.c -o RK4 -lm
This program generates a curve in four dimensions given values of curvature and torsion.
For each point, we define a local coordinate system given by the tangent, normal, etc.
Grouping the information into an array, I have a 5 by 4 array defined as follows.
state[0][0] = x coordinate
state[0][1] = y coordinate
state[0][2] = z coordinate
state[0][3] = t coordinate
state[1][0] = x unit tangent component
state[1][1] = y unit tangent component
state[1][2] = z unit tangent component
state[1][3] = t unit tangent component
state[2][0] = x unit normal component
state[2][1] = y unit normal component
state[2][2] = z unit normal component
state[2][3] = t unit normal component
state[3][0] = x unit binormal component
state[3][1] = y unit binormal component
state[3][2] = z unit binormal component
state[3][3] = t unit binormal component
state[4][0] = x unit trinormal component
state[4][1] = y unit trinormal component
state[4][2] = z unit trinormal component
state[4][3] = t unit trinormal component
Three curvatures are used.
K1 = standard curvature = 1/radius of curvature = measure of departure from line
K2 = standard torsion = measure of lift out of plane
K3 = next curvature = measure of lift out of volume
To generate a curve from the current state, we need to know the derivatives of our items.
Position: dr/ds = u; or D state[0][.] = state[1][.] // D is d/ds
Tangent: du/ds = K1*n; or D state[1][.] = K1*state[2][.]
Normal: dn/ds = K2*p - K1*u; or D state[2][.] = K2*state[3][.] - K1*state[1][.]
Binormal: db/ds = K3*q - K2*n; or D state[3][.] = K3*state[4][.] - K2*state[2][.]
Trinormal: dq/ds = - K3*b; or D state[4][.] = - K3*state[3][.]
To solve curves for three dimensions, as opposed to 4, set K3 to zero.
To solve for higher dimensions, extend the curvature formulas in the same patterns as above.
To start an integration, we will need to choose initial directions for our unit vectors.
We need an orthonormal set of initial directions.
In this program, I choose to set the state as follows:
position: 0 0 0 0 ;origin initially
tangent: 1 0 0 0 ;along x axis
normal: 0 1 0 0 ;along y axis
binormal: 0 0 1 0 ;along z axis
trinormal: 0 0 0 1 ;along t axis
*/
typedef double state[5][4];
#define npts 2500
double k1,k2,k3; // curvatures
double storage[4][npts];
void Derivative(double ds, state L, state LD)
// ordinary RK4 the first parameter is the position at which derivatives are calculated
// here, it is the distance away from the current position.
{
int i;
for (i=0;i<4;i++) {
LD[0][i] = L[1][i];
LD[1][i] = k1*L[2][i];
LD[2][i] = k2*L[3][i] - k1*L[1][i];
LD[3][i] = k3*L[4][i] - k2*L[2][i];
LD[4][i] = - k3*L[3][i];
}
}
void rk4(state L,state LD, int neqn, int ndim, double s, double ds, state NextValue)
{
state dym,dyt,yt;
double sh,hh,h6;
int i,j;
hh = 0.5*ds;
h6 = ds/6.0;
sh = s + hh;
for (i=0;i