3. Lengths and angles

probe measures its edges: their lengths, their directions, the angles between them, and the area of a face from a cross product.

program probe
!< Tutorial, chapter 3: lengths, directions and angles of the edges.
use penf, only : R8P
use vecfor, only : angle, ex, ey, ez, normL2, vector
implicit none
type(vector) :: o           ! the origin
type(vector) :: p(4)        ! the vertices
type(vector) :: a, b        ! two edges from p2
type(vector) :: u           ! a unit vector
type(vector) :: n           ! a normal
real(R8P)    :: deg         ! degrees per radian

p(1) = o ; p(2) = ex ; p(3) = ey ; p(4) = ez
deg = 180._R8P / acos(-1._R8P)
a = p(3) - p(2)
b = p(4) - p(2)
print '(A,F8.4)', 'length of p2->p3      ', a%normL2()
print '(A,F8.4)', 'squared length        ', a%sq_norm()
print '(A,F8.4)', 'length, free function ', normL2(b)
u = a%normalized()                    ! a unit vector along a; a is unchanged
call u%printf(prefix='direction of p2->p3 = ')
print '(A,F8.4)', 'its length            ', u%normL2()
call a%normalize                      ! a itself becomes a unit vector
call a%printf(prefix='a, normalized       = ')
a = p(3) - p(2)
print '(A,F8.4)', 'a . b                 ', a .dot. b
print '(A,F8.4)', 'p1->p2 . p1->p3       ', (p(2) - p(1)) .dot. (p(3) - p(1))
print '(A,F8.4)', 'angle at p2 [deg]     ', a%angle(b) * deg
print '(A,F8.4)', 'angle at p1 [deg]     ', angle(p(2) - p(1), p(3) - p(1)) * deg
n = a .cross. b                       ! orthogonal to a and b, as long as the area of their parallelogram
call n%printf(prefix='a x b = ')
print '(A,F8.4)', 'area of p2 p3 p4      ', 0.5_R8P * n%normL2()
endprogram probe

$ probe
length of p2->p3        1.4142
squared length          2.0000
length, free function   1.4142
direction of p2->p3 = -0.7071067811865475, +0.7071067811865475, +0.0
its length              1.0000
a, normalized       = -0.7071067811865475, +0.7071067811865475, +0.0
a . b                   1.0000
p1->p2 . p1->p3         0.0000
angle at p2 [deg]      60.0000
angle at p1 [deg]      90.0000
a x b = +1.0, +1.0, +1.0
area of p2 p3 p4        0.8660

Lengths

a%normL2() is the Euclidean length, $\sqrt{x^2+y^2+z^2}$; a%sq_norm() its square, $x^2+y^2+z^2$, cheaper (no square root) and exact for exact components: compare squared lengths when you only need to know which is longer.

Every method of a vector that returns a value is also a free function, with the vector as first argument: normL2(b) is b%normL2(). The free spelling works on any expression, the type-bound one only on a variable ((a - b)%normL2() is not Fortran); see A method or a free function.

Directions

a%normalized() returns a unit vector along a and leaves a alone; call a%normalize makes a itself a unit vector. The direction of the edge from p2 to p3 is $(-1,1,0)/\sqrt{2}$: -0.7071067811865475 is the double precision number nearest to $1/\sqrt{2}$, and the length of the result is 1 to the last digit. The zero vector has no direction: normalizing it gives the zero vector back, not a division by zero (see A unit vector). To be told when a vector is too short to have a direction, call normalized_checked(a, u, ok): ok is false and u is zero for such a vector (see A unit vector, or a report that there is none).

Dot product

a .dot. b is $a_x{b}_x+a_y{b}_y+a_z{b}_z$, a number: 1 for the two edges from p2, 0 for the edges from p1 along x and y, which are orthogonal.

Angles

a%angle(b), or angle(a, b), is the angle between the two vectors, in radians, between 0 and π: 60 degrees at p2, where the slanted face is equilateral, 90 degrees at p1. VecFor computes it with atan2, from the cross and the dot products of the unit vectors: accurate for nearly parallel and nearly opposite vectors too, where the textbook $\arccos (\hat{a}\cdot \hat{b})$ loses half of its digits.

Cross product

a .cross. b is the vector orthogonal to a and to b, oriented by the right-hand rule, as long as the area of the parallelogram they span: half of it is the area of the triangle p2 p3 p4, $\sqrt{3}/2=0.8660$.

.cross. and .dot. round each product before they subtract or add them. When the terms nearly cancel, as for nearly parallel vectors, the result loses digits; crossproduct_accurate(a, b) and dotproduct_accurate(a, b) keep them (see Products of nearly parallel vectors, accurately).

What you learned

normL2, sq_norm, normalized, normalize, .dot., angle, .cross.; every method that returns a value is also a free function. In passing: normalized_checked, crossproduct_accurate, dotproduct_accurate. Reference: Operators, Geometry.

Next: 4. Faces, area, volume.