#include #define NRANSI #include "nrutil.h" #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} void gaussj(double **a, int n, double **b, int m) { int *indxc,*indxr,*ipiv; int i,icol,irow,j,k,l,ll; double big,dum,pivinv,temp; indxc=ivector(1,n); indxr=ivector(1,n); ipiv=ivector(1,n); for (j=1;j<=n;j++) ipiv[j]=0; for (i=1;i<=n;i++) { big=0.0; for (j=1;j<=n;j++) if (ipiv[j] != 1) for (k=1;k<=n;k++) { if (ipiv[k] == 0) { if (fabs(a[j][k]) >= big) { big=fabs(a[j][k]); irow=j; icol=k; } } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1"); } ++(ipiv[icol]); if (irow != icol) { for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l]) for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l]) } indxr[i]=irow; indxc[i]=icol; if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2"); pivinv=1.0/a[icol][icol]; a[icol][icol]=1.0; for (l=1;l<=n;l++) a[icol][l] *= pivinv; for (l=1;l<=m;l++) b[icol][l] *= pivinv; for (ll=1;ll<=n;ll++) if (ll != icol) { dum=a[ll][icol]; a[ll][icol]=0.0; for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum; for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum; } } for (l=n;l>=1;l--) { if (indxr[l] != indxc[l]) for (k=1;k<=n;k++) SWAP(a[k][indxr[l]],a[k][indxc[l]]); } free_ivector(ipiv,1,n); free_ivector(indxr,1,n); free_ivector(indxc,1,n); } #undef SWAP #define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;} void covsrt(double **covar, int ma, int ia[], int mfit) { int i,j,k; double swap; for (i=mfit+1;i<=ma;i++) for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0; k=mfit; for (j=ma;j>=1;j--) { if (ia[j] == 1) { for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j]) for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i]) k--; } } } #define CON 1.4 #define CON2 (CON*CON) #define BIG 1.0e30 #define NTAB 10 #define SAFE 2.0 float dfridr(float (*func)(float), float x, float h, float *err) { int i,j; float errt,fac,hh,**a,ans; if (h == 0.0) nrerror("h must be nonzero in dfridr."); a=matrix(1,NTAB,1,NTAB); hh=h; a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); *err=BIG; for (i=2;i<=NTAB;i++) { hh /= CON; a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh); fac=CON2; for (j=2;j<=i;j++) { a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0); fac=CON2*fac; errt=FMAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1])); if (errt <= *err) { *err=errt; ans=a[j][i]; } } if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) break; } free_matrix(a,1,NTAB,1,NTAB); return ans; } #undef CON #undef CON2 #undef BIG #undef NTAB #undef SAFE float gammln(float xx) { double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j<=5;j++) ser += cof[j]/++y; return -tmp+log(2.5066282746310005*ser/x); } #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ a[k][l]=h+s*(g-h*tau); void jacobi(double **a, int n, double d[], double **v, int *nrot) { int j,iq,ip,i; double tresh,theta,tau,t,sm,s,h,g,c,*b,*z; b=dvector(1,n); z=dvector(1,n); for (ip=1;ip<=n;ip++) { for (iq=1;iq<=n;iq++) v[ip][iq]=0.0; v[ip][ip]=1.0; } for (ip=1;ip<=n;ip++) { b[ip]=d[ip]=a[ip][ip]; z[ip]=0.0; } *nrot=0; for (i=1;i<=50;i++) { sm=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]); } if (sm == 0.0) { free_dvector(z,1,n); free_dvector(b,1,n); return; } if (i < 4) tresh=0.2*sm/(n*n); else tresh=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) { g=100.0*fabs(a[ip][iq]); if (i > 4 && (fabs(d[ip])+g) == fabs(d[ip]) && (fabs(d[iq])+g) == fabs(d[iq])) a[ip][iq]=0.0; else if (fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; if ((fabs(h)+g) == fabs(h)) t=(a[ip][iq])/h; else { theta=0.5*h/(a[ip][iq]); t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } c=1.0/sqrt(1+t*t); s=t*c; tau=s/(1.0+c); h=t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq]=0.0; for (j=1;j<=ip-1;j++) { ROTATE(a,j,ip,j,iq) } for (j=ip+1;j<=iq-1;j++) { ROTATE(a,ip,j,j,iq) } for (j=iq+1;j<=n;j++) { ROTATE(a,ip,j,iq,j) } for (j=1;j<=n;j++) { ROTATE(v,j,ip,j,iq) } ++(*nrot); } } } for (ip=1;ip<=n;ip++) { b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.0; } } nrerror("Too many iterations in routine jacobi"); } #undef ROTATE #define ITMAX 100 #define EPS 3.0e-7 #define FPMIN 1.0e-30 void gcf(float *gammcf, float a, float x, float *gln) { float gammln(float xx); void nrerror(char error_text[]); int i; float an,b,c,d,del,h; *gln=gammln(a); b=x+1.0-a; c=1.0/FPMIN; d=1.0/b; h=d; for (i=1;i<=ITMAX;i++) { an = -i*(i-a); b += 2.0; d=an*d+b; if (fabs(d) < FPMIN) d=FPMIN; c=b+an/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; del=d*c; h *= del; if (fabs(del-1.0) < EPS) break; } if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf"); *gammcf=exp(-x+a*log(x)-(*gln))*h; } #undef ITMAX #undef EPS #undef FPMIN #define ITMAX 100 #define EPS 3.0e-7 void gser(float *gamser, float a, float x, float *gln) { float gammln(float xx); void nrerror(char error_text[]); int n; float sum,del,ap; *gln=gammln(a); if (x <= 0.0) { if (x < 0.0) nrerror("x less than 0 in routine gser"); *gamser=0.0; return; } else { ap=a; del=sum=1.0/a; for (n=1;n<=ITMAX;n++) { ++ap; del *= x/ap; sum += del; if (fabs(del) < fabs(sum)*EPS) { *gamser=sum*exp(-x+a*log(x)-(*gln)); return; } } nrerror("a too large, ITMAX too small in routine gser"); return; } } #undef ITMAX #undef EPS float gammp(float a, float x) { void gcf(float *gammcf, float a, float x, float *gln); void gser(float *gamser, float a, float x, float *gln); void nrerror(char error_text[]); float gamser,gammcf,gln; if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp"); if (x < (a+1.0)) { gser(&gamser,a,x,&gln); return gamser; } else { gcf(&gammcf,a,x,&gln); return 1.0-gammcf; } } void polint(float xa[], float ya[], int n, float x, float *y, float *dy) { int i,m,ns=1; float den,dif,dift,ho,hp,w; float *c,*d; dif=fabs(x-xa[1]); c=vector(1,n); d=vector(1,n); for (i=1;i<=n;i++) { if ( (dift=fabs(x-xa[i])) < dif) { ns=i; dif=dift; } c[i]=ya[i]; d[i]=ya[i]; } *y=ya[ns--]; for (m=1;m