猪冰龙

导航

c语言spline 

 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
spline.c

  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 <stdio.h>
 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_ */
nr.h

 

  1 #if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
  2 
  3 #include <stdio.h>
  4 #include <stddef.h>
  5 #include <stdlib.h>
  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 matrix
163 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
164 and ncol=nch-ncl+1. The routine should be called with the address
165 &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<nrow;i++,j++) m[j]=m[j-1]+ncol;
179     /* return pointer to array of pointers to rows */
180     return m;
181 }
182 
183 float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
184 /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
185 {
186     long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
187     float ***t;
188 
189     /* allocate pointers to pointers to rows */
190     t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
191     if (!t) nrerror("allocation failure 1 in f3tensor()");
192     t += NR_END;
193     t -= nrl;
194 
195     /* allocate pointers to rows and set pointers to them */
196     t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
197     if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
198     t[nrl] += NR_END;
199     t[nrl] -= ncl;
200 
201     /* allocate rows and set pointers to them */
202     t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
203     if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
204     t[nrl][ncl] += NR_END;
205     t[nrl][ncl] -= ndl;
206 
207     for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
208     for(i=nrl+1;i<=nrh;i++) {
209         t[i]=t[i-1]+ncol;
210         t[i][ncl]=t[i-1][ncl]+ncol*ndep;
211         for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
212     }
213 
214     /* return pointer to array of pointers to rows */
215     return t;
216 }
217 
218 void free_vector(float *v, long nl, long nh)
219 /* free a float vector allocated with vector() */
220 {
221     free((FREE_ARG) (v+nl-NR_END));
222 }
223 
224 void free_ivector(int *v, long nl, long nh)
225 /* free an int vector allocated with ivector() */
226 {
227     free((FREE_ARG) (v+nl-NR_END));
228 }
229 
230 void free_cvector(unsigned char *v, long nl, long nh)
231 /* free an unsigned char vector allocated with cvector() */
232 {
233     free((FREE_ARG) (v+nl-NR_END));
234 }
235 
236 void free_lvector(unsigned long *v, long nl, long nh)
237 /* free an unsigned long vector allocated with lvector() */
238 {
239     free((FREE_ARG) (v+nl-NR_END));
240 }
241 
242 void free_dvector(double *v, long nl, long nh)
243 /* free a double vector allocated with dvector() */
244 {
245     free((FREE_ARG) (v+nl-NR_END));
246 }
247 
248 void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
249 /* free a float matrix allocated by matrix() */
250 {
251     free((FREE_ARG) (m[nrl]+ncl-NR_END));
252     free((FREE_ARG) (m+nrl-NR_END));
253 }
254 
255 void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
256 /* free a double matrix allocated by dmatrix() */
257 {
258     free((FREE_ARG) (m[nrl]+ncl-NR_END));
259     free((FREE_ARG) (m+nrl-NR_END));
260 }
261 
262 void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
263 /* free an int matrix allocated by imatrix() */
264 {
265     free((FREE_ARG) (m[nrl]+ncl-NR_END));
266     free((FREE_ARG) (m+nrl-NR_END));
267 }
268 
269 void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
270 /* free a submatrix allocated by submatrix() */
271 {
272     free((FREE_ARG) (b+nrl-NR_END));
273 }
274 
275 void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
276 /* free a matrix allocated by convert_matrix() */
277 {
278     free((FREE_ARG) (b+nrl-NR_END));
279 }
280 
281 void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
282     long ndl, long ndh)
283 /* free a float f3tensor allocated by f3tensor() */
284 {
285     free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
286     free((FREE_ARG) (t[nrl]+ncl-NR_END));
287     free((FREE_ARG) (t+nrl-NR_END));
288 }
289 
290 #else /* ANSI */
291 /* traditional - K&R */
292 
293 #include <stdio.h>
294 #define NR_END 1
295 #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 matrix
465 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
466 and ncol=nch-ncl+1. The routine should be called with the address
467 &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<nrow;i++,j++) m[j]=m[j-1]+ncol;
481     /* return pointer to array of pointers to rows */
482     return m;
483 }
484 
485 float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
486 long nch,ncl,ndh,ndl,nrh,nrl;
487 /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
488 {
489     long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
490     float ***t;
491 
492     /* allocate pointers to pointers to rows */
493     t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
494     if (!t) nrerror("allocation failure 1 in f3tensor()");
495     t += NR_END;
496     t -= nrl;
497 
498     /* allocate pointers to rows and set pointers to them */
499     t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
500     if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
501     t[nrl] += NR_END;
502     t[nrl] -= ncl;
503 
504     /* allocate rows and set pointers to them */
505     t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
506     if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
507     t[nrl][ncl] += NR_END;
508     t[nrl][ncl] -= ndl;
509 
510     for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
511     for(i=nrl+1;i<=nrh;i++) {
512         t[i]=t[i-1]+ncol;
513         t[i][ncl]=t[i-1][ncl]+ncol*ndep;
514         for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
515     }
516 
517     /* return pointer to array of pointers to rows */
518     return t;
519 }
520 
521 void free_vector(v,nl,nh)
522 float *v;
523 long nh,nl;
524 /* free a float vector allocated with vector() */
525 {
526     free((FREE_ARG) (v+nl-NR_END));
527 }
528 
529 void free_ivector(v,nl,nh)
530 int *v;
531 long nh,nl;
532 /* free an int vector allocated with ivector() */
533 {
534     free((FREE_ARG) (v+nl-NR_END));
535 }
536 
537 void free_cvector(v,nl,nh)
538 long nh,nl;
539 unsigned char *v;
540 /* free an unsigned char vector allocated with cvector() */
541 {
542     free((FREE_ARG) (v+nl-NR_END));
543 }
544 
545 void free_lvector(v,nl,nh)
546 long nh,nl;
547 unsigned long *v;
548 /* free an unsigned long vector allocated with lvector() */
549 {
550     free((FREE_ARG) (v+nl-NR_END));
551 }
552 
553 void free_dvector(v,nl,nh)
554 double *v;
555 long nh,nl;
556 /* free a double vector allocated with dvector() */
557 {
558     free((FREE_ARG) (v+nl-NR_END));
559 }
560 
561 void free_matrix(m,nrl,nrh,ncl,nch)
562 float **m;
563 long nch,ncl,nrh,nrl;
564 /* free a float matrix allocated by matrix() */
565 {
566     free((FREE_ARG) (m[nrl]+ncl-NR_END));
567     free((FREE_ARG) (m+nrl-NR_END));
568 }
569 
570 void free_dmatrix(m,nrl,nrh,ncl,nch)
571 double **m;
572 long nch,ncl,nrh,nrl;
573 /* free a double matrix allocated by dmatrix() */
574 {
575     free((FREE_ARG) (m[nrl]+ncl-NR_END));
576     free((FREE_ARG) (m+nrl-NR_END));
577 }
578 
579 void free_imatrix(m,nrl,nrh,ncl,nch)
580 int **m;
581 long nch,ncl,nrh,nrl;
582 /* free an int matrix allocated by imatrix() */
583 {
584     free((FREE_ARG) (m[nrl]+ncl-NR_END));
585     free((FREE_ARG) (m+nrl-NR_END));
586 }
587 
588 void free_submatrix(b,nrl,nrh,ncl,nch)
589 float **b;
590 long nch,ncl,nrh,nrl;
591 /* free a submatrix allocated by submatrix() */
592 {
593     free((FREE_ARG) (b+nrl-NR_END));
594 }
595 
596 void free_convert_matrix(b,nrl,nrh,ncl,nch)
597 float **b;
598 long nch,ncl,nrh,nrl;
599 /* free a matrix allocated by convert_matrix() */
600 {
601     free((FREE_ARG) (b+nrl-NR_END));
602 }
603 
604 void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
605 float ***t;
606 long nch,ncl,ndh,ndl,nrh,nrl;
607 /* free a float f3tensor allocated by f3tensor() */
608 {
609     free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
610     free((FREE_ARG) (t[nrl]+ncl-NR_END));
611     free((FREE_ARG) (t+nrl-NR_END));
612 }
613 
614 #endif /* ANSI */
nrutil.c
 1 /* Driver for routine spline */
 2 
 3 #include <stdio.h>
 4 #include <math.h>
 5 #define NRANSI
 6 #include "nr.h"
 7 #include "nrutil.h"
 8 
 9 #define N 20
10 #define PI 3.1415926
11 
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
xspline.c

 

posted on 2018-03-12 22:27  猪冰龙  阅读(1302)  评论(0)  编辑  收藏  举报