C     ===================================================================   
      
C     ------- This subroutine initializes the fourier transforms and -----
C     --------     sets the coefficients for endpoint corrections---------
C     ------------         as well as the constants used -----------------
C      

      subroutine init(domega,tomega,dtime,z,time,omega,
     *                          ttime,acc,eta)
      
C set parameters

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

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



      integer j
      real*8 aif0(n/2), arf0(n/2), aif1(n/2), arf1(n/2),
     *                  aif2(n/2), arf2(n/2), aif3(n/2), arf3(n/2),
     *                aaf0(n/2), corfaf(n/2)
      real*8 ait0(n/2), art0(n/2), ait1(n/2), art1(n/2),
     *                  ait2(n/2), art2(n/2), ait3(n/2), art3(n/2),
     *                  aat0(n/2), corfat(n/2) 
      real*8 domega, dtime,tomega,time(n),omega(n),
     *                    ttime,acc,omega2(n), time2(n)
      real*8 corfac,eta
      real*8 a0i,a0r,a1i,a1r,a2i,a2r,a3i,a3r,aa0,s,c
      double complex z(n)
      real*8 wfftc(2*n+15),cpy(n)
        
      


      external  dftcorm
C     external iwkin
      common/corre1/ wfftc, cpy
      common/corre2/ aif0, arf0, aif1, arf1,
     *               aif2, arf2, aif3, arf3, aaf0, corfaf
      common/corre3/ ait0, art0, ait1, art1,
     *               ait2, art2, ait3, art3, aat0, corfat  
      common/corre4/ a0i,a0r,a1i,a1r,a2i,a2r,a3i,a3r,aa0,s,c
C      common/worksp/ rwksp

C ----------initialization of fft---------------------------------

C      call iwkin(786482)


C ----------set constants-----------------------------------------

      domega=tomega/dfloat(n)
      dtime=2.d0*pi/(tomega/2.d0)
      ttime=dtime*dfloat(n)
      do j=1,n
         omega(j)=domega*(dfloat(j)-dfloat(n)/2.d0)
         omega2(j)=domega*(dfloat(j-1))
         z(j)=omega(j)-zi*eta
         time(j)=dtime*(dfloat(j)-dfloat(n)/2.d0)
         time2(j)=dtime*(dfloat(j-1))
      enddo
      acc=acc*acc
      
C ---------Initialization of coefficients needed for endpoint corrections---
      
      do j=1,n/2
         call dftcorm(omega2(j), dtime,0.d0,ttime/2.d0,
     *                    corfac)
         corfaf(j)=corfac
         aif0(j)=a0i
         arf0(j)=a0r
         aif1(j)=a1i
         arf1(j)=a1r
         aif2(j)=a2i
         arf2(j)=a2r
         aif3(j)=a3i
         arf3(j)=a3r
         aaf0(j)=aa0
         
         call dftcorm(-time2(j), domega,0.d0, 
     *          tomega/2.d0,corfac)

         corfat(j)=corfac
         ait0(j)=a0i
         art0(j)=a0r
         ait1(j)=a1i
         art1(j)=a1r
         ait2(j)=a2i
         art2(j)=a2r
         ait3(j)=a3i
         art3(j)=a3r
         aat0(j)=aa0
      enddo
C-------- Initialization of datafile for output-------------------
c      write(50,*) 
c     *  'TitleText: Hubbard-DOS in 2nd Order Pertubation Theory'
c      write(50,*) 'XUnitText: energy (D)'
c      write(50,*) 'YUnitText: DOS (1/D)'
      return
      end

C     ===================================================================   

c subroutine dftcorm is the modidied imsl routine computing attenuation
c coefficients and coeff. for corrections to fourier integrals
c
      subroutine dftcorm(w,delta,a,b,corfac)
c     
      real*8 a, b, corfac, delta, w, s, c
      real*8 a0i, a0r, a1i, a1r, a2i, a2r, a3i, a3r, arg, aa0
      real*8 cth, ctth, spth2, sth, sth4i, stth, th, th2
      real*8 th4, th6, tmth2, tth4i
c
      common/corre4/  a0i, a0r, a1i, a1r, a2i, a2r, a3i, a3r, aa0, s,c
c
      th=w*delta
      th2=th*th
      th4=th2*th2
      th6=th2*th6
      if (a.ge.b.or.th.gt.2*3.14156d0) pause 'bad args for dftcorm'
      if (dabs(th).lt.5.d-2) then
         corfac=1.d0-(11.d0/720.d0)*th4+(23.d0/15120.d0)*th6
         a0r=(-2.d0/3.d0)+th2/45.d0+(103.d0/15120.d0)*th4
     $      -(169.d0/226800.d0)*th6
         a1r=(7.d0/24.d0)-(7.d0/180.d0)*th2+(5.d0/3456.d0)*th4
     $      -(7.d0/259200.d0)*th6
         a2r=(-1.d0/6.d0)+th2/45.d0-(5.d0/6048.d0)*th4+th6/64800.d0
         a3r=(1.d0/24.d0)-th2/180.d0+(5.d0/24192.d0)*th4-th6/259200.d0
         a0i=th*(2.d0/45.d0+(2.d0/105.d0)*th2-(8.d0/2835.d0)*th4
     $      +(86.d0/467775.d0)*th6)
         a1i=th*(7.d0/72.d0-th2/168.d0+(11.d0/72576.d0)*th4
     $      -(13.d0/5987520.d0)*th6)
         a2i=th*(-7.d0/90.d0+th2/210.d0-(11.d0/90720.d0)*th4
     $      +(13.d0/7484400.d0)*th6)
         a3i=th*(7.d0/360.d0-th2/840.d0+(11.d0/362880.d0)*th4
     $      -(13.d0/29937600.d0)*th6)
         aa0=(1.d0/3.d0)+(1.d0/45.d0)*th2-(8.d0/945.d0)*th4
     $      +(11.d0/14175.d0)*th6
      else
         cth=dcos(th)
         sth=dsin(th)
         ctth=cth**2-sth**2
         stth=2.d0*sth*cth
         tmth2=3.d0-th2
         spth2=6.d0+th2
         sth4i=1.d0/(6.d0*th4)
         tth4i=2.d0*sth4i
         corfac=tth4i*spth2*(3.d0-4.d0*cth+ctth)
         a0r=sth4i*(-42.d0+5.d0*th2+spth2*(8.d0*cth-ctth))
         a0i=sth4i*(th*(-12.d0+6.d0*th2)+spth2*stth)
         a1r=sth4i*(14.d0*tmth2-7.d0*spth2*cth)
         a1i=sth4i*(30.d0*th-5.d0*spth2*sth)
         a2r=tth4i*(-4.d0*tmth2+2.d0*spth2*cth)
         a2i=tth4i*(-12.d0*th+2.d0*spth2*sth)
         a3r=sth4i*(2.d0*tmth2-spth2*cth)
         a3i=sth4i*(6.d0*th-spth2*sth)
         aa0=sth4i*((-6.d0+11.d0*th2)+spth2*dcos(2.d0*th))
      endif
      arg=w*(b-a)
      c=dcos(arg)
      s=dsin(arg)

      return
      end

C     ===================================================================   
         
C ------The following subroutine fouriertansforms a discrete ---------
C--------   function from  frequency space to time space ------------
      subroutine fft_ftot(qqf,qqt,domega)
      integer n,i
      parameter(n=2048*4)
      real*8 pi
      parameter(pi=3.141592653589793d0)
      double complex zi
      parameter(zi=(0.d0,1.d0))
      double complex qqts(n/2),qqtsr(n/2),qqtsi(n/2),qqt(n),
     *               qqfs(n/2),qqf(n),qqfsr(n/2),qqfsi(n/2),
     *               qqta(n/2),qqtar(n/2),qqtai(n/2),
     *               qqfa(n/2),qqfar(n/2),qqfai(n/2)  
      real*8 domega
      external fftmin,fftplus
  
C     Decompose function in odd and even parts
      do i=1,(n/2)
         qqfs(i)=0.5d0*(qqf(n/2+i-1)+qqf(n/2-i+1))
         qqfa(i)=0.5d0*(qqf(n/2+i-1)-qqf(n/2-i+1))
      enddo
 

      do i=1,n/2
         qqfsr(i)=dble(qqfs(i))
         qqfsi(i)=dimag(qqfs(i))
         qqfar(i)=dble(qqfa(i))
         qqfai(i)=dimag(qqfa(i))
      enddo

C     Fouriertransform the even part
      call fftmin(qqfsr,qqtsr,domega)
      call fftmin(qqfsi,qqtsi,domega)
      do i=1,n/2
        qqts(i)=(dble(qqtsr(i))+zi*dble(qqtsi(i)))/pi  
      enddo
  
C     Fouriertransform the odd part
      call fftmin(qqfar,qqtar,domega)
      call fftmin(qqfai,qqtai,domega)
      do i=1,n/2
       qqta(i)=-2d0*(-zi*dimag(qqtar(i))
     *                 +dimag(qqtai(i)))/2.d0/pi
      enddo

C     Get the full time dependent function
      do i=1,n/2
        qqt(i)=qqts(n/2+1-i)-qqta(n/2+1-i)
        if (i.lt.n/2) then
            qqt(i+n/2)=qqts(i+1)+qqta(i+1)
        endif
      enddo
      qqt(n)=qqt(n-1)
     
      return
      end

C     ===================================================================

C ------The following subroutine fouriertansforms a discrete ---------
C--------   function from  time space to frequency space ------------
      subroutine fft_ttof(qqt,qqf,dtime)
      integer n,i
      real*8 pi
      parameter(pi=3.141592653589793d0)
      parameter(n=2048*4)
      double complex zi
      parameter(zi=(0.d0,1.d0))
      double complex qqfs(n/2),qqfsr(n/2),qqfsi(n/2),qqf(n),
     *               qqts(n/2),qqt(n),qqtsr(n/2),qqtsi(n/2),
     *               qqfa(n/2),qqfar(n/2),qqfai(n/2),
     *               qqta(n/2),qqtar(n/2),qqtai(n/2)  
      real*8 dtime
      external fftmin,fftplus
  
C     Decompose function in odd and even parts
      do i=1,n/2
         qqts(i)=0.5d0*(qqt(n/2+i-1)+qqt(n/2-i+1))
         qqta(i)=0.5d0*(qqt(n/2+i-1)-qqt(n/2-i+1))
      enddo

      do i=1,n/2
         qqtsr(i)=dble(qqts(i))
         qqtsi(i)=dimag(qqts(i))
         qqtar(i)=dble(qqta(i))
         qqtai(i)=dimag(qqta(i))
      enddo 
 
C     Fouriertransform the even part
      call fftplus(qqtsr,qqfsr,dtime)
      call fftplus(qqtsi,qqfsi,dtime)
      do i=1,n/2
        qqfs(i)=2.d0*(dble(qqfsr(i))+zi*dble(qqfsi(i)))
      enddo
  
C     Fouriertransform the odd part
      call fftplus(qqtar,qqfar,dtime)
      call fftplus(qqtai,qqfai,dtime)
      do i=1,n/2
        qqfa(i)=2.d0*(zi*dimag(qqfar(i))-dimag(qqfai(i)))
      enddo

C     Get the full frequency dependent function
      do i=1,n/2
        qqf(i)=qqfs(n/2+1-i)-qqfa(n/2+1-i)
        if (i.lt.n/2) then
         qqf(i+n/2)=qqfs(i+1)+qqfa(i+1)
        endif
      enddo
      qqf(n)=qqf(n-1)

      return
      end

C     ===================================================================

C ------- This subroutine fouriertransforms in (0, infinity) with the
C ---------------      integrant  f*exp(- omega*t*i)--------------------
      subroutine fftmin(qin,qout,delta)
C set parameters
      real*8 pi
      parameter(pi=3.141592653589793d0)
      
      integer n
      parameter(n=2048*4)
C     (n must be a power of 2.)
 
      double complex zi    
      parameter(zi=(0.d0, 1.d0))

      integer i
      real*8  delta, endpt(8)
      real*8 wfftc(2*n+15), cpy(n)
      real*8 ait0(n/2), art0(n/2),  ait1(n/2), art1(n/2),
     *                  ait2(n/2), art2(n/2),  ait3(n/2), art3(n/2),
     *                  aat0(n/2), corfat(n/2)
      double complex qin(n/2), qout(n/2),ca
       

      common / corre1 / wfftc, cpy
      common / corre2  / ait0, art0, ait1, art1,
     *                   ait2, art2, ait3, art3, aat0, corfat

C    -------- store endpoints needed for endpoint corrections ------
C              Set qin(n/2+1) = qin(n/2)
      do i=1,4
        endpt(i) = qin(i)
      enddo
      do i=1,3
         endpt(4+i)=qin(n/2-3+i)
      enddo
      endpt(8)= endpt(7)

C --------- call subroutine for dft -------------------
      do i=1,n/2
          qout(i)=qin(i)
      enddo

      call dfour1(qout,n/2,-1)

C --------compute endpoint corrections for attenuation coefficients
      do i=1,n/2
         ca= (art0(i)+ zi * ait0(i))*endpt(1)
     *      +(art1(i)+ zi * ait1(i))*endpt(2)
     *      +(art2(i)+ zi * ait2(i))*endpt(3)
     *      +(art3(i)+ zi * ait3(i))*endpt(4)
     *      +(aat0(i)- zi * ait0(i))*endpt(8)
     *      +(art1(i)- zi * ait1(i))*endpt(7)
     *      +(art2(i)- zi * ait2(i))*endpt(6)
     *      +(art3(i)- zi * ait3(i))*endpt(5)
         qout(i)=delta*(corfat(i)*qout(i)+ca)
      enddo
     

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

C ------The following subroutine fouriertransforms in (0, infinity)
C---------- with the integrant f*exp(+omega*t*i)--------------------

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

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

      integer  i
      real*8 delta,endpt(8)
      real*8 wfftc(2*n+15), cpy(n)
      real*8 aif0(n/2), arf0(n/2),  aif1(n/2), arf1(n/2),
     *                  aif2(n/2), arf2(n/2),  aif3(n/2), arf3(n/2),
     *                  aaf0(n/2), corfaf(n/2)
      double complex qin(n/2), qout(n/2),ca
       
      common / corre1 / wfftc, cpy
      common / corre3  / aif0, arf0, aif1, arf1,
     *                   aif2, arf2, aif3, arf3, aaf0, corfaf


C    -------- store endpoints needed for endpoint corrections ------
C              Set qin(n/2+1) = qin(n/2)
      do i=1,4
        endpt(i) = qin(i)
      enddo
      do i=1,3
         endpt(4+i)=qin(n/2-3+i)
      enddo
      endpt(8)= endpt(7)
      
C --------- call subroutine for dft -------------------
      
      do i=1,n/2
         qout(i)=qin(i)
      enddo

      call dfour1(qout,n/2,1)


C --------compute endpoint corrections for attenuation coefficients
      do i=1,n/2
         ca= (arf0(i)+ zi * aif0(i))*endpt(1)
     *      +(arf1(i)+ zi * aif1(i))*endpt(2)
     *      +(arf2(i)+ zi * aif2(i))*endpt(3)
     *      +(arf3(i)+ zi * aif3(i))*endpt(4)
     *      +(aaf0(i)- zi * aif0(i))*endpt(8)
     *      +(arf1(i)- zi * aif1(i))*endpt(7)
     *      +(arf2(i)- zi * aif2(i))*endpt(6)
     *      +(arf3(i)- zi * aif3(i))*endpt(5)
         qout(i)=delta*(corfaf(i)*qout(i)+ca)
      enddo

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

C If a matrix has to be fourier transformed, use the following subroutine!
      subroutine fftMttof(Mt,Mf,dtime)
C set parameters
      real*8 pi
      parameter(pi=3.141592653589793d0)
      
      integer n
      parameter(n=2048*4)
C     (n must be a power of 2.)

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

      integer i
      double complex Mf(2,2,n), Mt(2,2,n),
     *               Mt22(n),   Mt12(n),   Mt21(n),   Mt11(n),
     *               Mf22(n),   Mf12(n),   Mf21(n),   Mf11(n)
      real*8 dtime
       
      external fft_ttof

      do i=1,n
          Mt11(i)=Mt(1,1,i)
          Mt12(i)=Mt(1,2,i)
          Mt21(i)=Mt(2,1,i)
          Mt22(i)=Mt(2,2,i)
      enddo

      
      call fft_ttof(Mt11,Mf11,dtime)
      call fft_ttof(Mt12,Mf12,dtime)         
      call fft_ttof(Mt21,Mf21,dtime)
      call fft_ttof(Mt22,Mf22,dtime)


      do i=1,n
          Mf(1,1,i)= Mf11(i)
          Mf(1,2,i)= Mf12(i)
          Mf(2,1,i)= Mf21(i)
          Mf(2,2,i)= Mf22(i)
      enddo

      return
      end

C     ===================================================================

C If a matrix has to be fourier transformed, use the following subroutine!

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

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

      integer i
      double complex Mf(2,2,n), Mt(2,2,n),
     *               Mt22(n),   Mt12(n),   Mt21(n),   Mt11(n),
     *               Mf22(n),   Mf12(n),   Mf21(n),   Mf11(n)
      real*8 domega
      external fft_ftot
     
      do i=1,n
          Mf11(i)=Mf(1,1,i)
          Mf12(i)=Mf(1,2,i)
          Mf21(i)=Mf(2,1,i)
          Mf22(i)=Mf(2,2,i)
      enddo
      call fft_ftot(Mf11,Mt11,domega)
      call fft_ftot(Mf12,Mt12,domega)         
      call fft_ftot(Mf21,Mt21,domega)
      call fft_ftot(Mf22,Mt22,domega)
      do i=1,n
          Mt(1,1,i)= Mt11(i)
          Mt(1,2,i)= Mt12(i)
          Mt(2,1,i)= Mt21(i)
          Mt(2,2,i)= Mt22(i)
      enddo

      return
      end

C     ===================================================================

C--- subroutine for inverting a 2 by 2 matrix for all n frequencies ---

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

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

      integer i
      double complex A(2,2,n), B(2,2,n), det

      do i=1,n
         det=A(1,1,i)*A(2,2,i)-A(1,2,i)*A(2,1,i)
         if (cdabs(det).lt.1.d-80) then
            print*, 'Attention! Matrix might not be invertible:',i
         endif
      
         B(1,1,i)=A(2,2,i)/det
         B(2,2,i)=A(1,1,i)/det
         B(1,2,i)=-A(1,2,i)/det
         B(2,1,i)=-A(2,1,i)/det
      enddo
      
      return
      end
      
C     ===================================================================

      SUBROUTINE dfour1(data,nn,isign)
      INTEGER isign,nn
      REAL*8 data(2*nn)
      INTEGER i,istep,j,m,mmax,n
      REAL*8 tempi,tempr
      REAL*8 theta,wi,wpi,wpr,wr,wtemp
      n=2*nn
      j=1
      do 11 i=1,n,2
        if(j.gt.i)then
          tempr=data(j)
          tempi=data(j+1)
          data(j)=data(i)
          data(j+1)=data(i+1)
          data(i)=tempr
          data(i+1)=tempi
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
          j=j-m
          m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2
2     if (n.gt.mmax) then
        istep=2*mmax
        theta=6.28318530717959d0/(isign*mmax)
        wpr=-2.d0*sin(0.5d0*theta)**2
        wpi=sin(theta)
        wr=1.d0
        wi=0.d0
        do 13 m=1,mmax,2
          do 12 i=m,n,istep
            j=i+mmax
            tempr=wr*data(j)-wi*data(j+1)
            tempi=wr*data(j+1)+wi*data(j)
            data(j)=data(i)-tempr
            data(j+1)=data(i+1)-tempi
            data(i)=data(i)+tempr
            data(i+1)=data(i+1)+tempi
12        continue
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
      goto 2
      endif
      return
      END

C     ===================================================================

      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      real*8 yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      real*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END

C     ===================================================================

      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      real*8 x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END

C     ===================================================================
