zode.c

gridgen是一款强大的网格生成程序

/******************************************************************************
 *
 * File:           zode.c
 *
 * Created         4.4.2001
 *
 * Author:         Pavel Sakov
 *                 CSIRO Marine Laboratories
 *
 * Purpose:        Solving complex ODEs
 *
 * Description:    Wrapper to "ode.c" for solving complex ODEs
 *
 * Revisions:      none
 *	    
 *****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ode.h"
#include "zode.h"

typedef struct {
    int nz;
    zfunc zcalc;
    zdouble z0;
    zdouble z1;
    void* p;
    void* work;
} calcstruct;

static void calc(double x, double* y, double* y1, void* p)
{
    calcstruct* zs = (calcstruct*) p;
    zdouble* f = (zdouble*) zs->work;
    zdouble* f1 = &((zdouble*) zs->work)[zs->nz];
    zdouble z = (1.0 - x) * zs->z0 + x * zs->z1;
    int i, ii;

    for (i = 0, ii = 0; i < zs->nz; ++i, ii += 2) {
        f[i] = y[ii] + I * y[ii + 1];
        f1[i] = y1[ii] + I * y1[ii + 1];
    }

    zs->zcalc(z, f, f1, zs->p);

    for (i = 0, ii = 0; i < zs->nz; ++i, ii += 2) {
        y[ii] = creal(f[i]);
        y[ii + 1] = cimag(f[i]);
        y1[ii] = creal(f1[i]);
        y1[ii + 1] = cimag(f1[i]);
    }
}

/* Complex ODE integrator.
 * @param ode Real ODE integrator used
 * @param zcalc Derivative function
 * @param nz Dimension
 * @param z Initial argument value
 * @param f Function value (input/output)
 * @param zend Final argument value
 * @param Initial relative step guess (0 < h < 1)
 * @param eps Local precision
 * @param p Custom data; will be passed to `zcalc'
 * @return Passed from `ode'; 1 on success; 0 otherwhile
 */
int zode(integrator ode, zfunc zcalc, int nz, zdouble z, zdouble f[], zdouble zend, double* h, double eps, void* p)
{
    double* y = (double*) f;
    int status = -1;
    calcstruct cs;

    cs.nz = nz;
    cs.zcalc = zcalc;
    cs.z0 = z;
    cs.z1 = zend;
    cs.p = p;
    cs.work = malloc(sizeof(zdouble) * nz * 2);

    status = ode(calc, nz * 2, 0.0, y, 1.0, eps, 1.0, h, NULL, NULL, &cs);

    free(cs.work);

    return status;
}