program normperiod

c	*******************************************************************
c	THE PROGRAMME CALCULATES THE LOMB NORMALIZED PERIODOGRAM
c     -------------------------------------------------------------------
c     Input:
c	N: time series length
c	Tf: total time interval
c	fratio: oversampling parameter
c	file1: data file name
c	file2: time input file name
c     file3: output file name
c     conlev: confidence level
c
c	Output:
c     cconlev: value of the Lomb Periodogram at the specified conf. level
c	f: frequency
c     P: Lomb normalized periodogram value
c     -------------------------------------------------------------------
c	written by Stella Pytharouli on 12.10.2003, v2.1 on 12.01.2010
c	********************************************************************
	    
	implicit doubleprecision(a-h,o-z)
      dimension h(24600), t(24600)
      character*20 file1,file2,file3
		
      write(*,*) 'length of time series N'
	read(*,*) N
	
      write(*,*)'total time interval T:'
	read(*,*) Tf
	
	f0=1.d0/(Tf)
	fc=N/(2.d0*Tf)

	write(*,*)'fratio = ; gia fc= ', fc
	read(*,*) fratio

	
	fhi=fratio*fc
	Np=(fratio*4*N)/2
	fstep=(fhi-f0)/Np

      write(*,*)'Data file name (max 20 characters):'
	read(*,*) file1
	write(*,*)'Time input file name (max 20 characters):'
	read(*,*) file2
	write(*,*)'Output file name (max 20 characters):' 
	read(*,*) file3

c	 confidence level calculation

	write(*,*) 'confidence level (%):'
	read(*,*) conlev

	plevel=1-(conlev/100.d0)
	subroot=(1-plevel)**(1.d0/N)
	domln=1-subroot
	cconlev=-log(domln)

	write(*,*)
	write(*,*)'Lomb normalized periodogram amplitude'
	write(*,*)'for',conlev,'% confidence level:',cconlev 
	write(*,*)
	write(*,*)
	write(*,*)' processing....' 
	write(*,*)

	open(5,file=file1)
	open(15,file=file2)
	open(25,file=file3)

      sumx=0.d0
	do 10 i=1,N
	 read(5,*) h(i)
	 read(15,*) t(i)
	 sumx=sumx+h(i)
  10  continue	  
      hm=sumx/N

	sumsh=0.d0
	do 20 i=1,N
	 sumsh=sumsh+(h(i)-hm)**2.d0
  20  continue
      s2=sumsh/(N-1)
	
	p3=1.d0/(2.d0*s2)

	pi=3.1415926535897932384626433832795d0
	
	f=0.d0
	do while (f.le.fhi)
	 f=f+fstep
	 omega=2.d0*pi*f
	 sum1=0.d0
	 sum2=0.d0
	 do 30 i=1,N
	  sum1=sum1+sin(2.d0*omega*t(i))
	  sum2=sum2+cos(2.d0*omega*t(i))
 30    continue
       taf1=atan(sum1/sum2)
	 taf=(1.d0/(2*omega))*taf1
	 sum3=0.d0
	 sum4=0.d0
	 sum5=0.d0
	 sum6=0.d0
	 do 40 i=1,N
	  dh=h(i)-hm
	  dt=t(i)-taf
	  c1=cos(omega*dt) 
	  c2=sin(omega*dt)
        sum3=sum3+dh*c1
	  sum4=sum4+c1**2.d0
	  sum5=sum5+dh*c2
	  sum6=sum6+c2**2.d0
  40   continue
       p1a=sum3**2.d0
	 p1=p1a/sum4
	 p2a=sum5**2.d0
	 p2=p2a/sum6
	 P=p3*(p1+p2)
	 write(25,*) f,P
      end do
	

	close(5)
	close(15)
	close(25)

	stop
	end