code: def doeqdi(x, y, UP=False):
      
      xp, yp = y, x  # need to switch into geographic convention
      r = np.sqrt(xp**2+yp**2)
      z = 1.-r**2
      t = np.arcsin(z)
      if UP == 1:
          t = -t
      p = np.arctan2(yp, xp)
      dec, inc = np.degrees(p) % 360, np.degrees(t)
      return dec, inc
docs: Takes digitized x,y, data and returns the dec,inc, assuming an
      equal area projection
      Parameters
      __________________
          x : array of digitized x from point on equal area projection
          y : array of  igitized y from point on equal area projection
          UP : if True, is an upper hemisphere projection
      Output :
          dec : declination
          inc : inclination