I posted a VBA function to return The angle between two vectors, in 2D or 3D last year, and have just discovered that Python and Numpy are lacking this function. Since all the suggested code I found in a quick search used:
Cos θ = (a.b)/(|a||b|)
which gives inaccurate results for small angles, I have written my own, using the same procedure as the VBA version:
Tan θ = |(axb)|/ (a.b)
Here is the lengthy code:
import numpy as np
import numpy.linalg as la
@xl_func("numpy_row v1, numpy_row v2: float")
def py_ang(v1, v2):
""" Returns the angle in radians between vectors 'v1' and 'v2' """
cosang = np.dot(v1, v2)
sinang = la.norm(np.cross(v1, v2))
return np.arctan2(sinang, cosang)
For details on how to link to Python functions from Excel, using Pyxll, see: Installing Python, Scipy and Pyxll ; also an updated Glob_to_Loc function, using the py_ang function, will appear within the next few days.
Newton Excel Bach, not (just) an Excel Blog 
Doug,
you might want to further “harden” the floating-point calculation against inaccuracies introduced in case the norms of the two vectors are very different (what follows is symbolic code, not py!!) by adding the following lines at the beginning of def py_ang:
n_v1=la.norm(v1)
n_v2=la.norm(v2)
if (abs(log10(n_v1/nv2)) > 10) then
v1 = v1/n_v1
v2 = v2/n_v2
endif
This can be done very fast by defining a function that pulls out the floating-point binary exponents only of the variables n_v1 and n_v2 and then takes their absolute difference, because we do not need the accurate log here. Of course, then the limit should be set to 30 or something like that. Hopefully, la.norm uses Kahan’s compensated summation algorithm, because otherwise this step would be useless… I’ll check that tomorrow…
LikeLike
Thanks, Georg, I’ll have a look at that.
LikeLike
Edit 27 Dec 2016: The MPMath links provided by Georg below no longer work. The official site for MPMath is now: http://mpmath.org/. Also see http://sympy.org/ and http://sagemath.org/ for the Sympy and Sage systems, which include MPMath.
Doug,
it is probably much more straightforward to install and harness mpmath:
http://mpmath.googlecode.com/svn/trunk/doc/build/index.html
than to implement all these deliberate forks in PY code yourself. MPmath is a multi-precision PY module. It provides all the trig-functions you need:
http://mpmath.googlecode.com/svn/trunk/doc/build/functions/trigonometric.html
in arbitrary floating-point precision. So you might want to set mp.dps to 40 (or so) and then use the conversion functions described here:
http://mpmath.googlecode.com/svn/trunk/doc/build/general.html
in order to firstly convert your DP input into MP numbers, secondly do all the intermediate calculations of the simple acos formula in high precision, and finally convert the result back to DP. You have to pay attention whether the calculations are really done in MP or only in DP with many decimal digits in the output… I’m very pleased with version 0.17; yet, the recently relased 0.18 seems to add interval arithmetic, which (theoretically) allows for giving strict error bounds of floating-point results, and thus would be a remarkable improvement. There are much more valuable features which you might find interesting as well.
Nevertheless, even simple PY provides a decent compensated summation algorithm via math.fsum . You can see the difference between sum and fsum when you try to sum up 0.1 ten times. Of course, you could also try to port the XBLAS routines:
http://www.netlib.org/xblas/
you need. Unfortunately, these double-double or quad-double packages do not run reliably in VBA due to a double-rounding issue that occurs with an odds ratio of about 1/4096. So you would need to compile them in Fortran or C and make sure they use the SSE registers instead of the x87 stack by adding the corresponding compiler directions. This might cause a considarable loss of performance, because there is no hardware support for a LOG function in SSE, for example, whereas the x87 instruction set provides it.
LikeLike
Georg – thanks for the MPMath link. I will have a look, but for my purposes with structural analysis (where the materials properties will be accurate to 2 SF, if you are lucky, and the loads about 1 SF) it really looks like overkill; and good performance is important (particularly for solving large systems of simultaneous equations).
I’m surprised Numpy or Scipy don’t have an angle function. Perhaps I should look at getting them added
LikeLike
Interesting, but what would you use it for?
LikeLike
Finding the 3D orientation of beams and plates for computer analysis is the main use I have in mind, but I’m sure there are others.
LikeLike
Yes but that is a direct vector process. Better than trig. About I only works that out last week.
LikeLike
Bill – I get your point (I think). The rotation matrix used to convert between global and local axes uses only sin(Gamma) and cos(Gamma), so there is no need to find Gamma, since the sin and cos are defined by purely vector operations.
On the other hand, I use an angle to define the rotation of a 2 node beam, and even where I am defining a plate or beam using 3 nodes, I want to be able to define an additional rotation with an angle, so calculating the angle between the local y axis and the Global horizontal plane in all cases seems the simplest way to handle it.
LikeLike
Doug
I posted a note on Quartenaries without going far enough down the page and find it uses the sins and cosines. Will ook again. 2D is straightforward I think.
Bill
LikeLike
Bill – yes 2D is simple, my code passes the Sin and Cos without calculating the angle, using DX/L and DY/L.
I might do the same for 3D (using the vector relationships) when the additional rotation angle is 0 (which it is most of the time).
LikeLike