ARL-TN-0759 • MAY 2016

US Army Research Laboratory

Resolving the Orientation of Cylinders and

Cuboids from Projected Area Measurements

by Richard Saucier

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement

or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-TN-0759 • MAY 2016

US Army Research Laboratory

Resolving the Orientation of Cylinders and

Cuboids from Projected Area Measurements

by Richard Saucier

Survivability/Lethality Analysis Directorate, ARL

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704‐0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the

data needed, and completing and reviewing the collection information.  Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing

the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704‐0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202‐

4302.  Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a

currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD‐MM‐YYYY) 2. REPORT TYPE 3. DATES COVERED (From ‐ To)

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION REPORT

NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

13. SUPPLEMENTARY NOTES

14. ABSTRACT

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF:   17. LIMITATION

       OF  

       ABSTRACT

18. NUMBER

       OF  

        PAGES

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

May 2016 Final

Resolving the Orientation of Cylinders and Cuboids from Projected Area

Measurements

Richard Saucier

ARL-TN-0759

Approved for public release; distribution is unlimited.

November 2015–January 2016

AH80

US Army Research Laboratory

ATTN: RDRL-SLB-S

Aberdeen Proving Ground, MD 21005-5068

The FATEPEN model predicts the penetration of a mass striking a target plate for a variety of shapes, including cylinders and

cuboids—among others. Crucial to the use of the model is a good estimate of not only the mass and velocity but also the impact

orientation in terms of pitch, yaw, and roll. Yaw cards and orthogonal X-rays can provide estimates of the impact projected area,

but the model makes use of the impact angle, which is defined to be the minimum angle that a face makes with the target plate.

This report addresses this issue by calculating the pitch-yaw-roll rotation sequence that will bring about the orientation at impact

from orthogonal projected area measurements. It is shown that the impact angle is uniquely determined in the case of cylinders,

but that is not the case for cuboids. Furthermore, for cuboids multiple impact angles for the same projected area measurements can

lead to significantly different FATEPEN predictions. This leads to the conclusion that cylinders and not cuboids should be used for

FATEPEN validation.

FATEPEN, RCC, RPP, fragment impact, fragment orientation, projected area, orthogonal X-rays, pitch-yaw-roll

44

Richard Saucier

Unclassified Unclassified Unclassified UU 410-278-6721

ii
Approved for public release; distribution is unlimited.

Contents

List of Figures v

List of Tables vi

List of Listings vii

Acknowledgments viii

1. Summary 1

2. Introduction 1

3. Methods, Assumptions, and Procedures 2

3.1 Orientation of a Cylinder 3

3.1.1 Projected Areas from Cylinder Orientation 3

3.1.2 Cylinder Orientation from Projected Areas 4

3.2 Orientation of a Cuboid 6

3.2.1 Projected Areas from Cuboid Orientation 7

3.2.2 Cuboid Orientation from Projected Areas 9

4. Results and Discussion 14

4.1 Cylinder Orientation Sample Case 14

4.2 Cuboid Orientation Sample Case 16

4.2.1 Impact Angle and Effective Yaw Angle 17

5. Conclusions 19

6. References 21

Appendix A. Formula for an Euler Sequence of Rotations 23

Appendix B. Yaw Angle of a Cylinder as a Function of Shape Factor 29

List of Symbols, Abbreviations, and Acronyms 33

iii
Approved for public release; distribution is unlimited.

Distribution List 34

iv
Approved for public release; distribution is unlimited.

List of Figures

Fig. 1 Description of pitch, yaw, and roll. ................................................2

Fig. 2 An oriented cylinder with projected area on orthogonal planes ........... 16

Fig. 3 An oriented cuboid with projected area on orthogonal planes ............. 16

v
Approved for public release; distribution is unlimited.

List of Tables

Table 1 Pitch, yaw, and roll rotations of unit vectors ..................................3

Table 2 FATEPEN predictions for a 725-gr steel cuboid with a striking velocity of

1500 f/s impacting a 0.25-inch mild steel plate using version 3.3.10.3. In

each case, the projected areas on the orthogonal planes are the same. Note,

in particular, the sensitivity of the residual velocity, vr. ................... 17

Table 3 FATEPEN predictions for a 725-gr steel cuboid with a striking velocity of

1500 f/s impacting a 0.25-inch mild steel plate using version 3.3.18.0. In

each case, the projected areas on the orthogonal planes are the same. . 19

vi
Approved for public release; distribution is unlimited.

List of Listings

Listing 1 orient.cpp ............................................................ 12

Listing 2 rcc.cpp................................................................ 14

Listing A-1 reverse.cpp ........................................................... 27

Listing B-1 cyl.cpp................................................................ 31

vii
Approved for public release; distribution is unlimited.

Acknowledgments

I would like to thank Timothy Mallory for his interest in this problem, suggestions

that improved this report, and the numerous stand-alone FATEPEN runs he conducted

to check the results. I would also like to thank John Auten for catching some errors

and his technical review of this document. Of course, I accept full responsibility for

any errors that may remain.

viii
Approved for public release; distribution is unlimited.

1. Summary

The FATEPEN model predicts the penetration of a mass striking a target plate for

a variety of shapes, including cylinders and cuboids—among others. Crucial to

the use of the model is a good estimate of not only the mass and velocity but also

the impact orientation in terms of pitch, yaw, and roll. Yaw cards and orthogonal

X-rays can provide estimates of the impact projected area, but the model makes use

of the impact angle, which is defined to be the minimum angle that a penetrator

face makes with the target plate. This report addresses this issue by calculating the

pitch-yaw-roll rotation sequence that will bring about the orientation at impact from

orthogonal projected area measurements. It is shown that the impact angle is uniquely

determined in the case of cylinders, but that is not the case for cuboids. Furthermore,

for cuboids multiple impact angles for the same projected area measurements can

lead to significantly different FATEPEN predictions. This leads to the conclusion

that cylinders and not cuboids should be used for FATEPEN validation.

2. Introduction

The FATEPEN1 model is a Fortran code for simulating penetration of fragments,

long rods, and projectiles into target plates and predicting penetration, perforation,

residual mass, and residual velocity. It provides the following shapes: sphere, right

circular cylinder (RCC), round-nose cylinder, sharp-nose cylinder, tapered/truncated

cylinder, and rectangular parallelepiped (RPP). This report is limited to RCCs, which

we will simply call cylinders, and RPPs, which we will call cuboids. Unlike the

THOR2 model, which only requires an impact presented area, the FATEPEN model

requires a specific shape and orientation.

The impact presented area, or projected area, is a key parameter that affects the

penetration process since the penetration depth scales with the mass per unit presented area for a given striking velocity. But the FATEPEN model goes beyond

the presented area and also requires the angle that each flat face of the striking

mass makes with the target plate. This means that we need to account for both the

presented area and the orientation of the fragment to use FATEPEN properly.

1
Approved for public release; distribution is unlimited.

3. Methods, Assumptions, and Procedures

We first define a standard orientation of the cylinder or cuboid as the starting

orientation and then describe the operational procedure to give it a specific final

orientation. This standard orientation will be with the center at the origin of a righthanded Cartesian (x, y, z) coordinate system and the length of the cylinder or cuboid

aligned with the z-axis. (The target plate is taken to be parallel to the x-y plane, and

the fragment velocity vector is along the negative z axis.) The operational procedure

will be in terms pitch, yaw, and roll, where pitch is a counterclockwise rotation

about the x-axis, yaw is a counterclockwise rotation about the y-axis, and roll is a

counterclockwise rotation about the z-axis.∗ For convenience, we also define ˆı, ˆ,

and kˆ to be unit vectors along the x, y, and z axis, respectively, as shown in Fig. 1.

x

y

z

ˆı

ˆ

kˆ

yaw

pitch

roll

Fig. 1. Description of pitch, yaw, and roll

Pitch, yaw, and roll are then described by the operators Rˆı(φp), Rˆ(φy), and Rkˆ(φr),

where the unit vector subscript denotes the axis of rotation, and φp, φy, and φr are

pitch, yaw, and roll angles, respectively. A rotation sequence accounts for the fact

that the unit vectors are also transformed by prior rotations. Thus, a pitch-yaw-roll

rotation sequence is described by

R = Rkˆ00(φr)Rˆ0(φy)Rˆı(φp), (1)

where the operators are applied from right to left—first pitch, then yaw, then roll—

and whereˆ0is the transformed ˆ vector after the pitch rotation, and kˆ00 is the doubly

transformed kˆ vector after the pitch and yaw rotations. Now we make use of the

fact that we generate the same net rotation if we apply the sequence in reverse order

∗To avoid any confusion, note that FATEPEN uses a pitch-yaw-roll rotation sequence, which we

follow in this report, but much of the aerospace industry uses a yaw-pitch-roll rotation sequence.

2
Approved for public release; distribution is unlimited.

about the fixed axes.∗ This means that

R = Rkˆ00(φr)Rˆ0(φy)Rˆı(φp) = Rˆı(φp)Rˆ(φy)Rkˆ(φr) (2)

and we only have to concern ourselves with rotations about fixed axes. Table 1 lists

all 3 × 3 = 9 combinations of rotations applied to the unit vectors.

Table 1. Pitch, yaw, and roll rotations of unit vectors

Rotation ˆı ˆ kˆ

Pitch: Rˆı(φp) ˆı cos φpˆ + sin φpkˆ − sin φpˆ + cos φpkˆ

Yaw: Rˆ(φy) cos φyˆı − sin φykˆ ˆ sin φyˆı + cos φykˆ

Roll: Rkˆ(φr) cos φrˆı + sin φrˆ − sin φrˆı + cos φrˆ kˆ

3.1 Orientation of a Cylinder

The cylinder has a length L and diameter D. In standard orientation the center of the

cylinder is at the origin and the length is aligned with the z-axis.

3.1.1 Projected Areas from Cylinder Orientation

So let us apply a pitch-yaw-roll rotation sequence to the cylinder in standard orientation and work out the projected areas on the orthogonal y-z, x-z, and x-y planes.

The orientation of the cylinder is completely described by its axis of rotation, which

is along the unit vector kˆ. Let ˆu be the final orientation of the kˆ vector. Then, using

Table 1, we get

ˆu = Rˆı(φp)Rˆ(φy)kˆ

= Rˆı(φp)(sin φyˆı + cos φykˆ)

= sin φyˆı + cos φy(− sin φpˆ + cos φpkˆ). (3)

It is also convenient to specify ˆu in terms of the direction cosines:

ˆu = cos φ1ˆı + cos φ2ˆ + cos φ3kˆ, (4)

∗See Appendix A for a derivation of this result.

3
Approved for public release; distribution is unlimited.

where cos φ1 = ˆu ·ˆı, cos φ2 = ˆu ·ˆ, and cos φ3 = ˆu · kˆ. Comparing Eqs. 3 and 4

gives the correspondence

φ1 = cos−1(sin φy),

φ2 = cos−1(− sin φp cos φy), (5)

φ3 = cos−1(cos φp cos φy).

Let L be the length of the cylinder and D be its diameter. Then the projected areas

onto the orthogonal planes are

Ayz = LD| cos π

2

− φ1



| +

π

4

D

2

| cos φ1| = LD|sin φ1| +

π

4

D

2

| cos φ1|

Axz = LD| cos π

2

− φ2



| +

π

4

D

2

| cos φ2| = LD|sin φ2| +

π

4

D

2

| cos φ2| . (6)

Axy = LD| cos π

2

− φ3



| +

π

4

D

2

| cos φ3| = LD|sin φ3| +

π

4

D

2

| cos φ3|

Therefore, given the pitch and yaw of the cylinder, Eqs. 5 and 6 give us the projected

areas.

3.1.2 Cylinder Orientation from Projected Areas

Next, we work out the inverse problem of determining the pitch and yaw from the

measured projected areas. Thus, given the length L, diameter D, and the projected

areas of an RCC onto the orthogonal planes, we want to solve for the orientation of

the RCC. More specifically, we seek the FATEPEN pitch and yaw rotation sequence

that will bring about this orientation.∗It is convenient in what follows to specify the

final orientation in terms of the polar angle, θ, measured from the z-axis and φ, the

azimuthal angle measured from the x-axis in the x-y plane:

ˆu = sin θ cos φˆı + sin θ sin φˆ + cos θkˆ. (7)

There is no loss of generality if we restrict ˆu to lie in the first octant of the unit

sphere, so that

0 ≤ θ ≤ π/2 and 0 ≤ φ ≤ π/2. (8)

∗Ordinarily we would need pitch, yaw, and roll to completely specify orientation, but placing the

axis of the RCC along the z (roll) axis eliminates the need to consider roll.

4
Approved for public release; distribution is unlimited.

Comparing Eqs. 4 and 7 gives

θ = φ3 and φ = tan−1



cos φ2

cos φ1



. (9)

All the formulas in Eq. 6 are of the form

a sin α + b cos α = Ap, (10)

with α = φ1, φ2, or φ3, a = LD, b = πD2/4, and Ap the projected area, Ayz, Axz,

or Axy. We can solve this for the angle α by first introducing another angle β, defined

by

β ≡ tan−1(b/a), (11)

which gives a =

√

a

2 + b2 cos β and b =

√

a

2 + b2 sin β, and then Eq. 10 becomes

√

a

2 + b2 cos β sin α +

√

a

2 + b2 sin β cos α =

√

a

2 + b2 sin(α + β) = Ap. (12)

Now it is easy to show that √a

2 + b2 = Amax, the maximum presented area of the

RCC, so that 0 ≤ Ap/

√

a

2 + b2 ≤ 1, and therefore

sin(α + β) = Ap

√

a

2 + b2

. (13)

Two possible solutions to this equation are

α = sin−1



Ap

√

a

2 + b2



− tan−1



b

a



or

α = π − tan−1



b

a



− sin−1



Ap

√

a

2 + b2



. (14)

In practice we try both, so this gives 2

3 = 8 possible combinations for φ1, φ2, and

φ3. But by imposing the constraint that

cos2 φ1 + cos2 φ2 + cos2 φ3 = 1, (15)

we will find that only one combination will work. Once φ1, φ2, and φ3 have thus

been found, we return to Eq. 9 to get θ and φ.

5
Approved for public release; distribution is unlimited.

Next we need a rotation that will align kˆ with ˆu. It is easy to check that Rˆe(θ), where

ˆe =

kˆ × ˆu

|kˆ × ˆu|

= − sin φˆı + cos φˆ, (16)

is this rotation. That is, Rˆe(θ)kˆ = ˆu. The quaternion representation of this rotation

is

qˆe(θ) = cos(θ/2) + ˆe sin(θ/2) = cos(θ/2) + (− sin φˆı + cos φˆ) sin(θ/2). (17)

Finally, this quaternion can be factored into a pitch-yaw-roll rotation sequence.3 The

prescription is as follows:

1. Set p0 = cos(θ/2), p1 = − sin φ sin(θ/2), p2 = cos φ sin(θ/2), p3 = 0.

2. Set A = p1p2 − p0p3, B = p

2

1 − p

2

3

, D = p

2

0 − p

2

2

.

3. Then φr = tan−1[−2A/(B + D)].

4. Set c0 = cos(φr/2) and c3 = sin(φr/2).

5. Set q0 = p0c0 + p3c3, q1 = p1c0 − p2c3, q2 = p2c0 + p1c3, q3 = p3c0 − p0c3.

6. Then φp = 2 tan−1(q1/q0) and φy = 2 tan−1(q2/q0).

Thus, we accomplish what we set out to do: obtain the pitch and yaw rotation

sequence that will give the desired projected areas. For a cylinder there are actually

4 orientations that give the same projected areas, which are (φp, φy), (φp, −φy),

(−φp, φy), and (−φp, −φy).

3.2 Orientation of a Cuboid

The cuboid, or RPP, has a length L, width W, and thickness T, where L ≥ W ≥ T,

and is initially oriented so that the length is along the z-axis, the width is along the

x-axis and the thickness is along the y-axis. Thus, the initial areas along each of the

axes are Ax = LT, Ay = LW, and Az = W T. We takeˆı to be a unit vector along

the x-axis, ˆ to be a unit vector along the y-axis, and kˆ to be a unit vector along the

z-axis. Let

ˆu = xˆı + yˆ + zkˆ with x

2 + y2 + z2 = 1 (18)

6
Approved for public release; distribution is unlimited.

be a unit vector on the unit sphere. Then the projected area of the cuboid orthogonal

to ˆu is

Ap(x, y, z) = (Axˆı + Ayˆ + Azkˆ) · ˆu = Axx + Ayy + Azz. (19)

As we vary the direction of ˆu, the rate of change of the projected area is given by

DˆuAp = ∇Ap(x, y, z) · ˆu

= ||∇Ap(x, y, z)|| ||ˆu|| cos θ

= ||∇Ap(x, y, z)|| cos θ, (20)

where θ is the angle between the gradient and the unit vector. From this expression,

it is clear that the area is maximized when cos θ is a maximum, which occurs when

θ = 0. The gradient along this direction is

∇Ap(x, y, z) = ∂Ap

∂x ˆı +

∂Ap

∂y ˆ +

∂Ap

∂z kˆ = Axˆı + Ayˆ + Azkˆ, (21)

and therefore the magnitude of the maximum projected area is

Amax = ||∇Ap(x, y, z)|| =

q

A2

x + A2y + A2z

, (22)

and the view direction for this maximum is given by the unit vector

ˆumax =

Axˆı + Ayˆ + Azkˆ

Amax

. (23)

By construction, the minimum projected area is initially oriented along the z-axis,

Amin = Az, and is realized along the unit vector ˆumin = kˆ.

3.2.1 Projected Areas from Cuboid Orientation

Now let us work out the projected areas from a given pitch-yaw-roll rotation sequence, R. Again, making use of Eq. 2 and applying the rotations from Table 1

successively gives

Rˆ(φy)Rkˆ(φr)Axˆı = Ax cos φr(cos φyˆı − sin φykˆ) + Ax sin φrˆ,

Rˆ(φy)Rkˆ(φr)Ayˆ = −Ay sin φr(cos φyˆı − sin φykˆ) + Ay cos φrˆ, (24)

Rˆ(φy)Rkˆ(φr)Azkˆ = Az(sin φyˆı + cos φykˆ).

7
Approved for public release; distribution is unlimited.

Collecting terms we get

Rˆ(φy)Rkˆ(φr)Axˆı = Ax cos φr cos φyˆı + Ax sin φrˆ − Ax cos φr sin φykˆ,

Rˆ(φy)Rkˆ(φr)Ayˆ = −Ay sin φr cos φyˆı + Ay cos φrˆ + Ay sin φr sin φykˆ, (25)

Rˆ(φy)Rkˆ(φr)Azkˆ = Az sin φyˆı + Az cos φykˆ.

Finally, applying pitch, we get

RAxˆı = + Ax cos φr cos φyˆı

+ Ax sin φr(cos φpˆ + sin φpkˆ)

− Ax cos φr sin φy(− sin φpˆ + cos φpkˆ),

RAyˆ = − Ay sin φr cos φyˆı (26)

+ Ay cos φr(cos φpˆ + sin φpkˆ)

+ Ay sin φr sin φy(− sin φpˆ + cos φpkˆ),

RAzkˆ = + Az sin φyˆı

+ Az cos φy(− sin φpˆ + cos φpkˆ).

and collecting terms,

RAxˆı = + Ax cos φr cos φyˆı

+ Ax(sin φr cos φp + cos φr sin φy sin φp)ˆ

+ Ax(sin φr sin φp − cos φr sin φy cos φp)kˆ, (27)

RAyˆ = − Ay sin φr cos φyˆı

+ Ay(cos φr cos φp − sin φr sin φy sin φp)ˆ

+ Ay(cos φr sin φp + sin φr sin φy cos φp)kˆ, (28)

RAzkˆ = + Az sin φyˆı − Az cos φy sin φpˆ + Az cos φy cos φpkˆ. (29)

Now let Axy be the projection on the x-y plane, Axz be the projection on the x-z

plane, and Ayz be the projection on the y-z plane. Then we have

Axy = [R(Axˆı + Ayˆ + Azkˆ)] · kˆ, (30)

Axz = [R(Axˆı + Ayˆ + Azkˆ)] ·ˆ, (31)

Ayz = [R(Axˆı + Ayˆ + Azkˆ)] ·ˆı. (32)

8
Approved for public release; distribution is unlimited.

Three is the maximum number of sides that can project unto a plane for a cuboid, so

we can include all the contributions by using the absolute value. Thus, we get the

projected areas

Axy =Ax|sin φr sin φp − cos φr sin φy cos φp|+

Ay| cos φr sin φp + sin φr sin φy cos φp|+

Az| cos φy cos φp|, (33)

Axz =Ax|sin φr cos φp + cos φr sin φy sin φp|+

Ay| cos φr cos φp − sin φr sin φy sin φp|+

Az| cos φy sin φp|, (34)

Ayz =Ax| cos φr cos φy| + Ay|sin φr cos φy| + Az|sin φy|. (35)

3.2.2 Cuboid Orientation from Projected Areas

Next, given measured values for Axy, Axz, and Ayz, we want to solve Eqs. 33–35

for φp, φy, and φr. This is a system of 3 simultaneous nonlinear equations, which

can be solved numerically by Newton’s method in 3 dimensions. Define

f(φp, φy, φr) ≡Ax|sin φr sin φp − cos φr sin φy cos φp|+

Ay| cos φr sin φp + sin φr sin φy cos φp|+

Az| cos φy cos φp| − Axy, (36)

g(φp, φy, φr) ≡Ax|sin φr cos φp + cos φr sin φy sin φp|+

Ay| cos φr cos φp − sin φr sin φy sin φp|+

Az| cos φy sin φp| − Axz, (37)

h(φp, φy, φr) ≡ Ax| cos φr cos φy| + Ay|sin φr cos φy| + Az|sin φy| − Ayz. (38)

Then the problem is to find the roots of f, g, and h. Let A be the matrix of partial

derivatives

A =







fφp(φp, φy, φr) fφy(φp, φy, φr) fφr(φp, φy, φr)

gφp(φp, φy, φr) gφy(φp, φy, φr) gφr(φp, φy, φr)

hφp(φp, φy, φr) hφy(φp, φy, φr) hφr(φp, φy, φr)





 (39)

9
Approved for public release; distribution is unlimited.

where

fφp(φp, φy, φr) = ∂f(φp, φy, φr)

∂φp

(40)

and so on for the other partial derivatives. For this we need the derivative of the

absolute value:

d

dx|u| =

d

dx

√

u

2 =

1

2

(u

2

)

−1/2

2uu0 =

uu0

√

u

2

=

uu0

|u|

= sgn(u) u

0

, (41)

where sgn(x) is the sign function:

sgn(x) =







+1 if x > 0

0 if x = 0 .

−1 if x < 0

(42)

Then from Eqs. 36–38 we get

fφp(φp, φy, φr) = +Ax sgn(sin φr sin φp − cos φr sin φy cos φp)

(sin φr cos φp + cos φr sin φy sin φp)

+Ay sgn(cos φr sin φp + sin φr sin φy cos φp)

(cos φr cos φp − sin φr sin φy sin φp)

+Az sgn(cos φy cos φp)(− cos φy sin φp), (43)

fφy(φp, φy, φr) = +Ax sgn(sin φr sin φp − cos φr sin φy cos φp)

(− cos φr cos φy sin φp)

+Ay sgn(cos φr sin φp + sin φr sin φy cos φp)

(sin φr cos φy cos φp)

+Az sgn(cos φy cos φp)(− sin φy cos φp), (44)

fφr

(φp, φy, φr) = +Ax sgn(sin φr sin φp − cos φr sin φy cos φp)

(cos φr sin φp + sin φr sin φy cos φp)

+Ay sgn(cos φr sin φp + sin φr sin φy cos φp)

(− sin φr sin φp + cos φr sin φy cos φp), (45)

10
Approved for public release; distribution is unlimited.

gφp

(φp, φy, φr) = +Ax sgn(sin φr cos φp + cos φr sin φy sin φp)

(− sin φr sin φp + cos φr sin φy cos φp)

+Ay sgn(cos φr cos φp − sin φr sin φy sin φp)

(− cos φr sin φp − sin φr sin φy cos φp)

+Az sgn(cos φy sin φp)(cos φy cos φp), (46)

gφy(φp, φy, φr) = +Ax sgn(sin φr cos φp + cos φr sin φy sin φp)

(cos φr cos φy sin φp)

+Ay sgn(cos φr cos φp − sin φr sin φy sin φp)

(− sin φr cos φy sin φp)

+Az sgn(cos φy sin φp)(− sin φy sin φp), (47)

gφr(φp, φy, φr) = +Ax sgn(sin φr cos φp + cos φr sin φy sin φp)

(cos φr cos φp − sin φr sin φy sin φp)

+Ay sgn(cos φr cos φp − sin φr sin φy sin φp)

(− sin φr cos φp − cos φr sin φy sin φp), (48)

hφp(φp, φy, φr) = 0, (49)

hφy(φp, φy, φr) = +Ax sgn(cos φr cos φy)(− cos φr sin φy)+

+Ay sgn(sin φr cos φy)(− sin φr sin φy)+

+Az sgn(sin φy)(cos φy), (50)

hφr(φp, φy, φr) = +Ax sgn(cos φr cos φy)(− sin φr cos φy)

+Ay sgn(sin φr cos φy)(cos φr cos φy). (51)

Newton’s method for solving this 3-dimensional (3D) problem numerically for pitch,

yaw, and roll is expressed by the matrix iteration equation







δφp

δφy

δφr







≡







φp,n+1 − φp,n

φy,n+1 − φy,n

φr,n+1 − φr,n





 = A

−1







f(φp, φy, φr)

g(φp, φy, φr)

h(φp, φy, φr)





 , (52)

where the matrix A is given by Eq. 39. We pick a starting orientation and continue

iterating until either the δ’s fall below a preset tolerance value  or we exceed a

maximum number of iterations—in which case we try another starting orientation.

11
Approved for public release; distribution is unlimited.

If we write this coefficient matrix (Eq. 39) as

A =







a11 a12 a13

a21 a22 a23

a31 a32 a33





 , (53)

then the inverse is given by

A

−1 =

1

det A







a22a33 − a23a32 a13a32 − a12a33 a12a23 − a13a22

a23a31 − a21a33 a11a33 − a13a31 a13a21 − a11a23

a21a32 − a22a31 a12a31 − a11a32 a11a22 − a12a21





 . (54)

Thus, we now have all the ingredients for implementing the proposed solution with

the program in Listing 1.

Listing 1. orient.cpp

1 // orient.cpp: Given the dimensions of an RPP and projected areas onto three orthogonal planes,

2 // finds a pitch-yaw-roll rotation sequence that will orient the RPP for these projected areas.

3 // The orientation is not unique and will depend upon the initial starting point.

4 // Method of solution makes use of Newton’s method in three dimensions.

5 // R. Saucier, October 2015

6

7 #include <iostream>

8 #include <cstdlib>

9 #include <cmath>

10 #include <iomanip>

11 #include <chrono>

12 #include <random>

13

14 inline double sgn( double x ) {

15

16 if ( x > 0. ) return +1.;

17 else if ( x < 0. ) return -1.;

18 else return 0.;

19 }

20

21 int main( int argc, char* argv[] ) {

22

23 const double D2R = M_PI / 180.; // to convert deg to rad

24 const double R2D = 180. / M_PI; // to convert rad to deg

25 const double TOL = 1.e-9; // convergence criterion

26 const int N = 100; // max number of iterations

27

28 double L = 3.; // length

29 double W = 2.; // width

30 double T = 1.; // thickness

31

32 double Ax = L * T; // area of side (initially along i)

33 double Ay = L * W; // area of top (initially along j)

34 double Az = W * T; // area of front (initially along k)

35

36 double A_xy = 4.5; // projected area on x-y plane

37 double A_xz = 6.5; // projected area on x-z plane

38 double A_yz = 3.5; // projected area on y-z plane

39

40 double p, y, r, cp, sp, cy, sy, cr, sr, f, g, h;

41 double e1, e2, e3, e4, e5, e6, e7, e8, e9;

42 double s1, s2, s3, s4, s5, s6, s7, s8, s9;

43 double fp, fy, fr, gp, gy, gr, hp, hy, hr;

44 double a11, a12, a13, a21, a22, a23, a31, a32, a33, det;

45 double b11, b12, b13, b21, b22, b23, b31, b32, b33;

46 double del_p, del_y, del_r;

47

48 // initial estimate can’t be zero

49 unsigned int seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();

50 std::mt19937 rng( seed ); // Mersenne Twister engine

12
Approved for public release; distribution is unlimited.

51 std::uniform_real_distribution<double> uniform( 0., M_PI ); // uniform distribution

52 p = uniform( rng );

53 y = uniform( rng );

54 r = uniform( rng );

55

56 if ( argc == 4 ) { // optionally specify initial pitch, yaw, roll on commandline

57

58 p = atof( argv[1] ) * D2R;

59 y = atof( argv[2] ) * D2R;

60 r = atof( argv[3] ) * D2R;

61 }

62

63 for ( int i = 0; i < N; i++ ) {

64

65 cp = cos( p ); sp = sin( p );

66 cy = cos( y ); sy = sin( y );

67 cr = cos( r ); sr = sin( r );

68

69 e1 = sr * sp - cr * sy * cp;

70 e2 = cr * sp + sr * sy * cp;

71 e3 = cy * cp;

72 e4 = sr * cp + cr * sy * sp;

73 e5 = cr * cp - sr * sy * sp;

74 e6 = cy * sp;

75 e7 = cr * cy;

76 e8 = sr * cy;

77 e9 = sy;

78

79 s1 = sgn( e1 );

80 s2 = sgn( e2 );

81 s3 = sgn( e3 );

82 s4 = sgn( e4 );

83 s5 = sgn( e5 );

84 s6 = sgn( e6 );

85 s7 = sgn( e7 );

86 s8 = sgn( e8 );

87 s9 = sgn( e9 );

88

89 f = Ax * fabs( e1 ) + Ay * fabs( e2 ) + Az * fabs( e3 ) - A_xy;

90 g = Ax * fabs( e4 ) + Ay * fabs( e5 ) + Az * fabs( e6 ) - A_xz;

91 h = Ax * fabs( e7 ) + Ay * fabs( e8 ) + Az * fabs( e9 ) - A_yz;

92

93 fp = Ax * s1 * ( sr * cp + cr * sy * sp ) + Ay * s2 * ( cr * cp - sr * sy * sp ) + Az * s3 * ( -cy * sp );

94 fy = Ax * s1 * ( -cr * cy * sp ) + Ay * s2 * ( sr * cy * cp ) + Az * s3 * ( -sy * cp );

95 fr = Ax * s1 * ( cr * sp + sr * sy * cp ) + Ay * s2 * ( -sr * sp + cr * sy * cp );

96

97 gp = Ax * s4 * ( -sr * sp + cr * sy * cp ) + Ay * s5 * ( -cr * sp - sr * sy * cp ) + Az * s6 * ( cy * cp );

98 gy = Ax * s4 * ( cr * cy * sp ) + Ay * s5 * ( -sr * cy * sp ) + Az * s6 * ( -sy * sp );

99 gr = Ax * s4 * ( cr * cp - sr * sy * sp ) + Ay * s5 * ( -sr * cp - cr * sy * sp );

100

101 hp = 0.;

102 hy = Ax * s7 * ( -cr * sy ) + Ay * s8 * ( -sr * sy ) + Az * s9 * ( cy );

103 hr = Ax * s7 * ( -sr * cy ) + Ay * s8 * ( cr * cy );

104

105 a11 = fp; a12 = fy; a13 = fr;

106 a21 = gp; a22 = gy; a23 = gr;

107 a31 = hp; a32 = hy; a33 = hr;

108 det = a11 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * a32 - a22 * a31 );

109 if ( det == 0 ) {

110

111 std::cerr << "bad initial orientation: " << p * R2D << "\t" << y * R2D << "\t" << r * R2D << std::endl

112 << "program " << argv[0] << " stopped" << std::endl;

113 exit( EXIT_FAILURE );

114 }

115

116 b11 = ( a22 * a33 - a23 * a32 ) / det;

117 b12 = ( a13 * a32 - a12 * a33 ) / det;

118 b13 = ( a12 * a23 - a13 * a22 ) / det;

119 b21 = ( a23 * a31 - a21 * a33 ) / det;

120 b22 = ( a11 * a33 - a13 * a31 ) / det;

121 b23 = ( a13 * a21 - a11 * a23 ) / det;

122 b31 = ( a21 * a32 - a22 * a31 ) / det;

123 b32 = ( a12 * a31 - a11 * a32 ) / det;

124 b33 = ( a11 * a22 - a12 * a21 ) / det;

125

126 del_p = b11 * f + b12 * g + b13 * h;

127 del_y = b21 * f + b22 * g + b23 * h;

128 del_r = b31 * f + b32 * g + b33 * h;

129

130 p -= del_p;

131 y -= del_y;

132 r -= del_r;

133

134 if ( fabs( del_p ) < TOL && fabs( del_y ) < TOL && fabs( del_r ) < TOL ) break;

135

13
Approved for public release; distribution is unlimited.

136 if ( i == N-1 ) {

137

138 std::cerr << "failed to converge after " << N << " iterations: f = " << f << " g = " << g << " h = " << h << std

::endl

139 << "try another initial orientation" << std::endl

140 << "program " << argv[0] << " stopped" << std::endl;

141 exit( EXIT_FAILURE );

142 }

143 }

144 std::cout << std::setprecision(9) << std::fixed;

145 std::cout << "pitch = " << fmod( p * R2D, 360. ) << std::endl

146 << "yaw = " << fmod( y * R2D, 360. ) << std::endl

147 << "roll = " << fmod( r * R2D, 360. ) << std::endl;

148

149 return EXIT_SUCCESS;

150 }

The program may be compiled and run with the following commands:

1 c++ -O2 -Wall -std=c++11 -o orient orient.cpp -lm

2 ./orient

4. Results and Discussion

We now have

• a method for computing the projected areas on orthogonal planes, given the

cylinder or cuboid orientation, and

• a method for determining the pitch-yaw-roll rotation sequence to bring about

the orientation, given the projected area measurements.

This means that we have a way to check the results. We will show that projected area

measurements lead to a unique orientation for a cylinder, but that is not the case for

a cuboid.

4.1 Cylinder Orientation Sample Case

Consider an RCC with L = 1, D = 1, φp = 30◦, and φy = 15◦. The projected areas

are found to be Ayz = 1.16920, Axz = 1.25496, and Axy = 1.20494. Using these

values for the projected areas, the program in Listing 2 computes φp = ±30◦and

φy = ±15◦, with the 4 possible orientations displayed in Fig. 2.

Listing 2. rcc.cpp

1 // rcc.cpp: compute pitch and yaw from area projections for an RCC (sample case)

2 // R. Saucier, June 2006 (Revised October 2015)

3

4 #include "Rotation.h"

5 #include <iostream>

6 #include <cmath>

7 #include <cstdlib>

8 #include <cassert>

9 #include <iomanip>

10 using namespace std;

11

12 double angle1( double a, double b, double c ) { return asin( c / sqrt( a * a + b * b ) ) - atan( b / a ); }

13 double angle2( double a, double b, double c ) { return M_PI - asin( c / sqrt( a * a + b * b ) ) - atan( b / a ); }

14

15 int main( void ) {

14
Approved for public release; distribution is unlimited.

16

17 const double L = 1., D = 1.;

18 const double A = L * D;

19 const double B = 0.25 * M_PI * D * D;

20 const double A_MIN = A < B ? A : B;

21 const double A_MAX = sqrt( A * A + B * B );

22

23 cout << setprecision(6) << fixed;

24 cout << "L = " << L << endl;

25 cout << "D = " << D << endl;

26 cout << "A_MIN = " << A_MIN << endl;

27 cout << "A_MAX = " << A_MAX << endl;

28

29 double A_xy = 1.20494; // projected area on x-y plane

30 double A_xz = 1.25496; // projected area on x-z plane

31 double A_yz = 1.16920; // projected area on y-z plane

32

33 cout << "A_xy = " << A_xy << endl;

34 cout << "A_xz = " << A_xz << endl;

35 cout << "A_yz = " << A_yz << endl;

36

37 assert( A_MIN <= A_xy && A_xy <= A_MAX );

38 assert( A_MIN <= A_yz && A_yz <= A_MAX );

39 assert( A_MIN <= A_xz && A_xz <= A_MAX );

40

41 double phi_x[ 2 ], phi_y[ 2 ], phi_z[ 2 ];

42 phi_x[0] = angle1( A, B, A_yz );

43 phi_x[1] = angle2( A, B, A_yz );

44

45 phi_y[0] = angle1( A, B, A_xz );

46 phi_y[1] = angle2( A, B, A_xz );

47

48 phi_z[0] = angle1( A, B, A_xy );

49 phi_z[1] = angle2( A, B, A_xy );

50

51 double d, delta = 1.e36;

52 int ii = 0, jj = 0, kk = 0;

53 for ( int i = 0; i < 2; i++ ) {

54 for ( int j = 0; j < 2; j++ ) {

55 for ( int k = 0; k < 2 ; k++ ) {

56

57 d = cos( phi_x[i] ) * cos( phi_x[i] ) +

58 cos( phi_y[j] ) * cos( phi_y[j] ) +

59 cos( phi_z[k] ) * cos( phi_z[k] );

60 if ( fabs( d - 1. ) < delta ) {

61 ii = i;

62 jj = j;

63 kk = k;

64 delta = fabs( d - 1. );

65 }

66 }

67 }

68 }

69

70 double th = phi_z[kk];

71 double ph = atan( cos( phi_y[jj] ) / cos( phi_x[ii] ) );

72

73 cout << "Derived Angles (deg): " << endl

74 << " phi_x: " << phi_x[ii] * va::R2D << endl

75 << " phi_y: " << phi_y[jj] * va::R2D << endl

76 << " phi_z: " << phi_z[kk] * va::R2D << endl

77 << " theta: " << th * va::R2D << endl

78 << " phi: " << ph * va::R2D << endl;

79

80 va::Vector u = sin( th ) * cos( ph ) * va::Vector( 1., 0., 0. ) +

81 sin( th ) * sin( ph ) * va::Vector( 0., 1., 0. ) +

82 cos( th ) * va::Vector( 0., 0., 1. );

83

84 va::Rotation R( va::Vector( 0., 0., 1. ), u ); // find rotation that takes k to u

85 va::sequence s = factor( R, va::XYZ );

86

87 cout << "pyr = " << +s.first * va::R2D << "\t" << +s.second * va::R2D << "\t" << s.third * va::R2D << endl;

88 cout << "pyr = " << +s.first * va::R2D << "\t" << -s.second * va::R2D << "\t" << s.third * va::R2D << endl;

89 cout << "pyr = " << -s.first * va::R2D << "\t" << +s.second * va::R2D << "\t" << s.third * va::R2D << endl;

90 cout << "pyr = " << -s.first * va::R2D << "\t" << -s.second * va::R2D << "\t" << s.third * va::R2D << endl;

91

92 cout << "Total effective yaw = " << acos( cos( s.first ) * cos( s.second ) ) * va::R2D << endl;

93

94 return EXIT_SUCCESS;

95 }

15
Approved for public release; distribution is unlimited.

Fig. 2. An oriented cylinder for φp = ±30◦and φy = ±15◦is shown with its projected area in the

x-y plane in dark blue, the x-z plane in green, and the y-z plane in red. All of these have the

same projected areas on each of the orthogonal planes. The actual orientations are not quite

the same, but these are superficial differences since the effective yaw angle and the impact

angle as defined by FATEPEN are all the same (33.2259◦in this example) and thus give the

same results for plate impact.

4.2 Cuboid Orientation Sample Case

Consider an RPP with L = 3 cm, W = 2 cm, and T = 1 cm, and let the projected

areas be Axy = 4.5 cm2, Axz = 6.5 cm2, Ayz = 3.5 cm2. When this information is

given to the orient.cpp program (Listing 1), it finds 4 different orientations that will

work, as displayed in Fig. 3, demonstrating that the orientation is not unique.

Fig. 3. The projected area on the orthogonal planes does not uniquely determine the cuboid orientation. The oriented cuboid (L = 3 cm, W = 2 cm, T = 1 cm) is shown with its projected area

in the x-y plane in dark blue (4.5 cm2), the x-z plane in green (6.5 cm2), and the y-z plane in

red (3.5 cm2). All 4 have the same projected areas on each of the orthogonal planes, although

we see that the actual orientations are not the same. One consequence of this is that the impact

angle, μ, which is the minimum angle that a face makes with the x-y plane, along with the

effective yaw angle, φeff, have different values for the 4 orientations—and this will lead to

different predictions from FATEPEN.

Next, we want to see what difference, if any, these different orientations will make

in FATEPEN. The cuboid penetrator is taken to be steel 4140 with a Brinell hardness

16
Approved for public release; distribution is unlimited.

of 300. It has a length of 3 cm = 1.1811 inches, width of 2 cm = 0.7874 inch, and

thickness of 1 cm = 0.3937 inch. Using a steel density of 7.83 g/cm3, it has a mass of

725 gr. The plate target is taken to be 0.25-inch steel with a Brinell hardness of 150.

All impacts were at 0◦ obliquity. The cuboid is initially aligned so that its length

is along the z-axis, its width along the x-axis, and it thickness along the y-axis.

Then it is given the pitch-yaw-roll rotation sequence listed in Table 2 to give it the

orientation shown in Fig. 3, resulting in the following projected areas: Axy = 4.5

cm2, Axz = 6.5 cm2, and Ayz = 3.5 cm2. The target is taken to lie in the x-y plane.

This table also lists the residual mass mr, residual velocity vr, and limit velocity vL

for the 4 orientations. We see that there are significant differences in the residual

velocities. It is important to emphasize that all 4 orientations have the same projected

areas.

Table 2. FATEPEN predictions for a 725-gr steel cuboid with a striking velocity of 1500 f/s impacting

a 0.25-inch mild steel plate using version 3.3.10.3. In each case, the projected areas on the

orthogonal planes are the same. Note, in particular, the sensitivity of the residual velocity, vr.

Pitch (deg) Yaw (deg) Roll (deg) μ (deg) φeff (deg) mr (gr) vr (f/s) vL (f/s)

-160.576322 193.931472 1.054946 23.7 23.7 720.341 841.729 1194.12

-166.139458 299.729320 5.443168 61.2 61.2 723.672 547.712 1345.23

176.987194 67.658507 13.759713 67.7 67.7 723.349 376.671 1424.57

205.968278 176.788266 183.822078 26.2 26.2 721.834 830.670 1195.00

4.2.1 Impact Angle and Effective Yaw Angle

Let us also document how the impact angle and effective yaw angles are calculated

to check on the values output by FATEPEN. Let us consider the area vector for the

front, side, and top of the cuboid. Initially,

Af = W T kˆ, As = LT ˆı, and At = LW ˆ. (55)

And after the pitch-yaw-roll rotation sequence, these vectors can be calculated by

adapting Eqs. 30–32:

A0

f = W T[ sin φyˆı − cos φy sin φpˆ + cos φy cos φpkˆ], (56)

A0

s = LT[ (cos φr cos φy)ˆı+

(sin φr cos φp + cos φr sin φy sin φp)ˆ+

(sin φr sin φp − cos φr sin φy cos φp)kˆ], (57)

17
Approved for public release; distribution is unlimited.

A0

t = LW[(− sin φr cos φy)ˆı+

(cos φr cos φp − sin φr sin φy sin φp)ˆ+

(cos φr sin φp + sin φr sin φy cos φp)kˆ]. (58)

The effective yaw angle, φeff, is the angle between the length of the cuboid and the

target normal. It is obtained from the dot product between A0

f

and kˆ, so from Eq. 56,

φeff = cos−1(cos φy cos φp). However, if this angle turns out to be greater than 90◦,

then we should use the back face. Therefore, the correct angle is actually

φeff = min[cos−1(cos φy cos φp), π − cos−1(cos φy cos φp)]. (59)

The contact angle, μ, is the angle between each face of the cuboid and the target

normal, and the impact angle, μ, is the minimum of all the contact angles. Therefore,

from Eqs. 55, 56, and 57,

μ = min[ cos−1(cos φy cos φp),

π − cos−1(cos φy cos φp),

cos−1(sin φr sin φp − cos φr sin φy cos φp),

π − cos−1(sin φr sin φp − cos φr sin φy cos φp),

cos−1(cos φr sin φp + sin φr sin φy cos φp),

π − cos−1(cos φr sin φp + sin φr sin φy cos φp)]. (60)

Applying these formulas for the 4 orientations shown in Table 2, we find μ =

23.7454◦, μ = 30.4564◦, μ = 24.5429◦, and μ = 26.1524◦. Notice that 2 of these

angles disagree with the values output by FATEPEN as listed in Table 2. Even though

the FATEPEN documentation states that μ is the minimum contact angle, it seems to

always equal the effective yaw angle instead.

A newer version of FATEPEN (3.3.18.0) was also run,∗ which produced the results

listed in Table 3.

∗

I thank Timothy Mallory for making runs with both version 3.3.16.10 and 3.3.18.0, which

produced the same results shown in Table 3.

18
Approved for public release; distribution is unlimited.

Table 3. FATEPEN predictions for a 725-gr steel cuboid with a striking velocity of 1500 f/s impacting

a 0.25-inch mild steel plate using version 3.3.18.0. In each case, the projected areas on the

orthogonal planes are the same.

Pitch (deg) Yaw (deg) Roll (deg) μ (deg) φeff (deg) mr (gr) vr (f/s) vL (f/s)

-160.576322 193.931472 1.054946 23.7 23.7 720.338 757.72 1281.89

-166.139458 299.729320 5.443168 61.2 61.2 723.671 317.05 1456.00

176.987194 67.658507 13.759713 0.0 67.7 . . . 0.0 1549.32

205.968278 176.788266 183.822078 26.2 26.2 721.833 747.18 1282.89

So version 3.3.18.0 seems to say that the residual velocity can be half the striking

velocity, roughly one-fourth the striking velocity, or zero, depending upon angles

that we have no way of measuring during the test, even though the presented areas

on each of the orthogonal planes are all the same.

5. Conclusions

What are the implications from these results? In many cases, all we have are yaw

card estimates of the actual impact presented area, so that at best we only have the

projected area on one plane, the x-y plane. This may be adequate for a cylinder, but

it will not determine the orientation of a cuboid. What we have shown here is that

even if we have ideal measurements of the projected area on orthogonal planes, this

is still not enough to resolve the orientation. As far as FATEPEN is concerned, the

orientation is specified once the impact angle μ and the yaw angle φeff are known. I

am not aware of any measurement that can determine these during the tests. Thus, I

think we are justified in stating that this analysis shows that cuboids are not good

candidates for FATEPEN validation because the results are highly dependent upon

angles that we have no known way of measuring and no way of uniquely computing

from the measured projected areas.

What are the implications for trying to model natural fragments in FATEPEN? If we

expect to model natural fragments as multifaceted solids, then it seems we have a

real problem since we have no way of knowing the many contact angles, let alone

the impact angle. On the other hand, if we model natural fragments as the simplest

shape that has the same measured presented area, then yawed cylinders offer more

promise and can be checked against experiment.

Part of the problem here is that the FATEPEN code turns all shapes into “equivalent

19
Approved for public release; distribution is unlimited.

cylinders”.∗ So while we may be able to compute the contact angle of each face

of the cuboid, that is not necessarily the shape that the FATEPEN code is dealing

with at impact. We are much better off using cylinders since we can be assured that

FATEPEN is using the same shape. Further, it is more robust since we only need

to control the impact presented area and there is only one contact angle, which is

the angle the RCC face makes with the target normal. Once we know the impact

presented area and the cylinder dimensions, the impact angle is uniquely determined.

See Appendix B for an implementation of this procedure.

∗See Yatteau et al.,5 Section 2.1.4 on p. 2-5 and Section 4.1 on p. 4-1.

20
Approved for public release; distribution is unlimited.

6. References

1. Yatteau JD, Zernow RH, Recht GW, Edquist KT. FATEPEN. Version 3.0.0b.

Terminal Ballistic Penetration Model. Applied Research Associates (ARA)

Project 4714, prepared for Naval Surface Warfare Center, Dahlgren, VA; 1999

Jan.

2. Project THOR. The resistance of various metallic materials to perforation by

steel fragments: empirical relationships for fragment residual velocity and residual weight. Aberdeen Proving Ground (MD): Army Ballistic Research Laboratory (US); 1961 Apr. Report No.: TR-47.

3. Kuipers JB. Quaternions and rotation sequences: a primer with applications to

orbits, aerospace, and virtual reality. Princeton (NJ): Princeton University Press;

2002.

4. Mallory TD. FATEPEN steel-on-steel V50 estimates. Aberdeen Proving Ground

(MD): Army Research Laboratory (US); unpublished paper, 2015 Nov.

5. Yatteau JD, Zernow RH, Recht GW, Edquist KT. FATEPEN fast air target

encounter penetration (Ver. 3.0.0) terminal ballistic penetration model. Littleton (CO): Applied Research Associates, Inc.; 2001 Sep., revised 2005 Feb 2.

(Analyst’s manual; vol. 1).

21
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

22
Approved for public release; distribution is unlimited.

Appendix A. Formula for an Euler Sequence of Rotations

23
Approved for public release; distribution is unlimited.

Let the rotation be a Euler sequence about 3 (distinct or repeated) principal axes

unit-vectors ˆe1, ˆe2, ˆe3 as follows:

• First, a rotation of φ1 about ˆe1:

R1 = Rˆe1(φ1). (A-1)

• Second, a rotation of φ2 about the transformed ˆe2 unit vector:

R2 = Rˆe0

2

(φ2), (A-2)

where ˆe0

2 = R1ˆe2.

• Third, a rotation of φ3 about the (doubly) transformed ˆe3 unit vector:

R3 = Rˆe00

3

(φ3), (A-3)

where ˆe00

3 = R2ˆe03 = R2R1ˆe3.

The total combined rotation is

R = R3R2R1 = Rˆe00

3

(φ3)Rˆe0

2

(φ2)Rˆe1(φ1), (A-4)

where the rotations are applied successively from right to left.

Now, the rotation about ˆe0

2

can be obtained by undoing the effect of the first rotation,

performing the rotation about ˆe2, and then rotating back:

R2 = Rˆe0

2

(φ2)

= R1Rˆe2(φ2)R

−1

1

= Rˆe1(φ1)Rˆe2(φ2)R

−1

ˆe1

(φ1). (A-5)

24
Approved for public release; distribution is unlimited.

Similarly,

R3 = Rˆe00

3

(φ3)

= R2Rˆe0

3

(φ3)R

−1

2

= R2R1Rˆe3(φ3)R

−1

1 R

−1

2

(and using Eq. A-5)

= [Rˆe1(φ1)Rˆe2(φ2)R

−1

ˆe1

(φ1)]Rˆe1(φ1)Rˆe3(φ3)R

−1

ˆe1

(φ1)[Rˆe1(φ1)Rˆe2(φ2)R

−1

ˆe1

(φ1)]−1

= Rˆe1(φ1)Rˆe2(φ2)R

−1

ˆe1

(φ1)Rˆe1(φ1)Rˆe3(φ3)R

−1

ˆe1

(φ1)Rˆe1(φ1)R

−1

ˆe2

(φ2)R

−1

ˆe1

(φ1)

= Rˆe1(φ1)Rˆe2(φ2)Rˆe3(φ3)R

−1

ˆe2

(φ2)R

−1

ˆe1

(φ1). (A-6)

The combined rotation collapses into something very simple:

R = R3R2R1

= Rˆe1(φ1)Rˆe2(φ2)Rˆe3(φ3)R

−1

ˆe2

(φ2)R

−1

ˆe1

(φ1)Rˆe1(φ1)Rˆe2(φ2)R

−1

ˆe1

(φ1)Rˆe1(φ1)

= Rˆe1(φ1)Rˆe2(φ2)Rˆe3(φ3). (A-7)

Thus, we get the result that the combined rotation is equal to the successive rotations

about the original unit vectors but applied in reverse order. Explicitly, this means

Rˆe00

3

(φ3)Rˆe0

2

(φ2)Rˆe1(φ1) = Rˆe1(φ1)Rˆe2(φ2)Rˆe3(φ3) . (A-8)

What we have shown so far is more or less a plausibility argument. Here is a

derivation using quaternions. We use the notation

qˆu(φ) = cos φ

2

+ ˆu sin

φ

2

(A-9)

for the unit quaternion that represents a counterclockwise rotation of φ radians about

the unit vector ˆu. Then, referring to Eqs. A-1, A-2, and A-3,

R1 = qˆe1(φ1). (A-10)

R2 = qˆe0

2

(φ2) = cos φ2

2

+ ˆe0

2

sin

φ2

2

, (A-11)

where

ˆe0

2 = qˆe1

(φ1)ˆe2q

−1

ˆe1

(φ1), (A-12)

25
Approved for public release; distribution is unlimited.

so that

qˆe0

2

(φ2) = cos φ2

2

+ ˆe0

2

sin

φ2

2

= cos

φ2

2

+ qˆe1(φ1)ˆe2q

−1

ˆe1

(φ1) sin φ2

2

= qˆe1(φ1)



cos

φ2

2

+ ˆe2 sin

φ2

2



q

−1

ˆe1

(φ1)

= qˆe1(φ1)qˆe2(φ2)q

−1

ˆe1

(φ1), (A-13)

and

R3 = qˆe00

3

(φ3) = cos φ3

2

+ ˆe00

3

sin

φ3

2

, (A-14)

where

ˆe00

3 = qˆe0

2

(φ2)ˆe0

3

q

−1

ˆe0

2

(φ2)

= qˆe0

2

(φ2)qˆe1(φ1)ˆe3q

−1

ˆe1

(φ1)q

−1

ˆe0

2

(φ2)

= [qˆe1(φ1)qˆe2(φ2)q

−1

ˆe1

(φ1)]qˆe1(φ1)ˆe3q

−1

ˆe1

(φ1)[qˆe1(φ1)qˆe2(φ2)q

−1

ˆe1

(φ1)]−1

= qˆe1(φ1)qˆe2(φ2)q

−1

ˆe1

(φ1)qˆe1(φ1)ˆe3q

−1

ˆe1

(φ1)qˆe1(φ1)q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1)

= qˆe1(φ1)qˆe2(φ2)ˆe3q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1), (A-15)

so that

qˆe00

3

(φ3) = cos φ3

2

+ ˆe00

3

sin

φ3

2

= cos

φ3

2

+ qˆe1(φ1)qˆe2(φ2)ˆe3q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1) sin φ3

2

= qˆe1(φ1)qˆe2(φ2)



cos

φ3

2

+ ˆe3 sin

φ3

2



q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1)

= qˆe1(φ1)qˆe2(φ2)qˆe3(φ3)q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1). (A-16)

Therefore, the total combined rotation is

R = R3R2R1 = qˆe00

3

(φ3)qˆe0

2

(φ2)qˆe1(φ1)

= [qˆe1(φ1)qˆe2(φ2)qˆe3(φ3)q

−1

ˆe2

(φ2)q

−1

ˆe1

(φ1)][qˆe1(φ1)qˆe2(φ2)q

−1

ˆe1

(φ1)]qˆe1(φ1)

= qˆe1(φ1)qˆe2(φ2)qˆe3(φ3), (A-17)

as was to be shown.

26
Approved for public release; distribution is unlimited.

The program in Listing A-1 provides a way of checking this result.

Listing A-1. reverse.cpp

1 // reverse.cpp: program to check that a pitch-yaw-roll rotation sequence

2 // about transformed axes is equivalent to the reverse sequence about fixed axes

3

4 #include "Rotation.h"

5 #include <iostream>

6 #include <cstdlib>

7 #include <cmath>

8 #include <iomanip>

9 using namespace std;

10

11 int main( int argc, char* argv[] ) {

12

13 va::Vector i( 1., 0., 0. ), j( 0., 1., 0. ), k( 0., 0., 1. ); // three unit vectors

14 va::Vector i1, j1, k1, i2, j2, k2, i3, j3, k3; // transformed unit vectors

15 va::Rotation Rp, Ry, Rr; // rotations for pitch, yaw and roll

16 double p = 0., y = 0., r = 0.;

17 if ( argc == 4 ) {

18

19 p = atof( argv[1] ) * va::D2R;

20 y = atof( argv[2] ) * va::D2R;

21 r = atof( argv[3] ) * va::D2R;

22 }

23

24 // first, perform pitch about x-axis

25 Rp = va::Rotation( i, p );

26 i1 = Rp * i;

27 j1 = Rp * j;

28 k1 = Rp * k;

29

30 // second, perform yaw about transformed y-axis

31 Ry = va::Rotation( j1, y );

32 i2 = Ry * i1;

33 j2 = Ry * j1;

34 k2 = Ry * k1;

35

36 // third, perform roll about doubly transformed z-axis

37 Rr = va::Rotation( k2, r );

38 i3 = Rr * i2;

39 j3 = Rr * j2;

40 k3 = Rr * k2;

41

42 cout << "First done the conventional way:" << endl;

43 cout << setprecision(6) << fixed;

44 cout << "i3 = " << i3 << endl;

45 cout << "j3 = " << j3 << endl;

46 cout << "k3 = " << k3 << endl << endl;

47

48 // perform rotation sequence in reverse order about fixed axes

49 va::Rotation R = va::Rotation( i, p ) * va::Rotation( j, y ) * va::Rotation( k, r );

50

51 cout << "Now, the same rotations done in reverse order about fixed axes:" << endl;

52 cout << "i = " << R * i << endl;

53 cout << "j = " << R * j << endl;

54 cout << "j = " << R * k << endl << endl;

55

56 va::Vector v = 3.4 * i - 5.7 * j + 2.3 * k; // an arbitrary vector

57

58 cout << "original vector:" << endl;

59 cout << v << endl << endl;

60

61 cout << "transformed vector using the conventional procedure:" << endl;

62 cout << Rr * Ry * Rp * v << endl << endl;

63

64 cout << "transformed vector using the reversed order about fixed axes:" << endl;

65 cout << R * v << endl;

66

67 return 0;

68 }

27
Approved for public release; distribution is unlimited.

For example, the commands

1 c++ -O2 -Wall -std=c++11 -o reverse reverse.cpp -lm

2 ./reverse 65. -137. 54.

will print the results

1 First done the conventional way:

2 i3 = -0.429879 -0.021405 0.902633

3 j3 = 0.591678 0.748463 0.299535

4 k3 = -0.681998 0.662832 -0.309083

5

6 Now, the same rotations done in reverse order about fixed axes:

7 i = -0.429879 -0.021405 0.902633

8 j = 0.591678 0.748463 0.299535

9 j = -0.681998 0.662832 -0.309083

10

11 original vector:

12 3.400000 -5.700000 2.300000

13

14 transformed vector using the conventional procedure:

15 -6.402747 -2.814501 0.650707

16

17 transformed vector using the reversed order about fixed axes:

18 -6.402747 -2.814501 0.650707

28
Approved for public release; distribution is unlimited.

Appendix B. Yaw Angle of a Cylinder as a Function of Shape Factor

29
Approved for public release; distribution is unlimited.

The dimensionless shape factor γ is defined by the equation

Ap = γV 2/3, (B-1)

where Ap is the projected area and V is the volume. The formula for the dimensionless shape factor of a right-circular cylinder (RCC) as a function of yaw angle φy

is

γ(φy) = a sin φy + b cos φy, (B-2)

where

a ≡



π

4

L

D

−2/3

L

D

and b ≡



π

4

L

D

−2/3

π

4

. (B-3)

The yaw angle that gives the maximum shape factor is obtained by setting the

derivative with respect to φy equal to zero and solving for φy:



dγ

dφy



γ=γmax

= a cos φy − b sin φy = 0, (B-4)

which gives

φy γ=γmax = tan−1

a

b



. (B-5)

To simplify the notation, let φby denote this angle: φby ≡ φy γ=γmax. Then,

a = γmax sin φby and b = γmax cos φby, (B-6)

so that Eq. B-2 can be written as

γ(φy) =







γmax cos(φy − φby) if φy > φby

γmax cos(φby − φy) if φy < φby

, (B-7)

where φby = cos−1(b/γmax). Solving for the yaw angle, we get

φy =







cos−1(b/γmax) + cos−1(γ/γmax) if γ < b

cos−1(b/γmax) − cos−1(γ/γmax) if γ ≥ b

, (B-8)

which shows that we can easily get the orientation (yaw angle) of an RCC from just

the impact presented area.

30
Approved for public release; distribution is unlimited.

The code in Listing B-1 implements this procedure by sampling dimensionless shape

factors from a lognormal distribution to generate the appropriate yawed RCC. The

program may be compiled and run with the following commands:

1 c++ -O2 -Wall -std=c++11 -o cyl cyl.cpp -lm

2 ./cyl

Listing B-1. cyl.cpp

1 // cyl.cpp: Implementation of an algorithm for generating FATEPEN RCCs to represent a specified shape factor.

2 // Given a dimensionless shape factor, generates the L/D and yaw angle for the RCC to represent it.

3 // The RCCs are disk-like in that the L/D <= Pi/4, which is necessary since [0.5,4.5]

4 // is the range of shape factors for artillery and rod-like RCCs don’t span this range.

5 // R. Saucier, October 2011

6

7 #include <iostream>

8 #include <iomanip>

9 #include <cstdlib>

10 #include <cmath>

11 #include <random>

12 #include <chrono>

13 using namespace std;

14

15 int main( int argc, char* argv[] ) {

16

17 const int N = 1000; // number of samples

18 const double R2D = 180. / M_PI; // to convert from radians to degrees

19 const double SF_MIN = 0.5; // minimum shape factor (found from artillery fragments)

20 const double SF_MAX = 4.5; // maximum shape factor (found from artillery fragments)

21

22 // default values for the shape factor lognormal distribution from 122mm, 152mm and 155mm artillery

23 double mu = 0.590494; // these two parameters characterize the lognormal shape factor distribution

24 double sigma = 0.323433; // with mode = 1.63, median = 1.80 and mean = 1.90

25

26 // ability to override the default values from the command line by specifying min and max shape factor

27 // that is meant to capture 95% of the complete distribution (from 0.025 to 0.975)

28 if ( argc == 3 ) {

29

30 double sfmin = atof( argv[1] );

31 double sfmax = atof( argv[2] );

32 mu = 0.5 * log( sfmin * sfmax );

33 sigma = log( sfmax / sfmin ) / ( 2. * M_SQRT2 * 1.3859 );

34 }

35

36 unsigned int seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();

37 std::mt19937 rng( seed ); // Mersenne Twister engine

38 std::lognormal_distribution<double> lognormal( mu, sigma ); // lognormal shape factor distribution

39 std::uniform_real_distribution<double> uniform( 0.069, M_PI_4 ); // uniform distribution for L/D

40

41 double a, b, c, th, sf_min, sf_max, sf, l_d, l, d, yaw, V = 1.; // here we use a fragment with unit volume

42 std::cout << std::setprecision(6) <<std::fixed;

43

44 for ( int n = 0; n < N; n++ ) {

45

46 // normally, the shape factor would be provided, but here we get a shape factor within bounds [SF_MIN, SF_MAX]

47 do { sf = lognormal( rng ); } while ( sf < SF_MIN || sf > SF_MAX );

48

49 // now we want to realize this shape factor with a cylinder

50 do {

51 l_d = uniform( rng );

52 c = pow( M_PI_4 * l_d, -2./3. );

53 a = c * l_d;

54 b = c * M_PI_4;

55 sf_min = a;

56 sf_max = sqrt( a * a + b * b );

57 } while ( sf < sf_min || sf > sf_max );

58

59 if ( sf < b ) th = acos( b / sf_max ) + acos( sf / sf_max );

60 else th = acos( b / sf_max ) - acos( sf / sf_max );

61

62 d = pow( V / ( M_PI_4 * l_d ), 1./3. );

63 l = d * l_d;

64 yaw = th * R2D;

65

66 std::cout << l_d << "\t" << yaw << std::endl;

67 }

68

69 return EXIT_SUCCESS;

70 }

31
Approved for public release; distribution is unlimited.

INTENTIONALLY LEFT BLANK.

32
Approved for public release; distribution is unlimited.

List of Symbols, Abbreviations, and Acronyms

TERMS:

FATEPEN: Fast Air Target Encounter Penetration

RCC: right-circular cylinder

RPP: rectangular parallelepiped

MATHEMATICAL SYMBOLS:

Rˆe(θ): Rotation about the unit vector ˆe through the angle θ

33
Approved for public release; distribution is unlimited.

1

(PDF)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

2

(PDF)

DIRECTOR

US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS MGMT

1

(PDF)

GOVT PRINTG OFC

A MALHOTRA

1

(HC)

NVL SURFC WARFARE CTR

ATTN: D DICKINSON G24

6138 NORC AVE STE 313

DAHLGREN VA 22448-5157

1

(HC)

APPLIED RESEARCH ASSOCIATES

ATTN: R ZERNOW

10720 BRADFORD RD STE 110

LITTLETON CO 80127-4298

ABERDEEN PROVING GROUND

8

(PDF)

RDRL SLB D

J COLLINS

RDRL SLB G

D CARABETTA

T MALLORY

RDRL SLB S

J AUTEN

R DIBELKA

S MORRISON

N REED

R SAUCIER

1

(HC)

R SAUCIER

34