package numericalRecipes; import static java.lang.Math.*; /** Routines copied directly from numerical recipes in C. Only neccessary C -> java conversions made * Oh, and it's all in double's now * Also converted for arrays to be 0 based */ public class NumericalRecipes { final static double TINY = 1e-20; public static void dfour1(double data[], int nn, int isign) { int n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; double tempr,tempi; n=nn << 1; j=1; for (i=1;i i) { //SWAP(data[j-1],data[i-1]); double tmp = data[j-1]; data[j-1] = data[i-1]; data[i-1] = tmp; //SWAP(data[j],data[i]); tmp = data[j]; data[j] = data[i]; data[i] = tmp; } m=n >> 1; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax=2; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/mmax); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m>1); if (isign == 1) { c2 = -0.5; dfour1(data,n>>1,1); } else { c2=0.5; theta = -theta; } wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0+wpr; wi=wpi; np3=n+3; for (i=2;i<=(n>>2);i++) { i4=1+(i3=np3-(i2=1+(i1=i+i-1))); h1r=c1*(data[i1-1]+data[i3-1]); h1i=c1*(data[i2-1]-data[i4-1]); h2r = -c2*(data[i2-1]+data[i4-1]); h2i=c2*(data[i1-1]-data[i3-1]); data[i1-1]=h1r+wr*h2r-wi*h2i; data[i2-1]=h1i+wr*h2i+wi*h2r; data[i3-1]=h1r-wr*h2r+wi*h2i; data[i4-1] = -h1i+wr*h2i+wi*h2r; wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } if (isign == 1) { data[0] = (h1r=data[0])+data[1]; data[1] = h1r-data[1]; } else { data[0]=c1*((h1r=data[0])+data[1]); data[1]=c1*(h1r-data[1]); dfour1(data,n>>1,-1); } } public static void dtwofft(double data1[], double data2[], double fft1[], double fft2[],int n) { int nn3,nn2,jj,j; double rep,rem,aip,aim; nn3=1+(nn2=2+n+n); for (j=1,jj=2;j<=n;j++,jj+=2) { fft1[jj-2]=data1[j-1]; fft1[jj-1]=data2[j-1]; } dfour1(fft1,n,1); fft2[0]=fft1[1]; fft1[1]=fft2[1]=0.0; for (j=3;j<=n+1;j+=2) { rep=0.5*(fft1[j-1]+fft1[nn2-j-1]); rem=0.5*(fft1[j-1]-fft1[nn2-j-1]); aip=0.5*(fft1[j]+fft1[nn3-j-1]); aim=0.5*(fft1[j]-fft1[nn3-j-1]); fft1[j-1]=rep; fft1[j]=aim; fft1[nn2-j-1]=rep; fft1[nn3-j-1] = -aim; fft2[j-1]=aip; fft2[j] = -rem; fft2[nn2-j-1]=aip; fft2[nn3-j-1]=rem; } } public static void dconvlv(double data[], int n, double respns[], int m, int isign, double ans[]) { int i,no2; double dum,mag2,fft[]; //fft=vector(1,n<<1); fft = new double[n<<1]; for (i=1;i<=(m-1)/2;i++) respns[n-i]=respns[m-i]; for (i=(m+3)/2;i<=n-(m-1)/2;i++) respns[i-1]=0.0; dtwofft(data,respns,fft,ans,n); no2=n>>1; for (i=2;i<=n+2;i+=2) { if (isign == 1) { ans[i-2]=(fft[i-2]*(dum=ans[i-2])-fft[i-1]*ans[i-1])/no2; ans[i-1]=(fft[i-1]*dum+fft[i-2]*ans[i-1])/no2; } else if (isign == -1) { if ((mag2=SQR(ans[i-2])+SQR(ans[i-1])) == 0.0) throw new ArithmeticException("Deconvolving at response zero in convlv"); ans[i-2]=(fft[i-2]*(dum=ans[i-2])+fft[i-1]*ans[i-1])/mag2/no2; ans[i-1]=(fft[i-1]*dum-fft[i-2]*ans[i-1])/mag2/no2; } else throw new ArithmeticException("No meaning for isign in convlv"); } ans[1]=ans[n]; drealft(ans,n,-1); } public static double SQR(double a){ return a*a; } /** * Multiplies the complex fft arrays f and g together - This is for the convolution of two (at least real, maybe general) functions * @param f * @param g * @return */ public static double[] fftMul(double f[], double g[]){ int n = f.length; double ret[] = new double[n]; for(int i=2;i Math.abs(dum)) { dum=a[j-1][1-1]; i=j; } } indx[k-1]=i; if (dum == 0.0) a[k-1][1-1]=TINY; //Matrix is algorithmically singular, but proceed anyway with TINY pivot (desirable in some applications). //Interchange rows. if (i != k) { d = -d; for (int j=1;j<=mm;j++){ //SWAP(a[k][j],a[i][j]) double tmp = a[k-1][j-1]; a[k-1][j-1] = a[i-1][j-1]; a[i-1][j-1] = tmp; } } // Do the elimination. for (i=k+1;i<=l;i++) { dum=a[i-1][1-1]/a[k-1][1-1]; al[k-1][i-k-1]=dum; for (int j=2;j<=mm;j++) a[i-1][j-1-1]=a[i-1][j-1]-dum*a[k-1][j-1]; a[i-1][mm-1]=0.0; } } return d; } /** * Given the arrays a, al, and indx as returned from bandec, and given a right-hand side vector * b[1..n], solves the band diagonal linear equations A x = b. The solution vector x overwrites * b[1..n]. The other input arrays are not modi ed, and can be left in place for successive calls * with di erent right-hand sides. */ public static void banbks(double a[][], int n, int m1, int m2, double al[][], int indx[], double b[]){ int i,k,l; int mm; double dum; mm=m1+m2+1; l=m1; //Forward substitution, unscrambling the permuted rows for (k=1;k<=n;k++) { // as we go. i=indx[k-1]; if (i != k){ //SWAP(b[k],b[i]) double tmp = b[i-1]; b[i-1] = b[k-1]; b[k-1] = tmp; } if (l < n) l++; for (i=k+1;i<=l;i++) b[i-1] -= al[k-1][i-k-1]*b[k-1]; } l=1; //Backsubstitution. for (i=n;i>=1;i--) { dum=b[i-1]; for (k=2;k<=l;k++) dum -= a[i-1][k-1]*b[k+i-1-1]; b[i-1]=dum/a[i-1][1-1]; if (l < mm) l++; } } /**Matrix multiply b = A x, where A is band diagonal with m1 rows below the diagonal and m2 rows above. The input vector x and output vector b are stored as x[1..n] and b[1..n], respectively. The array a[1..n][1..m1+m2+1] stores A as follows: The diagonal elements are in a[1..n][m1+1]. Subdiagonal elements are in a[j..n][1..m1] (with j > 1 appropriate to the number of elements on each subdiagonal). Superdiagonal elements are in a[1..j][m1+2..m1+m2+1] with j < n appropriate to the number of elements on each superdiagonal.*/ public static double[] banmul(double a[][], int n, int m1, int m2, double x[]) { double b[] = new double[n]; for (int i=1;i<=n;i++) { b[i-1]=0.0; int k=i-m1-1; int j1=m1+m2+1; if( (n-k) < j1) j1 = n-k; int j0 = 1; if( (1-k) > j0) j0 = 1-k; for (int j=j0;j<=j1;j++) b[i-1] += a[i-1][j-1]*x[j+k-1]; } return b; } /** Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above; indx[1..n] is an output vector that records the row permutation effected by the partial pivoting; d is output as 1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix. */ public static double ludcmp(double a[][], int n, int indx[]) { int i,imax=0,j,k; double big,dum,sum,temp; double vv[]; // vv stores the implicit scaling of each row. double d; //for return vv = new double[n]; //vector(1,n); d = 1.0; //No row interchanges yet. for (i=1;i<=n;i++) { //Loop over rows to get the implicit scaling information. big = 0.0; for (j=1;j<=n;j++) if ((temp=Math.abs(a[i-1][j-1])) > big) big=temp; if (big == 0.0) throw new ArithmeticException("Singular matrix in routine ludcmp"); //No nonzero largest element. vv[i-1]=1.0/big; //Save the scaling. } for (j=1;j<=n;j++) { //This is the loop over columns of Crouts method. for (i=1;i= big) { //Is the figure of merit for the pivot better than the best so far? big=dum; imax=i; } } if (j != imax) { //Do we need to interchange rows? for (k=1;k<=n;k++) { //Yes, do so... dum=a[imax-1][k-1]; a[imax-1][k-1]=a[j-1][k-1]; a[j-1][k-1]=dum; } d = -d; //...and change the parity of d. vv[imax-1]=vv[j-1]; //Also interchange the scale factor. } indx[j-1] = imax; if (a[j-1][j-1] == 0.0) a[j-1][j-1] = TINY; //If the pivot element is zero the matrix is singular (at least to the precision of the //algorithm). For some applications on singular matrices, it is desirable to substitute //TINY for zero. if (j != n) { //Now, finally, divide by the pivot element. dum=1.0/(a[j-1][j-1]); for (i=j+1;i<=n;i++) a[i-1][j-1] *= dum; } } //Go back for the next column in the reduction. return d; } /**Solves the set of n linear equations A X = B. Here a[1..n][1..n] is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx[1..n] is input as the permutation vector returned by ludcmp. b[1..n] is input as the right-hand side vector B, and returns with the solution vector X. a, n, and indx are not modified by this routine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that b will begin with many zero elements, so it is efficient for use in matrix inversion.*/ public static void lubksb(double a[][], int n, int indx[], double b[]) { int i,ii=0,ip,j; double sum; for (i=1;i<=n;i++) { //When ii is set to a positive value, it will become the //index of the first nonvanishing element of b. Wenow //do the forward substitution, equation (2.3.6). The //only new wrinkle is to unscramble the permutation //as we go. ip=indx[i-1]; sum=b[ip-1]; b[ip-1]=b[i-1]; if(ii != 0) for (j=ii;j<=i-1;j++) sum -= a[i-1][j-1]*b[j-1]; else if (sum != 0) ii=i; //A nonzero element was encountered, so from now on we b[i-1]=sum; // will have to do the sums in the loop above. } for (i=n;i>=1;i--) { //Now we do the backsubstitution, equation (2.3.7). sum=b[i-1]; for (j=i+1;j<=n;j++) sum -= a[i-1][j-1]*b[j-1]; b[i-1] = sum / a[i-1][i-1]; //Store a component of the solution vector X. } //All done! } /** invert square matrix a, a is replaced with it's LU decomposition and the inverse is returned */ public static double[][] invertQuick(double a[][]){ int N = a.length; double y[][] = new double[N][N]; double col[] = new double[N]; int i,j; int indx[] = new int[N]; ludcmp(a,N,indx); //Decompose the matrix just once. for(j=1;j<=N;j++) { //Find inverse by columns. for(i=1;i<=N;i++) col[i-1]=0.0; col[j-1]=1.0; lubksb(a,N,indx,col); for(i=1;i<=N;i++) y[i-1][j-1]=col[i-1]; } return y; } /** inverse of matrix in, in is returned untouched and it's inverse if returned */ public static double[][] invertDontDestroy(double in[][]){ int N = in.length; double y[][] = new double[N][N]; double col[] = new double[N]; int indx[] = new int[N]; double a[][] = new double[N][N]; for(int i=0;iy[2-1]){ inhi=2; ihi=1; }else{ inhi=1; ihi=2; } for (i=1;i<=mpts;i++) { if (y[i-1] <= y[ilo-1]) ilo=i; if (y[i-1] > y[ihi-1]) { inhi=ihi; ihi=i; } else if (y[i-1] > y[inhi-1] && i != ihi) inhi=i; } rtol=2.0*StrictMath.abs(y[ihi-1]-y[ilo-1])/(StrictMath.abs(y[ihi-1])+StrictMath.abs(y[ilo-1])+TINY); //Compute the fractional range from highest to lowest and return if satisfactory. //If returning, put best point and value in slot 1. if (rtol < ftol || nfunk[0] >= nMax) { double tmp = y[1-1]; y[1-1] = y[ilo-1]; y[ilo-1] = tmp; for (i=1;i<=ndim;i++){ tmp = p[1-1][i-1]; p[1-1][i-1] = p[ilo-1][i-1]; p[ilo-1][i-1] = tmp; } break; } nfunk[0] += 2; //Begin a new iteration. First extrapolate by a factor 1 through the face of the simplex //across from the high point, i.e., reflect the simplex from the high point. ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0); if (ytry <= y[ilo-1]) //Gives a result better than the best point, so try an additional extrapolation by a //factor 2. ytry=amotry(p,y,psum,ndim,funk,ihi,2.0); else if (ytry >= y[inhi-1]) { //The reflected point is worse than the second-highest, so look for an intermediate //lower point, i.e., do a one-dimensional contraction. ysave=y[ihi-1]; ytry=amotry(p,y,psum,ndim,funk,ihi,0.5); // Cant seem to get rid of that high point. Better if (ytry >= ysave) { // contract around the lowest (best) point. for (i=1;i<=mpts;i++) { if (i != ilo) { for (j=1;j<=ndim;j++) p[i-1][j-1]=psum[j-1]=0.5*(p[i-1][j-1]+p[ilo-1][j-1]); //y[i]=(*funk)(psum); y[i-1] = funk.funk(psum); } } // Keep track of function evaluations. nfunk[0] += ndim; // Recompute psum. //GET_PSUM for (j=1;j<=ndim;j++) { for (sum=0.0,i=1;i<=mpts;i++) sum += p[i-1][j-1]; psum[j-1]=sum;} } // Correct the evaluation count. } else --nfunk[0]; //Go back for the test of doneness and the next iteration. } // free_vector(psum,1,ndim); } /** Extrapolates by a factor fac through the face of the simplex across from the high point, tries * it, and replaces the high point if the new point is better. */ private static double amotry(double p[][], double y[], double psum[], int ndim,NRFunctionProvider funk, int ihi, double fac) { int j; double fac1,fac2,ytry,ptry[]; //ptry=vector(1,ndim); ptry = new double[ndim]; fac1=(1.0-fac)/ndim; fac2=fac1-fac; for (j=1;j<=ndim;j++) ptry[j-1]=psum[j-1]*fac1-p[ihi-1][j-1]*fac2; // Evaluate the function at the trial point. //ytry=(*funk)(ptry); ytry = funk.funk(ptry); // If its better than the highest, then replace the highest. if (ytry < y[ihi-1]) { y[ihi-1]=ytry; for (j=1;j<=ndim;j++) { psum[j-1] += ptry[j-1]-p[ihi-1][j-1]; p[ihi-1][j-1]=ptry[j-1]; } } //free_vector(ptry,1,ndim); return ytry; } /**Minimization of a function func of n variables. Input consists of an initial starting point p[1..n]; an initial matrix xi[1..n][1..n], whose columns contain the initial set of di- rections (usually the n unit vectors); and ftol, the fractional tolerance in the function value such that failure to decrease by more than this amount on one iteration signals doneness. On output, p is set to the best point found, xi is the then-current direction set, fret is the returned function value at p, and iter is the number of iterations taken. The routine linmin is used.*/ /* Never finished this one //void powell(float p[], float **xi, int n, float ftol, int *iter, float *fret,float (*func)(float [])) public static int PowellDirectionSet(double p[], double xi[][], int n, double ftol, int iter[], double fret[],NRFunctionProvider func){ //void linmin(float p[], float xi[], int n, float *fret,float (*func)(float [])); int i,ibig,j; int iter; double del,fp,fptt,t,pt[],ptt[],xit[]; pt=new double[n]; ptt=new double[n]; xit=new double[n]; *fret=(*func)(p); // ?? for (j=1;j<=n;j++) pt[j]=p[j]; for (iter=1;;++(iter)) { fp=(*fret); ibig=0; // Will be the biggest function decrease. del=0.0; // In each iteration, loop over all directions in the set. for (i=1;i<=n;i++) { // Copy the direction, for (j=1;j<=n;j++) xit[j]=xi[j][i]; fptt=(*fret); // minimize along it, linmin(p,xit,n,fret,func); // and record it if it is the largest decrease if (fptt-(*fret) > del) { // so far. del=fptt-(*fret); ibig=i; } } if (2.0*(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))+TINY) { // Termination criterion. free_vector(xit,1,n); free_vector(ptt,1,n); free_vector(pt,1,n); return; } if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); // Construct the extrapolated point and the for (j=1;j<=n;j++) { // average direction moved. Save the ptt[j]=2.0*p[j]-pt[j]; // old starting point. xit[j]=p[j]-pt[j]; pt[j]=p[j]; } // Function value at extrapolated point. fptt=(*func)(ptt); if (fptt < fp) { t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); if (t < 0.0) { // Move to the minimum of the new direc- linmin(p,xit,n,fret,func); // tion, and save the new direction. for (j=1;j<=n;j++) { xi[j][ibig]=xi[j][n]; xi[j][n]=xit[j]; } } } // Back for another iteration. } } }*/ /** * Returns the modified Bessel function of the first kind, modified Bessel function of the second kind, and * its derivatives, all with fractional order. * @param nu The fractional order. * @param x The values at which the Bessel functions should be evaluated. * @return * ret[0] is the modified Bessel function of the first kind (Inu), * ret[1] the modified Bessel function of the second kind (Knu), * ret[2] the derivative of the modified Bessel function of the first kind (dInu/dx) * ret[3] the derivative of the modified Bessel function of the second kind (dKnu/dx) */ public static double[][] bessik(double nu, double[] x) { final double EPS = 1.0e-16; final double FPMIN = 1.0e-30; final int MAXIT = 10000; final double XMIN = 2.0; final double PI = 3.141592653589793; int i,l,nl; double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff, h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl, ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2; double[] ri = new double[x.length], rip = new double[x.length]; double[] rk = new double[x.length], rkp = new double[x.length]; double[] gam1 = new double[1], gam2 = new double[1], gampl = new double[1], gammi = new double[1]; for (int n=0;n=1;l--) { ritemp=fact*ril+ripl; fact -= xi; ripl=fact*ritemp+ril; ril=ritemp; } f=ripl/ril; if (x[n] < XMIN) { x2=0.5*x[n]; pimu=PI*xmu; fact = (Math.abs(pimu) < EPS ? 1.0 : pimu/Math.sin(pimu)); d = -Math.log(x2); e=xmu*d; fact2 = (Math.abs(e) < EPS ? 1.0 : Math.sinh(e)/e); beschb(xmu,gam1,gam2,gampl,gammi); ff=fact*(gam1[0]*Math.cosh(e)+gam2[0]*fact2*d); sum=ff; e=Math.exp(e); p=0.5*e/gampl[0]; q=0.5/(e*gammi[0]); c=1.0; d=x2*x2; sum1=p; for (i=1;i<=MAXIT;i++) { ff=(i*ff+p+q)/(i*i-xmu2); c *= (d/i); p /= (i-xmu); q /= (i+xmu); del=c*ff; sum += del; del1=c*(p-i*ff); sum1 += del1; if (Math.abs(del) < Math.abs(sum)*EPS) break; } if (i > MAXIT) throw new ArithmeticException("bessk series failed to converge"); rkmu=sum; rk1=sum1*xi2; } else { b=2.0*(1.0+x[n]); d=1.0/b; h=delh=d; q1=0.0; q2=1.0; a1=0.25-xmu2; q=c=a1; a = -a1; s=1.0+q*delh; for (i=2;i<=MAXIT;i++) { a -= 2*(i-1); c = -a*c/i; qnew=(q1-b*q2)/a; q1=q2; q2=qnew; q += c*qnew; b += 2.0; d=1.0/(b+a*d); delh=(b*d-1.0)*delh; h += delh; dels=q*delh; s += dels; if (Math.abs(dels/s) < EPS) break; } if (i > MAXIT) throw new ArithmeticException("bessik: failure to converge in cf2"); h=a1*h; rkmu=Math.sqrt(PI/(2.0*x[n]))*Math.exp(-x[n])/s; rk1=rkmu*(xmu+x[n]+0.5-h)*xi; } rkmup=xmu*xi*rkmu-rk1; rimu=xi/(f*rkmu-rkmup); ri[n]=(rimu*ril1)/ril; rip[n]=(rimu*rip1)/ril; for (i=1;i<=nl;i++) { rktemp=(xmu+i)*xi2*rk1+rkmu; rkmu=rk1; rk1=rktemp; } rk[n]=rkmu; rkp[n]=nu*xi*rkmu-rk1; } return new double[][] { ri, rk, rip, rkp }; } static final double beschb_c1[] = { -1.142022680371168e0,6.5165112670737e-3, 3.087090173086e-4,-3.4706269649e-6,6.9437664e-9, 3.67795e-11,-1.356e-13}; static final double beschb_c2[] = { 1.843740587300905e0,-7.68528408447867e-2, 1.2719271366546e-3,-4.9717367042e-6,-3.31261198e-8, 2.423096e-10,-1.702e-13,-1.49e-15}; /** * Evaluates Gamma1 and Gamma2 by Chebyshev expansion for abs(x) <= 0.5 Also returns * 1/Gamma(1+x) and 1/Gamma(1-x). * @param x */ public static void beschb(double x, double[] gam1, double[] gam2, double[] gampl, double[] gammi) { int NUSE1 = 7; int NUSE2 = 8; double xx; xx=8.0*x*x-1.0; gam1[0]=chebev(-1.0,1.0,beschb_c1,NUSE1,xx); gam2[0]=chebev(-1.0,1.0,beschb_c2,NUSE2,xx); gampl[0]= gam2[0]-x*(gam1[0]); gammi[0]= gam2[0]+x*(gam1[0]); } /** * Chebyshev evaluation: All arguments are input. c[0..m-1] is an array of Chebyshev coefficients, * the first m elements of c output from chebft (which must have been called with the * same a and b). The Chebyshev polynomial Pm-1 k=0 ckTk(y) − c0=2 is evaluated at a point * y = [x−(b+a)=2]=[(b− a)=2], and the result is returned as the function value. */ public static double chebev(double a, double b, double c[], int m, double x) { double d=0.0,dd=0.0,sv,y,y2; int j; if ((x-a)*(x-b) > 0.0) throw new ArithmeticException("x not in range in routine chebev"); y2=2.0*(y=(2.0*x-a-b)/(b-a)); for (j=m-1;j>=1;j--) { sv=d; d=y2*d-dd+c[j]; dd=sv; } return y*d-dd+0.5*c[0]; } }