CODEWITHSUNDEEP
pythonarraysnumpyscipyinterpolation
maelstromscientist
maelstromscientist

Reputation: 551

Get y coordinates given x and z for a plane in 3D

I've generated a 3D circular plane that has been rotated along the x-axis by 45 degrees:

enter image description here

I want to determine the y-coordinate of the plane, given the x and z coordinates, but I can't figure out how to do that. How do I interpolate the plane such that I can get the y-coordinate if I feed it an x- and z-coordinate?

Here is my code:

def coord_rotation(theta):
    # Convert to radians
    theta_1_rad = theta[0] * np.pi/180.0
    theta_2_rad = theta[1] * np.pi/180.0
    theta_3_rad = theta[2] * np.pi/180.0
    # The bicone and dust angles correspond to Euler angles which are 
    # (e1,e2,e3) -> (rotation about z, rotation about x, rotation about z again)
    theta_1,theta_2,theta_3 = theta_1_rad,theta_2_rad,theta_3_rad
    R_x = np.array([[1,         0,                  0                   ],
                    [0,         np.cos(theta_1),   np.sin(theta_1)   ],
                    [0,         -np.sin(theta_1),  np.cos(theta_1)    ]
                    ])
    R_y = np.array([[np.cos(theta_2),    0,        -np.sin(theta_2)    ],
                    [0,                   1,        0                    ],
                    [np.sin(theta_2),    0,         np.cos(theta_2)    ]
                    ])
    R_z = np.array([[np.cos(theta_3),       np.sin(theta_3),        0],
                    [-np.sin(theta_3),      np.cos(theta_3),       0],
                    [0,                      0,                      1]
                    ])             
    R = np.dot(R_z, np.dot( R_y, R_x ))
    return R

theta_D1_deg  = -45.0)
theta_D3_deg  = 0.0  
D = 2
sampling = 25
########################################################################################
phi       = 2*np.pi # rotation 
phi    = np.linspace(0,phi,sampling)
r      = np.linspace(-D,D,sampling)

ri,pi = np.ix_(r,phi) # get open grids            
X = ri*np.cos(pi)
Y = ri*np.sin(pi)
Z = np.zeros(np.shape(X))
# Rotate the dust plane in 3d
t = np.transpose(np.array([X,Y,Z]), (1,2,0))
R = coord_rotation((theta_D1_deg,0,theta_D3_deg))
xd,yd,zd = np.transpose(np.dot(t, R), (2,0,1))

# Make uniform grid
points = (xd.ravel(),yd.ravel())
values = zd.ravel()
xdgrid,ydgrid = np.meshgrid(np.linspace(-2,2,1000),np.linspace(-2,2,1000))
zdgrid = griddata(points, values, (xdgrid, ydgrid), method='linear')

# Plot
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(1,1,1, projection='3d')
ax1.view_init(elev=15, azim=35)
ax1.plot_wireframe(xdgrid,ydgrid,zdgrid,alpha=0.25,color='xkcd:orange',zorder=3)

fontsize = 12
# x-axis
ax1.set_xlim(-2,2)
ax1.set_xlabel(r'$x$',fontsize=fontsize)
xAxisLine = ((np.min(xd), np.max(xd)), (0, 0), (0,0))
ax1.plot(xAxisLine[0], xAxisLine[1], xAxisLine[2], color='black',zorder=1,alpha=0.5)
# y-axis
ax1.set_ylim(-2,2)
ax1.set_ylabel(r'$y$',fontsize=fontsize)
yAxisLine = ((0, 0), (np.min(yd), np.max(yd)), (0,0))
ax1.plot(yAxisLine[0], yAxisLine[1], yAxisLine[2], color='black',zorder=1,alpha=0.5)
# z-axis
ax1.set_zlim(-2,2)
ax1.set_zlabel(r'$z$',fontsize=fontsize)
zAxisLine = ((0, 0), (0,0), (np.min(xd), np.max(xd)))
ax1.plot(zAxisLine[0], zAxisLine[1], zAxisLine[2], color='black',zorder=1,alpha=0.5)
plt.tight_layout()

Upvotes: 0

Views: 793

Answers (2)

neutrino_logic
neutrino_logic

Reputation: 1299

First, generate the ellipse that the rotated disk projects onto the x-z plane:

https://math.stackexchange.com/questions/2388747/formula-for-ellipse-formed-from-projecting-a-tilted-circle-onto-the-xy-plane

It's a bit easier for this case because the disk is tilted 45 degrees. The major axis of the ellipse is just the diameter of the original disk, and the minor axis is related to the diameter by diameter**2 = 2 * (minor**2).

Once you have the major and minor axis of this ellipse, it's not too hard to determine if the x,z components of the test point lie inside that ellipse, using this equation:

https://math.stackexchange.com/questions/76457/check-if-a-point-is-within-an-ellipse

If the test point lies within the disk, then the y-component on the disk is as noted in above answer (for the 45 degree tilt only, but it is just a right triangle): y = -z You've now found the point in x-y-z space that lies on the rotated disk, so just subtract that y-value from the test point's y-value to determine the distance to the disk along the y-axis.

Upvotes: 0

Prune
Prune

Reputation: 77860

This is straightforward 3D analytic geometry. First, note that there is no such thing as a "circular plane"; you have describe a circle and its interior, which, by definition, are embedding in a particular plane.

The equation of that plane is y + z = 0; x is an unconstrained variable, except as used to define the boundaries of the circle.

Thus, your problem reduces to simply

y = -z

Upvotes: 1

Related Questions