当前位置：网站首页>PIC模拟（Particle-in-Cell Codes） （任务A和任务C）
PIC模拟（Particle-in-Cell Codes） （任务A和任务C）
2022-06-09 09:50:00 【小米米小life】
操作系统：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_ */
边栏推荐
- Kubernetes第七篇：Pod進階、Controller進階、Resource和Dashboard
- WPF 实现带明细的环形图表
- How many points can you get if the programmer's college entrance examination paper is exposed?
- Machine learning notes - breast cancer classification using keras and deep learning
- 用好条件访问，远程办公更安全高效
- 基于云的 LDAP 入门（下）
- 2289. 使数组按非递减顺序排列
- InfoQ 极客传媒 15 周年庆征文| 迁移 Eureka 到 Nacos 之双注册双订阅模式
- 如何实现合规的远程办公？
- 循环神经网络理论——【torch学习笔记】
猜你喜欢
“当你不再是程序员，很多事会脱离掌控”—— 对话全球最大独立开源公司SUSE CTO...
Openstack explanation (18) -- Nova service startup and service creation
Notes on the development of raspberry pie (15): Raspberry pie 4b+ compile and install MySQL database from the source code
How many points can you get if the programmer's college entrance examination paper is exposed?
Machine learning notes - Introduction to R language learning series I
【genius_platform软件平台开发】第三十五讲：UDP进行跨网段广播
失业潮？元宇宙开拓全新的就业机会
不加班的测试开发工程师不是好程序员？可能不是一只笨鸟，但一直在先飞......
Mingdao cloud on the list of office software in China's information and innovation industry in 2022
15个必知的Mysql索引失效场景，别再踩坑了！
随机推荐
Key configuration points of video fusion cloud service easycvr platform deployed in ECS
损失 3 亿美元后，IBM 宣布退出俄罗斯！
Sequence model - 【 torch learning notes 】
leetcode. 36 --- effective Sudoku
Go strconv package
Machine learning notes - Interpretation of u-net papers
Terrain learning summary (6) -- terrain practice based on Alibaba cloud platform
MicroNet：以极低的 FLOP 实现图像识别
1324. print word vertically - Li Kou Shuangbai code
MSF information collection based on FTP protocol
基于云的 LDAP 入门（上）
Machine learning notes - explore the keras dataset
什么样的数字藏品平台才是好平台？
循环神经网络理论——【torch学习笔记】
【genius_platform软件平台开发】第一万零一讲：电力项目dz产品windows环境vs2017编译遇到的报错汇总
Test development engineers who don't work overtime are not good programmers? It may not be a stupid bird, but it always flies first
Why do we want to be young when we grow up?
2289. arrange arrays in non decreasing order
构建词表与抽样——【torch学习笔记】
剑指 Offer 19. 正则表达式匹配-递归