C ---- subroutine for initializing the start GF ------------------------
      
      subroutine init_GF(Temp, gtin,gfin,gafin,t,omega,z,domega,
     *             mu,dtime,inguess,u,l)

C set parameters
      real*8 pi
      parameter(pi=3.141592653589793d0)
      
      integer n,mi
      parameter(n=2048*4)
C     (n must be a power of 2.)

      double complex zi
      parameter(zi=(0.d0, 1.d0))

      integer i,inguess,lat,l
      double complex gtin(2,2,n), gfin(2,2,n), gafin(n), Aint(n), 
     *            Atint(n),Ainf(n), Atinf(n),Attinf(n),Attint(n)
      double complex z(n),b,z2
      real*8  omega(n),t,u,tb,mu,domega,lh
      real*8 Temp,sgb,expon,nf,w,dtime,theta
      real*8 wfftc(2*n+15), cpy(n)
      integer kmu1,kmu2
      real*8 y,mu1,mu2

      external fft_ftot, fftMttof  
            
      common/corre1/ wfftc,cpy
     
      lh=dfloat(l)
      print*,'lh=',lh,l
      do i=1,n
         z2=z(i)+mu
         if (inguess.eq.1) then
            b=cdsqrt(z2*z2-4.d0*t*t)
            sgb=dabs(dimag(b))/dimag(b)
            gafin(i)=-(1.d0/2.d0/t/t)*(z2-sgb*b)
         else          
c          gafin(i)=0.5d0*(1.d0/(z2-u/2.d0-0.001*(0,1))+
c     *    1.d0/(z2+u/2.d0-0.001*(0,1)))
           gafin(i)=0
           do lat=-1,1,2
              z2=z(i)+mu+u/2*dfloat(lat)
              tb=u/8
              b=cdsqrt(z2*z2-4.d0*tb*tb)
              sgb=dabs(dimag(b))/dimag(b)
  
              if (lat.eq.1) then
                gafin(i)=-(1.d0/2.d0/tb/tb)*(z2-sgb*b)/2.d0/lh
     *                   +gafin(i)
              else
                gafin(i)=-(1.d0/2.d0/tb/tb)*(z2-sgb*b)/2.d0
     *              *(2.d0-1.d0/lh)
     *                   +gafin(i)
              endif  
           enddo
         endif
         if (inguess.eq.3) then
            z2=z(i)+dabs(mu)/mu*u/2
            w=u/2+1
            b=cdsqrt(z2*z2-4.d0*t*t*w*w)
            sgb=dabs(dimag(b))/dimag(b)
            gafin(i)=-(1.d0/2.d0/t/t/w/w)*(z2-sgb*b)
         endif
       enddo
         kmu1=mu/domega
         kmu2=kmu1+1
         mu1=kmu1*domega
         mu2=kmu2*domega
       do i=1,n
c         if (((i+kmu1).gt.0).and.((i+kmu2).lt.n)) then
c            y1=1/pi*dimag(gafin(i+kmu1))
c            y2=1/pi*dimag(gafin(i+kmu2))
c            y=(y1-y2)/(mu1-mu2)*(mu-mu1)+y1
c         else
c            y=0.d0
c         endif
         y=1/pi*dimag(gafin(i))

         Ainf(i)=y

         expon=omega(i)/Temp
         if (dabs(expon).gt.1.d2) then
             if (i.eq.n/2) then 
                 nf=0.5d0
             endif
             if (i.lt.n/2) then
                nf=1.d0
             endif
             if (i.gt.n/2) then
                nf=0.d0
             endif
         else
            nf=1.d0/(1.d0+dexp(expon))
         endif

         Atinf(i)=nf*Ainf(i)
         Attinf(i)=(1-nf)*Ainf(i)
      enddo
            
      call fft_ftot(Atinf, Atint,domega)
      call fft_ftot(Attinf, Attint, domega)
      call fft_ftot(Ainf,  Aint,domega)

      do i=1,n
         mi=1+n-i
         if (i.gt.n/2) then 
             theta=1.d0 
         else 
             theta=0.d0 
         endif
         if (i.eq.n/2) then
               theta=0.5
         endif
         gtin(2,2,i)=-2.d0*pi*zi*(theta*Attint(i)-
     *                        (1.-theta)*Atint(i))
C         gtin(1,1,mi)=-dconjg(gtin(2,2,i))
         gtin(2,1,i)=2.d0*pi*zi*Atint(i)
         gtin(1,2,i)=-2.d0*pi*zi*(Attint(i))
         gtin(1,1,i)=gtin(1,2,i)+gtin(2,1,i)-gtin(2,2,i)
      enddo

      call fftMttof(gtin,gfin,dtime)
      do i=1,n
          gafin(i)=gfin(2,2,i)-gfin(1,2,i)
      enddo

      return
      end
         
C====================================================================

C-----function for determining mu_tilde
      
      real*8 function f(dmu,gafinvhelp,nocc,domega,omega,Temp)

      real*8 pi
      parameter(pi=3.141592653589793d0)
      
      integer n
      parameter(n=2048*4)
C     (n must be a power of 2)

      double complex gafinvhelp(n),g0
      real*8 dmu,nocc,nhelp,omega(n),norm
      real*8 expon,Temp,domega,nf
      integer i

      nhelp=0.d0
      norm=0.d0
      do i=4,n-2,2
         g0=1./(gafinvhelp(i)+dmu)
         expon=omega(i)/Temp
         if (dabs(expon).ge.1.d2) then 
             if (i.eq.n/2) then 
                 nf=0.5d0
             endif
             if (i.lt.n/2) then
                nf=1.d0
             endif
             if (i.gt.n/2) then
                nf=0.d0
             endif
         else
             nf=1.d0/(1.d0+dexp(expon))
         endif
         nhelp=nhelp+1.d0/pi*dimag(g0)*nf*domega*2.d0/3.d0
         norm=norm+1.d0/pi*dimag(g0)*domega*2.d0/3.d0
       enddo
      do i=3,n-1,2
         g0=1./(gafinvhelp(i)+dmu)
         expon=omega(i)/Temp
         if (dabs(expon).ge.1.d2) then 
             if (i.eq.n/2) then 
                 nf=0.5d0
             endif
             if (i.lt.n/2) then
                nf=1.d0
             endif
             if (i.gt.n/2) then
                nf=0.d0
             endif
         else
             nf=1.d0/(1.d0+dexp(expon))
         endif
         nhelp=nhelp+1.d0/pi*dimag(g0)*nf*domega*4.d0/3.d0
         norm=norm+1.d0/pi*dimag(g0)*domega*4.d0/3.d0
       enddo
      do i=2,n,(n-2)
         g0=1./(gafinvhelp(i)+dmu)
         expon=omega(i)/Temp
         if (dabs(expon).ge.1.d2) then 
             if (i.eq.n/2) then 
                 nf=0.5d0
             endif
             if (i.lt.n/2) then
                nf=1.d0
             endif
             if (i.gt.n/2) then
                nf=0.d0
             endif
         else
             nf=1.d0/(1.d0+dexp(expon))
         endif
         nhelp=nhelp+1.d0/pi*dimag(g0)*nf*domega*1.d0/3.d0
         norm=norm+1.d0/pi*dimag(g0)*domega*1.d0/3.d0
       enddo

       nhelp=nhelp/norm
       f=nhelp-nocc
       return
       end



