c ###########################################################################
c Filtre passe Bande de Panteleev 
c
c Christian Bizouard - Novembrer 2014
c ###########################################################################
 
	implicit none
	
	integer m,l
	PARAMETER (m=100000,l=20000)
	real*8 TAB(80),T(m),X(m),Y(m),P(m),FREQ(l),FREQ1(l)	
	integer i,j,n,nfreq	
	real*8 a,sig_a,b,sig_b,cor_ab,pas,freq0,freqN,pi,amp,phase,f0,fc,tr,sig0,arg,t0,phi1,ti,tj
	COMPLEX*16 F(m),CI,DFT(l),DFT1(l),DFTI(m),filter,CONV(m)
	integer ncom,cold,col1,col2, opt
	CHARACTER*140 fichin, c_ncom, c_cold, c_col1, c_col2,com,c_fc, c_f0, c_t0 

	pi = 4.d0 * atan(1.d0)
	
	CI=(0.d0,1.d0)
	
      	call getarg(1,fichin     )  ! Entrer le nom fichier
        
	if (Len_trim(fichin).eq.0) then
	
	write(6,*) '<<     Filtrage passe bande Panteleev    >>'
	write(6,*) ' Usage : panteleva [filein][ncom][cold][col1][col2][fc][f0][t0]'
        write(6,*) ' [filein] ..... nom du fichier'
        write(6,*) ' [ncom] ....... nombre de lignes de commentaire       '
        write(6,*) ' [cold] ....... colonne des  dates                    '
        write(6,*) ' [col1] ....... 1ere colonne des valeurs              '
        write(6,*) ' [col2] ....... 2eme colonne des valeurs (0 si aucune)'
        write(6,*) ' [fc] ......... frequence centrale'		
        write(6,*) ' [f0] ......... Largeur fenetre f0'
	write(6,*) ' [t0] ......... origine temporelle de la phase de fc'	
	write(6,*) ' panteleva.res : t, amp, phase change, x , y'  
	stop     
	
	endif
	
	call getarg(2 ,c_ncom    )  ! Nombre de lignes de commentaires
        call getarg(3 ,c_cold    )  ! Numero de la colonne des dates
        call getarg(4 ,c_col1    )  ! Numero de la 1ere colonne des valeurs
        call getarg(5 ,c_col2    )  ! Numero de la 2eme colonne des valeurs
        call getarg(6 ,c_fc      )   
        call getarg(7 ,c_f0      )  
        call getarg(8 ,c_t0      )  
     
	read(c_ncom, *) ncom
	read(c_cold, *) cold
	read(c_col1, *) col1
	read(c_col2, *) col2
	read(c_fc, *)   fc
	read(c_f0, *)   f0
	read(c_t0, *)   t0
     	! write(6,*) 'largeur de bande =',1/f0

        open(10,file=fichin,access='sequential',status='old')

      	do i=1,ncom
        read(10,1000) com
      	enddo
1000   	format(a140)
 
       i = 0
 
10      read(10,*,end=20) (TAB(j),j=1, max(col1,col2,cold) )
        i = i + 1
        T(i) = TAB(cold) - t0
        X(i) = TAB(col1)	        
        if (col2.ne.0) then 
        Y(i) = TAB(col2) 
        else
        Y(i) = 0.d0 
        endif
       
        P(i) = 1.d0  ! poids des observations
        goto 10
        
 20     continue
        close(10)
        n = i



	opt=1
	
	
	! --------------------------
	!  Remove the mean and trend
	! --------------------------
       	!call DTREND(n,T,X,P,a,sig_a,b,sig_b,cor_ab,sig0)        
       	!call DTREND(n,T,Y,P,a,sig_a,b,sig_b,cor_ab,sig0)        

	! ----------------------------
	! formation du signal complexe
	! ----------------------------

        do i =1,n
        	if (col2.eq.0) F(i) = X(i)	        
        	if (col2.ne.0) F(i) = X(i) + CI * Y(i)
        enddo
	
	! #################
	! filtrage spectral
	! #################
	
	if (opt.eq.0) then 	

	! pas equidistant     
        pas = T(2) - T(1)   
         
	!frequence fondamentale      
        freq0 = 1.d0/(dble(n-1)*pas)
        !frequence de Nyquist
        freqN = 1.d0/(2.d0*pas)

	!boucle sur les frequences        
	i = 0
	FREQ(i) = 0.d0

	do while (FREQ(i).le.freqN)  
		i = i + 1
        		FREQ(i) = freq0 * dble(i-1) 		 		        
			DFT(i) = 0.d0 
		do j = 1,n
			DFT(i) = DFT(i) + F(j) * cdexp(-CI*2.d0*pi*FREQ(i)*T(j))
		enddo		 
	enddo
	NFREQ = i ! dimension du tableau des frequences
 
	if (col2.ne.0) then ! 2D
		do while(FREQ(i).ge.-freqN)  
			i = i + 1
        		FREQ(i) =-freq0 * dble(i-nfreq) 		 		        
			DFT(i) = 0.d0 
			do j = 1,n
				DFT(i) = DFT(i) + F(j) * cdexp(-CI*2.d0*pi*FREQ(i)*T(j))
			enddo		 
		enddo

		! tri par frequences croissantes  		
		do i = 1,nfreq-1
			FREQ1(i) = FREQ(2*nfreq-i)  
			DFT1(i)  =  DFT(2*nfreq-i)  
        	enddo
	
		do i = nfreq,2*nfreq-1
			FREQ1(i) = FREQ(i+1-nfreq) 
			DFT1(i)  = DFT(i+1-nfreq) 
       	 enddo
		
		NFREQ = nfreq + nfreq-1 

		do i = 1,nfreq
			FREQ(i) = FREQ1(i) 
			DFT(i)  = DFT1(i)   
		enddo

	endif

	!#####################
	!filtrage de Pentaleev
	!#####################
	do i=1,nfreq  		
		TR =  f0**4 / ( f0**4 + (FREQ(i) - fc)**4 ) 
		DFT(i) = DFT(i) * TR   
	enddo

	!#################################
	!Transformation de Fourier inverse
	!#################################
	DFTI(j) = 0 
	do j = 1, n 				 		        
		DFTI(j) = 0 
		do i = 1, nfreq 
			DFTI(j) = DFTI(j) + DFT(i) * cdexp( CI * 2.d0 * pi * FREQ(i) * T(j) )
		enddo
		DFTI(j) = DFTI(j)  / dble(n) 	
	enddo

	! ---------------------------------------------		
	! amplitude spectrale et ecriture des resultats 
	! ---------------------------------------------		

        open(11,file='panteleva.res',status='unknown')

	do i=1,n  
		if (  dreal(DFTI(i)).ne.0 ) arg= atan2 (dimag(DFTI(i)) , dreal(DFTI(i)) )  - dmod(2* pi * fc * T(i) ,2*pi)
		if ( (dreal(DFTI(i)).eq.0 ).and.(dimag(DFTI(i)).gt.0) ) arg = pi/2 - dmod(2* pi * fc * T(i) ,2*pi)
		if ( (dreal(DFTI(i)).eq.0 ).and.(dimag(DFTI(i)).lt.0) ) arg =-pi/2 - dmod(2* pi * fc * T(i) ,2*pi)
		arg = dmod(arg, 2 *pi) * 180/pi  	
		if (i.ge.2) then
			if ( arg - phi1.ge. 330.d0 ) arg = arg - 360.d0 
			if ( arg - phi1.le.-330.d0 ) arg = arg + 360.d0 
		endif
		phi1= arg
		write(11,1030) T(i)+t0,abs(DFTI(i)), arg, dreal(DFTI(i)), dimag(DFTI(i))
	enddo 
	
1030	format(5e20.8)
	close(11)
	
	else 
	
	! #################################################
	! produit convolution par la reponse impulsionnelle
	! #################################################
	
	
	do i = 1, n 
		ti = T(i)
		do j=1,n
		tj = T(j)
		CONV(i) = CONV(i) + F(j) * filter(ti -tj,f0,fc)     
		enddo
		CONV(i) = CONV(i) * (T(2)-T(1))  ! multiplication par le pas temporel				
		if (  dreal(CONV(i)).ne.0 ) arg= atan2 (dimag(CONV(i)) , dreal(CONV(i)) )  - dmod(2* pi * fc * T(i) ,2*pi)
		if ( (dreal(CONV(i)).eq.0 ).and.(dimag(CONV(i)).gt.0) ) arg = pi/2 - dmod(2* pi * fc * T(i) ,2*pi)
		if ( (dreal(CONV(i)).eq.0 ).and.(dimag(CONV(i)).lt.0) ) arg =-pi/2 - dmod(2* pi * fc * T(i) ,2*pi)
		arg = dmod(arg, 2 *pi) * 180/pi 
		
		if (i.ge.2) then
			if ( arg - phi1.ge. 330.d0 ) arg = arg - 360.d0 
			if ( arg - phi1.le.-330.d0 ) arg = arg + 360.d0 
		endif
		phi1= arg
		write(6,1030) T(i)+t0,abs(CONV(i)), arg , dreal(CONV(i)), dimag(CONV(i))
	enddo 	
	
	
	endif 
	
	end 
		
	include 'utils.f'

	
	! ######################
	! réponse impulsionnelle
	! ######################
	function filter(t,f0,fc)
	double precision t,f0,w0,fc
	complex*16 filter,CI
	pi = 4.d0 * atan(1.d0)	
	CI=(0.d0,1.d0)
	w0 = 2 * pi * f0
	filter = w0 / 2.d0 /dsqrt(2.d0) * cdexp(- (w0/dsqrt(2.d0) * dabs(t) - CI * 2*pi *fc * t) )
     & * ( dcos(w0/dsqrt(2.d0) *t) + dsin(w0/dsqrt(2.d0) *dabs(t) ) )
	end 
