#include <stdio.h>
/* #include <stddef.h> */
/* #include <stdlib.h> */
#include <math.h>
#include "nr.h"

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl+NR_END;
}

void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;

	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl+NR_END;
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate an integer matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	int **m;

	/* allocate pointers to rows */
	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an integer matrix allocated by imatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

char **cmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a char matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	char **m;

	/* allocate pointers to rows */
	m=(char **) malloc((size_t)((nrow+NR_END)*sizeof(char*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(char)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch)
/* free a char matrix allocated by dmatrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}

double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
	double ***t;

	/* allocate pointers to pointers to rows */
	t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
	if (!t) nrerror("allocation failure 1 in f3tensor()");
	t += NR_END;
	t -= nrl;

	/* allocate pointers to rows and set pointers to them */
	t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
	t[nrl] += NR_END;
	t[nrl] -= ncl;

	/* allocate rows and set pointers to them */
	t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
	t[nrl][ncl] += NR_END;
	t[nrl][ncl] -= ndl;

	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
	for(i=nrl+1;i<=nrh;i++) {
		t[i]=t[i-1]+ncol;
		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
	}

	/* return pointer to array of pointers to rows */
	return t;
}

void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
	long ndl, long ndh)
/* free a double f3tensor allocated by f3tensor() */
{
	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
	free((FREE_ARG) (t[nrl]+ncl-NR_END));
	free((FREE_ARG) (t+nrl-NR_END));
}

#define NRANSI

void convlv(double data[], unsigned long n, double respns[], unsigned long m,
	int isign, double ans[])
{
	unsigned long i,no2;
	double dum,mag2,*fft;

	fft=dvector(1,n<<1);
	for (i=1;i<=(m-1)/2;i++)
		respns[n+1-i]=respns[m+1-i];
	for (i=(m+3)/2;i<=n-(m-1)/2;i++)
		respns[i]=0.0;
	twofft(data,respns,fft,ans,n);
	no2=n>>1;
	for (i=2;i<=n+2;i+=2) {
		if (isign == 1) {
			ans[i-1]=(fft[i-1]*(dum=ans[i-1])-fft[i]*ans[i])/no2;
			ans[i]=(fft[i]*dum+fft[i-1]*ans[i])/no2;
		} else if (isign == -1) {
			if ((mag2=SQR(ans[i-1])+SQR(ans[i])) == 0.0)
				nrerror("Deconvolving at response zero in convlv");
			ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/mag2/no2;
			ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/mag2/no2;
		} else nrerror("No meaning for isign in convlv");
	}
	ans[2]=ans[n+1];
	realft(ans,n,-1);
	free_dvector(fft,1,n<<1);
}
#undef NRANSI

void correl(double data1[], double data2[], unsigned long n, double ans[])
{
        unsigned long no2,i;
        double dum,*fft;

        fft=dvector(1,n<<1);
        twofft(data1,data2,fft,ans,n);
        no2=n>>1;
        for (i=2;i<=n+2;i+=2) {
            ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/no2;
            ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/no2;
	}
        ans[2]=ans[n+1];
        realft(ans,n,-1);
        free_dvector(fft,1,n<<1);
}

void twofft(double data1[], double data2[], double fft1[], double fft2[],
	unsigned long n)
{
	unsigned long 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-1]=data1[j];
		fft1[jj]=data2[j];
	}
	four1(fft1,n,1);
	fft2[1]=fft1[2];
	fft1[2]=fft2[2]=0.0;
	for (j=3;j<=n+1;j+=2) {
		rep=0.5*(fft1[j]+fft1[nn2-j]);
		rem=0.5*(fft1[j]-fft1[nn2-j]);
		aip=0.5*(fft1[j+1]+fft1[nn3-j]);
		aim=0.5*(fft1[j+1]-fft1[nn3-j]);
		fft1[j]=rep;
		fft1[j+1]=aim;
		fft1[nn2-j]=rep;
		fft1[nn3-j] = -aim;
		fft2[j]=aip;
		fft2[j+1] = -rem;
		fft2[nn2-j]=aip;
		fft2[nn3-j]=rem;
	}
}

void realft(double data[], unsigned long n, int isign)
{
	void four1(double data[], unsigned long nn, int isign);
	unsigned long i,i1,i2,i3,i4,np3;
	double c1=0.5,c2,h1r,h1i,h2r,h2i;
	double wr,wi,wpr,wpi,wtemp,theta;

	theta=3.141592653589793/(double) (n>>1);
	if (isign == 1) {
		c2 = -0.5;
		four1(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]+data[i3]);
		h1i=c1*(data[i2]-data[i4]);
		h2r = -c2*(data[i2]+data[i4]);
		h2i=c2*(data[i1]-data[i3]);
		data[i1]=h1r+wr*h2r-wi*h2i;
		data[i2]=h1i+wr*h2i+wi*h2r;
		data[i3]=h1r-wr*h2r+wi*h2i;
		data[i4] = -h1i+wr*h2i+wi*h2r;
		wr=(wtemp=wr)*wpr-wi*wpi+wr;
		wi=wi*wpr+wtemp*wpi+wi;
	}
	if (isign == 1) {
		data[1] = (h1r=data[1])+data[2];
		data[2] = h1r-data[2];
	} else {
		data[1]=c1*((h1r=data[1])+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n>>1,-1);
	}
}

void four1(double data[], unsigned long nn, int isign)
{
	unsigned long n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	double tempr,tempi,temp;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		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<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

void mrqmin(double x[], double y[], double sig[], int ndata, double a[], 
int ia[], int ma, double **covar, double **alpha, double *chisq, 
void (*funcs)( double, double[], double *, double[], int), double *alamda) 
/* Levenberg-Marquardt method, attempting to reduce the value chi^2 of a fit */
/* between a set of data points x[1..ndata] , y[1..ndat a] with individual   */
/* standard deviations sig[1..ndata], and a nonlinear function dependent on  */
/* ma coeffcients a[1..ma]. The input array ia[1..ma] indicates by nonzero   */
/* entries those components of a that should be fitted for, and by zero      */
/* entries those components that should be held fixed at their input values. */
/*The program returns current best-fit values for the parameters a[1..ma] ,  */
/* and chi^2 = chisq. The arrays covar[1..ma][1. .ma], alpha[1..ma][1..ma]   */
/* are used as w rking space during most iterations. Supply a routine        */
/* funcs(x,a,yfit,dyda,ma) that evaluates the fitting function yfit, and its */
/* derivatives dyda[1..ma] with respect to the fitting parameters a at x.    */
/* On the first call p rovide an initial guess for the parameters a , and    */
/* set alamda<0 fo r initialization (which then sets alamda=.001). If a step */
/* succeeds chisq becomes smaller and alamda decreases by a facto r of 10.   */
/* If a step fails alamda grows by a factor of 10. You must call this        */
/* routine repeatedly until convergence is achieved. Then, make one final    */
/* call with alamda=0, so that covar[1..ma][1..ma] returns the covariance    */
/* matrix, and alpha the curvature matrix. (P a rameters held fixed will     */
/* return zero covariances.)                                                 */
{ 
  int j,k,l; 
  static int mfit; 
  static double ochisq,*atry ,*beta ,*da,**oneda; 

  if (*alamda < 0.0) {                                   /* Initialization.  */
    atry=dvector(1,ma); 
    beta=dvector(1,ma); 
    da=dvector(1,ma); 
      for (mfit=0,j=1;j<=ma;j++) 
        if (ia[j]) mfit++; 
      oneda=dmatrix(1,mfit,1,1); 
      *alamda=0.001; 
      mrqcof(x,y ,sig,ndata,a,ia,ma,alpha,beta,chisq,funcs); 
      ochisq=(*chisq); 
      for (j=1;j<=ma; j++) atry[j]=a [j] ; 
    } 
    /* Alter linearized fitting matrix, by augmenting diagonal elements.    */ 
    for (j=1;j<=mfit; j++) { 
      for (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k]; 
      covar[j][j]=alpha[j][j]*(1.0+(*alamda)); 
      oneda[j][1]=beta[j]; 
    } 
    gaussj(covar,mfit,oneda,1);                         /*Matrix solution.   */
    for (j=1;j<=mfit; j++) da[j]=oneda[j][1] ; 
    if (*alamda == 0.0) {    /* Once converged, evaluate covariance matrix.  */
      covsrt(covar,ma,ia,mfit); 
      free_dmatrix(oneda,1,mfit,1,1); 
      free_dvector(da,1,ma); 
      free_dvector(beta,1,ma); 
      free_dvector(atry,1,ma); 
    return; 
  } 
  for (j=0,l=1;l<=ma;l++)                         /* Did the trial succeed?  */
    if (ia[l]) atry[l]=a[l]+da[++j]; 
  mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq,funcs); 
  if (*chisq < ochisq) {               /* Success, accept the new solution.  */
    *alamda *= 0.1; 
    ochisq=(*chisq); 
    for (j=1;j<=mfit;j++) { 
      for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k] ; 
      beta[j]=da[j]; 
    } 
    for (l=1;l<=ma;l++) a[l]=atry[l]; 
  } 
  else {                            /* Failure, increase alamda and return.  */
    *alamda *= 10.0; 
    *chisq=ochisq; 
    } 
  } 

void mrqcof(double x[], double y[], double sig[], int ndata, double a[], 
int ia[], int ma, double **alpha, double beta[], double *chisq, 
void (*funcs)( double, double[], double *, double[], int)) 
/* Used by mrqmin to evaluate the linearized fitting matrix alpha, and       */
/* vector beta as in (15.5.8), and calculate chi^2.                          */
{ 
  int i,j,k,l,m,mfit=0; 
  double ymod,wt,sig2i,dy,*dyda; 

  dyda=dvector(1,ma); 
  for (j=1;j<=ma;j++) 
    if (ia[j]) mfit++; 
  for (j=1;j<=mfit; j++) {           /* Initialize (symmetric) alpha, beta.  */
    for (k=1;k<=j;k++) alpha[j][k]=0.0; 
    beta[j]=0.0; 
  } 
  *chisq=0.0; 
  for (i=1;i<=ndata ;i++) {               /* Summation loop over all data.  */
    (*funcs)(x[i],a,&ymod,dyda,ma); 
    sig2i=1.0/(sig[i]*sig[i]); 
    dy=y[i]-ymod; 
    for (j=0,l=1;l<=ma;l++) { 
      if (ia[l]) { 
        wt=dyda[l]*sig2i; 
        for (j++,k=0,m=1;m<=l;m++) 
        if (ia[m]) alpha[j][++k]+= wt*dyda[m]; 
        beta[j]+= dy*wt; 
      } 
    } 
    *chisq += dy*dy*sig2i;                               /* And find chi^2.  */
  } 
  for (j=2;j<=mfit; j++)                     /* Fill in the symmetric side.  */
  for (k=1;k<j;k++) alpha[k][j] =alpha[j][k]; 
  free_dvector(dyda,1,ma); 
}

void covsrt(double **covar, int ma, int ia[], int mfit) 
/* Expand in storage the covariance matrix covar , so as to take into    */
/* account parameters that are being held fixed. (For the latter, return */
/* zero covariances.)                                                    */
{ 
  int i,j,k; 
  double temp; 

  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]) { 
      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--; 
    } 
  } 
} 

void gaussj(double **a, int n, double **b, int m) 

/* Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) */
/* above. a[1..n][1..n ] is the input matrix. b[1..n][1..m ] is input     */
/* containing the m right-hand side vectors. On output, a is replaced by  */
/* its matrix inverse, and b is replaced by the corresponding set of      */
/* solution vectors.                                                      */

{ 
  int *indxc,*indxr,*ipiv; 
  int i,icol,irow,j,k,l,ll; 
  double big,dum,pivinv,temp; 
  /* The integer arrays ipiv, indxr, and indxc are used for bookkeeping   */
  /* on the pivoting.                                                     */
  indxc=ivector(1,n); 
  indxr=ivector(1,n); 
  ipiv=ivector(1,n); 
  for (j=1;j<=n; j++) ipiv[j]=0; 
  /* This is the main lo op over the columns to be reduced.               */
  for (i=1;i<=n; i++) { 
    big=0.0; 
    /* This is the outer loop of the search for a pivot element.          */
    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]); 
    /* We now have the pivot element, so we interchange rows, if needed,  */
    /* to put the pivot element on the diagonal. The columns are not      */
    /* physically interchanged, only relabeled: indxc[i], the column of   */
    /* the ith pivot element, is the ith column that is reduced, while    */
    /* indxr[i] is the row in which that pivot element was originally     */
    /* located. If indxr[i] != indxc[i] there is an implied column        */
    /* interchange. With this form of bookkeeping, the solution b's will  */
    /* end up in the correct order, and the inverse matrix will be        */
    /* scrambled by columns.                                              */
    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]) 
    } 
    /* We are now ready to divide the pivot row by the pivot element,     */
    /* located at irow and icol.                                          */
    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; 
    /* Next, we reduce the rows...                                        */
    /*...except for the pivot one, of course.                             */
    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; 
    } 
  } 
 /* This is the end of the main loop over columns of the reduction.       */
 /* It only remains to unscramble the solution in view of the column      */
 /* interchanges. We do this by interchanging pairs of columns in the     */
 /* reverse order that the permutation was built up.                      */
 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]]); 
 } 
   /* And we are done.                                                    */
   free_ivector(ipiv ,1,n); 
   free_ivector(indxr,1,n); 
   free_ivector(indxc,1,n); 
}

void ludcmp(double **a, int n, int *indx, double *d, int *singular) 
/* 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, a rranged 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.                                            */
{ 
  int i,imax,j,k ; 
  double big,dum,sum ,temp; 
  double *vv;               /* vv stores the implicit scaling of each row.   */

  vv=dvector(1 ,n) ; 
  *singular = 1;

  *d=1.0;                  /* No row interchanges yet.                      */
  for (i=1;i<=n; i++) {  
    /* Loop over ro ws to get the implicit scaling information.             */
    big=0.0; 
    for (j=1;j<=n;j ++) 
      if ((temp=fabs (a[ i] [j] )) > big) big=temp;
    if (big == 0.0)  /* Modification DG 27.5.2002 */
    {
      printf("Singular matrix in routine ludcmp\n");
      *singular = 0;
      return;
    }
    /* No nonzero largest element.                                          */
    vv[i]=1.0/ big ;   /* Save the scaling.                                 */
  } 
  for (j=1;j<=n; j++ ) {  
  /* This is the loop over columns of Crout's method.                       */
    for (i=1;i<j;i++) { 
    /* This is equation (2.3.12) except for i = j.                          */
      sum=a[i][j] ; 
      for (k=1;k<i;k ++) sum -= a[i][k]*a[k][j]; 
      a[i][j]=sum ; 
    } 
    big=0.0;   /* Initialize for the search for largest pivot element.      */
    for (i=j;i<=n;i ++) {   
    /* This is i = j of equation (2.3.12) and i = j + 1 ... N  of equation  */
    /* (2.3.13).                                                            */
      sum=a[i][j] ; 
      for (k=1;k<j;k ++) 
        sum -= a[i][k]*a[k][j]; 
        a[i][j]=sum ; 
        if ( (dum=vv[i]*fabs(sum) ) >= 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][k]; 
        a[imax][k]=a [j][k] ; 
        a[j][k]=dum; 
      } 
      *d = -(*d);        /* ...and change the parity of d .                 */
      vv[imax]=vv [j] ;  /* Also interchange the scale facto r.             */
    } 
    indx[j]=imax; 
    if (a[j][j] == 0.0) a[j][j]=TINY ; 
    /* If the pivot element is zero the matrix is singula r (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][j]) ; 
      for (i=j+1;i<= n;i++) a[i][j] *= dum; 
    } 
  } 
  /* Go back for the next column in the reduction.                          */
  free_dvector (vv ,1, n); 
} 

void lubksb(double **a, int n, int *indx, double b[]) 
/* Solves the set of n linea r 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.                                   */
{ 
  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. We now do the forward substitution,   */
  /* equation (2.3.6). The only new wrinkle is to unscramble the            */
  /* permutation as we go.                                                  */
    ip=indx[i]; 
    sum=b[ip]; 
    b[ip]=b[i]; 
    if (ii) 
      for (j=ii;j<=i-1; j++) sum -= a[i][j]*b[j]; 
    else if (sum) ii=i; 
    /* A nonzero element w as encountered, so from now on we will have to   */
    /* do the sums in the loop above.                                       */ 
    b[i]=sum; 
  } 
  for (i=n;i>=1; i--) { 
  /* Now we do the backsubstitution, equation (2.3.7).                      */
    sum=b[i]; 
    for (j=i+1;j<=n ;j++) sum -= a[i][j]*b[j]; 
    b[i]=sum/a[i][i]; 
    /* Store a component of the solution vector X .                         */
  }   /* All done!                                                          */
}

#define FREERETURN {free_dmatrix(fjac,1,n,1,n);free_dvector(fvec,1,n);\
    free_dvector(p,1,n);free_ivector(indx,1,n);return;} 

void mnewt(int ntrial, double x[], int n, double tolx, double tolf,
     int *singular,
     void (*usrfun)(double *x, int n, double *fvec, double **fjac, void *ptr), 
     void *ptr)
/* Given an initial guess x[1..n] for a root in n dimensions, take ntrial  */
/* Newton-Raphson steps to improve the root. Stop if the root converges in */
/* either summed absolute variable increments tolx or summed absolute      */
/* function values tolf.                                                   */
{ 
  int k,i,*indx; 
  double errx,errf,d,*fvec,**fjac,*p; 
  indx=ivector(1,n) ; 
  p=dvector(1, n);
  fvec=dvector(1, n); 
  fjac=dmatrix(1, n, 1 ,n) ; 
  for (k=1;k<=ntrial;k++) { 
    usrfun(x, n ,fvec ,fjac, ptr);
/* User function supplies function values at x in fvec and Jacobian       */
/* matrix in fjac.                                                        */
    errf=0.0;
    for (i=1;i<=n;i++) errf += fabs(fvec[i]);  
/*  Check function convergence.                                           */
    if (errf <= tolf) FREERETURN 
    for (i=1;i<=n;i ++) p[i] = -fvec[i]; 
/*  Right-hand side of linear equations.                                  */
    ludcmp(fjac,n ,indx ,&d, singular); 
/*  Solve linea r equations using LU decomp osition.                      */
    lubksb(fjac,n,indx,p); 
    errx=0.0; 
/*  Check root convergence.                                               */
    for (i=1;i<=n;i ++) {           /* Update solution.                   */
      errx += fabs(p[i]); 
      x[i] += p[i]; 
    } 
    if (errx <= tolx) FREERETURN 
  } 
  FREERETURN 
}

void fdjac(int n, double x[], double fvec[], double **df,
     void (*vecfunc)(int, double [], double [], void *ptr), void *ptr)
/* Computes forward-difference approximation to Jacobian. On input, x[1..n] */ 
/* is the point at which the Jacobian is to be evaluated, fvec[1..n] is the */
/* vector of function values at the point, and vecfunc(n,x,f) is a          */
/* user-supplied routine that returns the vector of functions at x.         */
/* On output, df[1..n][ 1.. n] is the Jacobian array.                       */
{ 
  int i,j;
  double h,temp,*f; 

  f=dvector(1, n); 
  for (j=1;j<=n; j++ ) { 
    temp=x[j]; 
    h=EPS*fabs(temp); 
    if (h==0.0) h=EPS; 
    x[j]=temp+h;         /* Trick to reduce finite precision error.       */ 
    h=x[j]-temp; 
    (*vecfunc) (n,x,f,ptr); 
    x[j]=temp; 
    for (i=1;i<=n;i++) df[i][j]=(f[i]-fvec[i])/h;  
/*  Forward difference formula.                                           */
  }
  free_dvector(f,1,n);
}

/* Here GOLD is the default ratio by which successive intervals are         */
/* magnified; GLIMIT is the maximum magnification allowed for a             */
/* parabolic-fit step.                                                      */

void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, 
     double *fc, double (*func)(double, void *), void *ptr) 
/* Given a function func, and given distinct initial points ax and bx, this */
/* routine searches in the downhill direction (defined by the function as   */
/* evaluated at the initial points) and returns new points ax, bx, cx that  */
/* bracket a minimum of the function. Also returned are the function values */
/* at the three points, fa, fb, and fc.                                     */
{ 
  double ulim,u,r,q,fu,dum;

  *fa=(*func)(*ax, ptr); 
  *fb=(*func)(*bx, ptr); 
printf("ax=%lf fa=%lf bx=%lf fb=%lf\n", *ax, *fa, *bx, *fb);
  if (*fb > *fa) { 
  /* Switch roles of a and b so that we can go  */
  /* downhill in the direction from a to b .    */
  SHFT(dum, *ax, *bx, dum) 
  SHFT(dum, *fb, *fa, dum) 
  } 
printf("ax=%lf fa=%lf bx=%lf fb=%lf\n", *ax, *fa, *bx, *fb);
  *cx=(*bx)+GOLD*(*bx - *ax);   /* First guess for c. */
printf("cx=%lf\n", *cx);
  *fc=(*func)(*cx, ptr); 
printf("ax=%lf fa=%lf bx=%lf fb=%lf cx=%lf fc=%lf\n", *ax, *fa, *bx, *fb,
*cx, *fc);
  while (*fb > *fc) {  /* Keep returning here until we bracket. */
  r=(*bx - *ax)*( *fb - *fc) ;  /* Compute u by parabolic extrapolation from */
  q=(*bx - *cx)*( *fb - *fa) ;  /* a, b, c . TINY is used to prevent any     */
  u=(*bx)-((*bx - *cx)*q-(*bx - *ax)*r)/    /* possible division by zero.    */
  (2.0*SIGN(DMAX(fabs(q-r),TINY),q-r));
  ulim=(*bx)+GLIMIT*(*cx - *bx); 
  /* We won't go farther than this. Test various possibilities:  */
  if ((*bx-u)*(u - *cx) > 0.0) { /* Parabolic u is between b and c : try it. */
    fu=(*func)(u, ptr); 
printf("u=%lf fu=%lf\n", u, fu);
    if (fu < *fc) {   /* Got a minimum between b and c . */ 
      *ax=(*bx); 
      *bx=u; 
      *fa=(*fb); 
      *fb=fu; 
      return; 
    } 
    else if (fu > *fb) { /* Got a minimum between a and u . */
      *cx=u; 
      *fc=fu; 
      return; 
    } 
    u=(*cx)+GOLD*(*cx - *bx);       /* Parabolic fit was no use. Use default */
    fu=(*func)(u, ptr);             /* magnification.                        */
printf("u=%lf fu=%lf\n", u, fu);
  } 
  else if ((*cx-u)*(u-ulim) > 0.0) {  /* Parabolic fit is between c and its  */
    fu=(*func)(u, ptr);               /* allowed limit.                      */
printf("u=%lf fu=%lf\n", u, fu);
    if (fu < *fc) { 
      SHFT(*bx, *cx ,u,  *cx+GOLD*(*cx - *bx)) 
      SHFT(*fb, *fc ,fu ,(*func)(u, ptr)) 
printf("u=%lf fu=%lf\n", u, fu);
    } 
  } else if ((u-ulim)*( ulim- *cx) >= 0.0) {  /* Limit parabolic u to   */
      u=ulim;                                 /* maximum allowed value. */
      fu=(*func)(u, ptr); 
printf("u=%lf fu=%lf\n", u, fu);
    } 
    else {   /* Reject parabolic u , use default magnification. */
      u=(*cx)+GOLD*( *cx - *bx); 
      fu=(*func)(u, ptr); 
printf("u=%lf fu=%lf\n", u, fu);
    } 
    SHFT(*ax, *bx, *cx, u)  /* Eliminate oldest point and continue. */
    SHFT(*fa,* fb, *fc, fu) 
  } 
} 

double golden(double ax, double bx, double cx, double (*f)(double, void *), 
       double tol, double *xmin, void *ptr) 
/* Given a function f, and given a bracketing triplet of abscissas ax, bx,  */
/* cx (such that bx is between ax and cx, and f(bx) is less than both f(ax) */
/* and f(cx)), this routine performs a  golden section search for the       */
/* minimum, isolating it to a fractional precision of about tol. The        */
/* abscissa of the minimum is returned as xmin , and the minimum function   */
/* value is returned as golden , the returned function value.               */
{ 
  double f1,f2,x0,x1,x2,x3;  

  /* At any given time we will keep track of four points, x0,x1,x2,x3.      */ 
  x3=cx; 
  if (fabs(cx-bx) > fabs(bx-ax)) {  /* Make x0 to x1 the smaller segment,   */
    x1=bx;                          /* and fill in the new point to be      */
    x2=bx+Cgolden*(cx-bx);                /* tried.                               */
  } 
  else { 
    x2=bx; 
    x1=bx-Cgolden*(bx-ax); 
  } 
  /* The initial function evaluations. Note that we never     */
  /* need to evaluate the function at the original endpoints. */
  f1=(*f)(x1, ptr);  
  f2=(*f)(x2, ptr);  
  while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) { 
    if (f2 < f1) {                      /* One possible outcome,            */ 
      SHFT3(x0,x1,x2,Rgolden*x1+Cgolden*x3)         /* its housekeeping,                */
      SHFT2(f1,f2,(*f)(x2, ptr))        /* and a new function evaluation.   */
    } else {                            /* The other outcome,               */
      SHFT3(x3,x2,x1,Rgolden*x2+Cgolden*x0) 
      SHFT2(f2,f1,(*f)(x1, ptr))        /* and its new function evaluation. */
    } 
  }                                     /* Back to see if we are done.      */
  if (f1 < f2) { /* We are done. Output the best of the two current values. */
    *xmin=x1; 
    return f1; 
  } else { 
    *xmin=x2; 
    return f2; 
  } 
} 


/* Here ITMAX is the maximum allowed number of iterations; CGOLD is the      */
/* golden ratio; ZEPS is a small number that protects against trying to      */
/* achieve fractional accuracy for a minimum that happens to be exactly zero.*/
double brent(double ax, double bx, double cx, double (*f)(double, void *), 
       double tol, double *xmin, void *ptr) 
/* Given a function f, and given a bracketing triplet of abscissas ax, bx,   */
/* cx (such that bx is between ax and cx, and f(bx) is less than both f(ax)  */
/* and f(cx)), this routine isolates the minimum to a fractional precision   */
/* of about tol using Brent's method. The abscissa of the minimum is         */
/* returned as xmin , and the minimum function value is returned as brent ,  */
/* the returned function value.                                              */
{ 
  int iter; 
  double a,b,d,etemp=0,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; 
  double e=0.0;  /* This will be the distance moved on the step before last. */
  a=(ax < cx ? ax : cx);  /* a and b must be in ascending order,             */
  b=(ax > cx ? ax : cx);  /* but input abscissas need not be.                */
  x=w=v=bx;               /* Initializations...                              */
  fw=fv=fx=(*f)(x, ptr); 
  for (iter=1; iter<=ITMAX; iter++) {               /* Main program loop.    */
    xm=0.5*(a+b); 
    tol2=2.0*(tol1=tol*fabs(x)+ZEPS); 
    if (fabs(x-xm) <= (tol2-0.5*(b-a))) {          /* Test for done here.    */
      *xmin=x; 
      return fx; 
    } 
    if (fabs(e) > tol1) {              /* Construct a trial parabolic fit.   */
      r=(x-w)*(fx-fv); 
      q=(x-v)*(fx-fw); 
      p=(x-v)*q-(x-w)*r; 
      q=2.0*(q-r) ; 
      if (q > 0.0) p = -p; 
      q=fabs(q); 
      if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) 
        d=CGOLD*(e=( x >= xm ? a-x : b-x)); 
      /* The above conditions determine the acceptability of the parabolic  */
      /* fit. Here we take the golden section step into the larger of the   */
      /* two segments.                                                      */
      else { 
        d=p/q;                                /* Take the parabolic step.   */
        u=x+d; 
        if (u-a < tol2 || b-u < tol2) 
          d=SIGN(tol1 ,xm-x ); 
        } 
    } 
    else { 
      d=CGOLD*(e= (x >= xm ? a-x : b-x)); 
    } 
    u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); 
    fu=(*f)(u, ptr); /* This is the one function evaluation per iteration.  */
    if (fu <= fx) {  /* Now decide what to do with our function evaluation. */
      if (u >= x) a=x; else b=x; 
      SHFT(v,w,x,u)                      /*Housekeepingfollows:             */ 
      SHFT(fv,fw,fx,fu) 
    } 
    else { 
      if (u < x) a=u; else b=u; 
      if (fu <= fw || w == x) { 
        v=w; 
        w=u; 
        fv=fw; 
        fw=fu; 
      } else if (fu <= fv || v == x || v == w) { 
        v=u; 
        fv=fu; 
      } 
    }                /* Done with housekeeping. Back for another iteration. */
  } 
  nrerror("Too many iterations in brent"); 

  return 0.0;
} 

double rtsafe(void (*funcd)(double, double *, double *, void *ptr), double x1,
        double x2, double xacc, void *ptr)
{
	void nrerror(char error_text[]);
	int j;
	double df,dx,dxold,f,fh,fl;
	double temp,xh,xl,rts;

	(*funcd)(x1,&fl,&df, ptr);
	(*funcd)(x2,&fh,&df, ptr);
	if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
		nrerror("Root must be bracketed in rtsafe");
	if (fl == 0.0) return x1;
	if (fh == 0.0) return x2;
	if (fl < 0.0) {
		xl=x1;
		xh=x2;
	} else {
		xh=x1;
		xl=x2;
	}
	rts=0.5*(x1+x2);
	dxold=fabs(x2-x1);
	dx=dxold;
	(*funcd)(rts,&f,&df, ptr);
	for (j=1;j<=MAXIT;j++) {
		if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
			|| (fabs(2.0*f) > fabs(dxold*df))) {
			dxold=dx;
			dx=0.5*(xh-xl);
			rts=xl+dx;
			if (xl == rts) return rts;
		} else {
			dxold=dx;
			dx=f/df;
			temp=rts;
			rts -= dx;
			if (temp == rts) return rts;
		}
		if (fabs(dx) < xacc) return rts;
		(*funcd)(rts,&f,&df, ptr);
		if (f < 0.0)
			xl=rts;
		else
			xh=rts;
	}
	nrerror("Maximum number of iterations exceeded in rtsafe");
	return 0.0;
}

double zbrent(double (*func)(double, void *ptr), double x1, double x2, 
       double tol, void *ptr)
{
	int iter;
	double a=x1,b=x2,c=x2,d,e,min1,min2;
	double fa=(*func)(a, ptr),fb=(*func)(b, ptr),fc,p,q,r,s,tol1,xm;

	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
		nrerror("Root must be bracketed in zbrent");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPSzbrent*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0) return b;
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0) q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += SIGN(tol1,xm);
		fb=(*func)(b, ptr);
	}
	nrerror("Maximum number of iterations exceeded in zbrent");
	return 0.0;
}

#define RADIX 2.0

void balanc(double **a, int n)
{
	int last,j,i;
	double s,r,g,f,c,sqrdx;

	sqrdx=RADIX*RADIX;
	last=0;
	while (last == 0) {
		last=1;
		for (i=1;i<=n;i++) {
			r=c=0.0;
			for (j=1;j<=n;j++)
				if (j != i) {
					c += fabs(a[j][i]);
					r += fabs(a[i][j]);
				}
			if (c && r) {
				g=r/RADIX;
				f=1.0;
				s=c+r;
				while (c<g) {
					f *= RADIX;
					c *= sqrdx;
				}
				g=r*RADIX;
				while (c>g) {
					f /= RADIX;
					c /= sqrdx;
				}
				if ((c+r)/f < 0.95*s) {
					last=0;
					g=1.0/f;
					for (j=1;j<=n;j++) a[i][j] *= g;
					for (j=1;j<=n;j++) a[j][i] *= f;
				}
			}
		}
	}
}
#undef RADIX

void elmhes(double **a, int n)
{
	int m,j,i;
	double y,x;
        double temp;

	for (m=2;m<n;m++) {
		x=0.0;
		i=m;
		for (j=m;j<=n;j++) {
			if (fabs(a[j][m-1]) > fabs(x)) {
				x=a[j][m-1];
				i=j;
			}
		}
		if (i != m) {
			for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j])
			for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m])
		}
		if (x) {
			for (i=m+1;i<=n;i++) {
				if ((y=a[i][m-1]) != 0.0) {
					y /= x;
					a[i][m-1]=y;
					for (j=m;j<=n;j++)
						a[i][j] -= y*a[m][j];
					for (j=1;j<=n;j++)
						a[j][m] += y*a[j][i];
				}
			}
		}
	}
}

#define NRANSI

void hqr(double **a, int n, double wr[], double wi[])
{
	int nn,m,l,k,j,its,i,mmin;
	double z,y,x,w,v,u,t,s,r,q,p,anorm;

	anorm=0.0;
	for (i=1;i<=n;i++)
		for (j=IMAX(i-1,1);j<=n;j++)
			anorm += fabs(a[i][j]);
	nn=n;
	t=0.0;
	while (nn >= 1) {
		its=0;
		do {
			for (l=nn;l>=2;l--) {
				s=fabs(a[l-1][l-1])+fabs(a[l][l]);
				if (s == 0.0) s=anorm;
				if ((double)(fabs(a[l][l-1]) + s) == s) break;
			}
			x=a[nn][nn];
			if (l == nn) {
				wr[nn]=x+t;
				wi[nn--]=0.0;
			} else {
				y=a[nn-1][nn-1];
				w=a[nn][nn-1]*a[nn-1][nn];
				if (l == (nn-1)) {
					p=0.5*(y-x);
					q=p*p+w;
					z=sqrt(fabs(q));
					x += t;
					if (q >= 0.0) {
						z=p+SIGN(z,p);
						wr[nn-1]=wr[nn]=x+z;
						if (z) wr[nn]=x-w/z;
						wi[nn-1]=wi[nn]=0.0;
					} else {
						wr[nn-1]=wr[nn]=x+p;
						wi[nn-1]= -(wi[nn]=z);
					}
					nn -= 2;
				} else {
					if (its == 30) nrerror("Too many iterations in hqr");
					if (its == 10 || its == 20) {
						t += x;
						for (i=1;i<=nn;i++) a[i][i] -= x;
						s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
						y=x=0.75*s;
						w = -0.4375*s*s;
					}
					++its;
					for (m=(nn-2);m>=l;m--) {
						z=a[m][m];
						r=x-z;
						s=y-z;
						p=(r*s-w)/a[m+1][m]+a[m][m+1];
						q=a[m+1][m+1]-z-r-s;
						r=a[m+2][m+1];
						s=fabs(p)+fabs(q)+fabs(r);
						p /= s;
						q /= s;
						r /= s;
						if (m == l) break;
						u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
						v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
						if ((double)(u+v) == v) break;
					}
					for (i=m+2;i<=nn;i++) {
						a[i][i-2]=0.0;
						if (i != (m+2)) a[i][i-3]=0.0;
					}
					for (k=m;k<=nn-1;k++) {
						if (k != m) {
							p=a[k][k-1];
							q=a[k+1][k-1];
							r=0.0;
							if (k != (nn-1)) r=a[k+2][k-1];
							if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) {
								p /= x;
								q /= x;
								r /= x;
							}
						}
						if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {
							if (k == m) {
								if (l != m)
								a[k][k-1] = -a[k][k-1];
							} else
								a[k][k-1] = -s*x;
							p += s;
							x=p/s;
							y=q/s;
							z=r/s;
							q /= p;
							r /= p;
							for (j=k;j<=nn;j++) {
								p=a[k][j]+q*a[k+1][j];
								if (k != (nn-1)) {
									p += r*a[k+2][j];
									a[k+2][j] -= p*z;
								}
								a[k+1][j] -= p*y;
								a[k][j] -= p*x;
							}
							mmin = nn<k+3 ? nn : k+3;
							for (i=l;i<=mmin;i++) {
								p=x*a[i][k]+y*a[i][k+1];
								if (k != (nn-1)) {
									p += z*a[i][k+2];
									a[i][k+2] -= p*r;
								}
								a[i][k+1] -= p*q;
								a[i][k] -= p;
							}
						}
					}
				}
			}
		} while (l < nn-1);
	}
}

#undef NRANSI

void tred2(double **a, int n, double d[], double e[])
{
	int l,k,j,i;
	double scale,hh,h,g,f;

	for (i=n;i>=2;i--) {
		l=i-1;
		h=scale=0.0;
		if (l > 1) {
			for (k=1;k<=l;k++)
				scale += fabs(a[i][k]);
			if (scale == 0.0)
				e[i]=a[i][l];
			else {
				for (k=1;k<=l;k++) {
					a[i][k] /= scale;
					h += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
				e[i]=scale*g;
				h -= f*g;
				a[i][l]=f-g;
				f=0.0;
				for (j=1;j<=l;j++) {
					a[j][i]=a[i][j]/h;
					g=0.0;
					for (k=1;k<=j;k++)
						g += a[j][k]*a[i][k];
					for (k=j+1;k<=l;k++)
						g += a[k][j]*a[i][k];
					e[j]=g/h;
					f += e[j]*a[i][j];
				}
				hh=f/(h+h);
				for (j=1;j<=l;j++) {
					f=a[i][j];
					e[j]=g=e[j]-hh*f;
					for (k=1;k<=j;k++)
						a[j][k] -= (f*e[k]+g*a[i][k]);
				}
			}
		} else
			e[i]=a[i][l];
		d[i]=h;
	}
	d[1]=0.0;
	e[1]=0.0;
	/* Contents of this loop can be omitted if eigenvectors not
			wanted except for statement d[i]=a[i][i]; */
	for (i=1;i<=n;i++) {
		l=i-1;
		if (d[i]) {
			for (j=1;j<=l;j++) {
				g=0.0;
				for (k=1;k<=l;k++)
					g += a[i][k]*a[k][j];
				for (k=1;k<=l;k++)
					a[k][j] -= g*a[k][i];
			}
		}
		d[i]=a[i][i];
		a[i][i]=1.0;
		for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;
	}
}

double pythag(double a, double b)
{
	double absa,absb;
	absa=fabs(a);
	absb=fabs(b);
	if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
	else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}

#define NRANSI

void tqli(double d[], double e[], int n, double **z)
{
	double pythag(double a, double b);
	int m,l,iter,i,k;
	double s,r,p,g,f,dd,c,b;

	for (i=2;i<=n;i++) e[i-1]=e[i];
	e[n]=0.0;
	for (l=1;l<=n;l++) {
		iter=0;
		do {
			for (m=l;m<=n-1;m++) {
				dd=fabs(d[m])+fabs(d[m+1]);
				if ((double)(fabs(e[m])+dd) == dd) break;
			}
			if (m != l) {
				if (iter++ == 30) nrerror("Too many iterations in tqli");
				g=(d[l+1]-d[l])/(2.0*e[l]);
				r=pythag(g,1.0);
				g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
				s=c=1.0;
				p=0.0;
				for (i=m-1;i>=l;i--) {
					f=s*e[i];
					b=c*e[i];
					e[i+1]=(r=pythag(f,g));
					if (r == 0.0) {
						d[i+1] -= p;
						e[m]=0.0;
						break;
					}
					s=f/r;
					c=g/r;
					g=d[i+1]-p;
					r=(d[i]-g)*s+2.0*c*b;
					d[i+1]=g+(p=s*r);
					g=c*r-b;
					for (k=1;k<=n;k++) {
						f=z[k][i+1];
						z[k][i+1]=s*z[k][i]+c*f;
						z[k][i]=c*z[k][i]-s*f;
					}
				}
				if (r == 0.0 && i >= l) continue;
				d[l] -= p;
				e[l]=g;
				e[m]=0.0;
			}
		} while (m != l);
	}
}
#undef NRANSI