cc	lsfdf.f

	program lsfdf

	implicit real*8 (a-h,o-z)

cc	================================================

cc	Set parameters:
cc	~~~~~~~~~~~~~~~

cc	nsim : number of simulated light curves
cc	(in Monte Carlo simulation)
	parameter(nsim = 100)

cc	asig : threshold value for sigma-clipping
	parameter(asig = 3.0d0)
cc	================================================

cc	Explicit declarations:
cc	~~~~~~~~~~~~~~~~~~~~~~

	character * ( 80 ) arg(100)
	character * ( 80 ) filefr
	character * ( 80 ) filets
      real*8 fff
	integer iargc
	integer numarg

	real*8 w(100)
	real*8 t(10000)
	real*8 x(10000)
	real*8 xw(10000)
	real*8 amp(1000,100), phs(1000,100)
	real*8 aa(100), pp(100)

	real*8 a(100),b(100)

	real*8 ampavr(100)
	real*8 ampdsp(100)
	real*8 phsavr(100)
	real*8 phsdsp(100)
cc	================================================

cc	Set basic logical parameters:
cc	~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

cc	iwrite = 1  -->  Write std. dev., epoch, Fourier
cc			     parameters, etc. on screen
	iwrite = 1

cc	iwhite = 1  -->  Write out the prewhitened
cc			     time series into 'lsfd.prw'
	iwhite = 1

cc	ierror = 1  -->  Calculate the amplitude and
cc			     phase errors by Monte Carlo simulation
	ierror = 1

	p7 = 8.0d0*datan(1.0d0)

cc	================================================

cc >	Read command line arguments:
	numarg = iargc ( )
	do ii=1,numarg
		call getarg ( ii, arg(ii) )
	end do
	
cc	================================================

cc >	Display help:

	if(arg(1).eq.'-help') then
	write(*,*) ''
	write(*,*) '_________________________________________________'
	write(*,*) ''
	write(*,*) '                    LSFD'
	write(*,*) '         Least Squares Fourier Decomposer'
	write(*,*) '_________________________________________________'
	write(*,*) ''
	write(*,*) 'Usage:      lsfd <freqfile> <magfile>'
	write(*,*) ' '
	write(*,*) '<freqfile> contains the frequencies to be fitted:'
	write(*,*) '            freq(1)'
	write(*,*) '            freq(2)'
	write(*,*) '            ...'
	write(*,*) '            freq(M)'
	write(*,*) ''
	write(*,*) '<magfile> contains the time series in the form:'
	write(*,*) '            HJD(1)  MAG(1)'
	write(*,*) '            HJD(1)  MAG(2)'
	write(*,*) '            ...'
	write(*,*) '            HJD(N)  MAG(N)'
	write(*,*) '_________________________________________________'
	write(*,*) ''
	write(*,*) 'Decomposition has the following form:'
	write(*,*) ' a0 + a1 * sin( 2*pi * f1 * (t-E) + phi1 ) + ...'
	write(*,*) ''
	write(*,*) 'Parameters are written to lsfd.fit in the form:'
	write(*,*) '    < Norder >  < a0 >  < epoch >'
	write(*,*) '        < f1 >  < a1 >   < phi1 >'
	write(*,*) '        < f1 >  < a1 >   < phi1 >'
	write(*,*) '          ..      ..        ..   '
	write(*,*) '        < f1 >  < a1 >   < phi1 >'
	write(*,*) ''
	write(*,*) '  Note: amplitudes are in mmag.'
	write(*,*) '_________________________________________________'
	write(*,*) ''
	go to 9999
	else if(numarg.ne.2) then
	write(*,*) '\n Usage:      lsfd <freqfile> <magfile>'
	write(*,*) 'More info:  lsfd -help\n'
	go to 9999
	end if

cc	================================================
cc	              WORK BEGINS HERE
cc	================================================


cc	File containing the frequencies:
	filefr = arg(1)
cc	File containing the time series:
	filets = arg(2)

cc >	Read in frequencies:
	open(1, file = filefr)
	nrows = 0
	ii = 1
 300	continue
		read(1,*,end = 3000) w(ii)
		nrows = nrows + 1
		ii = ii + 1
	go to 300
 3000	continue
	close(1)
cc	norder - number of frequencies
	norder = nrows


cc >	Read in time series:
	open(2, file = filets)
	nrows = 0
	jj = 1
 301	continue
		read(2,*,end = 3010) t(jj), x(jj)
		nrows = nrows + 1
		jj = jj + 1
	go to 301
 3010	continue
	close(2)
cc	ndata - number of data points in the time series
	ndata = nrows

	avrg=0.0d0
	do ndtpnt=1,ndata
c		write(*,*) x(ndtpnt)
		avrg = avrg + x(ndtpnt)
	end do
	avrg = avrg / ndata
	write(*,*) avrg


cc	================================================

	call fourdec(norder,w,ndata,t,x,a0,a,b,xw,disp,t0)

cc	================================================


	if(ierror.eq.1) then
		write(*,*) '\n Enter an INTEGER NUMBER for random number generator'
		write(*,*) ''
		read(*,*) iseed
c		iseed = 2000*iseed+1

		call errors(norder,w,ndata,t,x,a0,a,b,disp,t0,nsim,
     &	iseed,amp,phs)
		

		do 280 kk=1,norder
		do 281 it=1,nsim
			aa(it) = amp(it,kk)
			pp(it) = phs(it,kk)
 281		continue
c		call stats(aa,nsim,aav,adi)
c		call stats(pp,nsim,pav,pdi)
		call omit1(asig,aa,nsim,aav,adi,nclip)
		call omit1(asig,pp,nsim,pav,pdi,nclip)

		ampavr(kk) = aav
		ampdsp(kk) = adi
		
		phsavr(kk) = pav
		phsdsp(kk) = pdi

		write(*,*) aav, adi, pav, pdi, nclip

 280		continue

	end if

cc	================================================

cc >	Write out Fourier decomposition:
	open(5,file = 'lsfd.out')
	write(5,175) norder, a0, t0
	do jj=1,norder
		write(5,176) w(jj), 1000.0d0*a(jj), b(jj)
	end do
 175	format(1x,i8,2x,f7.3,2x,f6.1)
 176	format(1x,f8.5,2x,f6.2,2x,f6.4)
	close(5)


cc >	Write derived parameters on screen:
	if(iwrite.eq.1) then
	write(*,*) '_________________________________________________'
	write(*,*) ''
	write(*,*) '                 LSFD RESULTS'
	write(*,*) '_________________________________________________'
	write(*,*) ''
	write(*,801) 'NFREQ = ', norder
	write(*,802) 'EPOCH = ', t0
	write(*,802) 'STDEV = ', disp
	write(*,802) 'A0    = ', a0 
	write(*,*) ''
	write(*,*) '    Frequency        Amplitude        Phase'
	write(*,*) '   -----------------------------------------'
	write(*,*) ''
	do kk=1,norder
	write(*,803) ' f(',kk,')=',w(kk),'   a(',kk,')=',
     &a(kk)*1000.0d0,'   p(',kk,')=',b(kk)
	end do
	write(*,*) '_________________________________________________'
	write(*,*) ''
 801	format(3x,a,i10)
 802	format(3x,a,f14.5)
 803	format(1x,a,i2,a,f9.6,a,i2,a,f14.5,a,i2,a,f6.4)
	if(ierror.eq.1) then
	write(*,*) ''
	write(*,*) '    Errors ( NSIM = ',nsim,' ):'
	write(*,*) '   -----------------------------------------'
	write(*,*) ''
	do kk=1,norder
	write(*,803) ' f(',kk,')=',w(kk),'  da(',kk,')=',
     &ampdsp(kk),'  dp(',kk,')=',phsdsp(kk)
	end do
	write(*,*) '_________________________________________________'
	write(*,*) ''
	end if
	end if
	

cc >	Write prewhitened time series into 'lsfd.prw'
	if(iwhite.eq.1) then
	open(6,file = 'lsfd.prw')
	do ii=1,ndata
		write(6,804) t(ii), xw(ii)
	end do
	close(6)
	end if

 804  format(1x,f12.6,2x,f10.3)

 9999	continue


	return
	end

   
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


	subroutine errors(norder,w,ndata,t,x,a0,a,b,disp,t0,nsim,
     &	iseed,amp,phs)

	implicit real*8 (a-h,o-z)

	real*8 w(100)
	real*8 t(10000)
	real*8 x(10000)
	real*8 xwsim(10000)
	real*8 xsim(10000)
	
      real*8 a(100),b(100)
	real*8 asim(100), bsim(100)
	real*8 amp(1000,100), phs(1000,100)

cc	p7 = 2 * pi
	p7 = 8.0d0*datan(1.0d0)

	do 250 it=1,nsim

		do 260 ii=1,ndata
			s = a0
			do 2001 kk=1,norder
				phase = w(kk) * ( t(ii) - t0 )
				s = s + a(kk) * dsin( p7*phase + b(kk) )
 2001			continue
			xsim(ii) = s + disp*gauss(iseed)
			write(20,*) t(ii), xsim(ii)
 260		continue
            write(20,*) '\n\n\n'

		call fourdec(norder,w,ndata,t,xsim,a0sim,asim,bsim,
     &	xwsim,dispsim,t0)

		do 270 kk=1,norder
			amp(it,kk) = 1000.0d0 * asim(kk)
			phs(it,kk) = bsim(kk)
c			write(15,*) bsim(kk)
 270		continue
		write(*,*) it, dispsim
 250	continue


	return
	end


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


      subroutine fourdec(norder,w,ndata,t,x,a0,a,b,xw,disp,t0)

cc	This routine fits a linear least squares
cc	Fourier-sum to a time series.

cc	INPUT ARGUMENTS:
cc	~~~~~~~~~~~~~~~~

cc	norder - number of frequencies to be fitted
cc	ndata - number of data in the input time series
cc	ww(i) , i = 1...norder - frequencies to be fitted
cc	t(i) , i = 1...ndata - time
cc	x(i) , i = 1...ndata - time series

cc	OUTPUT ARGUMENTS:
cc	~~~~~~~~~~~~~~~~~

cc	a0 - zero frequency component
cc	a(i) , i = 1...norder - amplitude of the i-th component
cc	b(i) , i = 1...norder - phase of the j-th component
cc	disp - standard deviation of the fit
cc	t0 - epoch


cc	FOURIER DECOMPOSITION HAS THE FOLLOWING FORM:
cc	~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

cc	a0 + a1 * sin( 2*pi*w1*(t-t0) + b1 ) + ...
cc	t(i) , x(i) , i=1...ndata - input data


	implicit real*8 (a-h,o-z)

	real*8 w(100)
	real*8 t(10000)
	real*8 te(10000)
	real*8 x(10000)

      real*8 g(100,101)
      real*8 xs(10000), xw(10000)
      real*8 a(100),b(100)

      p7 = 8.0d0*datan(1.0d0)

cc	t0 - epoch
cc	Set epoch to the integer part
cc	of the first time of the time series:
	t0 = dble(int(t(1)))
c	t0 = 0.0d0

	do ii=1,ndata
		te(ii) = t(ii) - t0
	end do

	m2 = 2 * norder
	m3 = m2 + 1
	m4 = m3 + 1

	call norm(ndata, norder, te, x, w, g)
	call matinv(norder, g)

	a0 = 0.0d+00
	do j=1,m3
		a0 = a0 + g(1,j) * g(j,m4)
	end do

	do 32 j=2,m3
		s = 0.0d+00
		do 33 k=1,m3
			s = s + g(j,k) * g(k,m4)
 33		continue
		a(j) = s
 32	continue

	k = 0
	do 37 j=2,m3,2

		k = k + 1
		j1 = j + 1
		amp = dsqrt( a(j)*a(j) + a(j1)*a(j1) )
		ph = a(j1) / amp
		ph = dacos(ph)

		if(a(j) .lt. 0.0d+00) then
			ph = p7 - ph
		end if

		b(k) = ph
		a(k) = amp

 37	continue


cc >	Calculate dispersion:
	disp = 0.0d0
	dev = 0.0d0

	do 2000 i=1,ndata
		s = a0
		time = te(i)

		do 2001 k=1,norder
			phase = w(k) * time
			s = s + a(k) * dsin( p7*phase + b(k) )
 2001		continue

       	xs(i) = s
       	devi = x(i) - xs(i)

      	xw(i) = devi

      	disp = disp + devi*devi
      	dd = xs(i) - a0
      	dev = dev + dd*dd

 2000	continue

	disp = dsqrt(disp / dfloat(ndata - 2*norder - 1))
	dev = dsqrt(dev / dfloat(ndata))

	return
	end


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


       SUBROUTINE MATINV(M,G)
C      MATRIX INVERTATION
       IMPLICIT REAL*8 (A-H,O-Z)
C      INPUT: G(M3,M3), OUTPUT: G(M3,M3)
       DIMENSION G(100,101)
       M3=2*M+1
       DO 11 I=1,M3
       W=G(I,I)
       Q=1.0D+00/W
       DO 12 K=1,M3
       G(I,K)=Q*G(I,K)
       IF(K.NE.I) GO TO 12
       G(I,K)=Q
  12   CONTINUE
       DO 13 J=1,M3
       IF(J.EQ.I) GO TO 13
       TT=G(J,I)
       Q1=-TT/W
       DO 14 K=1,M3
       G(J,K)=G(J,K)-TT*G(I,K)
       IF(K.NE.I) GO TO 14
       G(J,K)=Q1
  14   CONTINUE
  13   CONTINUE
  11   CONTINUE
       RETURN
       END


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


       SUBROUTINE NORM(N,M,T,X,F,G)
C      CALCULATION OF THE ELEMENTS OF THE NORMAL MATRIX AND R.H.S.
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION G(100,101),T(10000),X(10000),F(100)
       P7=8.0D0*DATAN(1.0D0)
       M3=2*M+1
       M4=M3+1
       L=0
       G(1,1)=FLOAT(N)
       DO 1 J=2,M3,2
       S=0.0D+00
       C=0.0D+00
       L=L+1
       J1=J+1
       FL=F(L)
       DO 2 I=1,N
       PH=T(I)*FL
       PH=P7*PH
       S=S+DSIN(PH)
  2    C=C+DCOS(PH)
       G(1,J)=C
  1    G(1,J1)=S
       L1=0
       DO 3 J=2,M3,2
       L1=L1+1
       J1=J+1
       L2=L1-1
       FL1=F(L1)
       DO 4 K=J,M3,2
       K1=K+1
       L2=L2+1
       FL2=F(L2)
       S1=0.0D+00
       S2=0.0D+00
       C1=0.0D+00
       C2=0.0D+00
       DO 5 I=1,N
       PH1=T(I)*FL1
       PH1=P7*PH1
       PH2=T(I)*FL2
       PH2=P7*PH2
       SPH1=DSIN(PH1)
       CPH1=DCOS(PH1)
       SPH2=DSIN(PH2)
       CPH2=DCOS(PH2)
       C1=C1+CPH1*CPH2
       S1=S1+CPH2*SPH1
       C2=C2+SPH2*CPH1
  5    S2=S2+SPH1*SPH2
       G(J,K)=C1
       G(J,K1)=C2
       G(J1,K)=S1
  4    G(J1,K1)=S2
  3    CONTINUE
       S=0.0D+00
       DO 6 I=1,N
  6    S=S+X(I)
       G(1,M4)=S
       L=0
       DO 7 J=2,M3,2
       S=0.0D+00
       C=0.0D+00
       L=L+1
       J1=J+1
       FL=F(L)
       DO 8 I=1,N
       PH=T(I)*FL
       PH=P7*PH
       S=S+X(I)*DSIN(PH)
  8    C=C+X(I)*DCOS(PH)
       G(J,M4)=C
  7    G(J1,M4)=S
       DO 9 J=1,M3
       DO 10 K=J,M3
  10   G(K,J)=G(J,K)
   9   CONTINUE
       RETURN
       END


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


      FUNCTION GAUSS(ISEED)
C
C     this routine generates Gaussian random variables of unit variance
C
      implicit real*8 (a-h,o-z)
C
 1    CONTINUE
      P1=2.0D0*RANDOMNUMBER(ISEED)-1.0D0
      P2=2.0D0*RANDOMNUMBER(ISEED)-1.0D0
      RR=P1*P1+P2*P2
      IF(RR.GE.1.0D0) GO TO 1
      FAC=DSQRT(-2.0D0*DLOG(RR)/RR)
      GAUSS=P1*FAC
C
      RETURN
      END


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


      FUNCTION RANDOMNUMBER(IX)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K
      DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/
      XHI = IX/B16
      XALO = (IX - XHI * B16) * A
      LEFTLO = XALO / B16
      FHI = XHI * A + LEFTLO
      K = FHI / B15
      IX = ((( XALO - LEFTLO * B16 ) - P ) + (FHI - K * B15 ) * B16) + K
      IF (IX.LT.0) IX = IX + P
      RANDOMNUMBER = DFLOAT(IX) * 4.65661287D-10
      RETURN
      END


cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
cc	xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


	subroutine stats(xx,nn,aver,disper)

cc	This routine computes the average and the dispersion of xx.

	implicit real*8 (a-h,o-z)

	real*8 xx(100)

	s = 0.0d0

	do i=1,nn
		s = s + xx(i)
	end do
	aver = s / dfloat(nn)

	di = 0.0d0
	do i=1,nn
		d = xx(i) - aver
		di=di+d*d
	end do

	disper = 0.0d0
	if(nn.gt.1) then
		disper = dsqrt( di/dfloat(nn-1) )
	end if

	return
	end


        subroutine omit1(asig,x,n,xave,xsig,nclip)
c
c---------------------------------------------------------------------
c
c       This routine computes average and standard deviation 
c       of  { x(i) ; i = 1,2, ..., n } by iteratively omitting 
c       data which deviate more that  asig*sigma  
c
c       Input parameters/arrays:  asig, t(i), x(i), n 
c       ~~~~~~~~~~~~~~~~~~~~~~~~ 
c
c       Output parameters:    xave = average
c       ~~~~~~~~~~~~~~~~~~    xsig = standard deviation
c                             nclip = number of data points clipped
c
c---------------------------------------------------------------------       
c
        implicit real*8 (a-h,o-z)
c
        dimension x(100),y(100)
        dimension in(100)
c
        eps0 = 1.0d-15
        nn   = n
c
        do 21 i=1,nn
        y(i) = x(i)
 21     continue
c
 20     continue
c
        call stats(y,nn,aver,disper)
c
c...    Finding out if there is an outlier
c
        if(disper.le.eps0) return
c
        iout=0
        devm=asig*disper
        do 10 i=1,nn
        dev=dabs(y(i)-aver)
        in(i)=1
        if(dev.gt.devm) in(i)=0
 10     continue
c
c...    Omit outlier
c
        k=0
        do 1 i=1,nn
        if(in(i).eq.0) go to 1
        k=k+1
        y(k)=y(i)
 1      continue
        if(k.eq.nn) go to 22
        nn=k
        go to 20
 22     continue
        nclip=n-nn
        call stats(y,nn,xave,xsig)
c
        return
        end
c
c
