Electron Gas Susceptibility Code

c electron gas susceptibilities.
      parameter(n=400)
      dimension a(n)
      complex z,z1,z2
      complex g1,g2
      y0=5.0
      eta=0.001
      pi=3.1415927
c compute energy bin.
      y=n
      dy=y0/y
c prompt for x.
      write(*,1)
    1 format('Enter value for x.')
      read(*,*) x
      write(*,*) x
c open output file.
      open(7, file="output.dat", status="new")
c generate weights for transverse correction f.
      do 3 i=1,n
      y=(i-0.5)*dy
      a(i)=0.0
      yy=((y-x**2)/2.0/x)**2
      if(yy.lt.eta) yy=eta      
      if(yy.gt.1.0) go to 2
      a(i)=-yy*alog(yy)/x/2.0
    2 continue
      yy=((-y-x**2)/2.0/x)**2
      if(yy.lt.eta) yy=eta      
      if(yy.gt.1.0) go to 3
      a(i)=a(i)+yy*alog(yy)/x/2.0
    3 continue
c loop on frequencies.
      do 6 i=1,n
      y=(i-0.5)*dy
      z=cmplx(y,eta)
      z1=clog((z-x**2+2.0*x)/(z-x**2-2.0*x))
      z2=clog((-z-x**2+2.0*x)/(-z-x**2-2.0*x))
c compute longitudinal syceptibility.
      g1=-1
     *-((z-x**2)**2-(2.0*x)**2)/8.0/x**3*z1
     *-((-z-x**2)**2-(2.0*x)**2)/8.0/x**3*z2
c hilbert transform.
      r=0.0
      do 4 ii=1,n
      yy=(ii+i-1.0)
      r=r-a(ii)/yy
      if(i.eq.ii) go to 4
      yy=(i-ii)
      r=r+a(ii)/yy
    4 continue
c transverse correction.
      g2=cmplx(-r,a(i)*pi)
c compute transverse susceptibility.
      g2=g2+g1
c output results.
      write(7,5) y,g1,g2
    5 format(f10.4,4e18.7)
    6 continue
      stop
      end