def directp(node0, node1, node2, node3, hypocenter, reference, pp):
    """
    Get the Direct Point and the corresponding E-path as described in
    Spudich et al. (2013). This method also provides a logical variable
    stating if the DPP calculation must consider the neighbouring patch.
    To define the intersection point(Pd) of PpPh line segment and fault plane,
    we obtain the intersection points(Pd) with each side of fault plan, and
    check which intersection point(Pd) is the one fitting the definition in
    the Chiou and Spudich(2014) directivity model.
    Two possible locations for Pd, the first case, Pd locates on the side of
    the fault patch when Pp is not inside the fault patch. The second case is
    when Pp is inside the fault patch, then Pd=Pp.

    For the first case, it follows three conditions:
    1. the PpPh and PdPh line vector are the same,
    2. PpPh >= PdPh,
    3. Pd is not inside the fault patch.

    If we can not find solution for all the four possible intersection points
    for the first case, we check if the intersection point fit the second case
    by checking if Pp is inside the fault patch.

    Because of the coordinate system mapping(from geographic system to
    Catestian system), we allow an error when we check the location. The allow
    error will keep increasing after each loop when no solution in the two
    cases are found, until the solution get obtained.

    :param node0:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment.
    :param node1:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node2:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node3:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param hypocenter:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of floating hypocenter on each segment
        calculation. In the method, we take the direction point of the
        previous fault patch as hypocentre for the current fault patch.
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of reference point for projection
    :param pp:
        the projection of the site onto the plane containing the fault
        slipped area. A numpy array.
    :returns:
        Pd, a numpy array, representing the location of direction point
        E, the distance from direction point to hypocentre.
        go_next_patch, flag indicates if the calculation goes on the next
        fault patch. 1: yes, 0: no.
    """

    # Find the intersection point Pd, by checking if the PdPh share the
    # same vector with PpPh,  and PpPh >= PdPh
    # Transform to xyz coordinate

    node0_xyz = get_xyz_from_ll(node0, reference)
    node1_xyz = get_xyz_from_ll(node1, reference)
    node2_xyz = get_xyz_from_ll(node2, reference)
    node3_xyz = get_xyz_from_ll(node3, reference)
    hypocenter_xyz = get_xyz_from_ll(hypocenter, reference)
    hypocenter_xyz = np.array(hypocenter_xyz).flatten()

    pp_xyz = pp

    e = []

    # Loop each segments on the patch to find Pd
    segment_s = [node0_xyz, node1_xyz, node2_xyz, node3_xyz]
    segment_e = [node1_xyz, node2_xyz, node3_xyz, node0_xyz]

    # set buffering bu
    buf = 0.0001
    atol = 0.0001

    loop = True
    exit_flag = False
    looptime = 0.
    while loop:
        x_min = np.min(np.array([node0_xyz[0], node1_xyz[0], node2_xyz[0],
                                node3_xyz[0]])) - buf
        x_max = np.max(np.array([node0_xyz[0], node1_xyz[0], node2_xyz[0],
                                node3_xyz[0]])) + buf
        y_min = np.min(np.array([node0_xyz[1], node1_xyz[1], node2_xyz[1],
                                node3_xyz[1]])) - buf
        y_max = np.max(np.array([node0_xyz[1], node1_xyz[1], node2_xyz[1],
                                node3_xyz[1]])) + buf
        n_seg = 0
        exit_flag = False
        for (seg_s, seg_e) in zip(segment_s, segment_e):
            seg_s = np.array(seg_s).flatten()
            seg_e = np.array(seg_e).flatten()
            p_intersect, vector1, vector2, vector3, vector4 = _intersection(
                seg_s, seg_e, pp_xyz, hypocenter_xyz)

            ppph = dst.pdist([pp, hypocenter_xyz])
            pdph = dst.pdist([p_intersect.flatten(), hypocenter_xyz])
            n_seg = n_seg + 1

            # Check that the direction of the hyp-pp and hyp-pd vectors
            # have are the same.
            if (np.allclose(vector1.flatten(), vector2,
                            atol=atol, rtol=0.)):
                if ((np.allclose(vector3.flatten(), vector4, atol=atol,
                                 rtol=0.))):

                    # Check if ppph >= pdph.
                    if (ppph >= pdph):
                        if (p_intersect[0] >= x_min) & (p_intersect[0] <=
                                                        x_max):
                            if (p_intersect[1] >= y_min) & (p_intersect[1]
                                                            <= y_max):
                                e = pdph
                                pd = p_intersect
                                exit_flag = True
                                break
            # when the pp located within the fault rupture plane, e = ppph
            if not e:
                if (pp_xyz[0] >= x_min) & (pp_xyz[0] <= x_max):
                    if (pp_xyz[1] >= y_min) & (pp_xyz[1] <= y_max):
                        pd = pp_xyz
                        e = ppph
                        exit_flag = True
        if exit_flag:
            break

        if not e:
            looptime += 1
            atol = 0.0001 * looptime
            buf = 0.0001 * looptime
    # if pd is located at 2nd fault segment, then the DPP calculation will
    # keep going on the next fault patch
    if n_seg == 2:
        go_next_patch = True
    else:
        go_next_patch = False

    return pd, e, go_next_patch
def directp ( node0 , node1 , node2 , node3 , hypocenter , reference , pp ) : node0_xyz = get_xyz_from_ll ( node0 , reference ) node1_xyz = get_xyz_from_ll ( node1 , reference ) node2_xyz = get_xyz_from_ll ( node2 , reference ) node3_xyz = get_xyz_from_ll ( node3 , reference ) hypocenter_xyz = get_xyz_from_ll ( hypocenter , reference ) hypocenter_xyz = np . array ( hypocenter_xyz ) . flatten ( ) pp_xyz = pp e = [ ] segment_s = [ node0_xyz , node1_xyz , node2_xyz , node3_xyz ] segment_e = [ node1_xyz , node2_xyz , node3_xyz , node0_xyz ] buf = 0.0001 atol = 0.0001 loop = True exit_flag = False looptime = 0. while loop : x_min = np . min ( np . array ( [ node0_xyz [ 0 ] , node1_xyz [ 0 ] , node2_xyz [ 0 ] , node3_xyz [ 0 ] ] ) ) - buf x_max = np . max ( np . array ( [ node0_xyz [ 0 ] , node1_xyz [ 0 ] , node2_xyz [ 0 ] , node3_xyz [ 0 ] ] ) ) + buf y_min = np . min ( np . array ( [ node0_xyz [ 1 ] , node1_xyz [ 1 ] , node2_xyz [ 1 ] , node3_xyz [ 1 ] ] ) ) - buf y_max = np . max ( np . array ( [ node0_xyz [ 1 ] , node1_xyz [ 1 ] , node2_xyz [ 1 ] , node3_xyz [ 1 ] ] ) ) + buf n_seg = 0 exit_flag = False for ( seg_s , seg_e ) in zip ( segment_s , segment_e ) : seg_s = np . array ( seg_s ) . flatten ( ) seg_e = np . array ( seg_e ) . flatten ( ) p_intersect , vector1 , vector2 , vector3 , vector4 = _intersection ( seg_s , seg_e , pp_xyz , hypocenter_xyz ) ppph = dst . pdist ( [ pp , hypocenter_xyz ] ) pdph = dst . pdist ( [ p_intersect . flatten ( ) , hypocenter_xyz ] ) n_seg = n_seg + 1 if ( np . allclose ( vector1 . flatten ( ) , vector2 , atol = atol , rtol = 0. ) ) : if ( ( np . allclose ( vector3 . flatten ( ) , vector4 , atol = atol , rtol = 0. ) ) ) : if ( ppph >= pdph ) : if ( p_intersect [ 0 ] >= x_min ) & ( p_intersect [ 0 ] <= x_max ) : if ( p_intersect [ 1 ] >= y_min ) & ( p_intersect [ 1 ] <= y_max ) : e = pdph pd = p_intersect exit_flag = True break if not e : if ( pp_xyz [ 0 ] >= x_min ) & ( pp_xyz [ 0 ] <= x_max ) : if ( pp_xyz [ 1 ] >= y_min ) & ( pp_xyz [ 1 ] <= y_max ) : pd = pp_xyz e = ppph exit_flag = True if exit_flag : break if not e : looptime += 1 atol = 0.0001 * looptime buf = 0.0001 * looptime if n_seg == 2 : go_next_patch = True else : go_next_patch = False return pd , e , go_next_patch
Get the Direct Point and the corresponding E-path as described in
    Spudich et al. (2013). This method also provides a logical variable
    stating if the DPP calculation must consider the neighbouring patch.
    To define the intersection point(Pd) of PpPh line segment and fault plane,
    we obtain the intersection points(Pd) with each side of fault plan, and
    check which intersection point(Pd) is the one fitting the definition in
    the Chiou and Spudich(2014) directivity model.
    Two possible locations for Pd, the first case, Pd locates on the side of
    the fault patch when Pp is not inside the fault patch. The second case is
    when Pp is inside the fault patch, then Pd=Pp.

    For the first case, it follows three conditions:
    1. the PpPh and PdPh line vector are the same,
    2. PpPh >= PdPh,
    3. Pd is not inside the fault patch.

    If we can not find solution for all the four possible intersection points
    for the first case, we check if the intersection point fit the second case
    by checking if Pp is inside the fault patch.

    Because of the coordinate system mapping(from geographic system to
    Catestian system), we allow an error when we check the location. The allow
    error will keep increasing after each loop when no solution in the two
    cases are found, until the solution get obtained.

    :param node0:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment.
    :param node1:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node2:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node3:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param hypocenter:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of floating hypocenter on each segment
        calculation. In the method, we take the direction point of the
        previous fault patch as hypocentre for the current fault patch.
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of reference point for projection
    :param pp:
        the projection of the site onto the plane containing the fault
        slipped area. A numpy array.
    :returns:
        Pd, a numpy array, representing the location of direction point
        E, the distance from direction point to hypocentre.
        go_next_patch, flag indicates if the calculation goes on the next
        fault patch. 1: yes, 0: no.