#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define BIG 10000000
double d[BIG],dnew[BIG];

int main(int argc, char *argv[]){
FILE *fin,*fout,*ftmp;
int i=0,P,c,len,k=0,n;
double kappa,err,dav;

// input data file is first argument (1 col assumed) 
// output/appended data file is second argument
// output file contains error estimated from 64 bins

fin=fopen(argv[1],"r");
fout=fopen(argv[2],"a");
ftmp=fopen("err_analyze","w");

(void)printf("kappa coupling ?\n");
(void)scanf("%lg",&kappa);

while(fscanf(fin,"%lg",&d[i])!=EOF){i++;}
n=i;
if(n>BIG){(void)printf("overflow -- redimension arrays\n");exit(1);}


// round data file to nearest power of 2
P=(int)(log((double)n)/log(2.0));
n=(int)(pow(2.0,(double)P));
(void)printf("Using %d measurements\n",n);
if(n<64){(void)printf("Insufficient data for error analysis\n");exit(1);}

len=n;

// average 
dav=0.0;
for(i=0;i<len;i++){
dav+=d[i];}
dav/=len;

while(len>64){

// compute binned measurements in new array from last array
c=0;
for(i=0;i<len;i=i+2){
dnew[c]=(d[i]+d[i+1])/2;
c++;
}


// compute error estimate
len=len/2;
err=0.0;
for(i=0;i<len;i++){
err+=(dnew[i]-dav)*(dnew[i]-dav);}

err=sqrt(err)/len;

// write out error estimates
(void)fprintf(ftmp,"%d\t%lg\n",k,err);

// finally write/add 32 bin number to file 
if(len==64){(void)fprintf(fout,"%lg %lg %lg\n",kappa,dav,err);}

for(i=0;i<len;i++){
d[i]=dnew[i];}

// bin level
k++;
}

return(0);
}