module utils
   use nrtype
   implicit none
   real(dp),parameter::dw=0.05d0,gam=0.4d0,eps=1.0d0
   integer,parameter::L=512*2*4
   complex(dpc),parameter::I=cmplx(0.0d0,1.0d0)
end module

program hellofour
   use utils
   use nr
   !use nrtype
   implicit none
   real(dp)::wj,tk
   complex(dpc)::green,sigma,zk,zj
   complex(dpc),dimension(0:L-1)::g,sig
   integer::j,k

   call tabulate_g(g)
   do j=1,L-1,2; g(j)=-g(j); enddo
   call four1(g,-1)
   do j=1,L-1,2; g(j)=-g(j); enddo
   g=g*dw/twopi_d

   do k=0,L-1
      sig(k)=-g(k)**3
   enddo

   do j=1,L-1,2; sig(j)=-sig(j); enddo
   call four1(sig,+1)
   do j=1,L-1,2; sig(j)=-sig(j); enddo
   sig=sig*twopi_d/L/dw

   !testing
   open(10,file='sigw.dat')
   do j=0,L-1
      wj=dw*(j-L/2)
      zj=sigma(wj)
      write(10,"(5f10.5)") wj,real(sig(j)), imag(sig(j)), real(zj),imag(zj)
   enddo
   close(10)
end program

function green(w)
   use utils
   implicit none
   real(dp),intent(in)::w
   complex(dpc)::green
   green=1.0d0/(w-eps+I*gam)+1.0d0/(w+eps-I*gam)
end function


function sigma(w)
   use utils
   implicit none
   real(dp),intent(in)::w
   complex(dpc)::sigma
   sigma=-1.0d0/(w-3.0d0*eps+3.0d0*I*gam)-1.0d0/(w+3.0d0*eps-3.0d0*I*gam)
end function

subroutine tabulate_g(g)
   use utils
   implicit none
   complex(dpc),dimension(0:L-1)::g
   complex(dpc)::green
   integer::j
   real(dp)::wj
!   do j=L/2,L-1
!      wj=dw*(j-L/2); g(j-L/2)=green(wj)
!   enddo
!   do j=0,L/2-1
!      wj=dw*(j-L/2); g(j+L/2)=green(wj)
!   enddo
   do j=0,L-1
      wj=dw*(j-L/2); g(j)=green(wj)
   enddo
end subroutine


