code: def vgp_di(plat, plong, slat, slong):
      
      plong = plong % 360
      slong = slong % 360
      signdec = 1.
      delphi = abs(plong - slong)
      if delphi != 0:
          signdec = (plong - slong) / delphi
      if slat == 90.:
          slat = 89.99
      thetaS = np.radians(90. - slat)
      thetaP = np.radians(90. - plat)
      delphi = np.radians(delphi)
      cosp = np.cos(thetaS) * np.cos(thetaP) + np.sin(thetaS) * \
          np.sin(thetaP) * np.cos(delphi)
      thetaM = np.arccos(cosp)
      cosd = old_div((np.cos(thetaP) - np.cos(thetaM) *
                      np.cos(thetaS)), (np.sin(thetaM) * np.sin(thetaS)))
      C = abs(1. - cosd**2)
      if C != 0:
          dec = -np.arctan(cosd/np.sqrt(abs(C))) + (np.pi/2.)
      else:
          dec = np.arccos(cosd)
      if -np.pi < signdec * delphi and signdec < 0:
          dec = 2. * np.pi - dec  # checking quadrant
      if signdec * delphi > np.pi:
          dec = 2. * np.pi - dec
      dec = np.degrees(dec) % 360.
      inc = np.degrees(np.arctan2(2. * np.cos(thetaM), np.sin(thetaM)))
      return dec, inc
docs: Converts a pole position (pole latitude, pole longitude) to a direction
      (declination, inclination) at a given location (slat, slong) assuming a
      dipolar field.
  
      Parameters
      ----------
      plat : latitude of pole (vgp latitude)
      plong : longitude of pole (vgp longitude)
      slat : latitude of site
      slong : longitude of site
  
      Returns
      ----------
      dec,inc : tuple of declination and inclination