code: def calczeta(phi1, phi2, theta1, theta2):
      
      
      zeta = 0.0
      
      if phi1 == phi2 and theta1 == theta2:
          zeta = 0.0
      else:
          argument = sin(theta1)*sin(theta2)*cos(phi1-phi2) + \
            cos(theta1)*cos(theta2)
            
          if argument < -1:
              zeta = np.pi
          elif argument > 1:
              zeta = 0.0
          else:
              zeta = acos(argument)
      
      return zeta
docs: Calculate the angular separation between position (phi1, theta1) and
      (phi2, theta2)