/* 
	Multi-scale time series
	Ken Kiyono
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"
#include "intde2.c"
#define FUNC(x) ((*funk)(-log(x))/(x))
#define PI 3.1415926535897932384626433832795
#define SQR2PI 0.39894228040143267793994605993439
#define NUMAW 100000

void help(void);
int read_data(char filename[]);
void svdfit(double x[], double y[], int ndata, double a[], int ma, double *chisq);
void svdvar(double **v, int ma, double w[], double **cvm);
double polyfunc(double x, double p[], int np);
void lfit(double x[], double y[], int ndat, double a[], int ma, double *chisq);
void sort(unsigned long n, double arr[]);
double P_L(double s, double sgm);
double Gau(double x);
double G_r(double sg);
void intdeini(int lenaw, double tiny, double eps, double *aw);
void intdei(double (*f)(double), double a, 
        double *aw, double *nint, double *err);
int pdf(long n, double x[], int bin, double y[],double py[]);
/* Global variables */
char *pname;
double *data;
double extx;
double lambda;
double lambda2;
int nfunc;
double tiny, *aw, nint, err;
int lenaw;

main(int argc, char **argv)
{
    long nn,i,j,i0;
	long scale, ncg;
	long fsamp,scalef,nbin;
	int dp,pn;
	int methd,bck, cnt;
	double average, sd;
	double *px,*py,*pa,chisq;
	double *cgs, *iy;
	double logvar,logmean,lmd,mu3,tmpv;
	double *u,*Pu,*Put;
	double max,min,dx,shfty;
	char fname1[256];

	/* Read and interpret the command line. */
    pname = argv[0];
    scale = 1;
    dp = 0; /* 多項式の次数 */
    methd = 1;
    bck = 0;
    cnt = 0;
    fsamp = 1; /* サンプリング周波数 */
    nbin = 100;
    shfty = 1;
    for (i = 1; i < argc && *argv[i] == '-'; i++) {
      switch(argv[i][1]) {
  	case 's':
		scale = atol(argv[++i]);
		if(scale < 1) scale = 1;
		break;
 	case 'r':
		shfty = atof(argv[++i]);
		if(shfty == 0) shfty = 1;
		break;
  	case 'n':
		nbin = atol(argv[++i]);
		if(nbin < 10) scale = 10;
		break;
  	case 'f':
		fsamp = atol(argv[++i]);
		if(fsamp < 1) fsamp = 1;
		break;
	case 'd':
		dp = atoi(argv[++i]);
		break;
	case 'b':
		bck = 1;
		break;
	case 'c':
		cnt = 1;
		break;
	case 'm':
		methd = atoi(argv[++i]);
		if(methd < 0){
			methd = 1;
		}else if(methd > 2) methd = 2;
		break;
	case 'h':	/* print usage information and quit */
        default:
	  help();
	  exit(1);
      }
    }
    
    scalef = scale*fsamp;
	if(dp > scalef - 1) dp = scalef - 1;
	
	if(argc < 2){
	  help();
	  exit(1);
	}
	strcpy(fname1,argv[argc-1]);
	nn = read_data(fname1);
	dp = dp + 1; /* 多項式のパラメータ数 */
	
	average = 0;
	for(i=1;i<=nn;i++){
		average += data[i];
	}
	average = average/(double)nn;

	for(i=1;i<=nn;i++){
		data[i] = data[i] - average;
	}
	
	fprintf(stderr, "data length: %ld\n", nn);
	
	fprintf(stderr, "scale: %ld (sampling freq = %ld)\n", scale, fsamp);
	cgs = vector(1,2*(nn-scalef+1));

	if(dp <= 1){ /* detrendなし */
		fprintf(stderr, "WITHOUT detrending\n");
		ncg = nn - scalef + 1;

		for(i=1;i<=ncg;i++){
			cgs[i] = 0;
			for(j=i;j<=i+scalef-1;j++){
				cgs[i] += data[j];
			}
		}	
	
	}else{
		if(nn < 2*scalef) nrerror("The data length is not enough!\n");
		fprintf(stderr, "detrend using polynomial fits of order %d\n", dp-1);

		iy = vector(1,nn);
		pn = dp;
		px = vector(1,2*scalef);
		py = vector(1,2*scalef);
		pa = vector(1,pn);
		ncg = 0;

				
		iy[1] = data[1];

		for(i=2;i<=nn;i++){
			iy[i] = iy[i-1] + data[i];
		}

		for(i=1;i<=2*scalef;i++) px[i] = (double)i;
		    

		for(i=1; i<= nn - 2*scalef + 1; i+=scalef){
			for(j=1;j<=2*scalef;j++){
				py[j]=iy[j+i-1];
			}
			if(methd > 1){
				svdfit(px, py, 2*scalef, pa, pn, &chisq);
			}else{
				lfit(px, py, 2*scalef, pa, pn, &chisq);
			}
			
			for(j=1;j<=2*scalef;j++){
				py[j]=py[j] - polyfunc(px[j],pa,pn);
			}
		
			for(j=1;j<=scalef;j++){
				ncg++;
				cgs[ncg] = py[j+scalef] - py[j];
			}
		}

		if(bck > 0 && nn%scalef > 0){/* 後ろ向きに粗視化  */
			fprintf(stderr,"including backward computation\n");
			iy[1] = data[nn];
			for(i=2;i<=nn;i++){
				iy[i] = iy[i-1] + data[nn-i+1];
			}
			for(i=1; i<= nn - 2*scalef + 1; i+=scalef){
				for(j=1;j<=2*scalef;j++){
					py[j]=iy[j+i-1];
				}
				if(methd > 1){
					svdfit(px, py, 2*scalef, pa, pn, &chisq);
				}else{
					lfit(px, py, 2*scalef, pa, pn, &chisq);
				}
			
				for(j=1;j<=2*scalef;j++){
					py[j]=py[j] - polyfunc(px[j],pa,pn);
				}
		
				for(j=1;j<=scalef;j++){
					ncg++;
					cgs[ncg] = py[j+scalef] - py[j];
				}
			}
		}
		
		free_vector(px,1,2*scalef);
		free_vector(py,1,2*scalef);
		free_vector(pa,1,pn);
		free_vector(iy,1,nn);
	}
	
	if(cnt < 1){
		average = 0;
		for(i=1;i<=ncg;i++){
			average += cgs[i];
		}
		average = average/(double)ncg;
	}else{
		sort(ncg, cgs);
		if(ncg%2 > 0){
			average = cgs[(long)(ncg/2)+1];
		}else{
			average = (cgs[ncg/2]+cgs[ncg/2+1])/2;
		}
	}
		
	for(i=1;i<=ncg;i++){
		cgs[i] = cgs[i] - average; /* zero mean */
	}

	/* 対数振幅分散の計算 */
	logmean = 0;
	for(i = 1;i <= ncg; i++){
		logmean += log(fabs(cgs[i]));
	}
	logmean = logmean/(double)ncg;

	logvar = 0;
	mu3 = 0;
	for(i = 1;i <= ncg; i++){
		tmpv = (log(fabs(cgs[i])) - logmean)*(log(fabs(cgs[i])) - logmean);
		logvar += tmpv;
		mu3 += tmpv*(log(fabs(cgs[i])) - logmean);
	}
	
	logvar = logvar/(double)ncg - 1.2337005501361698273543113749845;
	mu3 = mu3/(double)ncg + 2.10359958052929;
	fprintf(stderr, "mu3 = %lf\n", mu3);

	if(logvar > 0){
		 lmd = sqrt(logvar);
	}else{
		lmd = 0;
	}
	fprintf(stderr, "lmd = %lf\n", lmd);

	/* 対数振幅分散の計算（ここまで） */
	
	/* 標準化された密度関数の推定 */
	sd = 0;
	for(i=1;i<=ncg;i++){
		sd += cgs[i]*cgs[i]; /* zero mean */
	}
	sd = sqrt(sd/(double)ncg);
	
	for(i=1;i<=ncg;i++){
		cgs[i]=cgs[i]/sd; /* unit variance */
	}

	
	u = vector(1,nbin);
	Pu = vector(1,nbin);
	Put = vector(1,nbin);
	
	i0 = pdf(ncg, cgs, nbin, u, Pu);
	
    lenaw = NUMAW;
	aw = vector(1,NUMAW);
    tiny = 1.0e-307;
	lambda = lmd;
	lambda2 = lambda*lambda;
	
    intdeiini(lenaw, tiny, 1.0e-16, aw);
	for(i = 1; i <= nbin; i++){
				
				extx = u[i];
				nfunc = 0;
				intdei(G_r, 0.0, aw, &nint, &err);
				Put[i] = nint;
	}	

	
	for(i=1;i<=nbin;i++){
		fprintf(stdout, "%lf\t%8.7lf\t%12.11lf\n", u[i], Pu[i]/shfty, Put[i]/shfty);
	}
	
		
	free_vector(data,1,nn);
	free_vector(cgs,1,2*(nn-scale+1));
	free_vector(u,1,nbin);
	free_vector(Pu,1,nbin);
	free_vector(Put,1,nbin);
	free_vector(aw,1,NUMAW);
    exit(0);
}

int read_data(char filename[]){
	/* データの入力 */
    long maxdat = 0L, N = 0L;
	double y;
	FILE *fp1;
		
    if((fp1 = fopen(filename,"r")) == NULL){
		fprintf(stderr,"File %s not found!\n",filename);
		return 0;
	}else{
		fprintf(stderr,"Input file: %s.\n",filename);
	}

	data = vector(1,10000);
	
	while (fscanf(fp1,"%lf", &y) == 1) {
        if (++N >= maxdat) {	    
			double *s;
			maxdat += 10000;	/* allow the input buffer to grow (the
				   increment is arbitrary) */
			if ((s = realloc(data, maxdat * sizeof(double))) == NULL) {
				fprintf(stderr,
					"%s: insufficient memory, truncating input at row %d\n",
					pname, N);
				break;
			}
			data = s;
		}
		data[N] = y;
    }
    
    fclose(fp1);
	if (N < 1){
		fprintf(stderr,"No data was read!\n");
	    return 0;
	}

	return N;
}


static char *help_strings[] = {
 "使用法: %s [OPTIONS ...] データファイル名\n",
 "設定できるOPTION:",
 " -n B             binの総数B",
 " -s W             粗視化スケール (サンプリング周波数を単位とする)",
 " -f n             サンプリング周波数 (デフォルトは1, 自然数のみ)",
 " -d N             多項式トレンドの次数N",
 " -c               中央値を中心に計算 (デフォルトでは平均値が中心)", 
 " -m M             最小２乗フィティングアルゴリズム(M=1:正規方程式，M=2:特異値分解)",
 " -b               後ろ向きの粗視化を追加",
 " -r Y             確率分布P(x)にYで割って出力",
 " -h               オプションの一覧を表示",
NULL
};

void help(void)
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}


#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
#define TOL 1.0e-10 /* 特異値分解 */

void svdfit(double x[], double y[], int ndata, double a[], int ma, double *chisq)
{
	void svbksb(double **u, double w[], double **v, int m, int n, double b[],
		double x[]);
	void fpoly(double x, double p[], int np);
	void svdcmp(double **a, int m, int n, double w[], double **v);
	int j,i;
	double wmax,thresh,sum,*b,*afunc,tmp;
	double **u, **v, *w;

	b=vector(1,ndata);
	afunc=vector(1,ma);
	u = matrix(1,ndata,1,ma);
	v = matrix(1,ma,1,ma);
	w = vector(1,ma);
	
	for (i=1;i<=ndata;i++) {
		fpoly(x[i],afunc,ma);
		for (j=1;j<=ma;j++) u[i][j]=afunc[j];
		b[i]=y[i];
	}
	svdcmp(u,ndata,ma,w,v);
	wmax=0.0;
	for (j=1;j<=ma;j++)
		if (w[j] > wmax) wmax=w[j];
	thresh=TOL*wmax;
	for (j=1;j<=ma;j++)
		if (w[j] < thresh) w[j]=0.0;
	svbksb(u,w,v,ndata,ma,b,a);
	*chisq=0.0;
	for (i=1;i<=ndata;i++) {
		fpoly(x[i],afunc,ma);
		for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];
		*chisq += (tmp=(y[i]-sum),tmp*tmp);
	}
	free_vector(afunc,1,ma);
	free_vector(b,1,ndata);
	free_matrix(u,1,ndata,1,ma);
	free_matrix(v, 1,ma,1,ma);
	free_vector(w,1,ma);

}


void svdvar(double **v, int ma, double w[], double **cvm)
{
	int k,j,i;
	double sum,*wti;

	wti=vector(1,ma);
	for (i=1;i<=ma;i++) {
		wti[i]=0.0;
		if (w[i]) wti[i]=1.0/(w[i]*w[i]);
	}
	for (i=1;i<=ma;i++) {
		for (j=1;j<=i;j++) {
			for (sum=0.0,k=1;k<=ma;k++) sum += v[i][k]*v[j][k]*wti[k];
			cvm[j][i]=cvm[i][j]=sum;
		}
	}
	free_vector(wti,1,ma);
}

void svbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
{
	int jj,j,i;
	double s,*tmp;

	tmp=vector(1,n);
	for (j=1;j<=n;j++) {
		s=0.0;
		if (w[j]) {
			for (i=1;i<=m;i++) s += u[i][j]*b[i];
			s /= w[j];
		}
		tmp[j]=s;
	}
	for (j=1;j<=n;j++) {
		s=0.0;
		for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
		x[j]=s;
	}
	free_vector(tmp,1,n);
}

void svdcmp(double **a, int m, int n, double w[], double **v)
{
	double pythag(double a, double b);
	int flag,i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z,*rv1;

	rv1=vector(1,n);
	g=scale=anorm=0.0;
	for (i=1;i<=n;i++) {
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i <= m) {
			for (k=i;k<=m;k++) scale += fabs(a[k][i]);
			if (scale) {
				for (k=i;k<=m;k++) {
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
				}
				for (k=i;k<=m;k++) a[k][i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i <= m && i != n) {
			for (k=l;k<=n;k++) scale += fabs(a[i][k]);
			if (scale) {
				for (k=l;k<=n;k++) {
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][l]=f-g;
				for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
				for (j=l;j<=m;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
					for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l;k<=n;k++) a[i][k] *= scale;
			}
		}
		anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n;i>=1;i--) {
		if (i < n) {
			if (g) {
				for (j=l;j<=n;j++)
					v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
		}
		v[i][i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=IMIN(m,n);i>=1;i--) {
		l=i+1;
		g=w[i];
		for (j=l;j<=n;j++) a[i][j]=0.0;
		if (g) {
			g=1.0/g;
			for (j=l;j<=n;j++) {
				for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
			}
			for (j=i;j<=m;j++) a[j][i] *= g;
		} else for (j=i;j<=m;j++) a[j][i]=0.0;
		++a[i][i];
	}
	for (k=n;k>=1;k--) {
		for (its=1;its<=30;its++) {
			flag=1;
			for (l=k;l>=1;l--) {
				nm=l-1;
				if ((double)(fabs(rv1[l])+anorm) == anorm) {
					flag=0;
					break;
				}
				if ((double)(fabs(w[nm])+anorm) == anorm) break;
			}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l;i<=k;i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((double)(fabs(f)+anorm) == anorm) break;
					g=w[i];
					h=pythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=1;j<=m;j++) {
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) {
				if (z < 0.0) {
					w[k] = -z;
					for (j=1;j<=n;j++) v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30) nrerror("no convergence in 30 svdcmp iterations");
			x=w[l];
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=pythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++) {
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=pythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g = g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=1;jj<=n;jj++) {
					x=v[jj][j];
					z=v[jj][i];
					v[jj][j]=x*c+z*s;
					v[jj][i]=z*c-x*s;
				}
				z=pythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=1;jj<=m;jj++) {
					y=a[jj][j];
					z=a[jj][i];
					a[jj][j]=y*c+z*s;
					a[jj][i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	free_vector(rv1,1,n);
}
#undef SWAP

double polyfunc(double x, double p[], int np)
{
	int j;
	double y;
	
	y = p[1];
	for (j=2;j<=np;j++) y += p[j]*pow(x,j-1);
	
	return y;
}

void fpoly(double x, double p[], int np)
{
	int j;

	p[1]=1.0;
	for (j=2;j<=np;j++) p[j]=p[j-1]*x;
}



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

void lfit(double x[], double y[], int ndat, double a[], int ma, double *chisq)
{
	void covsrt(double **covar, int ma);
	void gaussj(double **a, int n, double **b, int m);
	void fpoly(double x, double p[], int np);
	int i,j,k,l,m;
	double ym,wt,sum,**beta,*afunc;
	double **covar;

	beta=matrix(1,ma,1,1);
	afunc=vector(1,ma);
	covar = matrix(1,ma,1,ma);
	
	for (j=1;j<=ma;j++) {
		for (k=1;k<=ma;k++) covar[j][k]=0.0;
		beta[j][1]=0.0;
	}
	for (i=1;i<=ndat;i++) {
		fpoly(x[i],afunc,ma);
		ym=y[i];

		for (j=0,l=1;l<=ma;l++) {
			wt=afunc[l];
			for (j++,k=0,m=1;m<=l;m++)
				covar[j][++k] += wt*afunc[m];
			beta[j][1] += ym*wt;
		}
	}
	for (j=2;j<=ma;j++)
		for (k=1;k<j;k++)
			covar[k][j]=covar[j][k];
	gaussj(covar,ma,beta,1);
	for (j=0,l=1;l<=ma;l++)
		a[l]=beta[++j][1];
	*chisq=0.0;
	for (i=1;i<=ndat;i++) {
		fpoly(x[i],afunc,ma);
		for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];
		*chisq += SQR((y[i]-sum));
	}
	covsrt(covar,ma);
	free_vector(afunc,1,ma);
	free_matrix(beta,1,ma,1,1);
}

#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
void covsrt(double **covar, int ma)
{
	int i,j,k;
	double swap;

	k=ma;
	for (j=ma;j>=1;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--;
	}
}
#undef SWAP
#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) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

void sort(unsigned long n, double arr[])
{
	unsigned long i,ir=n,j,k,l=1;
	int jstack=0,*istack;
	double a,temp;

	istack=ivector(1,NSTACK);
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
				}
				arr[i+1]=a;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP(arr[l],arr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j]);
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			jstack += 2;
			if (jstack > NSTACK) nrerror("NSTACK too small in sort.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free_ivector(istack,1,NSTACK);
}
#undef M
#undef NSTACK
#undef SWAP

double G_r(double sg)
{
    extern int nfunc;
	double logs;
	logs = log(sg) + lambda2; // To standardize the PDF, this is necessary.
	logs = logs*logs;
    nfunc++;
	return (SQR2PI)* exp(- logs / (2*lambda2))*P_L(extx,sg)/ (sg*sg*lambda); 
}

double Gau(double x)
{
	return (SQR2PI) * exp(- (x * x) / 2); 
}

double P_L(double s, double sgm)
{
	return (SQR2PI) * exp(- (s*s)/(2*sgm*sgm));
}

int pdf(long n, double x[], int bin, double y[],double py[])
{
	long i;
	int bn,y0;
	double max,min,dy,d0;
	
	min = x[1];
	max = x[1];	
	for(i=2;i<= n;i++){
		if(x[i]<min) min = x[i]; /* unit variance */
		if(x[i]>max) max = x[i]; /* unit variance */
	}
	dy = (max-min)/(double)(bin-1);
	d0 = max-floor(max/dy)*dy;

	if(d0 > dy/2){
		min = min - d0 + dy/2;
	}else{
		min = min - d0 - dy/2;
	}

	for(i = 1;i<=bin;i++){
		y[i] = min + dy * ((double)i - 0.5);
		py[i] = 0;
	}
	
	for(i = 1;i<=n;i++){
		bn = (int)((x[i] - min)/dy+0.5);
		py[bn]++;
	}

	for(i = 1;i<=bin;i++){
		py[i]=py[i]/(double)n/dy;
	}

	y0 = (int)((0 - min)/dy+0.5);

	return y0;
}