myprojects-hg-reborn

Fork

#!/usr/bin/env python
import scipy as s
# import pylab as p
import numpy as n
import sys
import string
import scipy.linalg as sl
from mayavi import mlab


def random_rot():
    theta=s.arccos(1.-2.*s.random.uniform(0.,1.,1)[0])-s.pi/2.
    phi=s.random.uniform(-s.pi,s.pi,1)[0]
    psi=s.random.uniform(-s.pi,s.pi,1)[0]


    oneone=s.cos(theta)*s.cos(psi)
    onetwo=-s.cos(phi)*s.sin(psi)+s.sin(phi)*s.sin(theta)*s.cos(psi)
    onethree=s.sin(phi)*s.sin(psi)+s.cos(phi)*s.sin(theta)*s.cos(psi)

    twoone= s.cos(theta)*s.sin(psi)
    twotwo=s.cos(phi)*s.cos(psi)+s.sin(phi)*s.sin(theta)*s.sin(psi)
    twothree=-s.sin(phi)*s.cos(psi)+s.cos(phi)*s.sin(theta)*s.sin(psi)

    threeone=-s.sin(theta)
    threetwo=s.sin(phi)*s.cos(theta)
    threethree=s.cos(phi)*s.cos(theta)


    my_mat=s.zeros(9).reshape((3,3))

    my_mat[0,0]=oneone
    my_mat[0,1]=onetwo
    my_mat[0,2]=onethree

    my_mat[1,0]=twoone
    my_mat[1,1]=twotwo
    my_mat[1,2]=twothree

    my_mat[2,0]=threeone
    my_mat[2,1]=threetwo
    my_mat[2,2]=threethree

    return my_mat


def rotate_cluster(cluster):
    random_rot_mat= random_rot()
    n_row_col=s.shape(cluster)

    cluster_rot=s.zeros(s.prod(n_row_col)).reshape((n_row_col[0],\
                                                    n_row_col[1]))

    for i in s.arange(n_row_col[0]):
                
        cluster_rot[i,:]=s.dot(random_rot_mat, cluster[i,:])

    return cluster_rot

#############################################################

N_rot=4000

visualisation=1

ini_cluster=s.arange(3).reshape((1,3))*1.


ini_cluster[0,0]=1.
ini_cluster[0,1]=0.
ini_cluster[0,2]=0.


print "ini_cluster is, ", ini_cluster

rotated_clust=rotate_cluster(ini_cluster)

print "rotated_clust is, ", rotated_clust

# now iterate the rotations and save the position of the rotated vectors

rot_vec=s.copy(rotated_clust)

for i in xrange(N_rot):
    rotated_clust=rotate_cluster(ini_cluster)
    rot_vec=s.vstack((rot_vec,rotated_clust))
    
n.savetxt("rotated_vectors.dat", rot_vec)

if (visualisation==1):


    x=rot_vec[:,0]
    y=rot_vec[:,1]
    z=rot_vec[:,2]

    mlab.clf()
    pts=mlab.points3d(x, y, z,color=(0,0,1), scale_factor=.02)
    pts2=mlab.points3d(0, 0, 0,color=(1,0,0), scale_factor=2.,\
                       resolution=40)
    mlab.draw()
    mlab.show()




print "So far so good"