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

float cubinterpe();

float *bandpass1(data1,n1,dt1,w1,w2,w3,w4,n2,dt2,flag)
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 flag; /* 0: time -> time;
           * 1: time -> freq;
           * 2: freq -> time;
           * 3: freq -> freq; */
{
int i,m;
float w,wnq,length;
static int currentlen=0;
static float *data2;
float *work;
float x,temp;
double fac;
int N;
        if(flag<0 || flag>3) STOP(bandpass1: error -- illegal flag);
        if(flag==1 || flag==3) STOP(bandpass1: error -- to be fixed);
	
	wnq=PI/dt1;
	if(w3>wnq) w3=w4=wnq;   /* high pass */
	if(w4>wnq) STOP(bandpass1: error 1);
	if(w1>w2) STOP(bandpass1: error 2);
	if(w2>w3) STOP(bandpass1: error 3);
	if(w3>w4) STOP(bandpass1: 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));
	
	if(flag<2)
	{ for(i=0;i<n1;i++) work[i]=data1[i];
	  for(i=n1;i<m;i++) work[i]=0.0;
          cfftr(work,m);
        }
        else if(flag>1)
        { for(i=0;i<=m+1;i++) work[i]=data1[i];
        }
			
	length=dt1*(float)(n1-1);
	
	for(i=0;i<1+m/2;i++)
	{ complex aaa,Butterworth();
	
	  w=(F)(i)*wnq/(F)(m/2);
	  aaa=Butterworth(w,w2,0.15);
	  
	  temp=work[2*i];
	  work[2*i]=aaa.real*work[2*i]-aaa.imag*work[2*i+1];
	  work[2*i+1]=aaa.real*work[2*i+1]+aaa.imag*temp;
	  
	  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(bandpass1: 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);
}