Commit cc9031cc by Francois Gygi

rewritten to avoid dependencies


git-svn-id: http://qboxcode.org/svn/qb/trunk@494 cba15fb0-1239-40c8-b417-11db7ca47a34
parent bf1d090b
/*******************************************************************************
*
* sinft.c
*
******************************************************************************/
////////////////////////////////////////////////////////////////////////////////
//
// sinft.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: sinft.C,v 1.2 2007-08-13 21:23:51 fgygi Exp $
#include "sinft.h"
#include <math.h>
#include <assert.h>
#include "fftw.h"
#include <vector>
#include <complex>
using namespace std;
void sinft ( double *y, int n )
void sinft(int n, double *y)
{
/*
* Calculates the sine transform of a set of n-real-valued data points
* stored in array y[0:n-1]. The number n must be a power of 2. On exit
* y is replaced by its transform. This program, without change, also
* calculates the inverse sine transform, but in this case the output
* array should be multiplied by 2/n.
* Numerical Recipes, 2nd Edition.
*/
int j;
double sum,y1,y2,wi,wpi,wr,wpr,wtemp;
double theta = 3.141592653589793 /((double) n);
wr = 1.0;
wi = 0.0;
wtemp = sin ( 0.5 * theta );
wpr = -2.0 * wtemp * wtemp;
wpi = sin ( theta );
y[0] = 0.0;
for ( j = 0; j < n/2; j++ )
{
wtemp = wr;
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
y1 = wi * ( y[j+1] + y[n-j-1] );
y2 = 0.5 * ( y[j+1] - y[n-j-1] );
y[j+1] = y1 + y2;
y[n-j-1] = y1 - y2;
}
realft(y,n,1);
sum = 0.0;
y[0] *= 0.5;
y[1] = 0.0;
for ( j = 0; j < n-1; j++,j++ )
{
sum += y[j];
y[j] = y[j+1];
y[j+1] = sum;
}
}
void cosft1 ( double *y, int n )
{
/* Note: the array y contains n+1 elements */
int j;
double sum,y1,y2,wi,wpi,wr,wpr,wtemp;
double theta = 3.141592653589793 / ( (double) n );
wr = 1.0;
wi = 0.0;
wtemp = sin ( 0.5 * theta );
wpr = -2.0 * wtemp * wtemp;
wpi = sin ( theta );
sum = 0.5 * ( y[0] - y[n] );
y[0] = 0.5 * ( y[0] + y[n] );
for ( j = 0; j < n/2-1; j++ )
{
wtemp = wr;
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
y1 = 0.5 * ( y[j+1] + y[n-j-1] );
y2 = y[j+1] - y[n-j-1];
y[j+1] = y1 - wi * y2;
y[n-j-1] = y1 + wi * y2;
sum = sum + wr * y2;
}
realft(y,n,1);
y[n] = y[1];
y[1] = sum;
for ( j = 3; j < n; j++,j++ )
fftw_plan fwplan;
fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
vector<complex<double> > zin(2*n), zout(2*n);
zin[0] = 0.0;
for ( int i = 1; i < n; i++ )
{
sum += y[j];
y[j] = sum;
const double t = y[i];
zin[i] = t;
zin[2*n-i] = -t;
}
}
void realft ( double *data, int n, int isign )
{
int i,i1,i2,i3,i4,np3;
double c1,c2,h1i,h1r,h2i,h2r,wis,wrs,wr,wi,wpr,wpi,wtemp;
double theta = 3.141592653589793 / ( (double) n/2 );
c1 = 0.5;
if ( isign == 1 )
{
c2 = -0.5;
four1(data,n/2,1);
}
else
fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
for ( int i = 0; i < n; i++ )
{
c2 = 0.5;
theta = -theta;
}
wtemp = sin ( 0.5 * theta );
wpr = -2.0 * wtemp * wtemp;
wpi = sin ( theta );
wr = 1.0 + wpr;
wi = wpi;
np3 = n+3;
/*
Ftn pgm:
i1 = 3 5 7 9
i2 = 4 6 8 10
i3 = n-1 n-3 n-5 ...
i4 = n n-2 n-4 ...
C pgm:
i1 = 2 4 6 8
i2 = 3 5 7 9
i3 = n-2 n-4 n-6
i4 = n-1 n-3 n-5
*/
for ( i = 2; i <= n/4; i++ )
{
i1 = i+i-1;
i2 = i1+1;
i3 = np3-i2;
i4 = i3+1;
assert ( i1 >= 1 ); assert ( i1 <= n );
assert ( i2 >= 1 ); assert ( i2 <= n );
assert ( i3 >= 1 ); assert ( i3 <= n );
assert ( i4 >= 1 ); assert ( i4 <= n );
wrs = wr;
wis = wi;
h1r = c1 * ( data[i1-1] + data[i3-1] );
h1i = c1 * ( data[i2-1] - data[i4-1] );
h2r = -c2 * ( data[i2-1] + data[i4-1] );
h2i = c2 * ( data[i1-1] - data[i3-1] );
data[i1-1] = h1r + wrs * h2r - wis * h2i;
data[i2-1] = h1i + wrs * h2i + wis * h2r;
data[i3-1] = h1r - wrs * h2r + wis * h2i;
data[i4-1] = -h1i + wrs * h2i + wis * h2r;
wtemp = wr;
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
if ( isign == 1 )
{
h1r = data[0];
data[0] = h1r + data[1];
data[1] = h1r - data[1];
}
else
{
h1r = data[0];
data[0] = c1 * ( h1r + data[1] );
data[1] = c1 * ( h1r - data[1] );
four1(data,n/2,-1);
y[i] = -0.5 * imag(zout[i]);
}
fftw_destroy_plan(fwplan);
}
void four1 ( double *data, int nn, int isign )
void cosft1(int n, double *y)
{
int i,istep,j,m,mmax,n;
double tempr,tempi,wr,wi,wpr,wpi,wtemp,theta;
n = 2 * nn;
j = 1;
/* Note: the array y contains n+1 elements */
fftw_plan fwplan;
fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
vector<complex<double> > zin(2*n), zout(2*n);
for ( i = 1; i < n; i++,i++ )
zin[0] = y[0];
for ( int i = 1; i < n+1; i++ )
{
if ( j > i )
{
tempr = data[j-1];
tempi = data[j];
data[j-1] = data[i-1];
data[j] = data[i];
data[i-1] = tempr;
data[i] = tempi;
}
m = n/2;
while ( ( m >= 2 ) && ( j > m ) )
{
j -= m;
m /= 2;
}
j += m;
const double t = y[i];
zin[i] = t;
zin[2*n-i] = t;
}
mmax = 2;
while ( mmax < n )
fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
y[0] = 0.5 * real(zout[0]);
for ( int i = 1; i < n; i++ )
{
istep = 2 * mmax;
theta = isign * 2.0 * ( 3.141592653589793 / mmax );
wtemp = sin ( 0.5 * theta );
wpr = -2.0 * wtemp * wtemp;
wpi = sin ( theta );
wr = 1.0;
wi = 0.0;
for ( m = 1; m < mmax; m++,m++ )
{
for ( i = m; i <= n; i += istep )
{
j = i + mmax;
tempr = wr * data[j-1] - wi * data[j];
tempi = wr * data[j] + wi * data[j-1];
data[j-1] = data[i-1] - tempr;
data[j] = data[i] - tempi;
data[i-1] = data[i-1] + tempr;
data[i] = data[i] + tempi;
}
wtemp = wr;
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
y[i] = 0.5 * real(zout[i]);
}
fftw_destroy_plan(fwplan);
}
/*******************************************************************************
*
* sinft.h
*
******************************************************************************/
void sinft ( double *y, int n );
void cosft1 ( double *y, int n );
void realft ( double *data, int n, int isign );
void four1 ( double *data, int nn, int isign );
////////////////////////////////////////////////////////////////////////////////
//
// sinft.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: sinft.h,v 1.2 2007-08-13 21:23:51 fgygi Exp $
void sinft(int n, double *y);
void cosft1(int n, double *y);
/*******************************************************************************
*
* spline.c
*
******************************************************************************/
///////////////////////////////////////////////////////////////////////////////
//
// spline.C
//
///////////////////////////////////////////////////////////////////////////////
// $Id: spline.C,v 1.2 2007-08-13 21:23:51 fgygi Exp $
#include "spline.h"
#include <assert.h>
#include <cassert>
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
void tridsolve(int n, double* d, double* e, double* f, double* x)
{
// solve the tridiagonal system Ax=b
// d[i] = a(i,i)
// e[i] = a(i,i+1) (superdiagonal of A, e[n-1] not defined)
// f[i] = a(i,i-1) (subdiagonal of A, f[0] not defined)
// x[i] = right-hand side b as input
// x[i] = solution as output
for ( int i = 1; i < n; i++ )
{
f[i] /= d[i-1];
d[i] -= f[i]*e[i-1];
}
for ( int i = 1; i < n; i++ )
x[i] -= f[i]*x[i-1];
int i,k;
double p,qn,sig,un,*u = new double[n];
x[n-1] /= d[n-1];
if ( yp1 >= 1.e30 )
for ( int i = n-2; i >= 0; i-- )
x[i] = (x[i]-e[i]*x[i+1])/d[i];
}
void spline(int n, double *x, double *y, double yp_left, double yp_right,
int bcnat_left, int bcnat_right, double *y2)
{
const double third = 1.0/3.0;
const double sixth = 1.0/6.0;
double *d = new double[n];
double *e = new double[n];
double *f = new double[n];
if ( bcnat_left == 0 )
{
y2[0] = 0.0;
u[0] = 0.0;
// use derivative yp_left at x[0]
const double h = x[1]-x[0];
assert(h>0.0);
d[0] = third*h;
e[0] = sixth*h;
f[0] = 0.0;
y2[0] = (y[1]-y[0])/h - yp_left;
}
else
{
y2[0] = -0.5;
assert ( x[1] - x[0] > 0.0 );
u[0] = ( 3.0 / (x[1]-x[0]) ) * ( (y[1]-y[0]) / (x[1]-x[0]) - yp1 );
// use natural spline at x[0]
d[0] = 1.0;
e[0] = 0.0;
f[0] = 0.0;
y2[0] = 0.0;
}
for ( i = 1; i < n-1; i++ )
if ( bcnat_right == 0 )
{
assert ( x[i+1] > x[i] );
sig = ( x[i] - x[i-1] ) / ( x[i+1] - x[i-1] );
p = sig * y2[i-1] + 2.0;
y2[i] = ( sig - 1.0 ) / p;
u[i] = ( 6.0 * ( ( y[i+1] - y[i] ) / ( x[i+1] - x[i] ) -
( y[i] - y[i-1] ) / ( x[i] - x[i-1] ) ) /
( x[i+1] - x[i-1] ) - sig * u[i-1] ) / p;
}
if ( ypn >= 1.e30 )
{
qn = 0.0;
un = 0.0;
// use derivative yp_right at x[n-1]
const double h = x[n-1]-x[n-2];
assert(h>0.0);
d[n-1] = third*h;
e[n-1] = 0.0;
f[n-1] = sixth*h;
y2[n-1] = yp_right - (y[n-1]-y[n-2])/h;
}
else
{
qn = 0.5;
un = ( 3.0 / (x[n-1]-x[n-2]) ) *
( ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]) );
// use natural spline at x[n-1]
d[n-1] = 1.0;
e[n-1] = 0.0;
f[n-1] = 0.0;
y2[n-1] = 0.0;
}
y2[n-1] = ( un - qn * u[n-2] ) / ( qn * y2[n-2] + 1.0 );
for ( k = n-2; k >= 0; k-- )
// tridiagonal matrix
for ( int i = 1; i < n-1; i++ )
{
y2[k] = y2[k] * y2[k+1] + u[k];
const double hp = x[i+1]-x[i];
const double hm = x[i]-x[i-1];
assert(hp>0.0);
assert(hm>0.0);
d[i] = third * (hp+hm);
e[i] = sixth * hp;
f[i] = sixth * hm;
y2[i] = (y[i+1]-y[i])/hp - (y[i]-y[i-1])/hm;
}
delete [] u;
tridsolve(n,d,e,f,y2);
delete [] d;
delete [] e;
delete [] f;
}
void splint (double *xa, double *ya, double *y2a, int n, double x, double *y)
void splint (int n, double *xa, double *ya, double *y2a, double x, double *y)
{
int k,khi,klo;
int k;
double a,b,h;
klo = 0;
khi = n-1;
int kl = 0;
int kh = n-1;
while ( khi - klo > 1 )
while ( kh - kl > 1 )
{
k = ( khi + klo ) / 2;
k = ( kh + kl ) / 2;
if ( xa[k] > x )
khi = k;
kh = k;
else
klo = k;
kl = k;
}
h = xa[khi] - xa[klo];
h = xa[kh] - xa[kl];
assert ( h > 0.0 );
a = ( xa[khi] - x ) / h;
b = ( x - xa[klo] ) / h;
a = ( xa[kh] - x ) / h;
b = ( x - xa[kl] ) / h;
*y = a * ya[klo] + b * ya[khi] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[klo] + (b*b*b-b) * y2a[khi] );
*y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
}
void splintd (double *xa, double *ya, double *y2a,
int n, double x, double *y, double *dy)
void splintd (int n, double *xa, double *ya, double *y2a,
double x, double *y, double *dy)
{
int k,khi,klo;
int k;
double a,b,h;
klo = 0;
khi = n-1;
int kl = 0;
int kh = n-1;
while ( khi - klo > 1 )
while ( kh - kl > 1 )
{
k = ( khi + klo ) / 2;
k = ( kh + kl ) / 2;
if ( xa[k] > x )
khi = k;
kh = k;
else
klo = k;
kl = k;
}
h = xa[khi] - xa[klo];
h = xa[kh] - xa[kl];
assert ( h > 0.0 );
a = ( xa[khi] - x ) / h;
b = ( x - xa[klo] ) / h;
a = ( xa[kh] - x ) / h;
b = ( x - xa[kl] ) / h;
*y = a * ya[klo] + b * ya[khi] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[klo] + (b*b*b-b) * y2a[khi] );
*y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) *
( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] );
*dy = ( ya[khi] - ya[klo] ) / h +
h * ( ( (1.0/6.0) - 0.5 * a * a ) * y2a[klo] +
( 0.5 * b * b - (1.0/6.0) ) * y2a[khi] );
*dy = ( ya[kh] - ya[kl] ) / h +
h * ( ( (1.0/6.0) - 0.5 * a * a ) * y2a[kl] +
( 0.5 * b * b - (1.0/6.0) ) * y2a[kh] );
}
/*******************************************************************************
*
* spline.h
*
******************************************************************************/
#define SPLINE_FLAT_BC 0.0 /* Flat boundary condition (y'=0) */
#define SPLINE_NATURAL_BC 1.e31 /* Natural boundary condition (Y"=0) */
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2);
void splint (double *xa, double *ya, double *y2a, int n, double x, double *y);
void splintd (double *xa, double *ya, double *y2a,
int n, double x, double *y, double *dy);
///////////////////////////////////////////////////////////////////////////////
//
// spline.h
//
///////////////////////////////////////////////////////////////////////////////
// $Id: spline.h,v 1.2 2007-08-13 21:23:51 fgygi Exp $
void spline(int n, double *x, double *y, double yp_left, double yp_right,
int bcnat_left, int bcnat_right, double *y2);
void splint (int n, double *xa, double *ya, double *y2a, double x, double *y);
void splintd (int n, double *xa, double *ya, double *y2a,
double x, double *y, double *dy);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment