1 #define NRANSI 2 #include "nrutil.h" 3  4 void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]) 5 { 6     int i,k; 7     float p,qn,sig,un,*u; 8  9     u=vector(1,n-1);10     if (yp1 > 0.99e30)11         y2[1]=u[1]=0.0;12     else {13         y2[1] = -0.5;14         u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);15     }16     for (i=2;i<=n-1;i++) {17         sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);18         p=sig*y2[i-1]+2.0;19         y2[i]=(sig-1.0)/p;20         u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);21         u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;22     }23     if (ypn > 0.99e30)24         qn=un=0.0;25     else {26         qn=0.5;27         un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));28     }29     y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);30     for (k=n-1;k>=1;k--)31         y2[k]=y2[k]*y2[k+1]+u[k];32     free_vector(u,1,n-1);33 }34 #undef NRANSI

1 #ifndef _NR_H_  2 #define _NR_H_  3   4 #ifndef _FCOMPLEX_DECLARE_T_  5 typedef struct FCOMPLEX {
    float r,i;} fcomplex;  6 #define _FCOMPLEX_DECLARE_T_  7 #endif /* _FCOMPLEX_DECLARE_T_ */  8   9 #ifndef _ARITHCODE_DECLARE_T_ 10 typedef struct { 11     unsigned long *ilob,*iupb,*ncumfq,jdif,nc,minint,nch,ncum,nrad; 12 } arithcode; 13 #define _ARITHCODE_DECLARE_T_ 14 #endif /* _ARITHCODE_DECLARE_T_ */ 15  16 #ifndef _HUFFCODE_DECLARE_T_ 17 typedef struct { 18     unsigned long *icod,*ncod,*left,*right,nch,nodemax; 19 } huffcode; 20 #define _HUFFCODE_DECLARE_T_ 21 #endif /* _HUFFCODE_DECLARE_T_ */ 22  23 #include 
      
        24  25 #if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */ 26  27 void addint(double **uf, double **uc, double **res, int nf); 28 void airy(float x, float *ai, float *bi, float *aip, float *bip); 29 void amebsa(float **p, float y[], int ndim, float pb[],    float *yb, 30     float ftol, float (*funk)(float []), int *iter, float temptr); 31 void amoeba(float **p, float y[], int ndim, float ftol, 32     float (*funk)(float []), int *iter); 33 float amotry(float **p, float y[], float psum[], int ndim, 34     float (*funk)(float []), int ihi, float fac); 35 float amotsa(float **p, float y[], float psum[], int ndim, float pb[], 36     float *yb, float (*funk)(float []), int ihi, float *yhi, float fac); 37 void anneal(float x[], float y[], int iorder[], int ncity); 38 double anorm2(double **a, int n); 39 void arcmak(unsigned long nfreq[], unsigned long nchh, unsigned long nradd, 40     arithcode *acode); 41 void arcode(unsigned long *ich, unsigned char **codep, unsigned long *lcode, 42     unsigned long *lcd, int isign, arithcode *acode); 43 void arcsum(unsigned long iin[], unsigned long iout[], unsigned long ja, 44     int nwk, unsigned long nrad, unsigned long nc); 45 void asolve(unsigned long n, double b[], double x[], int itrnsp); 46 void atimes(unsigned long n, double x[], double r[], int itrnsp); 47 void avevar(float data[], unsigned long n, float *ave, float *var); 48 void balanc(float **a, int n); 49 void banbks(float **a, unsigned long n, int m1, int m2, float **al, 50     unsigned long indx[], float b[]); 51 void bandec(float **a, unsigned long n, int m1, int m2, float **al, 52     unsigned long indx[], float *d); 53 void banmul(float **a, unsigned long n, int m1, int m2, float x[], float b[]); 54 void bcucof(float y[], float y1[], float y2[], float y12[], float d1, 55     float d2, float **c); 56 void bcuint(float y[], float y1[], float y2[], float y12[], 57     float x1l, float x1u, float x2l, float x2u, float x1, 58     float x2, float *ansy, float *ansy1, float *ansy2); 59 void beschb(double x, double *gam1, double *gam2, double *gampl, 60     double *gammi); 61 float bessi(int n, float x); 62 float bessi0(float x); 63 float bessi1(float x); 64 void bessik(float x, float xnu, float *ri, float *rk, float *rip, 65     float *rkp); 66 float bessj(int n, float x); 67 float bessj0(float x); 68 float bessj1(float x); 69 void bessjy(float x, float xnu, float *rj, float *ry, float *rjp, 70     float *ryp); 71 float bessk(int n, float x); 72 float bessk0(float x); 73 float bessk1(float x); 74 float bessy(int n, float x); 75 float bessy0(float x); 76 float bessy1(float x); 77 float beta(float z, float w); 78 float betacf(float a, float b, float x); 79 float betai(float a, float b, float x); 80 float bico(int n, int k); 81 void bksub(int ne, int nb, int jf, int k1, int k2, float ***c); 82 float bnldev(float pp, int n, long *idum); 83 float brent(float ax, float bx, float cx, 84     float (*f)(float), float tol, float *xmin); 85 void broydn(float x[], int n, int *check, 86     void (*vecfunc)(int, float [], float [])); 87 void bsstep(float y[], float dydx[], int nv, float *xx, float htry, 88     float eps, float yscal[], float *hdid, float *hnext, 89     void (*derivs)(float, float [], float [])); 90 void caldat(long julian, int *mm, int *id, int *iyyy); 91 void chder(float a, float b, float c[], float cder[], int n); 92 float chebev(float a, float b, float c[], int m, float x); 93 void chebft(float a, float b, float c[], int n, float (*func)(float)); 94 void chebpc(float c[], float d[], int n); 95 void chint(float a, float b, float c[], float cint[], int n); 96 float chixy(float bang); 97 void choldc(float **a, int n, float p[]); 98 void cholsl(float **a, int n, float p[], float b[], float x[]); 99 void chsone(float bins[], float ebins[], int nbins, int knstrn,100     float *df, float *chsq, float *prob);101 void chstwo(float bins1[], float bins2[], int nbins, int knstrn,102     float *df, float *chsq, float *prob);103 void cisi(float x, float *ci, float *si);104 void cntab1(int **nn, int ni, int nj, float *chisq,105     float *df, float *prob, float *cramrv, float *ccc);106 void cntab2(int **nn, int ni, int nj, float *h, float *hx, float *hy,107     float *hygx, float *hxgy, float *uygx, float *uxgy, float *uxy);108 void convlv(float data[], unsigned long n, float respns[], unsigned long m,109     int isign, float ans[]);110 void copy(double **aout, double **ain, int n);111 void correl(float data1[], float data2[], unsigned long n, float ans[]);112 void cosft(float y[], int n, int isign);113 void cosft1(float y[], int n);114 void cosft2(float y[], int n, int isign);115 void covsrt(float **covar, int ma, int ia[], int mfit);116 void crank(unsigned long n, float w[], float *s);117 void cyclic(float a[], float b[], float c[], float alpha, float beta,118     float r[], float x[], unsigned long n);119 void daub4(float a[], unsigned long n, int isign);120 float dawson(float x);121 float dbrent(float ax, float bx, float cx,122     float (*f)(float), float (*df)(float), float tol, float *xmin);123 void ddpoly(float c[], int nc, float x, float pd[], int nd);124 int decchk(char string[], int n, char *ch);125 void derivs(float x, float y[], float dydx[]);126 float df1dim(float x);127 void dfour1(double data[], unsigned long nn, int isign);128 void dfpmin(float p[], int n, float gtol, int *iter, float *fret,129     float (*func)(float []), void (*dfunc)(float [], float []));130 float dfridr(float (*func)(float), float x, float h, float *err);131 void dftcor(float w, float delta, float a, float b, float endpts[],132     float *corre, float *corim, float *corfac);133 void dftint(float (*func)(float), float a, float b, float w,134     float *cosint, float *sinint);135 void difeq(int k, int k1, int k2, int jsf, int is1, int isf,136     int indexv[], int ne, float **s, float **y);137 void dlinmin(float p[], float xi[], int n, float *fret,138     float (*func)(float []), void (*dfunc)(float [], float[]));139 double dpythag(double a, double b);140 void drealft(double data[], unsigned long n, int isign);141 void dsprsax(double sa[], unsigned long ija[], double x[], double b[],142     unsigned long n);143 void dsprstx(double sa[], unsigned long ija[], double x[], double b[],144     unsigned long n);145 void dsvbksb(double **u, double w[], double **v, int m, int n, double b[],146     double x[]);147 void dsvdcmp(double **a, int m, int n, double w[], double **v);148 void eclass(int nf[], int n, int lista[], int listb[], int m);149 void eclazz(int nf[], int n, int (*equiv)(int, int));150 float ei(float x);151 void eigsrt(float d[], float **v, int n);152 float elle(float phi, float ak);153 float ellf(float phi, float ak);154 float ellpi(float phi, float en, float ak);155 void elmhes(float **a, int n);156 float erfcc(float x);157 float erff(float x);158 float erffc(float x);159 void eulsum(float *sum, float term, int jterm, float wksp[]);160 float evlmem(float fdt, float d[], int m, float xms);161 float expdev(long *idum);162 float expint(int n, float x);163 float f1(float x);164 float f1dim(float x);165 float f2(float y);166 float f3(float z);167 float factln(int n);168 float factrl(int n);169 void fasper(float x[], float y[], unsigned long n, float ofac, float hifac,170     float wk1[], float wk2[], unsigned long nwk, unsigned long *nout,171     unsigned long *jmax, float *prob);172 void fdjac(int n, float x[], float fvec[], float **df,173     void (*vecfunc)(int, float [], float []));174 void fgauss(float x, float a[], float *y, float dyda[], int na);175 void fill0(double **u, int n);176 void fit(float x[], float y[], int ndata, float sig[], int mwt,177     float *a, float *b, float *siga, float *sigb, float *chi2, float *q);178 void fitexy(float x[], float y[], int ndat, float sigx[], float sigy[],179     float *a, float *b, float *siga, float *sigb, float *chi2, float *q);180 void fixrts(float d[], int m);181 void fleg(float x, float pl[], int nl);182 void flmoon(int n, int nph, long *jd, float *frac);183 float fmin(float x[]);184 void four1(float data[], unsigned long nn, int isign);185 void fourew(FILE *file[5], int *na, int *nb, int *nc, int *nd);186 void fourfs(FILE *file[5], unsigned long nn[], int ndim, int isign);187 void fourn(float data[], unsigned long nn[], int ndim, int isign);188 void fpoly(float x, float p[], int np);189 void fred2(int n, float a, float b, float t[], float f[], float w[],190     float (*g)(float), float (*ak)(float, float));191 float fredin(float x, int n, float a, float b, float t[], float f[], float w[],192     float (*g)(float), float (*ak)(float, float));193 void frenel(float x, float *s, float *c);194 void frprmn(float p[], int n, float ftol, int *iter, float *fret,195     float (*func)(float []), void (*dfunc)(float [], float []));196 void ftest(float data1[], unsigned long n1, float data2[], unsigned long n2,197     float *f, float *prob);198 float gamdev(int ia, long *idum);199 float gammln(float xx);200 float gammp(float a, float x);201 float gammq(float a, float x);202 float gasdev(long *idum);203 void gaucof(int n, float a[], float b[], float amu0, float x[], float w[]);204 void gauher(float x[], float w[], int n);205 void gaujac(float x[], float w[], int n, float alf, float bet);206 void gaulag(float x[], float w[], int n, float alf);207 void gauleg(float x1, float x2, float x[], float w[], int n);208 void gaussj(float **a, int n, float **b, int m);209 void gcf(float *gammcf, float a, float x, float *gln);210 float golden(float ax, float bx, float cx, float (*f)(float), float tol,211     float *xmin);212 void gser(float *gamser, float a, float x, float *gln);213 void hpsel(unsigned long m, unsigned long n, float arr[], float heap[]);214 void hpsort(unsigned long n, float ra[]);215 void hqr(float **a, int n, float wr[], float wi[]);216 void hufapp(unsigned long index[], unsigned long nprob[], unsigned long n,217     unsigned long i);218 void hufdec(unsigned long *ich, unsigned char *code, unsigned long lcode,219     unsigned long *nb, huffcode *hcode);220 void hufenc(unsigned long ich, unsigned char **codep, unsigned long *lcode,221     unsigned long *nb, huffcode *hcode);222 void hufmak(unsigned long nfreq[], unsigned long nchin, unsigned long *ilong,223     unsigned long *nlong, huffcode *hcode);224 void hunt(float xx[], unsigned long n, float x, unsigned long *jlo);225 void hypdrv(float s, float yy[], float dyyds[]);226 fcomplex hypgeo(fcomplex a, fcomplex b, fcomplex c, fcomplex z);227 void hypser(fcomplex a, fcomplex b, fcomplex c, fcomplex z,228     fcomplex *series, fcomplex *deriv);229 unsigned short icrc(unsigned short crc, unsigned char *bufptr,230     unsigned long len, short jinit, int jrev);231 unsigned short icrc1(unsigned short crc, unsigned char onech);232 unsigned long igray(unsigned long n, int is);233 void iindexx(unsigned long n, long arr[], unsigned long indx[]);234 void indexx(unsigned long n, float arr[], unsigned long indx[]);235 void interp(double **uf, double **uc, int nf);236 int irbit1(unsigned long *iseed);237 int irbit2(unsigned long *iseed);238 void jacobi(float **a, int n, float d[], float **v, int *nrot);239 void jacobn(float x, float y[], float dfdx[], float **dfdy, int n);240 long julday(int mm, int id, int iyyy);241 void kendl1(float data1[], float data2[], unsigned long n, float *tau, float *z,242     float *prob);243 void kendl2(float **tab, int i, int j, float *tau, float *z, float *prob);244 void kermom(double w[], double y, int m);245 void ks2d1s(float x1[], float y1[], unsigned long n1,246     void (*quadvl)(float, float, float *, float *, float *, float *),247     float *d1, float *prob);248 void ks2d2s(float x1[], float y1[], unsigned long n1, float x2[], float y2[],249     unsigned long n2, float *d, float *prob);250 void ksone(float data[], unsigned long n, float (*func)(float), float *d,251     float *prob);252 void kstwo(float data1[], unsigned long n1, float data2[], unsigned long n2,253     float *d, float *prob);254 void laguer(fcomplex a[], int m, fcomplex *x, int *its);255 void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],256     int ma, float **covar, float *chisq, void (*funcs)(float, float [], int));257 void linbcg(unsigned long n, double b[], double x[], int itol, double tol,258      int itmax, int *iter, double *err);259 void linmin(float p[], float xi[], int n, float *fret,260     float (*func)(float []));261 void lnsrch(int n, float xold[], float fold, float g[], float p[], float x[],262      float *f, float stpmax, int *check, float (*func)(float []));263 void load(float x1, float v[], float y[]);264 void load1(float x1, float v1[], float y[]);265 void load2(float x2, float v2[], float y[]);266 void locate(float xx[], unsigned long n, float x, unsigned long *j);267 void lop(double **out, double **u, int n);268 void lubksb(float **a, int n, int *indx, float b[]);269 void ludcmp(float **a, int n, int *indx, float *d);270 void machar(int *ibeta, int *it, int *irnd, int *ngrd,271     int *machep, int *negep, int *iexp, int *minexp, int *maxexp,272     float *eps, float *epsneg, float *xmin, float *xmax);273 void matadd(double **a, double **b, double **c, int n);274 void matsub(double **a, double **b, double **c, int n);275 void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev);276 void memcof(float data[], int n, int m, float *xms, float d[]);277 int metrop(float de, float t);278 void mgfas(double **u, int n, int maxcyc);279 void mglin(double **u, int n, int ncycle);280 float midexp(float (*funk)(float), float aa, float bb, int n);281 float midinf(float (*funk)(float), float aa, float bb, int n);282 float midpnt(float (*func)(float), float a, float b, int n);283 float midsql(float (*funk)(float), float aa, float bb, int n);284 float midsqu(float (*funk)(float), float aa, float bb, int n);285 void miser(float (*func)(float []), float regn[], int ndim, unsigned long npts,286     float dith, float *ave, float *var);287 void mmid(float y[], float dydx[], int nvar, float xs, float htot,288     int nstep, float yout[], void (*derivs)(float, float[], float[]));289 void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,290     float *fc, float (*func)(float));291 void mnewt(int ntrial, float x[], int n, float tolx, float tolf);292 void moment(float data[], int n, float *ave, float *adev, float *sdev,293     float *var, float *skew, float *curt);294 void mp2dfr(unsigned char a[], unsigned char s[], int n, int *m);295 void mpadd(unsigned char w[], unsigned char u[], unsigned char v[], int n);296 void mpdiv(unsigned char q[], unsigned char r[], unsigned char u[],297     unsigned char v[], int n, int m);298 void mpinv(unsigned char u[], unsigned char v[], int n, int m);299 void mplsh(unsigned char u[], int n);300 void mpmov(unsigned char u[], unsigned char v[], int n);301 void mpmul(unsigned char w[], unsigned char u[], unsigned char v[], int n,302     int m);303 void mpneg(unsigned char u[], int n);304 void mppi(int n);305 void mprove(float **a, float **alud, int n, int indx[], float b[],306     float x[]);307 void mpsad(unsigned char w[], unsigned char u[], int n, int iv);308 void mpsdv(unsigned char w[], unsigned char u[], int n, int iv, int *ir);309 void mpsmu(unsigned char w[], unsigned char u[], int n, int iv);310 void mpsqrt(unsigned char w[], unsigned char u[], unsigned char v[], int n,311     int m);312 void mpsub(int *is, unsigned char w[], unsigned char u[], unsigned char v[],313     int n);314 void mrqcof(float x[], float y[], float sig[], int ndata, float a[],315     int ia[], int ma, float **alpha, float beta[], float *chisq,316     void (*funcs)(float, float [], float *, float [], int));317 void mrqmin(float x[], float y[], float sig[], int ndata, float a[],318     int ia[], int ma, float **covar, float **alpha, float *chisq,319     void (*funcs)(float, float [], float *, float [], int), float *alamda);320 void newt(float x[], int n, int *check,321     void (*vecfunc)(int, float [], float []));322 void odeint(float ystart[], int nvar, float x1, float x2,323     float eps, float h1, float hmin, int *nok, int *nbad,324     void (*derivs)(float, float [], float []),325     void (*rkqs)(float [], float [], int, float *, float, float,326     float [], float *, float *, void (*)(float, float [], float [])));327 void orthog(int n, float anu[], float alpha[], float beta[], float a[],328     float b[]);329 void pade(double cof[], int n, float *resid);330 void pccheb(float d[], float c[], int n);331 void pcshft(float a, float b, float d[], int n);332 void pearsn(float x[], float y[], unsigned long n, float *r, float *prob,333     float *z);334 void period(float x[], float y[], int n, float ofac, float hifac,335     float px[], float py[], int np, int *nout, int *jmax, float *prob);336 void piksr2(int n, float arr[], float brr[]);337 void piksrt(int n, float arr[]);338 void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k,339     float ***c, float **s);340 float plgndr(int l, int m, float x);341 float poidev(float xm, long *idum);342 void polcoe(float x[], float y[], int n, float cof[]);343 void polcof(float xa[], float ya[], int n, float cof[]);344 void poldiv(float u[], int n, float v[], int nv, float q[], float r[]);345 void polin2(float x1a[], float x2a[], float **ya, int m, int n,346     float x1, float x2, float *y, float *dy);347 void polint(float xa[], float ya[], int n, float x, float *y, float *dy);348 void powell(float p[], float **xi, int n, float ftol, int *iter, float *fret,349     float (*func)(float []));350 void predic(float data[], int ndata, float d[], int m, float future[], int nfut);351 float probks(float alam);352 void psdes(unsigned long *lword, unsigned long *irword);353 void pwt(float a[], unsigned long n, int isign);354 void pwtset(int n);355 float pythag(float a, float b);356 void pzextr(int iest, float xest, float yest[], float yz[], float dy[],357     int nv);358 float qgaus(float (*func)(float), float a, float b);359 void qrdcmp(float **a, int n, float *c, float *d, int *sing);360 float qromb(float (*func)(float), float a, float b);361 float qromo(float (*func)(float), float a, float b,362     float (*choose)(float (*)(float), float, float, int));363 void qroot(float p[], int n, float *b, float *c, float eps);364 void qrsolv(float **a, int n, float c[], float d[], float b[]);365 void qrupdt(float **r, float **qt, int n, float u[], float v[]);366 float qsimp(float (*func)(float), float a, float b);367 float qtrap(float (*func)(float), float a, float b);368 float quad3d(float (*func)(float, float, float), float x1, float x2);369 void quadct(float x, float y, float xx[], float yy[], unsigned long nn,370     float *fa, float *fb, float *fc, float *fd);371 void quadmx(float **a, int n);372 void quadvl(float x, float y, float *fa, float *fb, float *fc, float *fd);373 float ran0(long *idum);374 float ran1(long *idum);375 float ran2(long *idum);376 float ran3(long *idum);377 float ran4(long *idum);378 void rank(unsigned long n, unsigned long indx[], unsigned long irank[]);379 void ranpt(float pt[], float regn[], int n);380 void ratint(float xa[], float ya[], int n, float x, float *y, float *dy);381 void ratlsq(double (*fn)(double), double a, double b, int mm, int kk,382     double cof[], double *dev);383 double ratval(double x, double cof[], int mm, int kk);384 float rc(float x, float y);385 float rd(float x, float y, float z);386 void realft(float data[], unsigned long n, int isign);387 void rebin(float rc, int nd, float r[], float xin[], float xi[]);388 void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf,389     int ic1, int jc1, int jcf, int kc, float ***c, float **s);390 void relax(double **u, double **rhs, int n);391 void relax2(double **u, double **rhs, int n);392 void resid(double **res, double **u, double **rhs, int n);393 float revcst(float x[], float y[], int iorder[], int ncity, int n[]);394 void reverse(int iorder[], int ncity, int n[]);395 float rf(float x, float y, float z);396 float rj(float x, float y, float z, float p);397 void rk4(float y[], float dydx[], int n, float x, float h, float yout[],398     void (*derivs)(float, float [], float []));399 void rkck(float y[], float dydx[], int n, float x, float h,400     float yout[], float yerr[], void (*derivs)(float, float [], float []));401 void rkdumb(float vstart[], int nvar, float x1, float x2, int nstep,402     void (*derivs)(float, float [], float []));403 void rkqs(float y[], float dydx[], int n, float *x,404     float htry, float eps, float yscal[], float *hdid, float *hnext,405     void (*derivs)(float, float [], float []));406 void rlft3(float ***data, float **speq, unsigned long nn1,407     unsigned long nn2, unsigned long nn3, int isign);408 float rofunc(float b);409 void rotate(float **r, float **qt, int n, int i, float a, float b);410 void rsolv(float **a, int n, float d[], float b[]);411 void rstrct(double **uc, double **uf, int nc);412 float rtbis(float (*func)(float), float x1, float x2, float xacc);413 float rtflsp(float (*func)(float), float x1, float x2, float xacc);414 float rtnewt(void (*funcd)(float, float *, float *), float x1, float x2,415     float xacc);416 float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,417     float xacc);418 float rtsec(float (*func)(float), float x1, float x2, float xacc);419 void rzextr(int iest, float xest, float yest[], float yz[], float dy[], int nv);420 void savgol(float c[], int np, int nl, int nr, int ld, int m);421 void score(float xf, float y[], float f[]);422 void scrsho(float (*fx)(float));423 float select(unsigned long k, unsigned long n, float arr[]);424 float selip(unsigned long k, unsigned long n, float arr[]);425 void shell(unsigned long n, float a[]);426 void shoot(int n, float v[], float f[]);427 void shootf(int n, float v[], float f[]);428 void simp1(float **a, int mm, int ll[], int nll, int iabf, int *kp,429     float *bmax);430 void simp2(float **a, int m, int n, int *ip, int kp);431 void simp3(float **a, int i1, int k1, int ip, int kp);432 void simplx(float **a, int m, int n, int m1, int m2, int m3, int *icase,433     int izrov[], int iposv[]);434 void simpr(float y[], float dydx[], float dfdx[], float **dfdy,435     int n, float xs, float htot, int nstep, float yout[],436     void (*derivs)(float, float [], float []));437 void sinft(float y[], int n);438 void slvsm2(double **u, double **rhs);439 void slvsml(double **u, double **rhs);440 void sncndn(float uu, float emmc, float *sn, float *cn, float *dn);441 double snrm(unsigned long n, double sx[], int itol);442 void sobseq(int *n, float x[]);443 void solvde(int itmax, float conv, float slowc, float scalv[],444     int indexv[], int ne, int nb, int m, float **y, float ***c, float **s);445 void sor(double **a, double **b, double **c, double **d, double **e,446     double **f, double **u, int jmax, double rjac);447 void sort(unsigned long n, float arr[]);448 void sort2(unsigned long n, float arr[], float brr[]);449 void sort3(unsigned long n, float ra[], float rb[], float rc[]);450 void spctrm(FILE *fp, float p[], int m, int k, int ovrlap);451 void spear(float data1[], float data2[], unsigned long n, float *d, float *zd,452     float *probd, float *rs, float *probrs);453 void sphbes(int n, float x, float *sj, float *sy, float *sjp, float *syp);454 void splie2(float x1a[], float x2a[], float **ya, int m, int n, float **y2a);455 void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n,456     float x1, float x2, float *y);457 void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);458 void splint(float xa[], float ya[], float y2a[], int n, float x, float *y);459 void spread(float y, float yy[], unsigned long n, float x, int m);460 void sprsax(float sa[], unsigned long ija[], float x[], float b[],461     unsigned long n);462 void sprsin(float **a, int n, float thresh, unsigned long nmax, float sa[],463     unsigned long ija[]);464 void sprspm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],465     float sc[], unsigned long ijc[]);466 void sprstm(float sa[], unsigned long ija[], float sb[], unsigned long ijb[],467     float thresh, unsigned long nmax, float sc[], unsigned long ijc[]);468 void sprstp(float sa[], unsigned long ija[], float sb[], unsigned long ijb[]);469 void sprstx(float sa[], unsigned long ija[], float x[], float b[],470     unsigned long n);471 void stifbs(float y[], float dydx[], int nv, float *xx,472     float htry, float eps, float yscal[], float *hdid, float *hnext,473     void (*derivs)(float, float [], float []));474 void stiff(float y[], float dydx[], int n, float *x,475     float htry, float eps, float yscal[], float *hdid, float *hnext,476     void (*derivs)(float, float [], float []));477 void stoerm(float y[], float d2y[], int nv, float xs,478     float htot, int nstep, float yout[],479     void (*derivs)(float, float [], float []));480 void svbksb(float **u, float w[], float **v, int m, int n, float b[],481     float x[]);482 void svdcmp(float **a, int m, int n, float w[], float **v);483 void svdfit(float x[], float y[], float sig[], int ndata, float a[],484     int ma, float **u, float **v, float w[], float *chisq,485     void (*funcs)(float, float [], int));486 void svdvar(float **v, int ma, float w[], float **cvm);487 void toeplz(float r[], float x[], float y[], int n);488 void tptest(float data1[], float data2[], unsigned long n, float *t, float *prob);489 void tqli(float d[], float e[], int n, float **z);490 float trapzd(float (*func)(float), float a, float b, int n);491 void tred2(float **a, int n, float d[], float e[]);492 void tridag(float a[], float b[], float c[], float r[], float u[],493     unsigned long n);494 float trncst(float x[], float y[], int iorder[], int ncity, int n[]);495 void trnspt(int iorder[], int ncity, int n[]);496 void ttest(float data1[], unsigned long n1, float data2[], unsigned long n2,497     float *t, float *prob);498 void tutest(float data1[], unsigned long n1, float data2[], unsigned long n2,499     float *t, float *prob);500 void twofft(float data1[], float data2[], float fft1[], float fft2[],501     unsigned long n);502 void vander(double x[], double w[], double q[], int n);503 void vegas(float regn[], int ndim, float (*fxn)(float [], float), int init,504     unsigned long ncall, int itmx, int nprn, float *tgral, float *sd,505     float *chi2a);506 void voltra(int n, int m, float t0, float h, float *t, float **f,507     float (*g)(int, float), float (*ak)(int, int, float, float));508 void wt1(float a[], unsigned long n, int isign,509     void (*wtstep)(float [], unsigned long, int));510 void wtn(float a[], unsigned long nn[], int ndim, int isign,511     void (*wtstep)(float [], unsigned long, int));512 void wwghts(float wghts[], int n, float h,513     void (*kermom)(double [], double ,int));514 int zbrac(float (*func)(float), float *x1, float *x2);515 void zbrak(float (*fx)(float), float x1, float x2, int n, float xb1[],516     float xb2[], int *nb);517 float zbrent(float (*func)(float), float x1, float x2, float tol);518 void zrhqr(float a[], int m, float rtr[], float rti[]);519 float zriddr(float (*func)(float), float x1, float x2, float xacc);520 void zroots(fcomplex a[], int m, fcomplex roots[], int polish);521 522 #else /* ANSI */523 /* traditional - K&R */524 525 void addint();526 void airy();527 void amebsa();528 void amoeba();529 float amotry();530 float amotsa();531 void anneal();532 double anorm2();533 void arcmak();534 void arcode();535 void arcsum();536 void asolve();537 void atimes();538 void avevar();539 void balanc();540 void banbks();541 void bandec();542 void banmul();543 void bcucof();544 void bcuint();545 void beschb();546 float bessi();547 float bessi0();548 float bessi1();549 void bessik();550 float bessj();551 float bessj0();552 float bessj1();553 void bessjy();554 float bessk();555 float bessk0();556 float bessk1();557 float bessy();558 float bessy0();559 float bessy1();560 float beta();561 float betacf();562 float betai();563 float bico();564 void bksub();565 float bnldev();566 float brent();567 void broydn();568 void bsstep();569 void caldat();570 void chder();571 float chebev();572 void chebft();573 void chebpc();574 void chint();575 float chixy();576 void choldc();577 void cholsl();578 void chsone();579 void chstwo();580 void cisi();581 void cntab1();582 void cntab2();583 void convlv();584 void copy();585 void correl();586 void cosft();587 void cosft1();588 void cosft2();589 void covsrt();590 void crank();591 void cyclic();592 void daub4();593 float dawson();594 float dbrent();595 void ddpoly();596 int decchk();597 void derivs();598 float df1dim();599 void dfour1();600 void dfpmin();601 float dfridr();602 void dftcor();603 void dftint();604 void difeq();605 void dlinmin();606 double dpythag();607 void drealft();608 void dsprsax();609 void dsprstx();610 void dsvbksb();611 void dsvdcmp();612 void eclass();613 void eclazz();614 float ei();615 void eigsrt();616 float elle();617 float ellf();618 float ellpi();619 void elmhes();620 float erfcc();621 float erff();622 float erffc();623 void eulsum();624 float evlmem();625 float expdev();626 float expint();627 float f1();628 float f1dim();629 float f2();630 float f3();631 float factln();632 float factrl();633 void fasper();634 void fdjac();635 void fgauss();636 void fill0();637 void fit();638 void fitexy();639 void fixrts();640 void fleg();641 void flmoon();642 float fmin();643 void four1();644 void fourew();645 void fourfs();646 void fourn();647 void fpoly();648 void fred2();649 float fredin();650 void frenel();651 void frprmn();652 void ftest();653 float gamdev();654 float gammln();655 float gammp();656 float gammq();657 float gasdev();658 void gaucof();659 void gauher();660 void gaujac();661 void gaulag();662 void gauleg();663 void gaussj();664 void gcf();665 float golden();666 void gser();667 void hpsel();668 void hpsort();669 void hqr();670 void hufapp();671 void hufdec();672 void hufenc();673 void hufmak();674 void hunt();675 void hypdrv();676 fcomplex hypgeo();677 void hypser();678 unsigned short icrc();679 unsigned short icrc1();680 unsigned long igray();681 void iindexx();682 void indexx();683 void interp();684 int irbit1();685 int irbit2();686 void jacobi();687 void jacobn();688 long julday();689 void kendl1();690 void kendl2();691 void kermom();692 void ks2d1s();693 void ks2d2s();694 void ksone();695 void kstwo();696 void laguer();697 void lfit();698 void linbcg();699 void linmin();700 void lnsrch();701 void load();702 void load1();703 void load2();704 void locate();705 void lop();706 void lubksb();707 void ludcmp();708 void machar();709 void matadd();710 void matsub();711 void medfit();712 void memcof();713 int metrop();714 void mgfas();715 void mglin();716 float midexp();717 float midinf();718 float midpnt();719 float midsql();720 float midsqu();721 void miser();722 void mmid();723 void mnbrak();724 void mnewt();725 void moment();726 void mp2dfr();727 void mpadd();728 void mpdiv();729 void mpinv();730 void mplsh();731 void mpmov();732 void mpmul();733 void mpneg();734 void mppi();735 void mprove();736 void mpsad();737 void mpsdv();738 void mpsmu();739 void mpsqrt();740 void mpsub();741 void mrqcof();742 void mrqmin();743 void newt();744 void odeint();745 void orthog();746 void pade();747 void pccheb();748 void pcshft();749 void pearsn();750 void period();751 void piksr2();752 void piksrt();753 void pinvs();754 float plgndr();755 float poidev();756 void polcoe();757 void polcof();758 void poldiv();759 void polin2();760 void polint();761 void powell();762 void predic();763 float probks();764 void psdes();765 void pwt();766 void pwtset();767 float pythag();768 void pzextr();769 float qgaus();770 void qrdcmp();771 float qromb();772 float qromo();773 void qroot();774 void qrsolv();775 void qrupdt();776 float qsimp();777 float qtrap();778 float quad3d();779 void quadct();780 void quadmx();781 void quadvl();782 float ran0();783 float ran1();784 float ran2();785 float ran3();786 float ran4();787 void rank();788 void ranpt();789 void ratint();790 void ratlsq();791 double ratval();792 float rc();793 float rd();794 void realft();795 void rebin();796 void red();797 void relax();798 void relax2();799 void resid();800 float revcst();801 void reverse();802 float rf();803 float rj();804 void rk4();805 void rkck();806 void rkdumb();807 void rkqs();808 void rlft3();809 float rofunc();810 void rotate();811 void rsolv();812 void rstrct();813 float rtbis();814 float rtflsp();815 float rtnewt();816 float rtsafe();817 float rtsec();818 void rzextr();819 void savgol();820 void score();821 void scrsho();822 float select();823 float selip();824 void shell();825 void shoot();826 void shootf();827 void simp1();828 void simp2();829 void simp3();830 void simplx();831 void simpr();832 void sinft();833 void slvsm2();834 void slvsml();835 void sncndn();836 double snrm();837 void sobseq();838 void solvde();839 void sor();840 void sort();841 void sort2();842 void sort3();843 void spctrm();844 void spear();845 void sphbes();846 void splie2();847 void splin2();848 void spline();849 void splint();850 void spread();851 void sprsax();852 void sprsin();853 void sprspm();854 void sprstm();855 void sprstp();856 void sprstx();857 void stifbs();858 void stiff();859 void stoerm();860 void svbksb();861 void svdcmp();862 void svdfit();863 void svdvar();864 void toeplz();865 void tptest();866 void tqli();867 float trapzd();868 void tred2();869 void tridag();870 float trncst();871 void trnspt();872 void ttest();873 void tutest();874 void twofft();875 void vander();876 void vegas();877 void voltra();878 void wt1();879 void wtn();880 void wwghts();881 int zbrac();882 void zbrak();883 float zbrent();884 void zrhqr();885 float zriddr();886 void zroots();887 888 #endif /* ANSI */889 890 #endif /* _NR_H_ */
      

 

1 #if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */  2   3 #include 
      
         4 #include 
       
          5 #include 
        
           6 #define NR_END 1  7 #define FREE_ARG char*  8   9 void nrerror(char error_text[]) 10 /* Numerical Recipes standard error handler */ 11 { 12     fprintf(stderr,"Numerical Recipes run-time error...\n"); 13     fprintf(stderr,"%s\n",error_text); 14     fprintf(stderr,"...now exiting to system...\n"); 15     exit(1); 16 } 17  18 float *vector(long nl, long nh) 19 /* allocate a float vector with subscript range v[nl..nh] */ 20 { 21     float *v; 22  23     v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); 24     if (!v) nrerror("allocation failure in vector()"); 25     return v-nl+NR_END; 26 } 27  28 int *ivector(long nl, long nh) 29 /* allocate an int vector with subscript range v[nl..nh] */ 30 { 31     int *v; 32  33     v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); 34     if (!v) nrerror("allocation failure in ivector()"); 35     return v-nl+NR_END; 36 } 37  38 unsigned char *cvector(long nl, long nh) 39 /* allocate an unsigned char vector with subscript range v[nl..nh] */ 40 { 41     unsigned char *v; 42  43     v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); 44     if (!v) nrerror("allocation failure in cvector()"); 45     return v-nl+NR_END; 46 } 47  48 unsigned long *lvector(long nl, long nh) 49 /* allocate an unsigned long vector with subscript range v[nl..nh] */ 50 { 51     unsigned long *v; 52  53     v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); 54     if (!v) nrerror("allocation failure in lvector()"); 55     return v-nl+NR_END; 56 } 57  58 double *dvector(long nl, long nh) 59 /* allocate a double vector with subscript range v[nl..nh] */ 60 { 61     double *v; 62  63     v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); 64     if (!v) nrerror("allocation failure in dvector()"); 65     return v-nl+NR_END; 66 } 67  68 float **matrix(long nrl, long nrh, long ncl, long nch) 69 /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ 70 { 71     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 72     float **m; 73  74     /* allocate pointers to rows */ 75     m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); 76     if (!m) nrerror("allocation failure 1 in matrix()"); 77     m += NR_END; 78     m -= nrl; 79  80     /* allocate rows and set pointers to them */ 81     m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); 82     if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); 83     m[nrl] += NR_END; 84     m[nrl] -= ncl; 85  86     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; 87  88     /* return pointer to array of pointers to rows */ 89     return m; 90 } 91  92 double **dmatrix(long nrl, long nrh, long ncl, long nch) 93 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ 94 { 95     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; 96     double **m; 97  98     /* allocate pointers to rows */ 99     m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));100     if (!m) nrerror("allocation failure 1 in matrix()");101     m += NR_END;102     m -= nrl;103 104     /* allocate rows and set pointers to them */105     m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));106     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");107     m[nrl] += NR_END;108     m[nrl] -= ncl;109 110     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;111 112     /* return pointer to array of pointers to rows */113     return m;114 }115 116 int **imatrix(long nrl, long nrh, long ncl, long nch)117 /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */118 {119     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;120     int **m;121 122     /* allocate pointers to rows */123     m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));124     if (!m) nrerror("allocation failure 1 in matrix()");125     m += NR_END;126     m -= nrl;127 128 129     /* allocate rows and set pointers to them */130     m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));131     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");132     m[nrl] += NR_END;133     m[nrl] -= ncl;134 135     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;136 137     /* return pointer to array of pointers to rows */138     return m;139 }140 141 float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,142     long newrl, long newcl)143 /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */144 {145     long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;146     float **m;147 148     /* allocate array of pointers to rows */149     m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));150     if (!m) nrerror("allocation failure in submatrix()");151     m += NR_END;152     m -= newrl;153 154     /* set pointers to rows */155     for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;156 157     /* return pointer to array of pointers to rows */158     return m;159 }160 161 float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)162 /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix163 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1164 and ncol=nch-ncl+1. The routine should be called with the address165 &a[0][0] as the first argument. */166 {167     long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;168     float **m;169 170     /* allocate pointers to rows */171     m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));172     if (!m) nrerror("allocation failure in convert_matrix()");173     m += NR_END;174     m -= nrl;175 176     /* set pointers to rows */177     m[nrl]=a-ncl;178     for(i=1,j=nrl+1;i
         
          294 #define NR_END 1295 #define FREE_ARG char*296 297 void nrerror(error_text)298 char error_text[];299 /* Numerical Recipes standard error handler */300 {301     void exit();302 303     fprintf(stderr,"Numerical Recipes run-time error...\n");304     fprintf(stderr,"%s\n",error_text);305     fprintf(stderr,"...now exiting to system...\n");306     exit(1);307 }308 309 float *vector(nl,nh)310 long nh,nl;311 /* allocate a float vector with subscript range v[nl..nh] */312 {313     float *v;314 315     v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));316     if (!v) nrerror("allocation failure in vector()");317     return v-nl+NR_END;318 }319 320 int *ivector(nl,nh)321 long nh,nl;322 /* allocate an int vector with subscript range v[nl..nh] */323 {324     int *v;325 326     v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));327     if (!v) nrerror("allocation failure in ivector()");328     return v-nl+NR_END;329 }330 331 unsigned char *cvector(nl,nh)332 long nh,nl;333 /* allocate an unsigned char vector with subscript range v[nl..nh] */334 {335     unsigned char *v;336 337     v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));338     if (!v) nrerror("allocation failure in cvector()");339     return v-nl+NR_END;340 }341 342 unsigned long *lvector(nl,nh)343 long nh,nl;344 /* allocate an unsigned long vector with subscript range v[nl..nh] */345 {346     unsigned long *v;347 348     v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));349     if (!v) nrerror("allocation failure in lvector()");350     return v-nl+NR_END;351 }352 353 double *dvector(nl,nh)354 long nh,nl;355 /* allocate a double vector with subscript range v[nl..nh] */356 {357     double *v;358 359     v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));360     if (!v) nrerror("allocation failure in dvector()");361     return v-nl+NR_END;362 }363 364 float **matrix(nrl,nrh,ncl,nch)365 long nch,ncl,nrh,nrl;366 /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */367 {368     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;369     float **m;370 371     /* allocate pointers to rows */372     m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));373     if (!m) nrerror("allocation failure 1 in matrix()");374     m += NR_END;375     m -= nrl;376 377     /* allocate rows and set pointers to them */378     m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));379     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");380     m[nrl] += NR_END;381     m[nrl] -= ncl;382 383     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;384 385     /* return pointer to array of pointers to rows */386     return m;387 }388 389 double **dmatrix(nrl,nrh,ncl,nch)390 long nch,ncl,nrh,nrl;391 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */392 {393     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;394     double **m;395 396     /* allocate pointers to rows */397     m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));398     if (!m) nrerror("allocation failure 1 in matrix()");399     m += NR_END;400     m -= nrl;401 402     /* allocate rows and set pointers to them */403     m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));404     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");405     m[nrl] += NR_END;406     m[nrl] -= ncl;407 408     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;409 410     /* return pointer to array of pointers to rows */411     return m;412 }413 414 int **imatrix(nrl,nrh,ncl,nch)415 long nch,ncl,nrh,nrl;416 /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */417 {418     long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;419     int **m;420 421     /* allocate pointers to rows */422     m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));423     if (!m) nrerror("allocation failure 1 in matrix()");424     m += NR_END;425     m -= nrl;426 427 428     /* allocate rows and set pointers to them */429     m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));430     if (!m[nrl]) nrerror("allocation failure 2 in matrix()");431     m[nrl] += NR_END;432     m[nrl] -= ncl;433 434     for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;435 436     /* return pointer to array of pointers to rows */437     return m;438 }439 440 float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)441 float **a;442 long newcl,newrl,oldch,oldcl,oldrh,oldrl;443 /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */444 {445     long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;446     float **m;447 448     /* allocate array of pointers to rows */449     m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));450     if (!m) nrerror("allocation failure in submatrix()");451     m += NR_END;452     m -= newrl;453 454     /* set pointers to rows */455     for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;456 457     /* return pointer to array of pointers to rows */458     return m;459 }460 461 float **convert_matrix(a,nrl,nrh,ncl,nch)462 float *a;463 long nch,ncl,nrh,nrl;464 /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix465 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1466 and ncol=nch-ncl+1. The routine should be called with the address467 &a[0][0] as the first argument. */468 {469     long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;470     float **m;471 472     /* allocate pointers to rows */473     m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));474     if (!m)    nrerror("allocation failure in convert_matrix()");475     m += NR_END;476     m -= nrl;477 478     /* set pointers to rows */479     m[nrl]=a-ncl;480     for(i=1,j=nrl+1;i
         
        
       
      

1 /* Driver for routine spline */ 2  3 #include 
      
        4 #include 
       
         5 #define NRANSI 6 #include "nr.h" 7 #include "nrutil.h" 8  9 #define N 2010 #define PI 3.141592611 12 int main(void)13 {14     int i;15     float yp1,ypn,*x,*y,*y2;16 17     x=vector(1,N);18     y=vector(1,N);19     y2=vector(1,N);20     printf("\nsecond-derivatives for sin(x) from 0 to pi\n");21     /* Generate array for interpolation */22     for (i=1;i<=20;i++) {23         x[i]=i*PI/N;24         y[i]=sin(x[i]);25     }26     /* calculate 2nd derivative with spline */27     yp1=cos(x[1]);28     ypn=cos(x[N]);29     spline(x,y,N,yp1,ypn,y2);30     /* test result */31     printf("%23s %16s\n","spline","actual");32     printf("%11s %14s %16s\n","angle","2nd deriv","2nd deriv");33     for (i=1;i<=N;i++)34         printf("%10.2f %16.6f %16.6f\n",x[i],y2[i],-sin(x[i]));35     free_vector(y2,1,N);36     free_vector(y,1,N);37     free_vector(x,1,N);38     return 0;39 }40 #undef NRANSI