anchor.diff:  program susy_qpi
     implicit none
  -  integer, parameter :: dp=kind(0.d0)
  -  integer :: i,j
  -
  -  do i=1,100
  -    do j=1,10
  -      call write_data(i*0.5_dp,j*0.2_dp)
  -    end do
  -  end do
  -
  -end program susy_qpi
  -
  -subroutine write_data(om, del)
  -  implicit none
  -  external sG0
   
   ! Double precision
     integer, parameter :: dp=kind(0.d0)
  -
  -! Parameters
  -!    omega = om + i * del
  -  real(dp), intent(in) :: om, del
   
   ! Global variables
   !   epsq = -2 * t ( cos(kx) + cos(ky) ) - mu
  
positive.code: !     epsf   real.
  !   Inv[G(k,omega)] = (omega)*I - H_k
  !     Gkinv  real array of dimension 16 (4x4 matrix)
  !            Inverse of greens function
  !     Gkqinv real array of dimension 16 (4x4 matrix)
  !            Inverse of greens function at k+q
  !     omega  complex.
  !            Frequency (energy) used for G
  
    common /gvars/      t, mu, x0, epsf, qx, qy, Pi, Gkinv, Gkqinv
    real(dp)         :: t, mu, x0, epsf, qx, qy, Pi, Uc, Uf, V
    complex(dp)      :: Gkinv(16), Gkqinv(16)
    common /freq/       omega
    complex(dp)      :: omega
  
  ! Local variables
  !     steps  integer.
  
negative.code:   abserr = 1e-6_dp
    relerr = 1e-3_dp
  
  ! frequency
    omega = dcmplx(om,del)
  
  ! Step through 1/8 triangle of BZ
    steps=100
  