Calculates the GC2 coordinates for the nodes of the upper edge of the
        fault
507021
def _get_gc2_coordinates_for_rupture(self, edge_sets):
        

        # Establish GC2 length - for use with Ry0
        rup_gc2t, rup_gc2u = self.get_generalised_coordinates(
            edge_sets[:, 0], edge_sets[:, 1])
        # GC2 length should be the largest positive GC2 value of the edges
        self.gc_length = numpy.max(rup_gc2u)
281977
def nodes_touching_edge(edge_coord):
    
    a, b = hex_digit(edge_coord, 1), hex_digit(edge_coord, 2)
    if a % 2 == 0 and b % 2 == 0:
        return [coord_from_hex_digits(a, b + 1),
                coord_from_hex_digits(a + 1, b)]
    else:
        return [coord_from_hex_digits(a, b),
                coord_from_hex_digits(a + 1, b + 1)]
980332
function positionToGeographic(position) {
    var x = position[0];
    var y = position[1];
    return [radToDeg(x / EarthRadius) - (Math.floor((radToDeg(x / EarthRadius) + 180) / 360) * 360), radToDeg((Math.PI / 2) - (2 * Math.atan(Math.exp(-1.0 * y / EarthRadius))))];
  }
1181745
def zone_calc(surf, cortex, src):
    

    cortex_vertices, cortex_triangles = surf_keep_cortex(surf, cortex)

    dist_vals = np.zeros((len(source_nodes), len(cortex_vertices)))

    for x in range(len(source_nodes)):

        translated_source_nodes = translate_src(source_nodes[x], cortex)
        dist_vals[x, :] = gdist.compute_gdist(cortex_vertices, cortex_triangles, source_indices = translated_source_nodes)

    data = np.argsort(dist_vals, axis=0)[0, :] + 1

    zone = recort(data, surf, cortex)

    del data

    return zone
947792
def _get_ln_y_ref(self, rup, dists, C):
        
        # reverse faulting flag
        Frv = 1 if 30 <= rup.rake <= 150 else 0
        # normal faulting flag
        Fnm = 1 if -120 <= rup.rake <= -60 else 0
        # hanging wall flag
        Fhw = (dists.rx >= 0)
        # aftershock flag. always zero since we only consider main shock
        AS = 0

        ln_y_ref = (
            # first line of eq. 13a
            C['c1']
            + (C['c1a'] * Frv
               + C['c1b'] * Fnm
               + C['c7'] * (rup.ztor - 4))
            * (1 - AS)
            + (C['c10'] + C['c7a'] * (rup.ztor - 4)) * AS
            # second line
            + C['c2'] * (rup.mag - 6)
            + ((C['c2'] - C['c3']) / C['cn'])
            * np.log(1 + np.exp(C['cn'] * (C['cm'] - rup.mag)))
            # third line
            + C['c4']
            * np.log(dists.rrup
                     + C['c5']
                     * np.cosh(C['c6'] * max(rup.mag - C['chm'], 0)))
            # fourth line
            + (C['c4a'] - C['c4'])
            * np.log(np.sqrt(dists.rrup ** 2 + C['crb'] ** 2))
            # fifth line
            + (C['cg1'] + C['cg2'] / (np.cosh(max(rup.mag - C['cg3'], 0))))
            * dists.rrup
            # sixth line
            + C['c9'] * Fhw
            * np.tanh(dists.rx
                      * (np.cos(np.radians(rup.dip)) ** 2)
                      / C['c9a'])
            * (1 - np.sqrt(dists.rjb ** 2 + rup.ztor ** 2)
               / (dists.rrup + 0.001))
        )
        return ln_y_ref
989325
def B(u, v, dfs_data):
    """Bu(v) = {wx | w is in S*(v)}"""
    if a(v, dfs_data) != u:
        return None

    return list(set([edge_id for w in S_star(v, dfs_data) for edge_id in dfs_data['graph'].get_node(w)['edges']]))
1179117
def compute_coordinate(self, index):
        ''''''
        j = index%self.board_size
        i = (index - j) // self.board_size
        return (i, j)
680935
def _compute_imt1100(self, C, sites, rup, dists, get_pga_site=False):
        
        # Calculates simple site response term assuming all sites 1100 m/s
        fsite = (C['c10'] + (C['k2'] * C['n'])) * log(1100. / C['k1'])
        # Calculates the PGA on rock
        pga1100 = np.exp(self._compute_magnitude_term(C, rup.mag) +
                         self._compute_distance_term(C, rup, dists) +
                         self._compute_style_of_faulting_term(C, rup) +
                         self._compute_hanging_wall_term(C, rup, dists) +
                         self._compute_basin_response_term(C, sites.z2pt5) +
                         fsite)
        # If PGA at the site is needed then remove factor for rock and
        # re-calculate on correct site condition
        if get_pga_site:
            pga_site = np.exp(np.log(pga1100) - fsite)
            fsite = self._compute_shallow_site_response(C, sites, pga1100)
            pga_site = np.exp(np.log(pga_site) + fsite)
        else:
            pga_site = None
        return pga1100, pga_site
121192
def fix_double_nodes(coor, ngroups, conns, eps):
    
    n_nod, dim = coor.shape
    cmap = find_map( coor, nm.zeros( (0,dim) ), eps = eps, allow_double = True )
    if cmap.size:
        output('double nodes in input mesh!')
        output('trying to fix...')

        while cmap.size:
            print cmap.size

            # Just like in Variable.equation_mapping()...
            ii = nm.argsort( cmap[:,1] )
            scmap = cmap[ii]

            eq = nm.arange( n_nod )
            eq[scmap[:,1]] = -1
            eqi = eq[eq >= 0]
            eq[eqi] = nm.arange( eqi.shape[0] )
            remap = eq.copy()
            remap[scmap[:,1]] = eq[scmap[:,0]]
            print coor.shape
            coor = coor[eqi]
            ngroups = ngroups[eqi]
            print coor.shape
            ccs = []
            for conn in conns:
                ccs.append( remap[conn] )
            conns = ccs
            cmap = find_map( coor, nm.zeros( (0,dim) ), eps = eps,
                            allow_double = True )
        output('...done')
    return coor, ngroups, conns
1177902
def dfs_edges(G, start, depth_limit=1, get_only=True, get_path=False):
    """"""
    depth_limit = depth_limit - 1

    # creates unsigned int array (2 Byte)
    output_nodes = array('L')
    output_depth = array('I')
    # creates float array (4 Byte)
    output_weights = array('f')
    apath = []

    if G.node.get(start) is None:
        # raise KeyError('Start node not found')
        print('Start node not found')
        return output_nodes, output_weights, output_depth, apath

    visited = set()
    visited.add(start)

    # Save the start node with its data to the stack
    stack = [(start, G.edges_iter(start, data=True), 1.0)]
    visited.add(start)
    while stack:
        if len(output_nodes) > 80100100:
            print("To many nodes for: {}".format(start))
            del output_nodes
            del output_weights
            del output_depth
            output_nodes = array('L')
            output_depth = array('I')
            # creates float array (4 Byte)
            output_weights = array('f')
            gc.collect()
            break
        parent, children, weight = stack[-1]
        try:
            parent_, child, child_keys = next(children)
            # print "child: {}, parent_data: {}".format(child, parent_data)
            if child not in visited:
                weight = child_keys.get('weight', 1.0) * weight
                visited.add(child)

                if len(stack) >= depth_limit or weight <= 0.00001:
                    visited.remove(child)
                else:
                    stack.append((child, G.edges_iter(child, data=True),
                                 weight))

                # if its not and user.
                if get_only and child > 100000000000:
                    # if get_only and G.node[child].get('Type') != get_only:
                    continue

                output_nodes.append(child)
                output_weights.append(weight)
                output_depth.append(len(stack))
                if get_path:
                    apath.append([step[0] for step in stack])
        except StopIteration:
            stack.pop()
            visited.remove(parent)
            # if data.get('Type') == "Node":

    return output_nodes, output_weights, output_depth, apath
749203
def _validate_simple_fault_geometry(self, node, _float_re):
        
        try:
            # Parse the geometry
            coords = split_coords_2d(~node.LineString.posList)
            trace = geo.Line([geo.Point(*p) for p in coords])
        except ValueError:
            # If the geometry cannot be created then use the LogicTreeError
            # to point the user to the incorrect node. Hence, if trace is
            # compiled successfully then len(trace) is True, otherwise it is
            # False
            trace = []
        if len(trace):
            return
        raise LogicTreeError(
            node, self.filename,
            "'simpleFaultGeometry' node is not valid")
422627
public static void computePositions(Node node, int depth, double aoff, double awid, int maxdepth) {
    double r = depth / (maxdepth - 1.);
    node.x = FastMath.sin(aoff + awid * .5) * r;
    node.y = FastMath.cos(aoff + awid * .5) * r;

    double cpos = aoff;
    double cwid = awid / node.weight;
    for(Node c : node.children) {
      computePositions(c, depth + 1, cpos, cwid * c.weight, maxdepth);
      cpos += cwid * c.weight;
    }
  }
544424
def setxy(self):
        
        self._edge_inverter()
        self._detect_alignment_conflicts()
        inf = float('infinity')
        # initialize vertex coordinates attributes:
        for l in self.layers:
            for v in l:
                self.grx[v].root  = v
                self.grx[v].align = v
                self.grx[v].sink  = v
                self.grx[v].shift = inf
                self.grx[v].X     = None
                self.grx[v].x     = [0.0]*4
        curvh = self.dirvh # save current dirvh value
        for dirvh in xrange(4):
            self.dirvh = dirvh
            self._coord_vertical_alignment()
            self._coord_horizontal_compact()
        self.dirvh = curvh # restore it
        # vertical coordinate assigment of all nodes:
        Y = 0
        for l in self.layers:
            dY = max([v.view.h/2. for v in l])
            for v in l:
                vx = sorted(self.grx[v].x)
                # mean of the 2 medians out of the 4 x-coord computed above:
                avgm = (vx[1]+vx[2])/2.
                # final xy-coordinates :
                v.view.xy = (avgm,Y+dY)
            Y += 2*dY+self.yspace
        self._edge_inverter()
61979
def chk_edges(self):
        """"""
        goids = set(self.go2obj)
        self.chk_edges_nodes(self.edges, goids, "is_a")
        for reltype, edges in self.edges_rel.items():
            self.chk_edges_nodes(edges, goids, reltype)
1363003
def build_multi_point_source_node(multi_point_source):
    
    # parse geometry
    pos = []
    for p in multi_point_source.mesh:
        pos.append(p.x)
        pos.append(p.y)
    mesh_node = Node('gml:posList', text=pos)
    upper_depth_node = Node(
        "upperSeismoDepth", text=multi_point_source.upper_seismogenic_depth)
    lower_depth_node = Node(
        "lowerSeismoDepth", text=multi_point_source.lower_seismogenic_depth)
    source_nodes = [Node(
        "multiPointGeometry",
        nodes=[mesh_node, upper_depth_node, lower_depth_node])]
    # parse common distributed attributes
    source_nodes.extend(get_distributed_seismicity_source_nodes(
        multi_point_source))
    return Node("multiPointSource",
                get_source_attributes(multi_point_source),
                nodes=source_nodes)
599254
def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej):
    
    # The input argument are the dot products
    #
    #   <e1, e2>
    #   <e2, e0>
    #   <e0, e1>
    #
    # of the edges
    #
    #   e0: x1->x2,
    #   e1: x2->x0,
    #   e2: x0->x1.
    #
    # Note that edge e_i is opposite of node i and the edges add up to 0.

    # The trilinear coordinates of the circumcenter are
    #
    #   cos(alpha0) : cos(alpha1) : cos(alpha2)
    #
    # where alpha_k is the angle at point k, opposite of edge k. The Cartesian
    # coordinates are (see
    # <https://en.wikipedia.org/wiki/Trilinear_coordinates#Between_Cartesian_and_trilinear_coordinates>)
    #
    #     C = sum_i ||e_i|| * cos(alpha_i)/beta * P_i
    #
    # with
    #
    #     beta = sum ||e_i||*cos(alpha_i)
    #
    # Incidentally, the cosines are
    #
    #    cos(alpha0) = <e1, e2> / ||e1|| / ||e2||,
    #
    # so in total
    #
    #    C = <e_0, e0> <e1, e2> / sum_i (<e_i, e_i> <e{i+1}, e{i+2}>) P0
    #      + ... P1
    #      + ... P2.
    #
    # An even nicer formula is given on
    # <https://en.wikipedia.org/wiki/Circumscribed_circle#Barycentric_coordinates>: The
    # barycentric coordinates of the circumcenter are
    #
    #   a^2 (b^2 + c^2 - a^2) : b^2 (c^2 + a^2 - b^2) : c^2 (a^2 + b^2 - c^2).
    #
    # This is only using the squared edge lengths, too!
    #
    alpha = ei_dot_ei * ei_dot_ej
    alpha_sum = alpha[0] + alpha[1] + alpha[2]
    beta = alpha / alpha_sum[None]
    a = X * beta[..., None]
    cc = a[0] + a[1] + a[2]

    # alpha = numpy.array([
    #     ei_dot_ei[0] * (ei_dot_ei[1] + ei_dot_ei[2] - ei_dot_ei[0]),
    #     ei_dot_ei[1] * (ei_dot_ei[2] + ei_dot_ei[0] - ei_dot_ei[1]),
    #     ei_dot_ei[2] * (ei_dot_ei[0] + ei_dot_ei[1] - ei_dot_ei[2]),
    # ])
    # alpha /= numpy.sum(alpha, axis=0)
    # cc = (X[0].T * alpha[0] + X[1].T * alpha[1] + X[2].T * alpha[2]).T
    return cc
993797
def sites_at_edges( self ):
        
        min_x = min( [ s.r[0] for s in self.sites ] )
        max_x = max( [ s.r[0] for s in self.sites ] )
        min_y = min( [ s.r[1] for s in self.sites ] )
        max_y = max( [ s.r[1] for s in self.sites ] )
        min_z = min( [ s.r[2] for s in self.sites ] )
        max_z = max( [ s.r[2] for s in self.sites ] )
        x_max = [ s for s in self.sites if s.r[0] == min_x ]
        x_min = [ s for s in self.sites if s.r[0] == max_x ]
        y_max = [ s for s in self.sites if s.r[1] == min_y ]
        y_min = [ s for s in self.sites if s.r[1] == max_y ]
        z_max = [ s for s in self.sites if s.r[2] == min_z ]
        z_min = [ s for s in self.sites if s.r[2] == max_z ]
        return ( x_max, x_min, y_max, y_min, z_max, z_min )