c语言spline

时间：2022-12-23 16:16:35
 #define NRANSI

 #include "nrutil.h"

 void spline(float x[], float y[], int n, float yp1, float ypn, float y2[])

 {

     int i,k;

     float p,qn,sig,un,*u;

     u=vector(,n-);

     if (yp1 > 0.99e30)

         y2[]=u[]=0.0;

     else {

         y2[] = -0.5;

         u[]=(3.0/(x[]-x[]))*((y[]-y[])/(x[]-x[])-yp1);

     }

     for (i=;i<=n-;i++) {

         sig=(x[i]-x[i-])/(x[i+]-x[i-]);

         p=sig*y2[i-]+2.0;

         y2[i]=(sig-1.0)/p;

         u[i]=(y[i+]-y[i])/(x[i+]-x[i]) - (y[i]-y[i-])/(x[i]-x[i-]);

         u[i]=(6.0*u[i]/(x[i+]-x[i-])-sig*u[i-])/p;

     }

     if (ypn > 0.99e30)

         qn=un=0.0;

     else {

         qn=0.5;

         un=(3.0/(x[n]-x[n-]))*(ypn-(y[n]-y[n-])/(x[n]-x[n-]));

     }

     y2[n]=(un-qn*u[n-])/(qn*y2[n-]+1.0);

     for (k=n-;k>=;k--)

         y2[k]=y2[k]*y2[k+]+u[k];

     free_vector(u,,n-);

 }

 #undef NRANSI

spline.c

c语言spline

 #ifndef _NR_H_

 #define _NR_H_

 #ifndef _FCOMPLEX_DECLARE_T_

 typedef struct FCOMPLEX {float r,i;} fcomplex;

 #define _FCOMPLEX_DECLARE_T_

 #endif /* _FCOMPLEX_DECLARE_T_ */

 #ifndef _ARITHCODE_DECLARE_T_

 typedef struct {

     unsigned long *ilob,*iupb,*ncumfq,jdif,nc,minint,nch,ncum,nrad;

 } arithcode;

 #define _ARITHCODE_DECLARE_T_

 #endif /* _ARITHCODE_DECLARE_T_ */

 #ifndef _HUFFCODE_DECLARE_T_

 typedef struct {

     unsigned long *icod,*ncod,*left,*right,nch,nodemax;

 } huffcode;

 #define _HUFFCODE_DECLARE_T_

 #endif /* _HUFFCODE_DECLARE_T_ */

 #include <stdio.h>

 #if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */

 void addint(double **uf, double **uc, double **res, int nf);

 void airy(float x, float *ai, float *bi, float *aip, float *bip);

 void amebsa(float **p, float y[], int ndim, float pb[],    float *yb,

     float ftol, float (*funk)(float []), int *iter, float temptr);

 void amoeba(float **p, float y[], int ndim, float ftol,

     float (*funk)(float []), int *iter);

 float amotry(float **p, float y[], float psum[], int ndim,

     float (*funk)(float []), int ihi, float fac);

 float amotsa(float **p, float y[], float psum[], int ndim, float pb[],

     float *yb, float (*funk)(float []), int ihi, float *yhi, float fac);

 void anneal(float x[], float y[], int iorder[], int ncity);

 double anorm2(double **a, int n);

 void arcmak(unsigned long nfreq[], unsigned long nchh, unsigned long nradd,

     arithcode *acode);

 void arcode(unsigned long *ich, unsigned char **codep, unsigned long *lcode,

     unsigned long *lcd, int isign, arithcode *acode);

 void arcsum(unsigned long iin[], unsigned long iout[], unsigned long ja,

     int nwk, unsigned long nrad, unsigned long nc);

 void asolve(unsigned long n, double b[], double x[], int itrnsp);

 void atimes(unsigned long n, double x[], double r[], int itrnsp);

 void avevar(float data[], unsigned long n, float *ave, float *var);

 void balanc(float **a, int n);

 void banbks(float **a, unsigned long n, int m1, int m2, float **al,

     unsigned long indx[], float b[]);

 void bandec(float **a, unsigned long n, int m1, int m2, float **al,

     unsigned long indx[], float *d);

 void banmul(float **a, unsigned long n, int m1, int m2, float x[], float b[]);

 void bcucof(float y[], float y1[], float y2[], float y12[], float d1,

     float d2, float **c);

 void bcuint(float y[], float y1[], float y2[], float y12[],

     float x1l, float x1u, float x2l, float x2u, float x1,

     float x2, float *ansy, float *ansy1, float *ansy2);

 void beschb(double x, double *gam1, double *gam2, double *gampl,

     double *gammi);

 float bessi(int n, float x);

 float bessi0(float x);

 float bessi1(float x);

 void bessik(float x, float xnu, float *ri, float *rk, float *rip,

     float *rkp);

 float bessj(int n, float x);

 float bessj0(float x);

 float bessj1(float x);

 void bessjy(float x, float xnu, float *rj, float *ry, float *rjp,

     float *ryp);

 float bessk(int n, float x);

 float bessk0(float x);

 float bessk1(float x);

 float bessy(int n, float x);

 float bessy0(float x);

 float bessy1(float x);

 float beta(float z, float w);

 float betacf(float a, float b, float x);

 float betai(float a, float b, float x);

 float bico(int n, int k);

 void bksub(int ne, int nb, int jf, int k1, int k2, float ***c);

 float bnldev(float pp, int n, long *idum);

 float brent(float ax, float bx, float cx,

     float (*f)(float), float tol, float *xmin);

 void broydn(float x[], int n, int *check,

     void (*vecfunc)(int, float [], float []));

 void bsstep(float y[], float dydx[], int nv, float *xx, float htry,

     float eps, float yscal[], float *hdid, float *hnext,

     void (*derivs)(float, float [], float []));

 void caldat(long julian, int *mm, int *id, int *iyyy);

 void chder(float a, float b, float c[], float cder[], int n);

 float chebev(float a, float b, float c[], int m, float x);

 void chebft(float a, float b, float c[], int n, float (*func)(float));

 void chebpc(float c[], float d[], int n);

 void chint(float a, float b, float c[], float cint[], int n);

 float chixy(float bang);

 void choldc(float **a, int n, float p[]);

 void cholsl(float **a, int n, float p[], float b[], float x[]);

 void chsone(float bins[], float ebins[], int nbins, int knstrn,

     float *df, float *chsq, float *prob);

 void chstwo(float bins1[], float bins2[], int nbins, int knstrn,

     float *df, float *chsq, float *prob);

 void cisi(float x, float *ci, float *si);

 void cntab1(int **nn, int ni, int nj, float *chisq,

     float *df, float *prob, float *cramrv, float *ccc);

 void cntab2(int **nn, int ni, int nj, float *h, float *hx, float *hy,

     float *hygx, float *hxgy, float *uygx, float *uxgy, float *uxy);

 void convlv(float data[], unsigned long n, float respns[], unsigned long m,

     int isign, float ans[]);

 void copy(double **aout, double **ain, int n);

 void correl(float data1[], float data2[], unsigned long n, float ans[]);

 void cosft(float y[], int n, int isign);

 void cosft1(float y[], int n);

 void cosft2(float y[], int n, int isign);

 void covsrt(float **covar, int ma, int ia[], int mfit);

 void crank(unsigned long n, float w[], float *s);

 void cyclic(float a[], float b[], float c[], float alpha, float beta,

     float r[], float x[], unsigned long n);

 void daub4(float a[], unsigned long n, int isign);

 float dawson(float x);

 float dbrent(float ax, float bx, float cx,

     float (*f)(float), float (*df)(float), float tol, float *xmin);

 void ddpoly(float c[], int nc, float x, float pd[], int nd);

 int decchk(char string[], int n, char *ch);

 void derivs(float x, float y[], float dydx[]);

 float df1dim(float x);

 void dfour1(double data[], unsigned long nn, int isign);

 void dfpmin(float p[], int n, float gtol, int *iter, float *fret,

     float (*func)(float []), void (*dfunc)(float [], float []));

 float dfridr(float (*func)(float), float x, float h, float *err);

 void dftcor(float w, float delta, float a, float b, float endpts[],

     float *corre, float *corim, float *corfac);

 void dftint(float (*func)(float), float a, float b, float w,

     float *cosint, float *sinint);

 void difeq(int k, int k1, int k2, int jsf, int is1, int isf,

     int indexv[], int ne, float **s, float **y);

 void dlinmin(float p[], float xi[], int n, float *fret,

     float (*func)(float []), void (*dfunc)(float [], float[]));

 double dpythag(double a, double b);

 void drealft(double data[], unsigned long n, int isign);

 void dsprsax(double sa[], unsigned long ija[], double x[], double b[],

     unsigned long n);

 void dsprstx(double sa[], unsigned long ija[], double x[], double b[],

     unsigned long n);

 void dsvbksb(double **u, double w[], double **v, int m, int n, double b[],

     double x[]);

 void dsvdcmp(double **a, int m, int n, double w[], double **v);

 void eclass(int nf[], int n, int lista[], int listb[], int m);

 void eclazz(int nf[], int n, int (*equiv)(int, int));

 float ei(float x);

 void eigsrt(float d[], float **v, int n);

 float elle(float phi, float ak);

 float ellf(float phi, float ak);

 float ellpi(float phi, float en, float ak);

 void elmhes(float **a, int n);

 float erfcc(float x);

 float erff(float x);

 float erffc(float x);

 void eulsum(float *sum, float term, int jterm, float wksp[]);

 float evlmem(float fdt, float d[], int m, float xms);

 float expdev(long *idum);

 float expint(int n, float x);

 float f1(float x);

 float f1dim(float x);

 float f2(float y);

 float f3(float z);

 float factln(int n);

 float factrl(int n);

 void fasper(float x[], float y[], unsigned long n, float ofac, float hifac,

     float wk1[], float wk2[], unsigned long nwk, unsigned long *nout,

     unsigned long *jmax, float *prob);

 void fdjac(int n, float x[], float fvec[], float **df,

     void (*vecfunc)(int, float [], float []));

 void fgauss(float x, float a[], float *y, float dyda[], int na);

 void fill0(double **u, int n);

 void fit(float x[], float y[], int ndata, float sig[], int mwt,

     float *a, float *b, float *siga, float *sigb, float *chi2, float *q);

 void fitexy(float x[], float y[], int ndat, float sigx[], float sigy[],

     float *a, float *b, float *siga, float *sigb, float *chi2, float *q);

 void fixrts(float d[], int m);

 void fleg(float x, float pl[], int nl);

 void flmoon(int n, int nph, long *jd, float *frac);

 float fmin(float x[]);

 void four1(float data[], unsigned long nn, int isign);

 void fourew(FILE *file[], int *na, int *nb, int *nc, int *nd);

 void fourfs(FILE *file[], unsigned long nn[], int ndim, int isign);

 void fourn(float data[], unsigned long nn[], int ndim, int isign);

 void fpoly(float x, float p[], int np);

 void fred2(int n, float a, float b, float t[], float f[], float w[],

     float (*g)(float), float (*ak)(float, float));

 float fredin(float x, int n, float a, float b, float t[], float f[], float w[],

     float (*g)(float), float (*ak)(float, float));

 void frenel(float x, float *s, float *c);

 void frprmn(float p[], int n, float ftol, int *iter, float *fret,

     float (*func)(float []), void (*dfunc)(float [], float []));

 void ftest(float data1[], unsigned long n1, float data2[], unsigned long n2,

     float *f, float *prob);

 float gamdev(int ia, long *idum);

 float gammln(float xx);

 float gammp(float a, float x);

 float gammq(float a, float x);

 float gasdev(long *idum);

 void gaucof(int n, float a[], float b[], float amu0, float x[], float w[]);

 void gauher(float x[], float w[], int n);

 void gaujac(float x[], float w[], int n, float alf, float bet);

 void gaulag(float x[], float w[], int n, float alf);

 void gauleg(float x1, float x2, float x[], float w[], int n);

 void gaussj(float **a, int n, float **b, int m);

 void gcf(float *gammcf, float a, float x, float *gln);

 float golden(float ax, float bx, float cx, float (*f)(float), float tol,

     float *xmin);

 void gser(float *gamser, float a, float x, float *gln);

 void hpsel(unsigned long m, unsigned long n, float arr[], float heap[]);

 void hpsort(unsigned long n, float ra[]);

 void hqr(float **a, int n, float wr[], float wi[]);

 void hufapp(unsigned long index[], unsigned long nprob[], unsigned long n,

     unsigned long i);

 void hufdec(unsigned long *ich, unsigned char *code, unsigned long lcode,

     unsigned long *nb, huffcode *hcode);

 void hufenc(unsigned long ich, unsigned char **codep, unsigned long *lcode,

     unsigned long *nb, huffcode *hcode);

 void hufmak(unsigned long nfreq[], unsigned long nchin, unsigned long *ilong,

     unsigned long *nlong, huffcode *hcode);

 void hunt(float xx[], unsigned long n, float x, unsigned long *jlo);

 void hypdrv(float s, float yy[], float dyyds[]);

 fcomplex hypgeo(fcomplex a, fcomplex b, fcomplex c, fcomplex z);

 void hypser(fcomplex a, fcomplex b, fcomplex c, fcomplex z,

     fcomplex *series, fcomplex *deriv);

 unsigned short icrc(unsigned short crc, unsigned char *bufptr,

     unsigned long len, short jinit, int jrev);

 unsigned short icrc1(unsigned short crc, unsigned char onech);

 unsigned long igray(unsigned long n, int is);

 void iindexx(unsigned long n, long arr[], unsigned long indx[]);

 void indexx(unsigned long n, float arr[], unsigned long indx[]);

 void interp(double **uf, double **uc, int nf);

 int irbit1(unsigned long *iseed);

 int irbit2(unsigned long *iseed);

 void jacobi(float **a, int n, float d[], float **v, int *nrot);

 void jacobn(float x, float y[], float dfdx[], float **dfdy, int n);

 long julday(int mm, int id, int iyyy);

 void kendl1(float data1[], float data2[], unsigned long n, float *tau, float *z,

     float *prob);

 void kendl2(float **tab, int i, int j, float *tau, float *z, float *prob);

 void kermom(double w[], double y, int m);

 void ks2d1s(float x1[], float y1[], unsigned long n1,

     void (*quadvl)(float, float, float *, float *, float *, float *),

     float *d1, float *prob);

 void ks2d2s(float x1[], float y1[], unsigned long n1, float x2[], float y2[],

     unsigned long n2, float *d, float *prob);

 void ksone(float data[], unsigned long n, float (*func)(float), float *d,

     float *prob);

 void kstwo(float data1[], unsigned long n1, float data2[], unsigned long n2,

     float *d, float *prob);

 void laguer(fcomplex a[], int m, fcomplex *x, int *its);

 void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],

     int ma, float **covar, float *chisq, void (*funcs)(float, float [], int));

 void linbcg(unsigned long n, double b[], double x[], int itol, double tol,

      int itmax, int *iter, double *err);

 void linmin(float p[], float xi[], int n, float *fret,

     float (*func)(float []));

 void lnsrch(int n, float xold[], float fold, float g[], float p[], float x[],

      float *f, float stpmax, int *check, float (*func)(float []));

 void load(float x1, float v[], float y[]);

 void load1(float x1, float v1[], float y[]);

 void load2(float x2, float v2[], float y[]);

 void locate(float xx[], unsigned long n, float x, unsigned long *j);

 void lop(double **out, double **u, int n);

 void lubksb(float **a, int n, int *indx, float b[]);

 void ludcmp(float **a, int n, int *indx, float *d);

 void machar(int *ibeta, int *it, int *irnd, int *ngrd,

     int *machep, int *negep, int *iexp, int *minexp, int *maxexp,

     float *eps, float *epsneg, float *xmin, float *xmax);

 void matadd(double **a, double **b, double **c, int n);

 void matsub(double **a, double **b, double **c, int n);

 void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev);

 void memcof(float data[], int n, int m, float *xms, float d[]);

 int metrop(float de, float t);

 void mgfas(double **u, int n, int maxcyc);

 void mglin(double **u, int n, int ncycle);

 float midexp(float (*funk)(float), float aa, float bb, int n);

 float midinf(float (*funk)(float), float aa, float bb, int n);

 float midpnt(float (*func)(float), float a, float b, int n);

 float midsql(float (*funk)(float), float aa, float bb, int n);

 float midsqu(float (*funk)(float), float aa, float bb, int n);

 void miser(float (*func)(float []), float regn[], int ndim, unsigned long npts,

     float dith, float *ave, float *var);

 void mmid(float y[], float dydx[], int nvar, float xs, float htot,

     int nstep, float yout[], void (*derivs)(float, float[], float[]));

 void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,

     float *fc, float (*func)(float));

 void mnewt(int ntrial, float x[], int n, float tolx, float tolf);

 void moment(float data[], int n, float *ave, float *adev, float *sdev,

     float *var, float *skew, float *curt);

 void mp2dfr(unsigned char a[], unsigned char s[], int n, int *m);

 void mpadd(unsigned char w[], unsigned char u[], unsigned char v[], int n);

 void mpdiv(unsigned char q[], unsigned char r[], unsigned char u[],

     unsigned char v[], int n, int m);

 void mpinv(unsigned char u[], unsigned char v[], int n, int m);

 void mplsh(unsigned char u[], int n);

 void mpmov(unsigned char u[], unsigned char v[], int n);

 void mpmul(unsigned char w[], unsigned char u[], unsigned char v[], int n,

     int m);

 void mpneg(unsigned char u[], int n);

 void mppi(int n);

 void mprove(float **a, float **alud, int n, int indx[], float b[],

     float x[]);

 void mpsad(unsigned char w[], unsigned char u[], int n, int iv);

 void mpsdv(unsigned char w[], unsigned char u[], int n, int iv, int *ir);

 void mpsmu(unsigned char w[], unsigned char u[], int n, int iv);

 void mpsqrt(unsigned char w[], unsigned char u[], unsigned char v[], int n,

     int m);

 void mpsub(int *is, unsigned char w[], unsigned char u[], unsigned char v[],

     int n);

 void mrqcof(float x[], float y[], float sig[], int ndata, float a[],

     int ia[], int ma, float **alpha, float beta[], float *chisq,

     void (*funcs)(float, float [], float *, float [], int));

 void mrqmin(float x[], float y[], float sig[], int ndata, float a[],

     int ia[], int ma, float **covar, float **alpha, float *chisq,

     void (*funcs)(float, float [], float *, float [], int), float *alamda);

 void newt(float x[], int n, int *check,

     void (*vecfunc)(int, float [], float []));

 void odeint(float ystart[], int nvar, float x1, float x2,

     float eps, float h1, float hmin, int *nok, int *nbad,

     void (*derivs)(float, float [], float []),

     void (*rkqs)(float [], float [], int, float *, float, float,

     float [], float *, float *, void (*)(float, float [], float [])));

 void orthog(int n, float anu[], float alpha[], float beta[], float a[],

     float b[]);

 void pade(double cof[], int n, float *resid);

 void pccheb(float d[], float c[], int n);

 void pcshft(float a, float b, float d[], int n);

 void pearsn(float x[], float y[], unsigned long n, float *r, float *prob,

     float *z);

 void period(float x[], float y[], int n, float ofac, float hifac,

     float px[], float py[], int np, int *nout, int *jmax, float *prob);

 void piksr2(int n, float arr[], float brr[]);

 void piksrt(int n, float arr[]);

 void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k,

     float ***c, float **s);

 float plgndr(int l, int m, float x);

 float poidev(float xm, long *idum);

 void polcoe(float x[], float y[], int n, float cof[]);

 void polcof(float xa[], float ya[], int n, float cof[]);

 void poldiv(float u[], int n, float v[], int nv, float q[], float r[]);

 void polin2(float x1a[], float x2a[], float **ya, int m, int n,

     float x1, float x2, float *y, float *dy);

 void polint(float xa[], float ya[], int n, float x, float *y, float *dy);

 void powell(float p[], float **xi, int n, float ftol, int *iter, float *fret,

     float (*func)(float []));

 void predic(float data[], int ndata, float d[], int m, float future[], int nfut);

 float probks(float alam);

 void psdes(unsigned long *lword, unsigned long *irword);

 void pwt(float a[], unsigned long n, int isign);

 void pwtset(int n);

 float pythag(float a, float b);

 void pzextr(int iest, float xest, float yest[], float yz[], float dy[],

     int nv);

 float qgaus(float (*func)(float), float a, float b);

 void qrdcmp(float **a, int n, float *c, float *d, int *sing);

 float qromb(float (*func)(float), float a, float b);

 float qromo(float (*func)(float), float a, float b,

     float (*choose)(float (*)(float), float, float, int));

 void qroot(float p[], int n, float *b, float *c, float eps);

 void qrsolv(float **a, int n, float c[], float d[], float b[]);

 void qrupdt(float **r, float **qt, int n, float u[], float v[]);

 float qsimp(float (*func)(float), float a, float b);

 float qtrap(float (*func)(float), float a, float b);

 float quad3d(float (*func)(float, float, float), float x1, float x2);

 void quadct(float x, float y, float xx[], float yy[], unsigned long nn,

     float *fa, float *fb, float *fc, float *fd);

 void quadmx(float **a, int n);

 void quadvl(float x, float y, float *fa, float *fb, float *fc, float *fd);

 float ran0(long *idum);

 float ran1(long *idum);

 float ran2(long *idum);

 float ran3(long *idum);

 float ran4(long *idum);

 void rank(unsigned long n, unsigned long indx[], unsigned long irank[]);

 void ranpt(float pt[], float regn[], int n);

 void ratint(float xa[], float ya[], int n, float x, float *y, float *dy);

 void ratlsq(double (*fn)(double), double a, double b, int mm, int kk,

     double cof[], double *dev);

 double ratval(double x, double cof[], int mm, int kk);

 float rc(float x, float y);

 float rd(float x, float y, float z);

 void realft(float data[], unsigned long n, int isign);

 void rebin(float rc, int nd, float r[], float xin[], float xi[]);

 void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf,

     int ic1, int jc1, int jcf, int kc, float ***c, float **s);

 void relax(double **u, double **rhs, int n);

 void relax2(double **u, double **rhs, int n);

 void resid(double **res, double **u, double **rhs, int n);

 float revcst(float x[], float y[], int iorder[], int ncity, int n[]);

 void reverse(int iorder[], int ncity, int n[]);

 float rf(float x, float y, float z);

 float rj(float x, float y, float z, float p);

 void rk4(float y[], float dydx[], int n, float x, float h, float yout[],

     void (*derivs)(float, float [], float []));

 void rkck(float y[], float dydx[], int n, float x, float h,

     float yout[], float yerr[], void (*derivs)(float, float [], float []));

 void rkdumb(float vstart[], int nvar, float x1, float x2, int nstep,

     void (*derivs)(float, float [], float []));

 void rkqs(float y[], float dydx[], int n, float *x,

     float htry, float eps, float yscal[], float *hdid, float *hnext,

     void (*derivs)(float, float [], float []));

 void rlft3(float ***data, float **speq, unsigned long nn1,

     unsigned long nn2, unsigned long nn3, int isign);

 float rofunc(float b);

 void rotate(float **r, float **qt, int n, int i, float a, float b);

 void rsolv(float **a, int n, float d[], float b[]);

 void rstrct(double **uc, double **uf, int nc);

 float rtbis(float (*func)(float), float x1, float x2, float xacc);

 float rtflsp(float (*func)(float), float x1, float x2, float xacc);

 float rtnewt(void (*funcd)(float, float *, float *), float x1, float x2,

     float xacc);

 float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,

     float xacc);

 float rtsec(float (*func)(float), float x1, float x2, float xacc);

 void rzextr(int iest, float xest, float yest[], float yz[], float dy[], int nv);

 void savgol(float c[], int np, int nl, int nr, int ld, int m);

 void score(float xf, float y[], float f[]);

 void scrsho(float (*fx)(float));

 float select(unsigned long k, unsigned long n, float arr[]);

 float selip(unsigned long k, unsigned long n, float arr[]);

 void shell(unsigned long n, float a[]);

 void shoot(int n, float v[], float f[]);

 void shootf(int n, float v[], float f[]);

 void simp1(float **a, int mm, int ll[], int nll, int iabf, int *kp,

     float *bmax);

 void simp2(float **a, int m, int n, int *ip, int kp);

 void simp3(float **a, int i1, int k1, int ip, int kp);

 void simplx(float **a, int m, int n, int m1, int m2, int m3, int *icase,

     int izrov[], int iposv[]);

 void simpr(float y[], float dydx[], float dfdx[], float **dfdy,

     int n, float xs, float htot, int nstep, float yout[],

     void (*derivs)(float, float [], float []));

 void sinft(float y[], int n);

 void slvsm2(double **u, double **rhs);

 void slvsml(double **u, double **rhs);

 void sncndn(float uu, float emmc, float *sn, float *cn, float *dn);

 double snrm(unsigned long n, double sx[], int itol);

 void sobseq(int *n, float x[]);

 void solvde(int itmax, float conv, float slowc, float scalv[],

     int indexv[], int ne, int nb, int m, float **y, float ***c, float **s);

 void sor(double **a, double **b, double **c, double **d, double **e,

     double **f, double **u, int jmax, double rjac);

 void sort(unsigned long n, float arr[]);

 void sort2(unsigned long n, float arr[], float brr[]);

 void sort3(unsigned long n, float ra[], float rb[], float rc[]);

 void spctrm(FILE *fp, float p[], int m, int k, int ovrlap);

 void spear(float data1[], float data2[], unsigned long n, float *d, float *zd,

     float *probd, float *rs, float *probrs);

 void sphbes(int n, float x, float *sj, float *sy, float *sjp, float *syp);

 void splie2(float x1a[], float x2a[], float **ya, int m, int n, float **y2a);

 void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n,

     float x1, float x2, float *y);

 void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);

 void splint(float xa[], float ya[], float y2a[], int n, float x, float *y);

 void spread(float y, float yy[], unsigned long n, float x, int m);

 void sprsax(float sa[], unsigned long ija[], float x[], float b[],

     unsigned long n);

 void sprsin(float **a, int n, float thresh, unsigned long nmax, float sa[],

     unsigned long ija[]);

 void sprspm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],

     float sc[], unsigned long ijc[]);

 void sprstm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],

     float thresh, unsigned long nmax, float sc[], unsigned long ijc[]);

 void sprstp(float sa[], unsigned long ija[], float sb[], unsigned long ijb[]);

 void sprstx(float sa[], unsigned long ija[], float x[], float b[],

     unsigned long n);

 void stifbs(float y[], float dydx[], int nv, float *xx,

     float htry, float eps, float yscal[], float *hdid, float *hnext,

     void (*derivs)(float, float [], float []));

 void stiff(float y[], float dydx[], int n, float *x,

     float htry, float eps, float yscal[], float *hdid, float *hnext,

     void (*derivs)(float, float [], float []));

 void stoerm(float y[], float d2y[], int nv, float xs,

     float htot, int nstep, float yout[],

     void (*derivs)(float, float [], float []));

 void svbksb(float **u, float w[], float **v, int m, int n, float b[],

     float x[]);

 void svdcmp(float **a, int m, int n, float w[], float **v);

 void svdfit(float x[], float y[], float sig[], int ndata, float a[],

     int ma, float **u, float **v, float w[], float *chisq,

     void (*funcs)(float, float [], int));

 void svdvar(float **v, int ma, float w[], float **cvm);

 void toeplz(float r[], float x[], float y[], int n);

 void tptest(float data1[], float data2[], unsigned long n, float *t, float *prob);

 void tqli(float d[], float e[], int n, float **z);

 float trapzd(float (*func)(float), float a, float b, int n);

 void tred2(float **a, int n, float d[], float e[]);

 void tridag(float a[], float b[], float c[], float r[], float u[],

     unsigned long n);

 float trncst(float x[], float y[], int iorder[], int ncity, int n[]);

 void trnspt(int iorder[], int ncity, int n[]);

 void ttest(float data1[], unsigned long n1, float data2[], unsigned long n2,

     float *t, float *prob);

 void tutest(float data1[], unsigned long n1, float data2[], unsigned long n2,

     float *t, float *prob);

 void twofft(float data1[], float data2[], float fft1[], float fft2[],

     unsigned long n);

 void vander(double x[], double w[], double q[], int n);

 void vegas(float regn[], int ndim, float (*fxn)(float [], float), int init,

     unsigned long ncall, int itmx, int nprn, float *tgral, float *sd,

     float *chi2a);

 void voltra(int n, int m, float t0, float h, float *t, float **f,

     float (*g)(int, float), float (*ak)(int, int, float, float));

 void wt1(float a[], unsigned long n, int isign,

     void (*wtstep)(float [], unsigned long, int));

 void wtn(float a[], unsigned long nn[], int ndim, int isign,

     void (*wtstep)(float [], unsigned long, int));

 void wwghts(float wghts[], int n, float h,

     void (*kermom)(double [], double ,int));

 int zbrac(float (*func)(float), float *x1, float *x2);

 void zbrak(float (*fx)(float), float x1, float x2, int n, float xb1[],

     float xb2[], int *nb);

 float zbrent(float (*func)(float), float x1, float x2, float tol);

 void zrhqr(float a[], int m, float rtr[], float rti[]);

 float zriddr(float (*func)(float), float x1, float x2, float xacc);

 void zroots(fcomplex a[], int m, fcomplex roots[], int polish);

 #else /* ANSI */

 /* traditional - K&R */

 void addint();

 void airy();

 void amebsa();

 void amoeba();

 float amotry();

 float amotsa();

 void anneal();

 double anorm2();

 void arcmak();

 void arcode();

 void arcsum();

 void asolve();

 void atimes();

 void avevar();

 void balanc();

 void banbks();

 void bandec();

 void banmul();

 void bcucof();

 void bcuint();

 void beschb();

 float bessi();

 float bessi0();

 float bessi1();

 void bessik();

 float bessj();

 float bessj0();

 float bessj1();

 void bessjy();

 float bessk();

 float bessk0();

 float bessk1();

 float bessy();

 float bessy0();

 float bessy1();

 float beta();

 float betacf();

 float betai();

 float bico();

 void bksub();

 float bnldev();

 float brent();

 void broydn();

 void bsstep();

 void caldat();

 void chder();

 float chebev();

 void chebft();

 void chebpc();

 void chint();

 float chixy();

 void choldc();

 void cholsl();

 void chsone();

 void chstwo();

 void cisi();

 void cntab1();

 void cntab2();

 void convlv();

 void copy();

 void correl();

 void cosft();

 void cosft1();

 void cosft2();

 void covsrt();

 void crank();

 void cyclic();

 void daub4();

 float dawson();

 float dbrent();

 void ddpoly();

 int decchk();

 void derivs();

 float df1dim();

 void dfour1();

 void dfpmin();

 float dfridr();

 void dftcor();

 void dftint();

 void difeq();

 void dlinmin();

 double dpythag();

 void drealft();

 void dsprsax();

 void dsprstx();

 void dsvbksb();

 void dsvdcmp();

 void eclass();

 void eclazz();

 float ei();

 void eigsrt();

 float elle();

 float ellf();

 float ellpi();

 void elmhes();

 float erfcc();

 float erff();

 float erffc();

 void eulsum();

 float evlmem();

 float expdev();

 float expint();

 float f1();

 float f1dim();

 float f2();

 float f3();

 float factln();

 float factrl();

 void fasper();

 void fdjac();

 void fgauss();

 void fill0();

 void fit();

 void fitexy();

 void fixrts();

 void fleg();

 void flmoon();

 float fmin();

 void four1();

 void fourew();

 void fourfs();

 void fourn();

 void fpoly();

 void fred2();

 float fredin();

 void frenel();

 void frprmn();

 void ftest();

 float gamdev();

 float gammln();

 float gammp();

 float gammq();

 float gasdev();

 void gaucof();

 void gauher();

 void gaujac();

 void gaulag();

 void gauleg();

 void gaussj();

 void gcf();

 float golden();

 void gser();

 void hpsel();

 void hpsort();

 void hqr();

 void hufapp();

 void hufdec();

 void hufenc();

 void hufmak();

 void hunt();

 void hypdrv();

 fcomplex hypgeo();

 void hypser();

 unsigned short icrc();

 unsigned short icrc1();

 unsigned long igray();

 void iindexx();

 void indexx();

 void interp();

 int irbit1();

 int irbit2();

 void jacobi();

 void jacobn();

 long julday();

 void kendl1();

 void kendl2();

 void kermom();

 void ks2d1s();

 void ks2d2s();

 void ksone();

 void kstwo();

 void laguer();

 void lfit();

 void linbcg();

 void linmin();

 void lnsrch();

 void load();

 void load1();

 void load2();

 void locate();

 void lop();

 void lubksb();

 void ludcmp();

 void machar();

 void matadd();

 void matsub();

 void medfit();

 void memcof();

 int metrop();

 void mgfas();

 void mglin();

 float midexp();

 float midinf();

 float midpnt();

 float midsql();

 float midsqu();

 void miser();

 void mmid();

 void mnbrak();

 void mnewt();

 void moment();

 void mp2dfr();

 void mpadd();

 void mpdiv();

 void mpinv();

 void mplsh();

 void mpmov();

 void mpmul();

 void mpneg();

 void mppi();

 void mprove();

 void mpsad();

 void mpsdv();

 void mpsmu();

 void mpsqrt();

 void mpsub();

 void mrqcof();

 void mrqmin();

 void newt();

 void odeint();

 void orthog();

 void pade();

 void pccheb();

 void pcshft();

 void pearsn();

 void period();

 void piksr2();

 void piksrt();

 void pinvs();

 float plgndr();

 float poidev();

 void polcoe();

 void polcof();

 void poldiv();

 void polin2();

 void polint();

 void powell();

 void predic();

 float probks();

 void psdes();

 void pwt();

 void pwtset();

 float pythag();

 void pzextr();

 float qgaus();

 void qrdcmp();

 float qromb();

 float qromo();

 void qroot();

 void qrsolv();

 void qrupdt();

 float qsimp();

 float qtrap();

 float quad3d();

 void quadct();

 void quadmx();

 void quadvl();

 float ran0();

 float ran1();

 float ran2();

 float ran3();

 float ran4();

 void rank();

 void ranpt();

 void ratint();

 void ratlsq();

 double ratval();

 float rc();

 float rd();

 void realft();

 void rebin();

 void red();

 void relax();

 void relax2();

 void resid();

 float revcst();

 void reverse();

 float rf();

 float rj();

 void rk4();

 void rkck();

 void rkdumb();

 void rkqs();

 void rlft3();

 float rofunc();

 void rotate();

 void rsolv();

 void rstrct();

 float rtbis();

 float rtflsp();

 float rtnewt();

 float rtsafe();

 float rtsec();

 void rzextr();

 void savgol();

 void score();

 void scrsho();

 float select();

 float selip();

 void shell();

 void shoot();

 void shootf();

 void simp1();

 void simp2();

 void simp3();

 void simplx();

 void simpr();

 void sinft();

 void slvsm2();

 void slvsml();

 void sncndn();

 double snrm();

 void sobseq();

 void solvde();

 void sor();

 void sort();

 void sort2();

 void sort3();

 void spctrm();

 void spear();

 void sphbes();

 void splie2();

 void splin2();

 void spline();

 void splint();

 void spread();

 void sprsax();

 void sprsin();

 void sprspm();

 void sprstm();

 void sprstp();

 void sprstx();

 void stifbs();

 void stiff();

 void stoerm();

 void svbksb();

 void svdcmp();

 void svdfit();

 void svdvar();

 void toeplz();

 void tptest();

 void tqli();

 float trapzd();

 void tred2();

 void tridag();

 float trncst();

 void trnspt();

 void ttest();

 void tutest();

 void twofft();

 void vander();

 void vegas();

 void voltra();

 void wt1();

 void wtn();

 void wwghts();

 int zbrac();

 void zbrak();

 float zbrent();

 void zrhqr();

 float zriddr();

 void zroots();

 #endif /* ANSI */

 #endif /* _NR_H_ */

nr.h

 #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();

 }

 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++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++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++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++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++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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,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+,ncol=nch-ncl+;

     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=,j=nrl+;i<nrow;i++,j++) m[j]=m[j-]+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+,ncol=nch-ncl+,ndep=ndh-ndl+;

     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+;j<=nch;j++) t[nrl][j]=t[nrl][j-]+ndep;

     for(i=nrl+;i<=nrh;i++) {

         t[i]=t[i-]+ncol;

         t[i][ncl]=t[i-][ncl]+ncol*ndep;

         for(j=ncl+;j<=nch;j++) t[i][j]=t[i][j-]+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();

 }

 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++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++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++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++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++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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,ncol=nch-ncl+;

     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+;i<=nrh;i++) m[i]=m[i-]+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+,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+,ncol=nch-ncl+;

     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=,j=nrl+;i<nrow;i++,j++) m[j]=m[j-]+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+,ncol=nch-ncl+,ndep=ndh-ndl+;

     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+;j<=nch;j++) t[nrl][j]=t[nrl][j-]+ndep;

     for(i=nrl+;i<=nrh;i++) {

         t[i]=t[i-]+ncol;

         t[i][ncl]=t[i-][ncl]+ncol*ndep;

         for(j=ncl+;j<=nch;j++) t[i][j]=t[i][j-]+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.c

 /* Driver for routine spline */

 #include <stdio.h>

 #include <math.h>

 #define NRANSI

 #include "nr.h"

 #include "nrutil.h"

 #define N 20

 #define PI 3.1415926

 int main(void)

 {

     int i;

     float yp1,ypn,*x,*y,*y2;

     x=vector(,N);

     y=vector(,N);

     y2=vector(,N);

     printf("\nsecond-derivatives for sin(x) from 0 to pi\n");

     /* Generate array for interpolation */

     for (i=;i<=;i++) {

         x[i]=i*PI/N;

         y[i]=sin(x[i]);

     }

     /* calculate 2nd derivative with spline */

     yp1=cos(x[]);

     ypn=cos(x[N]);

     spline(x,y,N,yp1,ypn,y2);

     /* test result */

     printf("%23s %16s\n","spline","actual");

     printf("%11s %14s %16s\n","angle","2nd deriv","2nd deriv");

     for (i=;i<=N;i++)

         printf("%10.2f %16.6f %16.6f\n",x[i],y2[i],-sin(x[i]));

     free_vector(y2,,N);

     free_vector(y,,N);

     free_vector(x,,N);

     return ;

 }

 #undef NRANSI

xspline.c

c语言spline的更多相关文章

  1. R语言︱常用统计方法包&plus;机器学习包（名称、简介）

    一.一些函数包大汇总 转载于:http://www.dataguru.cn/thread-116761-1-1.html 时间上有点过期,下面的资料供大家参考基本的R包已经实现了传统多元统计的很多功能 ...

  2. D01-R语言基础学习

    R语言基础学习——D01 20190410内容纲要: 1.R的下载与安装 2.R包的安装与使用方法 (1)查看已安装的包 (2)查看是否安装过包 (3)安装包 (4)更新包 3.结果的重用 4.R处理 ...

  3. 【R】R语言常用函数

    R语言常用函数 基本 一.数据管理vector:向量 numeric:数值型向量 logical:逻辑型向量character:字符型向量 list:列表 data.frame:数据框c:连接为向量或 ...

  4. R语言画图小结

    本文以1950年到2010年期间我国的火灾统计数据为例,数据如下所示: (0)加载数据 data<-read.csv("E:\\MyDocument\\p\\Data\\1950~20 ...

  5. R语言常用包汇总

    转载于:https://blog.csdn.net/sinat_26917383/article/details/50651464?locationNum=2&fps=1 一.一些函数包大汇总 ...

  6. C语言 &&num;183&semi; 高精度加法

    问题描述 输入两个整数a和b,输出这两个整数的和.a和b都不超过100位. 算法描述 由于a和b都比较大,所以不能直接使用语言中的标准数据类型来存储.对于这种问题,一般使用数组来处理. 定义一个数组A ...

  7. Windows server 2012 添加中文语言包&lpar;英文转为中文&rpar;（离线）

    Windows server 2012 添加中文语言包(英文转为中文)(离线) 相关资料: 公司环境:亚马孙aws虚拟机 英文版Windows2012 中文SQL Server2012安装包,需要安装 ...

  8. iOS开发系列--Swift语言

    概述 Swift是苹果2014年推出的全新的编程语言,它继承了C语言.ObjC的特性,且克服了C语言的兼容性问题.Swift发展过程中不仅保留了ObjC很多语法特性,它也借鉴了多种现代化语言的特点,在 ...

  9. C语言 &&num;183&semi; Anagrams问题

    问题描述 Anagrams指的是具有如下特性的两个单词:在这两个单词当中,每一个英文字母(不区分大小写)所出现的次数都是相同的.例如,"Unclear"和"Nuclear ...

随机推荐

  1. 2&period;ASP&period;NET MVC 中使用Crystal Report水晶报表

    上一篇,介绍了怎么导出Excel文件,这篇文章介绍在ASP.NET MVC中使用水晶报表. 项目源码下载:https://github.com/caofangsheng93/CrystalReport ...

  2. Apache Kylin 部署之不完全指南

    1. 引言 Apache Kylin(麒麟)是由eBay开源的分布式分析引擎,提供Hadoop之上的SQL查询接口及多维分析(OLAP)能力以支持超大规模数据.底层存储用的是HBase,数据输入与cu ...

  3. 【python】判断字符串日期是否有效

    来源: http://www.jb51.net/article/66014.htm http://www.runoob.com/python/att-time-strptime.html 用time模 ...

  4. Java集合---LinkedList源码解析

    一.源码解析1. LinkedList类定义2.LinkedList数据结构原理3.私有属性4.构造方法5.元素添加add()及原理6.删除数据remove()7.数据获取get()8.数据复制clo ...

  5. Http的请求对象

    Servlet 客户端 HTTP 请求 当浏览器请求网页时,它会向 Web 服务器发送特定信息,这些信息不能被直接读取,因为这些信息是作为 HTTP 请求的头的一部分进行传输的.您可以查看 HTTP ...

  6. oracle中区分audit&lowbar;file&lowbar;dest&comma; background&lowbar;dump&lowbar;dest&comma; core&lowbar;dump&lowbar;dest&comma; user&lowbar;dump&lowbar;dest

    一般在$ORACLE_HOME\admin\{SID}目录下: audit_file_dest                 = /u01/app/oracle/admin/{SID}/adump ...

  7. Java中常见数据结构Map之HashMap

    之前很早就在博客中写过HashMap的一些东西: 彻底搞懂HashMap,HashTableConcurrentHashMap关联: http://www.cnblogs.com/wang-meng/ ...

  8. UNIX网络编程——TCP带外数据小结

    带外数据概念实际上时向接收端传送三个不同的信息:(1)发送端进入紧急模式这个事实.接收进程得以通知这个事实的手段不外乎SIGURG信号或select调用.本通知在发送进程发送带外字节后由发送端TCP立 ...

  9. spring boot中的约定优于配置

    Spring Boot并不是一个全新的框架,而是将已有的Spring组件整合起来. Spring Boot可以说是遵循约定优于配置这个理念产生的.它的特点是简单.快速和便捷. 既然遵循约定优于配置,则 ...

  10. Stall Reservations POJ - 3190（贪心）

    Oh those picky N (1 <= N <= 50,000) cows! They are so picky that each one will only be milked ...

相关文章