操作系统：windows
编译器：vscode（Dev C++可能会编译错误，因为不支持for循环开头括号里有定义）
（关于fftw在linux上跑的完整代码后续也会开源出来哦）
homework.cpp
#include <bits/stdc++.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
//#include <blitz/array.h>
#include "nrutil.h"
#include "nrutil.c"
//#include <fftw.h>
using namespace std;
//using namespace blitz;
void Output(char *fn1, char *fn2, double t,
double* r, double* v);
void Density(double* r, double* n);
void Electric(double* phi, double* E);
void Poisson1D(double* &u, double* v, double kappa);
void rk4_fixed(double &x, double* y,
void (*rhs_eval)(double, double*, double* ),
double h,int n);
void rhs_eval(double t, double* y, double* dydt);
void Load(double* r, double* v, double* y);
void UnLoad(double* y, double* r, double* v);
double distribution(double vb);
double L;
int N, J;
int main()
{
// Parameters
L; // Domain of solution 0 <= x <= L (in Debye lengths)
N; // Number of electrons
J; // Number of grid points
double vb; // Beam velocity
double dt; // Time-step (in inverse plasma frequencies)
double tmax; // Simulation run from t = 0. to t = tmax
// // Get parameters
// printf("Please input N: ");
// scanf("%d", &N);
// printf("Please input vb: ");
// scanf("%lf", &vb);
// printf("Please input L: ");
// scanf("%lf", &L);
// printf("Please input J: ");
// scanf("%d", &J);
// printf("Please input dt: ");
// scanf("%lf", &dt);
N=20000;
vb=3;
L=100;
J=1000;
dt=0.1;
printf("Please input tmax:(默认为17.5) ");
scanf("%lf", &tmax);
int skip = int(tmax / dt) / 10;
if ((N < 1) || (J < 2) || (L <= 0.) || (vb <= 0.) || (dt <= 0.) || (tmax <= 0.) || (skip < 1))
{
printf("Error - invalid input parameters\n");
exit(1);
}
// Set names of output files
char *phase[11];
char *data[11];
phase[0] = "phase0.out";
phase[1] = "phase1.out";
phase[2] = "phase2.out";
phase[3] = "phase3.out";
phase[4] = "phase4.out";
phase[5] = "phase5.out";
phase[6] = "phase6.out";
phase[7] = "phase7.out";
phase[8] = "phase8.out";
phase[9] = "phase9.out";
phase[10] = "phase10.out";
data[0] = "data0.out";
data[1] = "data1.out";
data[2] = "data2.out";
data[3] = "data3.out";
data[4] = "data4.out";
data[5] = "data5.out";
data[6] = "data6.out";
data[7] = "data7.out";
data[8] = "data8.out";
data[9] = "data9.out";
data[10] = "data10.out";
// Initialize solution
double t = 0.;
int seed = time(NULL);
srand(seed);
double* r=dvector(0,N), *v=dvector(0,N);
for (int i = 0; i < N; i++)
{
r[i] = L * double(rand()) / double(RAND_MAX);
v[i] = distribution(vb);
//v[i]=0;
}
Output(phase[0], data[0], t, r, v);
// Evolve solution
double* y=dvector(0,2 * N);
Load(r, v, y);
for (int k = 1; k <= 10; k++)
{
for (int kk = 0; kk < skip; kk++)
{
// Take time-step
rk4_fixed(t, y, rhs_eval, dt,2*N);
// Make sure all coordinates in range 0 to L.
for (int i = 0; i < N; i++)
{
if (y[i] < 0.)
y[i] += L;
if (y[i] > L)
y[i] -= L;
}
printf("t = %11.4e\n", t);
}
printf("Plot %3d\n", k);
// Output data
UnLoad(y, r, v);
Output(phase[k], data[k], t, r, v);
}
return 0;
}
void rk4_fixed(double &x, double* y,
void (*rhs_eval)(double, double*, double* ),
double h,int n)
{
// Array y assumed to be of extent n, where n is no. of coupled
// equations
//int n = y.extent(0);
// Declare local arrays
double* k1=dvector(0,n), *k2=dvector(0,n), *k3=dvector(0,n),
*k4=dvector(0,n), *f=dvector(0,n), *dydx=dvector(0,n);
// Zeroth intermediate step
(*rhs_eval)(x, y, dydx);
for (int j = 0; j < n; j++)
{
k1[j] = h * dydx[j];
f[j] = y[j] + k1[j] / 2.;
}
// First intermediate step
(*rhs_eval)(x + h / 2., f, dydx);
for (int j = 0; j < n; j++)
{
k2[j] = h * dydx[j];
f[j] = y[j] + k2[j] / 2.;
}
// Second intermediate step
(*rhs_eval)(x + h / 2., f, dydx);
for (int j = 0; j < n; j++)
{
k3[j] = h * dydx[j];
f[j] = y[j] + k3[j];
}
// Third intermediate step
(*rhs_eval)(x + h, f, dydx);
for (int j = 0; j < n; j++)
{
k4[j] = h * dydx[j];
}
// Actual step
for (int j = 0; j < n; j++)
{
y[j] += k1[j] / 6. + k2[j] / 3. + k3[j] / 3. + k4[j] / 6.;
}
x += h;
return;
}
void Output(char *fn1, char *fn2, double t,
double* r, double* v)
{
// Write phase-space data
FILE *file = fopen(fn1, "w");
for (int i = 0; i < N; i++)
fprintf(file, "%e,%e\n", r[i], v[i]);
fclose(file);
// Write electric field data
double* ne=dvector(0,J), *n=dvector(0,J), *phi=dvector(0,J), *E=dvector(0,J);
Density(r, ne);
for (int j = 0; j < J; j++)
n[j] = double(J) * ne[j] / double(N) - 1.;
double kappa = 2. * M_PI / L;
//Poisson1D(phi, n, kappa);
for(int i=0;i<J;i++){
phi[i]=sin(M_PI*r[i]/L);
}
Electric(phi, E);
file = fopen(fn2, "w");
for (int j = 0; j < J; j++)
{
double x = double(j) * L / double(J);
fprintf(file, "%e,%e,%e,%e\n", x, ne[j], n[j], E[j]);
}
double x = L;
fprintf(file, "%e,%e,%e,%e\n", x, ne[0], n[0], E[0]);
fclose(file);
}
double distribution(double vb)
{
// Initialize random number generator
static int flag = 0;
if (flag == 0)
{
int seed = time(NULL);
srand(seed);
flag = 1;
}
// Generate random v value
double fmax = 0.5 * (1. + exp(-2. * vb * vb));
double vmin = -5. * vb;
double vmax = +5. * vb;
double v = vmin + (vmax - vmin) * double(rand()) / double(RAND_MAX);
// Accept/reject value
double f = 0.5 * (exp(-(v - vb) * (v - vb) / 2.) +
exp(-(v + vb) * (v + vb) / 2.));
double x = fmax * double(rand()) / double(RAND_MAX);
if (x > f)
return distribution(vb);
else
return v;
}
void Density(double* r, double * n)
{
// Initialize
double dx = L / double(J);
n[0] = 0.;
// Evaluate number density.
for (int i = 0; i < N; i++)
{
int j = int(r[i] / dx);
double y = r[i] / dx - double(j);
n[j] += (1. - y) / dx;
if (j + 1 == J)
n[0] += y / dx;
else
n[j + 1] += y / dx;
}
}
/* void fft_forward(Array<double, 1> f, Array<double, 1> &Fr,
Array<double, 1> &Fi)
{
fftw_complex ff[J], FF[J];
// Load data
for (int j = 0; j < J; j++)
{
c_re(ff[j]) = f(j);
c_im(ff[j]) = 0.;
}
// Call fftw routine
fftw_plan p = fftw_create_plan(J, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_one(p, ff, FF);
fftw_destroy_plan(p);
// Unload data
for (int j = 0; j < J; j++)
{
Fr(j) = c_re(FF[j]);
Fi(j) = c_im(FF[j]);
}
// Normalize data
Fr /= double(J);
Fi /= double(J);
}
// Calculates inverse Fourier transform of arrays Fr and Fi in array f
void fft_backward(Array<double, 1> Fr, Array<double, 1> Fi,
Array<double, 1> &f)
{
fftw_complex ff[J], FF[J];
// Load data
for (int j = 0; j < J; j++)
{
c_re(FF[j]) = Fr(j);
c_im(FF[j]) = Fi(j);
}
// Call fftw routine
fftw_plan p = fftw_create_plan(J, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_one(p, FF, ff);
fftw_destroy_plan(p);
// Unload data
for (int j = 0; j < J; j++)
f(j) = c_re(ff[j]);
} */
// The following routine solves Poisson’s equation in 1 - D to find the instantaneous electric potential on a uniform grid.
// Solves 1-d Poisson equation:
// d^u / dx^2 = v for 0 <= x <= L
// Periodic boundary conditions:
// u(x + L) = u(x), v(x + L) = v(x)
// Arrays u and v assumed to be of length J.
// Now, jth grid point corresponds to
// x_j = j dx for j = 0,J-1
// where dx = L / J.
// Also,
// kappa = 2 pi / L
void Poisson1D(double* &u, double* v, double kappa)
{
// Declare local arrays.
double * Vr=dvector(0,J), *Vi=dvector(0,J), *Ur=dvector(0,J), *Ui=dvector(0,J);
// Fourier transform source term
/* fft_forward(v, Vr, Vi);
// Calculate Fourier transform of u
Ur(0) = Ui(0) = 0.;
for (int j = 1; j <= J / 2; j++)
{
Ur(j) = -Vr(j) / double(j * j) / kappa / kappa;
Ui(j) = -Vi(j) / double(j * j) / kappa / kappa;
}
for (int j = J / 2; j < J; j++)
{
Ur(j) = Ur(J - j);
Ui(j) = -Ui(J - j);
}
// Inverse Fourier transform to obtain u
fft_backward(Ur, Ui, u); */
}
// Calculate electric field from potential
void Electric(double* phi, double* E)
{
double dx = L / double(J);
for (int j = 1; j < J - 1; j++)
E[j] = (phi[j - 1] - phi[j + 1]) / 2. / dx;
E[0] = (phi[J - 1] - phi[1]) / 2. / dx;
E[J - 1] = (phi[J - 2] - phi[0]) / 2. / dx;
}
void rhs_eval(double t, double* y, double* dydt)
{
// Declare local arrays
double* r=dvector(0,N), *v=dvector(0,N), *rdot=dvector(0,N), *vdot=dvector(0,N), *r0=dvector(0,N);
double* ne=dvector(0,J), *rho=dvector(0,J), *phi=dvector(0,J), *E=dvector(0,J);
// Unload data from y
UnLoad(y, r, v);
// Make sure all coordinates in range 0 to L
r0 = r;
for (int i = 0; i < N; i++)
{
if (r0[i] < 0.)
r0[i] += L;
if (r0[i] > L)
r0[i] -= L;
}
// Calculate electron number density
Density(r0, ne);
// Solve Poisson’s equation
double n0 = double(N) / L;
for (int j = 0; j < J; j++)
rho[j] = ne[j] / n0 - 1.;
double kappa = 2. * M_PI / L;
//Poisson1D(phi, rho, kappa);
//Poisson1D(phi, rho, kappa);
for(int i=0;i<J;i++){
phi[i]=sin(M_PI*r0[i]/L);
}
// Calculate electric field
Electric(phi, E);
// Equations of motion
for (int i = 0; i < N; i++)
{
double dx = L / double(J);
int j = int(r0[i] / dx);
double y = r0[i] / dx - double(j);
double Efield;
if (j + 1 == J)
Efield = E[j] * (1. - y) + E[0] * y;
else
Efield = E[j] * (1. - y) + E[j + 1] * y;
rdot[i] = v[i];
vdot[i] = -Efield;
}
// Load data into dydt
Load(rdot, vdot, dydt);
}
// Load particle coordinates into solution vector
void Load(double* r, double* v, double* y)
{
for (int i = 0; i < N; i++)
{
y[i] = r[i];
y[N + i] = v[i];
}
}
// Unload particle coordinates from solution vector
// Load particle coordinates into solution vector
void UnLoad(double* y, double* r, double* v)
{
for (int i = 0; i < N; i++)
{
r[i] = y[i];
v[i] = y[N + i];
}
}nrutil.c
/* CAUTION: This is the combined ANSI and traditional K&R C version
of the Numerical Recipes utility file nrutil.c. Do not confuse
this file with the same-named file nrutil.c that is supplied in
the 'recipes' subdirectory. *That* file contains only ANSI.
*This* file contains both ANSI and traditional K&R versions,
along with #ifdef macros to select the correct version. */
#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(long nl, long nh)
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
#else /* ANSI */
/* traditional - K&R */
#include <stdio.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(error_text)
char error_text[];
/* Numerical Recipes standard error handler */
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(nl,nh)
long nh,nl;
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(nl,nh)
long nh,nl;
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(nl,nh)
long nh,nl;
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(nl,nh)
long nh,nl;
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(nl,nh)
long nh,nl;
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(nrl,nrh,ncl,nch)
long nch,ncl,nrh,nrl;
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
float **a;
long newcl,newrl,oldch,oldcl,oldrh,oldrl;
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(a,nrl,nrh,ncl,nch)
float *a;
long nch,ncl,nrh,nrl;
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
long nch,ncl,ndh,ndl,nrh,nrl;
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
void free_vector(v,nl,nh)
float *v;
long nh,nl;
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(v,nl,nh)
int *v;
long nh,nl;
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(v,nl,nh)
long nh,nl;
unsigned char *v;
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(v,nl,nh)
long nh,nl;
unsigned long *v;
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(v,nl,nh)
double *v;
long nh,nl;
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
long nch,ncl,nrh,nrl;
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
long nch,ncl,nrh,nrl;
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
long nch,ncl,nrh,nrl;
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(b,nrl,nrh,ncl,nch)
float **b;
long nch,ncl,nrh,nrl;
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
float ***t;
long nch,ncl,ndh,ndl,nrh,nrl;
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
#endif /* ANSI */
nrutil.h
/* CAUTION: This is the ANSI C (only) version of the Numerical Recipes
utility file nrutil.h. Do not confuse this file with the same-named
file nrutil.h that is supplied in the 'misc' subdirectory.
*That* file is the one from the book, and contains both ANSI and
traditional K&R versions, along with #ifdef macros to select the
correct version. *This* file contains only ANSI C. */
#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_
static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
(dmaxarg1) : (dmaxarg2))
static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
(dminarg1) : (dminarg2))
static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
(minarg1) : (minarg2))
static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))
static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
(imaxarg1) : (imaxarg2))
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
(iminarg1) : (iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void nrerror(char error_text[]);
float *vector(long nl, long nh); //�����ȸ��������� x[nl...nh]
int *ivector(long nl, long nh); //��������
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
double *dvector(long nl, long nh); //˫�������� x[nl..nh]
float **matrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl);
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_cvector(unsigned char *v, long nl, long nh);
void free_lvector(unsigned long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
#endif /* _NR_UTILS_H_ */
1062
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
