CODEWITHSUNDEEP
pythonnumpymathlinear-algebra
user115657
user115657

Reputation: 19

Find planes orthogonal to a known plane

I was hoping to get some help here. I'm a geologist so my math background isn't as strong as it should be. I have a geological fault that I've digitized in a CAD program. I've placed three representative points on the fault plane and extracted those coordinate points. I've found the unit vector of that plane using the following code:

p1 = np.array([503562, 2811504, 1771], dtype='float64')
p2 = np.array([504122, 2810766, 1820], dtype='float64')
p3 = np.array([504325, 2811311, 1042], dtype='float64')

#  Calculate vectors
v1 = p3 - p1
v2 = p2 - p1

#  Find the cross product (vector normal to the plane)
cp = np.cross(v1, v2)

#  Assign to a, b and c for components of plane equation
a, b, c = cp

#  Calculate d component with a dot product and the following equation:
#    a*x3 + b*y3 + c*z3 = d
d = np.dot(cp, p3)

#  Find the unit vector
uv = cp / (cp**2).sum()**0.5

This vector represents a plane which is my maximum direction of continuity (a geological term). I need to get the two planes which are at 90 degree angles to the plane defined above. I know this is probably quite simple in terms of linear algebra but I'm at a loss right now. I eventually need this to get the strike and dip of the three planes of continuity defined above.

Edit: added a picture for clarity

The blue feature is a fault. The green rugby ball is my search ellipse. I need to figure out the planes which are defined by that shape. I used those green points digitized on the fault to calculate the first plane. I don't know how to get the other two

fault and search ellipse

Upvotes: 1

Views: 747

Answers (2)

jberrio
jberrio

Reputation: 1124

From your description, this is what I believe you want:

  • For the strike, you want the intersection of your plane and a horizontal plane (so cross product of your unit vector and a (0,0,1) vector.

  • For the dip, you want the cross product of your plane's unit vector and the strike's unit vector.

I agree with other commenters - maybe sort out the math on your own... then you can come to stackoverflow if you have problems implementing it in Python :)

Upvotes: 2

Julien
Julien

Reputation: 15249

So assuming you are happy with any such planes:

You already have cp defining your "continuity plane". ortho1 = np.cross(v1, cp) will give you one orthogonal direction to cp i.e one of your possible requested planes. Then ortho2 = np.cross(cp, ortho1) will give you the second corresponding one.

Upvotes: 3

Related Questions