/* cosine tapers are applied; n2,dt2 returned;
 */
 
#include <math.h>
#include "const.h"
#include "message.h"

float cubinterpe();

float *bandpass(data1,n1,dt1,w1,w2,w3,w4,n2,dt2)
float *data1,dt1,*dt2; /* if *dt2!=0.0 data will be resampled */
float w1;          /* cut-off omega for hi-pass */
float w2;          /* corner  omega for hi-pass */
float w3;          /* corner  omega for low-pass */
float w4;          /* cut-off omega for low-pass */
int n1,*n2;
{
int i,m;
float w,wnq,length;
static int currentlen=0;
static float *data2;
float *work;
float x,temp;
double fac;
int N;
	
	STOP(bandpass: use newer version -- bandpass 1);
		
	length=dt1*(float)(n1-1);
	wnq=PI/dt1;
	if(w3>wnq) w3=w4=wnq;   /* high pass */
	if(w4>wnq) STOP(bandpass: error 1);
	if(w1>w2) STOP(bandpass: error 2);
	if(w2>w3) STOP(bandpass: error 3);
	if(w3>w4) STOP(bandpass: error 4);
	
	m=min2n(n1);
	
	N=1;
	if(*dt2!=0.0)
	{ temp=wnq;	/* make sure at lease 8 samples per cycle */
	  while(temp<4.*w4)
	  { temp*=2.0;
	    N*=2;
	  }
	}
		
	work=(float *)malloc(sizeof(float)*(N*m+2));
	
	for(i=0;i<n1;i++) work[i]=data1[i];
	for(i=n1;i<m;i++) work[i]=0.0;
	
        cfftr(work,m);
			
	for(i=0;i<1+m/2;i++)
	{ w=(F)(i)*wnq/(F)(m/2);
	  if(w>=w4 || w<= w1) 
	  { work[2*i]=work[2*i+1]=0.0;
	  }
	  else if(w<w4 && w>w3)
	  { fac=PI*(w-w3)/(w4-w3);
	    fac=0.5*(1.+cos(fac));
	    work[2*i]*=fac;
	    work[2*i+1]*=fac;
	  }
	  else if(w>w1 && w<w2)
	  { fac=PI*(w2-w)/(w2-w1);
	    fac=0.5*(1.+cos(fac));
	    work[2*i]*=fac;
	    work[2*i+1]*=fac;
	  }
	}
	
	for(i=m+2;i<N*m+2;i++) work[i]=0.0;
	
	cfftri(work,N*m);
	
	for(i=0;i<N*m;i++) work[i]*=2./(F)m;
	
	if(*dt2!=0.0)
	{ *dt2=PI/w4;
	  if(*dt2>=10.&&*dt2<100.) *dt2=(int)(*dt2);/*round-off to 2-d*/
	  else if(*dt2>=1.0 && *dt2<10.) *dt2=0.1*(int)(*dt2*10.);
	  else STOP(bandpass: to be fixed);
	
	  *n2=1 + (int)(length/(*dt2));
	  
	  if(*n2>currentlen)
	  { if(currentlen) free(data2);
	    currentlen=*n2;
	    data2=(float *)malloc(sizeof(float)*(*n2));
	  }
		
/* interpolating */
	
	  data2[0]=cubinterpe(work,N*m,dt1/N,0.0,1);
	  for(i=1;i<*n2;i++)
	  { x=*dt2*i;
	    data2[i]=cubinterpe(work,N*m,dt1/N,x,0);
	  }
	}
	else 
	{ *n2=n1;
	  *dt2=dt1;
	  
	  if(*n2>currentlen)
	  { if(currentlen) free(data2);
	    currentlen=*n2;
	    data2=(float *)malloc(sizeof(float)*(*n2));
	  }
	  
	  for(i=0;i<n1;i++) data2[i]=work[i];
	}
		
	free(work);
	
	return(data2);
}