Source Code for Module lib.geometry.vectors

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004-2005,2008-2010,2013-2014 Edward d'Auvergne               # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Collection of functions for vector operations.""" 
 24   
 25  # Python module imports. 
 26  from math import acos, atan2, cos, pi, sin 
 27  from numpy import array, cross, dot, float64, sqrt 
 28  from numpy.linalg import norm 
 29  from random import uniform 
 30   
 31   
 32 -def complex_inner_product(v1=None, v2_conj=None): 
 33      """Calculate the inner product <v1|v2> for the two complex vectors v1 and v2. 
 34   
 35      This is calculated as:: 
 36                     ___ 
 37                     \ 
 38          <v1|v2> =   >   v1i . v2i* , 
 39                     /__ 
 40                      i 
 41   
 42      where * is the complex conjugate. 
 43   
 44   
 45      @keyword v1:        The first vector. 
 46      @type v1:           numpy rank-1 complex array 
 47      @keyword v2_conj:   The conjugate of the second vector.  This is already in conjugate form to allow for non-standard definitions of the conjugate (for example Sm* = (-1)^m S-m). 
 48      @type v2_conj:      numpy rank-1 complex array 
 49      @return:            The value of the inner product <v1|v2>. 
 50      @rtype:             float 
 51      """ 
 52   
 53      # Return the inner product. 
 54      return dot(v1, v2_conj) 
 55   
 56   
 57 -def random_unit_vector(vector): 
 58      """Generate a random rotation axis. 
 59   
 60      Uniform point sampling on a unit sphere is used to generate a random axis orientation. 
 61   
 62      @param vector:  The 3D rotation axis. 
 63      @type vector:   numpy 3D, rank-1 array 
 64      """ 
 65   
 66      # Random azimuthal angle. 
 67      u = uniform(0, 1) 
 68      theta = 2*pi*u 
 69   
 70      # Random polar angle. 
 71      v = uniform(0, 1) 
 72      phi = acos(2.0*v - 1) 
 73   
 74      # Random unit vector. 
 75      vector[0] = cos(theta) * sin(phi) 
 76      vector[1] = sin(theta) * sin(phi) 
 77      vector[2] = cos(phi) 
 78   
 79   
 80 -def unit_vector_from_2point(point1, point2): 
 81      """Generate the unit vector connecting point 1 to point 2. 
 82   
 83      @param point1:  The first point. 
 84      @type point1:   list of float or numpy array 
 85      @param point2:  The second point. 
 86      @type point2:   list of float or numpy array 
 87      @return:        The unit vector. 
 88      @rtype:         numpy float64 array 
 89      """ 
 90   
 91      # Convert to numpy data structures. 
 92      point1 = array(point1, float64) 
 93      point2 = array(point2, float64) 
 94   
 95      # The vector. 
 96      vect = point2 - point1 
 97   
 98      # Return the unit vector. 
 99      return vect / norm(vect) 
100   
101   
102 -def vector_angle_acos(vector1, vector2): 
103      """Calculate the angle between two N-dimensional vectors using the acos formula. 
104   
105      The formula is:: 
106   
107          angle = acos(dot(a / norm(a), b / norm(b))). 
108   
109   
110      @param vector1:     The first vector. 
111      @type vector1:      numpy rank-1 array 
112      @param vector2:     The second vector. 
113      @type vector2:      numpy rank-1 array 
114      @return:            The angle between 0 and pi. 
115      @rtype:             float 
116      """ 
117   
118      # Calculate and return the angle. 
119      return acos(dot(vector1 / norm(vector1), vector2 / norm(vector2))) 
120   
121   
122 -def vector_angle_atan2(vector1, vector2): 
123      """Calculate the angle between two N-dimensional vectors using the atan2 formula. 
124   
125      The formula is:: 
126   
127          angle = atan2(norm(cross(a, b)), dot(a, b)). 
128   
129      This is more numerically stable for angles close to 0 or pi than the acos() formula. 
130   
131   
132      @param vector1:     The first vector. 
133      @type vector1:      numpy rank-1 array 
134      @param vector2:     The second vector. 
135      @type vector2:      numpy rank-1 array 
136      @return:            The angle between 0 and pi. 
137      @rtype:             float 
138      """ 
139   
140      # Calculate and return the angle. 
141      return atan2(norm(cross(vector1, vector2)), dot(vector1, vector2)) 
142   
143   
144 -def vector_angle_complex_conjugate(v1=None, v2=None, v1_conj=None, v2_conj=None): 
145      """Calculate the inter-vector angle between two complex vectors using the arccos formula. 
146   
147      The formula is:: 
148   
149          theta = arccos(Re(<v1|v2>) / (|v1|.|v2|)) , 
150   
151      where:: 
152                     ___ 
153                     \  
154          <v1|v2> =   >   v1i . v2i* , 
155                     /__ 
156                      i 
157   
158      and:: 
159   
160          |v1| = Re(<v1|v1>) . 
161   
162   
163      @keyword v1:        The first vector. 
164      @type v1:           numpy rank-1 complex array 
165      @keyword v2:        The second vector. 
166      @type v2:           numpy rank-1 complex array 
167      @keyword v1_conj:   The conjugate of the first vector.  This is already in conjugate form to allow for non-standard definitions of the conjugate (for example Sm* = (-1)^m S-m). 
168      @type v1_conj:      numpy rank-1 complex array 
169      @keyword v2_conj:   The conjugate of the second vector.  This is already in conjugate form to allow for non-standard definitions of the conjugate (for example Sm* = (-1)^m S-m). 
170      @type v2_conj:      numpy rank-1 complex array 
171      @return:            The angle between 0 and pi. 
172      @rtype:             float 
173      """ 
174   
175      # The inner products. 
176      inner_v1v2 = complex_inner_product(v1=v1, v2_conj=v2_conj) 
177      inner_v1v1 = complex_inner_product(v1=v1, v2_conj=v1_conj) 
178      inner_v2v2 = complex_inner_product(v1=v2, v2_conj=v2_conj) 
179   
180      # The normalised inner product, correcting for truncation errors creating ratios slightly above 1.0. 
181      ratio = inner_v1v2.real / (sqrt(inner_v1v1).real * sqrt(inner_v2v2).real) 
182      if ratio > 1.0: 
183          ratio = 1.0 
184   
185      # Calculate and return the angle. 
186      return acos(ratio) 
187   
188   
189 -def vector_angle_normal(vector1, vector2, normal): 
190      """Calculate the directional angle between two N-dimensional vectors. 
191   
192      @param vector1:     The first vector. 
193      @type vector1:      numpy rank-1 array 
194      @param vector2:     The second vector. 
195      @type vector2:      numpy rank-1 array 
196      @param normal:      The vector defining the plane, to determine the sign. 
197      @type normal:       numpy rank-1 array 
198      @return:            The angle between -pi and pi. 
199      @rtype:             float 
200      """ 
201   
202      # Normalise the vectors (without changing the original vectors). 
203      vector1 = vector1 / norm(vector1) 
204      vector2 = vector2 / norm(vector2) 
205   
206      # The cross product. 
207      cp = cross(vector1, vector2) 
208   
209      # The angle. 
210      angle = acos(dot(vector1, vector2)) 
211      if dot(cp, normal) < 0.0: 
212          angle = -angle 
213   
214      # Return the signed angle. 
215      return angle 
216