The next step is to rotate a point in space using matrices, and then to explore the concept of gimbal lock. I will probably tackle gimbal lock in a future posting. Before delving into 3-D rotations, I needed to figure out how to extend my plotting knowledge.

I cleaned up the code from the previous posting, moving the spherical to Cartesian coordinates conversion code to a separate module -

**sphcoords**. sphcoords also includes a function that takes spherical coordinates and returns a 1x3 NumPy matrix that represents a point in space.

import numpy as np

from numpy import matrix

def x_cart(incl_sph_coords, azim_sph_coords):

return np.sin(incl_sph_coords) * np.cos(azim_sph_coords)

def y_cart(incl_sph_coords, azim_sph_coords):

return np.sin(incl_sph_coords) * np.sin(azim_sph_coords)

def z_cart(incl_sph_coords):

return np.cos(incl_sph_coords)

def sph2cart(incl_sph_coord, azim_sph_coord):

return np.matrix([

x_cart(incl_sph_coord, azim_sph_coord),

y_cart(incl_sph_coord, azim_sph_coord),

z_cart(incl_sph_coord)])

Below is an example of

**sphcoords**function,

**sph2cart**which generates a 1x3 matrix containing Cartesian coordinates. In this case, spherical coordinates of 30 degrees inclination and 180 degrees azimuth are converted.

`poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))`

By using matplotlib's subplot capability along with spherical to Cartesian coordinate conversions, I was able to depict a 3-D spherical model of space. A single red "point of interest" will be rotated using rotational matrices.

The initial position of the red point is arbitrary. The red point will be initially rotated about the z axis. The rotated position of the red point will remain along the point's original circle of inclination. This is what I was hoping to illustrate with the cyan-colored annotations.

Before going further, let me apologize for my inability to generate circles. The elliptical, squashed figures in my plots will be circular after I figure out how to set fixed aspect ratios.

The Initial Position of the Red Point of Interest |

**rotmat**using the NumPy Matrix.

```
import numpy as np
from numpy import matrix
def x_rot(rads):
return np.matrix([
[1,0,0],
[0, np.cos(rads), -np.sin(rads)],
[0, np.sin(rads), np.cos(rads)]])
def y_rot(rads):
return np.matrix([
[np.cos(rads), 0, np.sin(rads)],
[0, 1, 0],
[-np.sin(rads), 0, np.cos(rads)]])
def z_rot(rads):
return np.matrix([
[np.cos(rads), -np.sin(rads), 0],
[np.sin(rads), np.cos(rads), 0],
[0,0,1]])
```

` `

` `

The 1x3 matrix containing the point of interest's (POI) Cartesian coordinates is rotated about the z axis thusly:`poi = poi*rotmat.z_rot(-5*azi_incr)`

The plot below shows the POI rotated 90 degrees about the z axis. The cyan annotations show the range of the next rotation which is about the x axis. The POI will be rotated 30 degrees about the x axis. The most interesting update will be on the upper right plot that will show the the POI jumping from the yellow 30 degrees circle to the black circle at 60 degrees.

After Rotating the POI about the Z axis. |

After Rotating the POI about the X axis. |

` `

The final position of the POI after rotating about each Euler Axis. |

` `

` `

import numpy as np

import matplotlib.pylab as plt

from numpy import matrix

import rotmat

import sphcoords

def generate_plot(inclin, azi, poi):

# plot x-y plane

plt.subplot(221)

plt.plot(

sphcoords.x_cart(inclin[0],azi[0]),

sphcoords.y_cart(inclin[0],azi[0]), "g-o",

sphcoords.x_cart(inclin[1],azi[1]),

sphcoords.y_cart(inclin[1],azi[1]), "y-o",

sphcoords.x_cart(inclin[2],azi[2]),

sphcoords.y_cart(inclin[2],azi[2]), "y-o",

sphcoords.x_cart(inclin[3],azi[3]),

sphcoords.y_cart(inclin[3],azi[3]), "k-o",

sphcoords.x_cart(inclin[4],azi[4]),

sphcoords.y_cart(inclin[4],azi[4]), "k-o",

sphcoords.x_cart(inclin[5],azi[5]),

sphcoords.y_cart(inclin[5],azi[5]), "b-o",

sphcoords.x_cart(inclin[6],azi[6]),

sphcoords.y_cart(inclin[6],azi[6]), "b-o",

poi[0,0], poi[0,1], "r-o")

plt.grid(True)

plt.ylabel('Y Axis')

plt.xlabel('X Axis')

# plot y-z plane

plt.subplot(222)

plt.plot(

sphcoords.z_cart(inclin[0]),

sphcoords.y_cart(inclin[0],azi[0]), "g-o",

sphcoords.z_cart(inclin[1]),

sphcoords.y_cart(inclin[1],azi[1]), "y-o",

sphcoords.z_cart(inclin[2]),

sphcoords.y_cart(inclin[2],azi[2]), "y-o",

sphcoords.z_cart(inclin[3]),

sphcoords.y_cart(inclin[3],azi[3]), "k-o",

sphcoords.z_cart(inclin[4]),

sphcoords.y_cart(inclin[4],azi[4]), "k-o",

sphcoords.z_cart(inclin[5]),

sphcoords.y_cart(inclin[5],azi[6]), "b-o",

sphcoords.z_cart(inclin[6]),

sphcoords.y_cart(inclin[6],azi[6]), "b-o",

poi[0,2], poi[0,1], "r-o")

plt.xlabel('Z Axis')

plt.grid(True)

#plot x-z plane

plt.subplot(223)

plt.plot(

sphcoords.x_cart(inclin[0],azi[0]),

sphcoords.z_cart(inclin[0]),"g-o",

sphcoords.x_cart(inclin[1],azi[1]),

sphcoords.z_cart(inclin[1]),"y-o",

sphcoords.x_cart(inclin[2],azi[2]),

sphcoords.z_cart(inclin[2]),"y-o",

sphcoords.x_cart(inclin[3],azi[3]),

sphcoords.z_cart(inclin[3]),"k-o",

sphcoords.x_cart(inclin[4],azi[4]),

sphcoords.z_cart(inclin[4]),"k-o",

sphcoords.x_cart(inclin[5],azi[5]),

sphcoords.z_cart(inclin[5]),"b-o",

sphcoords.x_cart(inclin[6],azi[6]),

sphcoords.z_cart(inclin[6]),"b-o",

poi[0,0], poi[0,2], "r-o")

plt.ylabel('Z Axis')

plt.grid(True)

plt.show()

if __name__ == '__main__':

inclin = []

azi = []

inclin_points = range(0,20)

azi_incr = (2 * np.pi)/float(20)

azi_points = np.arange(-np.pi, np.pi, azi_incr)

#90 deg inclination

inclin.append([np.pi/2 for i in inclin_points])

#30,150 deg inclinations

inclin.append([np.pi/6 for i in inclin_points])

inclin.append([np.pi-np.pi/6 for i in inclin_points])

#60,120 deg inclinations

inclin.append([np.pi/3 for i in inclin_points])

inclin.append([np.pi-np.pi/3 for i in inclin_points])

# poles

inclin.append([0])

inclin.append([np.pi])

azi.append(azi_points)

azi.append(azi_points)

azi.append(azi_points)

azi.append(azi_points)

azi.append(azi_points)

azi.append([np.pi])

azi.append([np.pi])

poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))

poi = poi*rotmat.z_rot(-5*azi_incr)

poi = poi*rotmat.x_rot(np.pi/6)

poi = poi*rotmat.y_rot(np.pi/2)

generate_plot(inclin, azi, poi)

` `

## 5 comments:

Is this just fun and games or are you trying to achieve something with this? No offence, just curious.

I think this stuff is interesting enough to experiment with on my on time. The graphical nature of this work helps to clarify this issue at hand and makes it worthwhile to share with others. It is good to update the blog after a long hiatus. Over all, I think it may be fun to come up with Python, NumPy, matplotlib hacks.

It would be interesting to create some abstract spherical coordinate animations.

I was looking for something that was informative but less work than animation. In fact, if just looking at the numerical results providing meaningful information, I would have stopped there. The difference between raw numbers and the 2-D plots is night and day.

Also, showing each axis rotation discretely illustrates how the rotation of Euler angles/matrices works. Real world rotations occur about each axis simultaneously, this is where matrix rotation can fail due to gimbal lock.

Finally, I really would like to increase my Python, NumPy, and matplotlib skills. Most engineers coming out of college today know matlab and get a lot of use of it. I didn't get exposed to matlab but want that same type of power available to me. I hope to leverage and build upon my knowledge of Python. Thanks for stopping by.

You'll probably be interested in <a href="http://www.mathworks.com/discovery/quaternion.html>quaternion</a> math with MATLAB (or it's alternatives) if you start exploring gimbal lock.

Post a Comment