def vgp_di(plat, plong, slat, slong):
    """
    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
    """
    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
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 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
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