double precision function datan2 (sn, cs)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision sn, cs, pi, abssn, abscs, sml, big, d1mach,
     1  datan
      external d1mach, datan
      data pi / 3.1415926535 8979323846 2643383279 50 d0 /
      data sml, big /2*0.d0/
c
      if (sml.ne.0.d0) go to 10
      sml = d1mach(1)
      big = d1mach(2)
c
c we now make sure sn can be divided by cs.  it is painful.
c
 10   abssn = dabs(sn)
      abscs = dabs(cs)
      if (abscs.le.abssn) go to 20
c
      if (abscs.le.1.0d0) go to 30
      if (abssn.gt.abscs*sml) go to 30
      datan2 = 0.0d0
      go to 40
c
 20   if (abscs.ge.1.0d0) go to 30
      if (abssn.lt.abscs*big) go to 30
c
      if (sn.eq.0.0d0) call seteru (
     1  31hdatan2  both arguments are zero, 31, 1, 2)
      datan2 = dsign (0.5d0*pi, sn)
      return
c
 30   datan2 = datan (sn/cs)
 40   if (cs.lt.0.d0) datan2 = datan2 + pi
      if (datan2.gt.pi) datan2 = datan2 - 2.0d0*pi
      return
c
      end