program datdinf
c       data analysis of the time dependent green's functions
c       from the D=\infty-program

        integer nbin,l,nmax,lmax
c ie: band index
        integer il,id,ldummy, ie   
        parameter (nmax=150,lmax=301)
        real*8  p(2)
        real*8  G1(2,lmax), G1q(2,lmax), G2(2,lmax), G2q(2,lmax)
        real*8  G12(2,lmax), G12q(2,lmax)
        real*8  Gsum(lmax), Gsumq(lmax)
        real*8  GA(2,lmax), GAq(2,lmax), GB(2,lmax), GBq(2,lmax)
        real*8  dG,dGA,dGB
        real*8  Gtau, GAtau, GBtau, rbin
        real*8  G1M, G2M, G1qM, G2qM
        character*25  dummy25
        character*2  dummy2
        integer time(lmax)
c        real*8   time(lmax)


        iemax = 1

        read(5,'(a25,i5)') dummy25, Nbin
        read(5,'(a1,i4)') dummy1, l
        write(6,*) Nbin,l,iemax
        rbin=dfloat(Nbin)
        if (Nbin.eq.1) rbin=rbin+1.0D-6

c ---------------- set to zero ---------------------
        do ie=1,iemax
        do il=1,l+1
          G1(ie,il)  = 0.0
          G1q(ie,il) = 0.0
          G2(ie,il)  = 0.0
          G2q(ie,il) = 0.0
          G12(ie,il)  = 0.0
          G12q(ie,il) = 0.0
        enddo
        enddo
c --------------------------------------------------

c ------------------ read data ---------------------
c       Loop over the data-files (nbin)
cttttttttttttttttttt
c        if ((id.eq.1).or.(id.eq.2).or.(id.eq.3).or.
c     + (id.eq.6).or.(id.eq.8)) then
        do 10 id=1,nbin

c       Loop over the sublattices
          do ie=1,iemax

c       Loop over the imaginary time slices (l)
            do il=1,l+1
              read(5,*) time(il),G1tau,G2tau,dummy
              G1(ie,il)  = G1(ie,il)+G1tau
              G1q(ie,il) = G1q(ie,il)+G1tau**2
              G2(ie,il)  = G2(ie,il)+G2tau
              G2q(ie,il) = G2q(ie,il)+G2tau**2
              G12(ie,il)  = G12(ie,il)+(G1tau+G2tau)
              G12q(ie,il) = G12q(ie,il)+G1tau**2+G2tau**2
            enddo
            if (id.lt.nbin) read(5,'(a5)') dummy25
          enddo
 
 10     continue

cttttttttttttttttttt
c        endif 
c -----------------------------------------------

c ---------------- output results ---------------
c       Loop over the sublattices
       do ie=1,iemax
         OPEN (UNIT=9,FILE='g1.dat',STATUS='unknown')
         OPEN(UNIT=19,FILE='g2.dat',STATUS='unknown')
         OPEN(UNIT=49,FILE='g12.dat',STATUS='unknown')
         OPEN(UNIT=59,FILE='g1M.dat',STATUS='unknown')
         OPEN(UNIT=69,FILE='g2M.dat',STATUS='unknown')
         write (9,*) 0,0,0
         write (19,*)  0,0,0
         write (49,*)  0,0,0
         if (ie.eq.2) then
         write (59,*)  0,0,0
         write (69,*)  0,0,0
         endif

ctttttttt           
         do il=1,l+1,1
ctttttttt
         dG1  = sqrt((G1q(ie,il)/rbin-(G1(ie,il)/rbin)**2)/(rbin-1.d0))
         dG2  = sqrt((G2q(ie,il)/rbin-(G2(ie,il)/rbin)**2)/(rbin-1.d0))
         dG12=sqrt((G12q(ie,il)/rbin/2.-(G12(ie,il)/rbin/2.)**2)
     c                                 /(2.*rbin-1.d0))
c         dG  = sqrt(( Gq(1,il)/rbin-( G(1,il)/rbin)**2)/(rbin-1.d0))

         write (9,1000) time(il),G1(ie,il)/rbin,dG1
         write (19,1000) time(il),G2(ie,il)/rbin,dG2
         write (49,1000) time(il),G12(ie,il)/rbin/2.,dG12
         if (ie.eq.2) then
            G1qM=.5D0*(G1q(1,il)+G1q(2,il))
            G2qM=.5D0*(G2q(1,il)+G2q(2,il))
            G1M=.5D0*(G1(1,il)+G1(2,il))
            G2M=.5D0*(G2(1,il)+G2(2,il))
            dG1M  = sqrt((G1qM/rbin-(G1M/rbin)**2)/(rbin-1.d0))
            dG2M  = sqrt((G2qM/rbin-(G2M/rbin)**2)/(rbin-1.d0))           
            write (59,1000) time(il),G1M/rbin,dG1M
            write (69,1000) time(il),G2M/rbin,dG2M
         endif
         enddo

       enddo

       OPEN(UNIT=59,FILE='gsum.dat',STATUS='unknown')
ctttttttt
c       do il=1,l+1
c         Gsum(il) = 0.d0
c         Gsumq(il) = 0.d0
c         do ie=1,iemax
c           Gsum(il) = Gsum(il)+p(ie)*G(ie,il)
c           dG = (Gq(ie,il)/rbin-( G(ie,il)/rbin)**2)/(rbin-1.d0)
c           Gsumq(il) = Gsumq(il)+p(ie)*dG
c         enddo
 

c          write(59,1000) time(il),Gsum(il)/rbin,sqrt(Gsumq(il))c
c
c       enddo

 1000  format(1x,i6,4f12.6)
c 1000  format(1x,5f12.6)

       write(6,*)
       write(6,*) '----------- ok ----------'

       stop
       end