/* RRIを等時間間隔にリサンプリング  */
/* 16/2/2009 Ken KIYONO */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nr.h"
#include "nrutil.h"
#include "nrutil.c"


void help(void);
void splint(double xa[], double ya[], double y2a[], int n, double x, double *y);
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]);
int read_data(char filename[], double x[]);

/* Global variables. */
char *pname;	/* this program's name (for use in error messages) */
double *data;

main(int argc, char **argv)
{
    long nn,i;
	long hrlength,ttmp;
    double y, yp = 0.0, ysum = 0.0;
	double *hr,tt;
	double rint;
	double unit;
	double totalt;
	double intt;
	double slope;
	double yp1,ypn,*xs,*ys,*ys2,ysp;
	int linear;
	int j,k;
	int outc;
	char fname1[256];

	/* Read and interpret the command line. */
    pname = argv[0];
    rint = 0.25;
	linear = 0;
	outc = 0;
    for (i = 1; i < argc && *argv[i] == '-'; i++) {
      switch(argv[i][1]) {
  	case 'd':
		linear = 0;
		break;
  	case 's':
		linear = 2;
		break;
  	case 'l':
		linear = 1;
		break;
  	case 'i':
		outc = 1;
		break;
  	case 't':
		rint = atof(argv[++i]);
		if(rint < 0.01) rint = 0.01;
		break;
	case 'h':	/* print usage information and quit */
        default:
	  help();
	  exit(1);
      }
    }
	
	if(argc < 2){
	  help();
	  exit(1);
	}
	strcpy(fname1,argv[argc-1]);
	nn = read_data(fname1, data);
	
	ysum = 0;
	for(i=1;i<=nn;i++){
		ysum += data[i];
	}

	if(data[1] > 100) unit = 1000; /* 単位　msec　と　sec の判別 */
	else unit = 1;

	totalt = (ysum-data[nn])/unit;

	hrlength = (long)((totalt)/rint);
	hr = vector(1,hrlength);
	ysum = ysum/(double)(nn);

	if(unit > 1){
		for(j=1;j<=nn;j++){
			data[j] = data[j]/unit;
		}
	}

	fprintf(stderr,"number of input data: %d\n", nn);
	fprintf(stderr,"total time: %lf\n", totalt);
	fprintf(stderr,"resampling interval (frequency): %lf (%lf)\n", rint, 1/rint);

	ttmp = 1;
	tt = data[ttmp];

	if(linear < 1){ /* ステップ補間 */
		fprintf(stderr, "Resamping as a step function\n");	
		for(k = 1;k <= hrlength; k++){
			intt = rint*(double)(k-1);
//			if(tt < intt){
			while(tt < intt){
				ttmp++;
				tt = tt + data[ttmp];
			}
			hr[k] = data[ttmp];
		}
	}else if(linear == 1){ /* 線形補完 */
		fprintf(stderr, "Resamping with linear interpolation.\n");	
		for(k = 1;k <= hrlength; k++){
			intt = rint*(double)(k-1);

			while(tt < intt){
				ttmp++;
				tt = tt + data[ttmp];
			}
			slope = (data[ttmp+1]-data[ttmp])/data[ttmp];
			hr[k] = data[ttmp+1]+slope*(intt - tt);	
		}

	}else{
		
		fprintf(stderr, "Resamping with cubic spline interpolation.\n");	
			
//		ttmp = 3;
		xs=vector(1,nn);
		ys=vector(1,nn);
		ys2=vector(1,nn);

		tt = 0;
		for(i = 1;i <= nn;i++){
			xs[i] = tt;
			ys[i] = data[i];
			tt += data[i];
		}
		
		yp1 = 1e30; /* 自然スプラインの境界条件 */
		ypn = yp1; /* 自然スプラインの境界条件 */
		spline(xs,ys,nn,yp1,ypn,ys2);
		
		for(k = 1;k<= hrlength; k++){
			tt = rint*(double)(k-1);
			splint(xs,ys,ys2,nn,tt,&ysp);
			hr[k] = ysp;
		}
		free_vector(ys2,1,nn);
		free_vector(ys,1,nn);
		free_vector(xs,1,nn);

	}


	for(k = 1;k<= hrlength; k++){
		if(outc < 1) fprintf(stdout,"%lf\t%lf\n",rint*(double)(k-1),1000*hr[k]);
		else fprintf(stdout,"%lf\n",1000*hr[k]);
	}
    /* Output the results. */

    /* Release allocated memory. */

	free_vector(hr,1,hrlength);
    exit(0);
}

int read_data(char filename[], double x[]){
	/* データの入力 */
    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:",
 " -d               ステップ補間",
 " -s               スプライン補間",
 " -l               線形補間",
 " -i               y列のみの出力",
 " -t T             リサンプリング間隔をTに設定",
 " -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 NRANSI
#include "nrutil.h"

void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
	int i,k;
	double p,qn,sig,un,*u;

	u=vector(1,n-1);
	if (yp1 > 0.99e30)
		y2[1]=u[1]=0.0;
	else {
		y2[1] = -0.5;
		u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
	}
	for (i=2;i<=n-1;i++) {
		sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
		p=sig*y2[i-1]+2.0;
		y2[i]=(sig-1.0)/p;
		u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
		u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
	}
	if (ypn > 0.99e30)
		qn=un=0.0;
	else {
		qn=0.5;
		un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
	}
	y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
	for (k=n-1;k>=1;k--)
		y2[k]=y2[k]*y2[k+1]+u[k];
	free_vector(u,1,n-1);
}
#undef NRANSI

void splint(double xa[], double ya[], double y2a[], int n, double x, double *y)
{
	void nrerror(char error_text[]);
	int klo,khi,k;
	double h,b,a;

	klo=1;
	khi=n;
	while (khi-klo > 1) {
		k=(khi+klo) >> 1;
		if (xa[k] > x) khi=k;
		else klo=k;
	}
	h=xa[khi]-xa[klo];
	if (h == 0.0) nrerror("Bad xa input to routine splint");
	a=(xa[khi]-x)/h;
	b=(x-xa[klo])/h;
	*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}