code: def gha(julian_day, f):
      
      rad = old_div(np.pi, 180.)
      d = julian_day - 2451545.0 + f
      L = 280.460 + 0.9856474 * d
      g = 357.528 + 0.9856003 * d
      L = L % 360.
      g = g % 360.
  # ecliptic longitude
      lamb = L + 1.915 * np.sin(g * rad) + .02 * np.sin(2 * g * rad)
  # obliquity of ecliptic
      epsilon = 23.439 - 0.0000004 * d
  # right ascension (in same quadrant as lambda)
      t = (np.tan(old_div((epsilon * rad), 2)))**2
      r = old_div(1, rad)
      rl = lamb * rad
      alpha = lamb - r * t * np.sin(2 * rl) + \
          (old_div(r, 2)) * t * t * np.sin(4 * rl)
  #       alpha=mod(alpha,360.0)
  # declination
      delta = np.sin(epsilon * rad) * np.sin(lamb * rad)
      delta = old_div(np.arcsin(delta), rad)
  # equation of time
      eqt = (L - alpha)
  #
      utm = f * 24 * 60
      H = old_div(utm, 4) + eqt + 180
      H = H % 360.0
      return H, delta
docs: returns greenwich hour angle