program gr
  implicit none
  real(8) :: p, dp, sig, a, b, a1
  real(8) :: Cr1, Cr2, Ci
  integer :: i
  integer, parameter :: N=1000000
  dp=4./dble(N)

  open(17, file='eady_velo.dat', status='unknown')
  do i = 1,N
     p=dp*dble(i)
     a = 1. / tanh(p/2.) - p/2.     
     b = p / 2. - tanh(p/2.)     
     if ( a .lt. 0 ) then
     a1 = abs(a)    
     a = 0.
     sig = sqrt(a1*b)
     Cr1 = 0.5 + sig/p
     Cr2 = 0.5 - sig/p     
     Ci = 0.
     write (17,'(4e12.4)') p, Cr1, Cr2, Ci
  else
     sig = sqrt(a*b)
     Ci = sig / p
     Cr1 = 0.5
     Cr2 = 0.5
     write (17,'(4e12.4)') p, Cr1, Cr2, Ci
  end if
  end do
  close(17)
end program gr
