	program datfit
c.. Purpose:
c    to fit data time series by 6 time functions
c    time window is allowed to process selected data within the window
c    post-fit unit weight sigma is also calculated
c
c
c.. Author:
c    Xiaoping Frank Wu
c.. Date:
c    06/23/2004 at JPL
c...parameters:
c    parameter input
c      maxy ---- maximum number of data for each series
c      maxx ---- number of fit parameters 
c      tbeg ---- beginging time for fitting
c      tend ---- ending time for fitting

	implicit none
      character*120 rcsID
      data rcsID
     &/'$Id: datfit.f,v 1.1 2004/09/28 21:36:08 xpw Exp xpw $'/
        include 'stdunits.h'

        integer i,j,k,l,m,n,ldim,ndata
        integer imon,itmon,iyr,ifile
        integer maxx,maxy,npar
        parameter (maxx=6,maxy=5660,l=0,ldim=1)
        double precision pi
        double precision t,omega1,omega2,tbeg,tend,epoch
        double precision hx(maxy,maxx),y(maxy),dcov(maxy)
        double precision dcom(maxx,maxx),sol(maxx)
        double precision hc(maxy,ldim),dc(ldim,ldim),pxc(maxx,ldim)
        double precision sxc(maxx,ldim),dcon(maxx,maxx)
        double precision hxt(maxx,maxy),p(maxy),b1(maxy+1,maxx+1)
        double precision b2(maxx,ldim),b3(ldim,ldim)
	double precision xap(maxx),Cap(maxx,maxx)
	double precision data, sigma, rms
	double precision ym(maxy),time(maxy)
        Namelist /input/npar,tbeg,tend

	pi=4.d0*datan(1.d0)
	epoch=2004.d0
        Open ( Unit = 1, Status = 'OLD',
     &         FILE= 'datfit.nml')
        read(unit=1, nml=input)
	close (1)

        omega1=2.d0*pi/365.25d0
        omega2=2.d0*pi/365.25d0*2.d0

	k=0
	do i=1,maxy
	read(5,*,end=100)t,data,sigma
        if(t.ge.tbeg .and. t.le.tend)then
	k=k+1
	time(k)=t
	t=(t-epoch)*365.25d0
	y(k)=data
	dcov(k)=sigma*sigma
	if(dcov(k).eq.0.d0)dcov(k)=1.d0
        hx(k,1)=1.d0
        hx(k,2)=t/365.25d0
        hx(k,3)=cos(omega1*t)
        hx(k,4)=sin(omega1*t)
        hx(k,5)=cos(omega2*t)
        hx(k,6)=sin(omega2*t)
	end if
	end do

 100	continue
	ndata=k
        do n=1,maxx
        do m=1,maxx
        dcom(n,m)=0.d0
        end do
        end do

        call perconuw2(ndata,npar,l,maxy,hx,dcov,y,hc,ldim,dc,maxx,
     &                 pxc,sxc,dcom,sol,dcon,hxt,p,b1,b2,b3)

	call amppha(npar/2,maxx,sol,dcom,xap,Cap)
	rms=0.d0
        open(2,file='sdensigma',access='sequential',
     &       form='formatted')
	do i=1,ndata
	ym(i)=0.d0
	do j=1,npar
	ym(i)=ym(i)+hx(i,j)*sol(j)
	end do ! j
	rms=rms+(y(i)-ym(i))*(y(i)-ym(i))/dcov(i)
	end do !  i
	rms=sqrt(rms/dble(ndata-npar))
	write(stdout,*) ' ndata, posteriori unit weight sigma = ',
     &                  ndata, rms
	write(stdout,*)' bias a = ',sol(1),sqrt(dcom(1,1))*rms
	write(stdout,*)' rate b = ',sol(2),sqrt(dcom(2,2))*rms
	write(stdout,*)' annual in = ',sol(3),sqrt(dcom(3,3))*rms
	write(stdout,*)' annual out = ',sol(4),sqrt(dcom(4,4))*rms
	write(stdout,*)' annual amp = ',xap(3),sqrt(Cap(3,3))*rms
        if(xap(4).lt.0.d0)xap(4)=xap(4)+2.d0*pi
	write(stdout,*)' annual max time = ',xap(4)/omega1,
     &                 sqrt(Cap(4,4))/omega1*rms
	write(stdout,*)' semi-annual in = ',sol(5),
     &                   sqrt(dcom(5,5))*rms
	write(stdout,*)' semi-annual out = ',sol(6),
     &                   sqrt(dcom(6,6))*rms
	write(stdout,*)' semi-annual amp = ',xap(5),
     &                   sqrt(Cap(5,5))*rms
        if(xap(6).lt.0.d0)xap(6)=xap(6)+2.d0*pi
	write(stdout,*)' semi-annual max time = ',xap(6)/omega2,
     &                 sqrt(Cap(6,6))/omega2*rms
	do i=1,ndata
        write(2,'(f20.10, 4e15.6)')time(i), y(i),ym(i), y(i)-ym(i),
     &                             sqrt(dcov(i))
	end do ! i
	end

