function atan2 (sn, cs)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      external atan, r1mach
      data pi /3.1415926535 8979323846e0/
      data sml, big / 2*0.0 /
c
      if (sml.ne.0.0) go to 10
      sml = r1mach(1)
      big = r1mach(2)
c
c we now make sure sn can be divided by cs.  it is painful.
c
 10   abssn = abs(sn)
      abscs = abs(cs)
      if (abscs.le.abssn) go to 20
c
      if (abscs.le.1.0) go to 30
      if (abssn.gt.abscs*sml) go to 30
      atan2 = 0.0
      go to 40
c
 20   if (abscs.ge.1.0) go to 30
      if (abssn.lt.abscs*big) go to 30
c
      if (sn.eq.0.0) call seteru (
     1  31hatan2   both arguments are zero, 31, 1, 2)
      atan2 = sign (0.5*pi, sn)
      return
c
 30   atan2 = atan (sn/cs)
 40   if (cs.lt.0.0) atan2 = atan2 + pi
      if (atan2.gt.pi) atan2 = atan2 - 2.0*pi
      return
c
      end

.