Linear Theory of a Projectile With a Rotating Internal Part 

in Atmospheric Flight 

by Geoffrey Frost and Mark Costello 

ARL-CR-530 July 2003 

prepared by 

Geoffrey Frost and Mark Costello 

Oregon State University 

Corvallis, OR 97331 

under contract 

DAAD17-02-P-1104 

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. 
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5066 

ARL-CR-530 July 2003 

Linear Theory of a Projectile With a Rotating Internal Part 

in Atmospheric Flight 

Geoffrey Frost and Mark Costello 

Oregon State University 

prepared by 

Geoffrey Frost and Mark Costello 

Oregon State University 

Corvallis, OR 97331 

under contract 

DAAD17-02-P-1104 

Approved for public release; distribution is unlimited. 
ii

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)

July 2003 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

May 2002–May 2003 

5a. CONTRACT NUMBER 

DAAD17-02-P-1104 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Linear Theory of a Projectile With a Rotating Internal Part in Atmospheric Flight 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

Geoffrey Frost and Mark Costello 

5f. WORK UNIT NUMBER 

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

Oregon State University 

Department of Mechanical Engineering 

Corvallis, OR 97331 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

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

U.S. Army Research Laboratory 

ATTN: AMSRL-WM-BC 

Aberdeen Proving Ground, MD 21005-5066 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

ARL-CR-530 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

Dynamic modeling of the atmospheric flight mechanics of a projectile equipped with an internal rotating disk is investigated 

and a modified projectile linear theory is established for this configuration. To model this type of projectile requires alteration 

of several of the coefficients of the epicyclic dynamics leading to changes in the fast and slow epicyclic modes. A study of 

the frequency and damping properties of the epicyclic modes is conducted by systematically varying disk mass, rotational 

speed, and location. It is shown that the presence of an internal rotating disk can cause substantial changes in the epicyclic

dynamics, including instability, in some configurations. 

15. SUBJECT TERMS

flight dynamics, projectile linear theory 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

Peter Plostins 

a. REPORT 

UNCLASSIFIED 

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

52 

19b. TELEPHONE NUMBER (Include area code) 

410-278-8878 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18
iii

Contents 

List of Figures iv

1. Introduction 1

2. Rotating Internal Part Projectile Dynamic Model 2

3. Aerodynamic Forces and Moments 9

4. Rotating Internal Part Projectile Linear Theory 11

5. Example Results 14

6. Conclusions 24

7. References 25

Appendix A. Coefficients of the Epicyclic Dynamic Equations in C-Form 27

Appendix B. Coefficients of the Epicyclic Dynamic Equations in Symbolic-Form for a 

Special Case 41

List of Abbreviations and Symbols 43
iv

List of Figures 

Figure 1. Position coordinates schematic of a rotating internal part projectile. .............................3 

Figure 2. Attitude coordinates schematic of a rotating internal part projectile. .............................3 

Figure 3. Rotating internal part projectile geometry.......................................................................4 

Figure 4. Root-locus for mass ratio = 1/100 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................15 

Figure 5. Root-locus for mass ratio = 1/10 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................15 

Figure 6. Root-locus for mass ratio = 1/5 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................16 

Figure 7. Root-locus for spin ratio = 1/2 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................16 

Figure 8. Root-locus for spin ratio = 1 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................17 

Figure 9. Root-locus for spin ratio = 10 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile)..................................................18 

Figure 10. Root-locus for a disk located 3/4 in off-axis of symmetry in the waterline 

direction (spin ratio = 10, mass ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk angle; 

□ = 180o disk angle; x = rigid projectile).................................................................................19 

Figure 11. Fast mode root-locus for a disk located 1/4 in off-axis of symmetry in the 

waterline and buttline directions (spin ratio = 10; mass ratio = 1/10). ....................................21 

Figure 12. Slow mode root-locus for a disk located 1/4 in off-axis of symmetry in the 

waterline and buttline directions (spin ratio = 10; mass ratio = 1/10). ....................................22 

Figure 13. Root-locus for a disk located at various distances off-axis of symmetry in a 

direction equidistant from the waterline and buttline directions (spin ratio = 10; mass 

ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk angle; □ = 180o disk angle; x = rigid 

projectile). ................................................................................................................................23 

Figure 14. Root-locus for a disk located 3/4 in off-axis of symmetry in a direction 

equidistant from the waterline and buttline directions for various spin ratios (mass 

ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk angle; □ = 180o disk angle; x = rigid 

projectile). ................................................................................................................................23 
1

1. Introduction 

Many conventional projectile configurations contain internal parts that slightly move in flight in 

some way, shape, or form. Fuze mechanisms used on some indirect fire ammunition employ a 

rotor that is permitted to move slightly with respect to the main projectile body. Submunitions 

deployed from indirect fire projectiles are keyed into place inside the round; however, small 

relative motion between parts occurs. These configurations can experience dynamic instability 

typified by large loss in range and large spin decay (1). Soper (2) evaluated the stability of a 

spinning projectile that contains a cylindrical mass fitted loosely into a cylindrical cavity. The 

cylinder is constrained to spin with the main body projectile. It is shown that an unstable coning 

motion exists in which spin decay and cone angle grow proportional to friction coefficient 

between the mass and cavity and the maximum cant angle between the mass and the projectile. 

Using a similar geometric configuration, Murphy (3) developed a quasi-linear solution for a 

projectile with an internal moving part. It provided an explanation of the unusual flight behavior 

exhibited by four projectiles that all contained parts with slight relative motion between 

components. Later, D’Amico (4) performed a detailed series of experiments where a projectile 

with a loose internal part was driven by the rotor of a freely gimbaled gyroscope. The gyroscope 

yaw history and the orbital motion of the loose part were measured and used to predict the 

moment and resulting yaw growth caused by the loose part. Hodapp (5) expanded the work of 

Soper and Murphy (2, 3) by considering a projectile configuration with a partially restrained 

internal member with a mass center offset. Results of this study indicate that small mass center 

offset of the partially restrained internal member can reduce the instability caused by the loose 

internal part. 

New projectiles configurations have emerged that contain multiple moving parts which are 

fundamental to the basic design and operation of the projectile. The gimbal nose projectile 

configuration is an example of one such configuration. It consists of a standard projectile shape 

with a nose section that is free to rotate with respect to the main body. Goddard (6) originally 

conceived this device for aircraft control, and later Barrett and Stutts (7) considered this 

mechanism for active control of munitions. Schmidt and Donovan (8) as well as Costello and 

Agarwalla (9) showed that dispersion of a fin-stabilized direct fire projectile could be reduced by 

more than 50% using a passive gimbal nose to significantly reduce aerodynamic jump. Another 

multiple component configuration is the dual-spin projectile, which consists of forward and aft 

sections connected through a bearing, allowing different spin rates for each section. Smith et al. 

(10) used a dual-spin projectile for active control of an artillery shell by mounting canards on the 

forward section of the projectile. The forward section was roll stabilized to aid the functionality 

of the canards while the aft section provided spin stability. Costello and Peterson (11) developed 

a linear theory for dual-spin projectiles that predicts stability for this configuration while 

Burchett et al. (12) predicted swerve of a dual-spin projectile caused by lateral pulse jets. 
2

An internal rotating disk is an important dynamic component of some new projectile 

configurations. For example, a new concept for generating real-time bomb damage information 

relies on releasing a relatively small sensor projectile that is tethered to a parent bomb. As the 

two projectiles separate, a reel on the parent munition spins. In another application, active 

trajectory control is achieved by controlling the spin rate of the external projectile body, thus 

predictably changing the aerodynamic loads. To implement this control concept requires an 

internal rotating disk. In these cases, weapon system designers require guidance on the effect of 

the rotating internal part as well as guidance on how to optimally configure such a system. The 

work reported here sheds light on these matters by first developing a projectile linear theory 

specific to projectiles with an internal rotating disk. This is followed by a study of the effect of 

the physical parameters such as the placement, orientation, mass, and speed of the rotating disk 

on the epicyclic modes of vibration. 

2. Rotating Internal Part Projectile Dynamic Model 

A projectile containing an axisymmetric rotating internal part that spins at a constant rate, Ω, is 

considered as shown in Figure 1. The mathematical model describing the motion of the 

projectile allows for three translational and three rotational rigid-body degrees of freedom. In 

order to develop the dynamic equations of motion for these six degrees of freedom, three 

separate reference frames are used as shown in Figure 1. The ground surface is used as an 

inertial reference frame with KI

v positive down. A body frame is fixed on the projectile at the 

mass center of the two-body system with BI

v positive out the nose of the projectile. Due to the 

axisymmetric nature of the rotating internal part, the mass center of the internal part is fixed with 

respect to the body frame. The disk is considered to have a known constant spin rate, and its axis 

of spin is specified in the body frame by a set of direction cosine elements. 

The three translational degrees of freedom are the three components of the two-body, mass 

center position vector. 

O C I I KI r xI yJ z v v v v→ = + + . (1) 

A sequence of rotations from the inertial frame to the projectile frame is defined by a set of 

body-fixed rotations that are ordered in the conventional manner as shown in Figure 2. The three 

rotational degrees of freedom are the Euler roll angle (φ ), pitch angle (θ ), and yaw angle (ψ ). 

In the normal process of simplifying the equations of motion using projectile linear theory, an 

intermediate frame is utilized, namely, the no-roll frame, defined as an intermediate frame before 

roll angle rotation. Figure 3 shows the relative locations of the projectile and disk centers of 

gravity and the projectile body centers of pressure. 
3

I J

v

φ

ψ D

KB

v

B J

v

y 

x 

-z 

KI

v

I I

v

φ B J

v

N J

v

KB

v

KN

v

Ω

B I

v

Figure 1. Position coordinates schematic of a rotating internal part projectile. 

KI

v I J

v

I I

v

KN

v

KB

v

B J

v

N J

v

φ

ψ

θ

φ

φ

Figure 2. Attitude coordinates schematic of a rotating internal part projectile. 
4

P D r →

v

P A r →

v

Dr

v

P M r →

v

Pr

v

Figure 3. Rotating internal part projectile geometry. 

The kinematic differential equations define 6 of the 12 dynamic equations needed to describe the 

motion of the states: x, y, z, φ, θ, and ψ. The no-roll variables, , ~, ~, ~, ~, ~u v w p q and r

~ , are chosen 

for the remaining 12 state variables. The transformation from the no-roll frame (N) to the inertial 

frame (I ) is 

[ ] 



















=









































−

−

=





















N

N

N

N

N

N

N

I

I

I

K

J

I

T

K

J

I

s c

c s c s s

c c s s c

K

J

I

v

v

v

v

v

v

v

v

v

θ θ

θ ψ ψ θ ψ

θ ψ ψ θ ψ

0

, (2) 

while the transformation from the projectile body frame, (B), to the no-roll frame, (N), is 

[ ] 



















=









































= −





















B

B

B

B

B

B

N

N

N

K

J

I

T

K

J

I

s c

c s

K

J

I

v

v

v

v

v

v

v

v

v

φ

φ φ

φ φ

0

0

1 0 0

, (3) 

leading to the transformation from the body frame (B) to the inertial frame (I) described by 

[ ] 



















=









































−

+ −

− +

=





















B

B

B

B

B

B

B

I

I

I

K

J

I

T

K

J

I

s s c c c

c s s s s c c c s s s c

c c s s c c s c s c s s

K

J

I

v

v

v

v

v

v

v

v

v

θ φ θ φ θ

θ ψ φ θ ψ φ ψ φ θ ψ φ ψ

θ ψ φ θ ψ φ ψ φ θ ψ φ ψ

. (4) 

In the above equations and the equations shown below, the standard shorthand notation for 

trigonometric functions is used: α α sin( , ) ≡ s α α cos( , ) ≡ c α α tan( . ) ≡ t
5

The mass center velocity vector of the two-body system is defined in each of the reference 

frames previously discussed, 

C I N N N B B B I I KI v uI vJ wK uI vJ wK xI yJ z

v

& v

& v& v v v v v v v = + + = + + = + + ~ ~ ~ / , (5) 

as is the angular velocity vector of the projectile body, 

B I N N N B B KB pI qJ rK pI qJ r v v v v v v v = + + = + + ~ ~ ~ ω / , (6) 

and the angular velocity vector of the no-roll reference frame, 

N I N N KN B rt I qJ r v v v v ~ ~ ~ ω / = − θ + + . (7) 

The angular velocity vector of the internal rotating part with respect to the inertial frame is found 

by summing the angular velocity vector of the projectile body with respect to the inertial frame 

with the angular velocity vector of the internal rotating part with respect to the projectile body 

frame: 

ω D / I ωB / I ω D / B

v v v = + , (8) 

where the angular velocity vector of the internal rotating part with respect to the projectile body 

frame is 

D P x B y B z KB n I n J n v v v v

ω / = Ω + Ω + Ω . (9) 

Applying the transformation given in equation 2 to the mass center velocity components 

expressed in the no-roll reference frame yields 









































−

−

=





















w

v

u

s c

c s c s s

c c s s c

z

y

x

~

~

~

0 θ θ

θ ψ ψ θ ψ

θ ψ ψ θ ψ

&

&

&

. (10) 

Equating the projectile angular velocity vectors described using no-roll frame components and 

using Euler angle rates generates 









































=





















r

q

p

c

t

~

~

~

0 0 1

0 1 0

1 0

θ

θ

ψ

θ

φ

&

&

&

. (11) 

The kinetic differential equations are derived by separating the two-body system at the disk axle 

connection points and considering the reaction forces and moments associated with each 

individual part. A constraint force FR

v and a constraint moment MRv , applied at the disk center 

of gravity, couple the disk and projectile bodies. The translational kinetic differential equations 

for each body are given by equations 12 and 13: 
6

D D I FR WD m a v v v = + ; (12) 

P P I FR WP FA m a v v v v = − + + . (13) 

Summing equations 12 and 13 yields the expression for the translational dynamic equation of 

motion for the two-body system, 

C C I FA WC m a v v v = + , (14) 

where 

mC aC I mDaD / I mPaP/ I

v v v = + (15) 

and 

WC WP WD

v v v = + ; (16) 

N I C I

C I

N

C I V

dt

dV

a / /

/

/

v v v v = +ω × . (17) 

The constraint force is obtained by subtracting equation 13 from equation 12: 

















+ 













 = − +

D P

D P

P

A

R D I P I

m m

m m

m

F F a a

v v v v . (18) 

The acceleration of the mass center of the disk, aD I

v , can be expressed in terms of the 

acceleration of the projectile mass center by using the formula for two points fixed on a rigid 

body. 

 ) ( D / I P / I P / I P D P / I P / I P D a = a + × r → + × × r →

v v v v v v v α ω ω . (19) 

After substituting this expression into equation 18, the constraint force is expressed in the 

following manner: 

















+ 













 = × → + × × → +

D P

D P

P

A

R P I P D P I P I P D

m m

m m

m

F F r r

v v v v v v v ( ) α / ω / ω / . (20) 

The rotational kinetic equations of motion for the projectile and disk bodies are given by 

R

D I

I

M

dt

dH v v

= / ; (21) 

R P D R A

P I

I

M r F M

dt

dH v v v v v

= − − → × + / . (22) 

Summing equations 21 and 22 eliminates the reaction moment and forms the rotational kinetic 

equation for the two-body system expressed in the body reference frame. 
7

P D R A

D I

I

P I

I

r F M

dt

dH

dt

dH v v v v v

+ = − → × + / / . (23) 

Unlike the kinetic translational equation of motion, it is easier to form the rotational equation of 

motion in the body frame and later convert it to the no-roll reference frame when expressing in 

component form. Therefore, the angular momentum derivatives are expressed as 

B I D I

D I

B

D I

I

H

dt

dH

dt

dH

/ /

/ / v v v

= +ω × ; (24) 

B I P I

P I

B

P I

I

H

dt

dH

dt

dH

/ /

/ / v v v

= +ω × . (25) 

The translational dynamic equation given in equation 26 is expressed in the no-roll frame. 









































− −

−

−





















=





















w

v

u

q rt

r rt

r q

Z

Y

X

m

w

v

u

C ~

~

~

0 ~ ~

~ 0 ~

~ ~ 0

~

~

~

1

~

~

~

θ

θ

&

&

&

, (26) 

where the weight force and aerodynamic loads expressed in the no-roll frame are written as 

A A N A N Z AKN F X I Y J v ~ v ~ v ~ v = + + ; (27) 

C W N W N ZW KN W X I Y J v ~ v ~ v ~ v = + + , (28) 

so that 





















+





















=





















W

W

W

A

A

A

Z

Y

X

Z

Y

X

Z

Y

X

~

~

~

~

~

~

~

~

~

, (29) 

where 



















−

=





















θ

θ

c

s

m g

Z

Y

X

C

W

W

W

0 ~

~

~

. (30) 

The aerodynamic forces are described in section 3. 

The body frame components of the rotational dynamics equation of motion given by equation 23 

can be written as 
8

[ ] RD { } BRD

r

q

p

A =





















&

&

&

, (31) 

where 

D D S RPD RPD

T ARD = I P + TD I T − m S S ; (32) 





















+ Ω

+ Ω

+ Ω

−





















−





















−





















−





















=

z

y

x

D D

T

P D

A

A

A

RPD

P

S

PD

PD

PD

s RPD

A

A

A

RD

r n

q n

p n

S T I T

r

q

p

S I

Z

Y

X

S

m

m

z

y

x

m S S S

N

M

L

B ω ω ω ω ; (33) 

and 





















−

−

−

=

0

0

0

PD PD

PD PD

PD PD

RPD

y x

z x

z y

S ; (34) 

[ ]





















−

−

−

=

0

0

0

q p

r p

r q

Sω . (35) 

The unit vectors, nx, ny, and nz, are the set of direction cosine elements that specify the axis of 

spin of the disk in the body frame. 

To be consistent with projectile linear theory, equation 31 is converted to the no-roll frame 

through multiplication of the equation by T Tφ . Also, a change of variables is introduced from 

body frame angular velocity components to no-roll frame components. With these conversions, 

equation 31 is expressed as 

[ ] RD { } BRD

r

q

p

A ~

~

~

~

~ =





















&

&

&

, (36) 

where 

( ) φ Tφ A T I T I T m S S D D S RPD RPD

T

P D

T

RD = + − ~ , (37) 
9





















−





















−





















−





















=

r

q

p

T S I T

Z

Y

X

T S

m

m

z

y

x

m T S S S

N

M

L

B P

T

A

A

A

RPD

T

P

S

PD

PD

PD

RPD

T

S

A

A

A

RD

~

~

~

~

~

~

~

~

~

~ φ ω ω φ φ ω φ 

( ) 



















− + −





















Ω

Ω

Ω

−





















−

r

q

p

T I T I T m S S T

n

n

n

T S T I T

r

q

p

T S T I T T D D S RPD RPD

T

P D

T

z

y

x

D D

T

D

T

D D

T

D

T

~

~

~

~

~

~

φ ω φ φ ω φ φ

& , (38) 

and 





















−

=

φ φ

φ φ φ

s c

T c s

0

0

1 0 0

; (39) 

( )





















− −

= + −

φ φ

φ θ φ φ

c s

T p t r s c

0

0

0 0 0

& ~ ~ . (40)

In the cross product operator matrix, Sω , in equation 38, the body frame angular velocity 

components p, q, and r are replaced by p

~ , c q s r ~ ~ φ + φ , and s q c r ~ ~ − φ + φ , respectively. 

Equations 10, 11, 26, and 36 provide 12 nonlinear differential equations that govern atmospheric 

flight of a projectile equipped with an axisymmetric rotating internal component. With a given 

set of initial conditions, these equations can be numerically integrated forward in time. 

3. Aerodynamic Forces and Moments 

The equations of motion previously discussed are largely driven by the aerodynamic forces and 

moments exerted on the projectile body. The aerodynamic loads consist of steady aerodynamic 

forces and linear Magnus forces and are formulated separately. 





















+





















=





















M

M

M

S

S

S

A

A

A

Z

Y

X

Z

Y

X

Z

Y

X

~

~

~

~

~

~

~

~

~

. (41) 

The steady aerodynamic forces act at the center of pressure of the projectile body and are 

provided in equation 42: 
10





















+

+

+ +

= −





















α

β

α β

α

0 1

0 1

2

2

2

0 2

~

~

~

Z ZB

Y YB

X X X

S

S

S

C C

C C

C C C

q

Z

Y

X

. (42) 

The Magnus force acts at the Magnus force center of pressure, which is different from the center 

of pressure of the steady aerodynamic forces. Figure 3 shows the relative locations of the 

projectile body centers of pressure. 





































−

=





















V

pDC

V

pDC

q

Z

Y

X

YPA

YPA

M

M

M

2

~ 2

~ 0

~

~

~

β

α α . (43) 

The longitudinal and lateral aerodynamic angles of attack used in equations 42 and 43 are given 

in equation 44: 

w u w u v u v u

~/ ~ ) ~/ ~ tan ( ~/ ~ ) ~/ ~ tan (

1 1 = ≅ = ≅ − − α β ; (44) 

) ~ ~ ~ ( 8

2 2 2 2

q = ρD u + v + w

π α . (45) 

Aerodynamic coefficients in equations 42 and 43 depend on the local Mach number at the 

projectile mass center. 

The externally applied moments on the projectile body found on the right-hand side of the 

rotational kinetic equations contain contributions from steady and unsteady aerodynamics and 

Magnus moments. 





















+





















+





















=





















M

M

M

U

U

U

S

S

S

A

A

A

N

M

L

N

M

L

N

M

L

N

M

L

~

~

~

~

~

~

~

~

~

~

~

~

. (46) 

The steady aerodynamic moments are computed for the projectile body with a cross product 

between the steady body aerodynamic force vector and the distance vector from the center of 

gravity to the center of pressure. Magnus moments on the body are computed in a similar way, 

with a cross product between the Magnus force vector and the distance vector from the center of 

gravity to the Magnus center of pressure. The unsteady body aerodynamic moments provide a 

damping source for projectile angular motion and are given by equation 47: 
11



































 +

=





















V

rDC

V

qDC

V

pDC C

D

N

M

L

NR

MQ

LP

DD

U

U

U

2

~ 2

~ 2

~

~

~

~

~

α γ . (47) 

Air density is computed using the standard atmosphere (13). 

4. Rotating Internal Part Projectile Linear Theory 

The preceding equations of motion in their current state are highly nonlinear and while a 

solution, given an initial set of conditions, may be obtained numerically, it is desirable to solve 

them with a closed form solution for increased understanding of the dynamic behavior. Linear 

theory for symmetric rigid projectiles introduces a series of assumptions that yield a refined set 

of linear differential equations that can be solved in closed form. These equations form the basis 

of classic projectile stability theory. This same set of assumptions can be used to establish a 

linear theory for projectiles containing an axisymmetric rotating internal part in atmospheric 

flight. The necessary assumptions are as follows: 

• The variable is changed from no-roll, station line velocity u, to total velocity V. Equations 

48 and 49 relate V and u and their derivatives: 

~2 ~2 ~2 V = u + v + w ; (48) 

V (uu v v ww)/V & ~~& ~~& ~ ~& = + + . (49) 

• The variables are changed from time, t, to dimensionless arc length s. The dimensionless 

arc length, as defined by Murphy (14) is given in equation 50 and is measured in units of 

distance traveled: 

∫ =

t

Vd

D

s

0

1

τ . (50) 

Time and arc length derivatives of a dummy variable ζ are related by 

ζ = (V / D)ζ ′ & ; (51) 

 ) ( / ) ( / 2 ζ

&& = V D ζ ′′ + ζ ′

V′ V . (52) 

• Euler yaw and pitch angles are small. 

sθ ≈ θ 1 cθ ≈ sψ ≈ψ 1 cψ ≈ . (53) 
12

• Aerodynamic angles of attack are small. 

w/V ~ α ≈ v /V ~ β ≈ . (54) 

• The effect of mass and inertia changes of a projectile on stability are well known. To 

properly evaluate the effects on stability solely due to the rotating internal part, the total 

mass of the two-component system is held constant and individual component mass and 

intertia properties are appropriately modified. This is accomplished by reevaluating the 

inertial properties of the projectile after the addition of the internal part based on its 

nominal characteristics prior to the addition of the internal part. Equation 55 provides the 

shift in the projectile’s center of gravity due to the addition of the internal part, and 

equation 56 is used to populate the projectile’s inertia matrix for various disk 

configurations. 





















−

−

−

=





















E E D D

E E D D

E E D D

P

P

P

P

m z m z

m y m y

m x m x

m

z

y

x

1 ; (55) 





















+

+

+

− −





















+

+

+

= −

2 2

2 2

2 2

2 2

2 2

2 2

D D D D D D

D D D D D D

D D D D D D

D D D

T

D

P P P P P P

P P P P P P

P P P P P P

P E P

x z y z x y

x y x z y z

y z x y x z

T I T m

x z y z x y

x y x z y z

y z x y x z

I I m , (56) 

where 

mP = mE − mD . (57) 

• The projectile is aerodynamically symmetric. 

CNR = CMQ . (58) 

 0 CY 0 = CZ 0 = . (59) 

CYB1 = CZB1 = CNA . (60) 

• A flat fire trajectory is assumed, and the force of gravity is neglected for stability analysis. 

• The quantities V and φ are large compared to , ~, ~, ~ θ ,ψ, q r v and w

~ , such that products of 

small quantities and their derivatives are negligible. 

After application of these assumptions, the differential equations can be solved in closed form as 

stated previously. 

x′ = D . (61) 

y′ = (D V )v +ψD ~ . (62) 
13

z′ = (D V )w −θD ~ . (63) 

(D V )p

~ φ′ = . (64) 

(D V )q

~ θ′ = . (65) 

(D V )r

~ ψ′ = . (66) 

( ) V SDCX V 2mC / = − ρ 0 ′ . (67) 

( )

( ) 















+

+ ′ =

DXX PXX

LP DD

I I

SD DC p C V

p 4

2 2 ~ ρ . (68) 





















































−

−

−

=





























′

′

′

′

r

q

w

v

C B F E

B C E F

A D

A D

r

q

w

v

~

~

~

~

0 0

0 0

. (69) 

Equations 61–69 are linear, except for the total velocity, V, which is retained in several of the 

equations. It is assumed that V changes slowly with respect to the other state variables, and is 

considered to be constant where it appears in other dynamic equations. With this assumption, 

the total velocity, the angle of attack dynamics, and the roll dynamics all become uncoupled, 

linear-time invariant equations of motion. 

Equation 69 is the matrix form of the epicyclic dynamic equations. The equations largely 

determine the angle-of-attack stability of the projectile. It should be noted that the basic 

structure of the epicyclic dynamic equations of a rigid projectile are the same as previously 

shown. Differences in the epicyclic dynamics of both configurations are contained in the 

coefficients B, C, E, and F. In the most general case, where the disk is located off the axis of 

symmetry and canted at an arbitrary angle, these coefficients are algebraically lengthy. While 

straightforward and computationally trivial to compute, space limitations here prevent the most 

general form of these coefficients to be listed. However, Appendix A provides these 

coefficients, as they would appear in the C-programming language. The symbolic form for the 

special case of the disk located on and aligned with the projectile axis of symmetry is also 

included in Appendix B. 

The four roots of the characteristic equation shown in equation 69 are given in equation 70: 

() ( )

() ( ) 























 



 + + ± − + − + + −





 



 + − ± − + − + − −

=

A E iF A E CD F i EF BD AF

A E iF A E CD F i AF BD EF

s

4 2 2

4 2 2

2 2

2

1

2 2

2

1

. (70) 
14

These results are identical to conventional rigid projectile analysis. Consequently, rotating 

internal part projectile analysis can be approached in essentially the same manner that rigid 

projectiles are analyzed. 

5. Example Results 

In order to examine the changes an internal rotating disk induces on the dynamic behavior of a 

projectile, the following analysis documents how the fast and slow epicyclic modes change for 

various rotating disk arrangements. Results are shown for a typical 155-mm spin-stabilized 

artillery shell having a nominal weight of 94.88 lbf with a reference area of S = 0.20 ft2

 and a 

reference diameter of D = 0.51 ft. The nominal stationline, buttline, and waterline distances to the 

projectile’s center of gravity are xE = 1.06 ft, yE = 0.00 ft, and zE = 0.00 ft. The standard 

aerodynamic center of pressure and the aerodynamic center of Magnus are located at the following 

distances along the stationline, buttline, and waterline: xA = 1.76 ft, yA = 0.00 ft, 

zA = 0.00 ft, and xM = 1.76 ft, yM = 0.00 ft, zM = 0.00 ft. The moments of inertia about the bodyaxis are: IEXX = 0.11 slug – ft2

, and IEYY = 1.40 slug − ft2, and IEZZ = 1.40 slug – ft2. The nominal 

disk is 0.33 ft in diameter, 0.06 ft thick, and weighs 9.5 lbf. It is nominally located on the 

projectile center of gravity and has the following inertia properties: IDXX = 0.0041 slug – ft2

, IDYY = 

0.0021 slug – ft2

, and IDZZ = 0.0021 slug – ft2. Euler angles of the projectile are φ = 0.00 rad, θ = 

0.22 rad, and ψ = 0.00 rad. The projectile has a forward velocity of u

~ = 2710 ft/s and a spin rate 

of p

~ = 1674.10 rad/s. The pitch and yaw rates and the side velocities are all equal to zero. The 

disk spin rate and angle are varied. 

In the analysis to follow, the ratio of the disk mass to projectile mass is dubbed the mass ratio, 

MR, the ratio of the disk spin rate to the projectile spin rate is called the spin ratio, SR, and the 

angle that the disk spin axis is precessed from the projectile axis is given as the disk angle, ψD. 

The disk angle is defined in the inset of Figure 4 for a disk located off the projectile axis of 

symmetry in the buttline direction ( B J

v axis). A disk located on the projectile axis of symmetry 

can be visualized by setting the buttline offset distance equal to zero. 

The effect on the epicyclic dynamics of the system of varying the disk angle for three different 

mass ratios with a spin ratio of 5 is demonstrated in the root-locus plots shown in Figures 4–6. 

The disk angle is varied from 0° (represented by the diamond) to 360°. The circle represents the 

eigenvalues for a disk angle of 90°, the square for a disk angle of 180°, and the bold X denotes 

the location of the eigenvalues determined from a similar rigid artillery round without a rotating 

internal part. These conventions are used for all root-locus plots. In Figures 4–8, the epicyclic 

modes travel from the diamond to the square as the disk angle is increased from 0° to 180°, and 

march back up the same path towards the diamond as the disk angle is further increased from 

180° to 360°. In Figure 4 it can be seen that a disk of mass ratio 1/100 slightly affects the 

projectile fast and slow modes at angles other than 90°. Comparison of the three root-locus plots 
15

-0.22 -0.2 -0.18 -0.16 -0.14 -0.12 -0.1

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

Re(s) x 1000

Im(s)

Figure 4. Root-locus for mass ratio = 1/100 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile). 

-0.25 -0.2 -0.15 -0.1 -0.05

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

Re(s) x 1000

Im(s)

Figure 5. Root-locus for mass ratio = 1/10 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile).
16

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

Re(s) x 1000

Im(s)

Figure 6. Root-locus for mass ratio = 1/5 and spin ratio = 5 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile). 

-0.22 -0.2 -0.18 -0.16 -0.14 -0.12 -0.1

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

Re(s) x 1000

Im(s)

Figure 7. Root-locus for spin ratio = 1/2 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile).
17

-0.22 -0.2 -0.18 -0.16 -0.14 -0.12 -0.1

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

Re(s) x 1000

Im(s)

Figure 8. Root-locus for spin ratio = 1 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile). 

reveals that increasing the mass ratio increases the effect that the spinning disk has on the 

epicyclic dynamics. As shown in Figure 6, when the mass ratio is sufficiently large, the fast 

mode can become unstable for large disk angles greater than 90°. In this case, spin stability of 

the complete round is adversely impacted by the component of disk angular velocity in the 

opposite direction of projectile spin. Also note that the mass ratio has very little effect on the 

stability of the system for a disk angle of 90°. For disk angles >180°, the eigenvalues march 

back along the same path toward 0°. This is expected because for a disk located on the projectile 

axis of symmetry the disk angles from 0° to 180° relative to the projectile are equivalent to the 

angles from 180° to 360° with the direction of spin reversed. 

Figures 7–9 are the root-locus plots obtained by varying the disk angle for three different spin 

ratios with a disk to projectile mass ratio of 1/10. By comparing Figures 7–9, it is shown that 

decreasing spin ratio diminishes the effect that a spinning disk of a given mass has on the 

epicyclic modes. The opposite holds true as well. Increasing the spin ratio increases the 

dynamic effects and is capable of driving the system unstable for certain disk angles. 

Comparison of Figures 6 and 9 shows that the same dynamic effects can be achieved with either 

mass ratio or spin ratio. Figures 4–9 also demonstrate that a disk of any mass ratio, spinning in 

the same direction as the projectile (i.e., disk angle <90°), tends to stabilize the fast mode of the 

system while destabilizing the slow mode. Disk angles of >90° up to 180° have the opposite 

effect. However, for large spin ratios and large mass ratios, the movement of the modes is much 
18

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

Re(s) x 1000

Im(s)

Figure 9. Root-locus for spin ratio = 10 and mass ratio = 1/10 (◊ = 0o, 360o disk angle; 

○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile). 

greater per increase in disk angle. While not explicitly shown previously, stationline position of 

the disk has no effect on the fast and slow modes when the disk mass center is located on the 

projectile axis of symmetry. Moreover, in this case, if the product MRSR is held constant, the fast 

and slow modes do not change with disk location along the projectile axis of symmetry. 

Figures 9 and 10 examine the fast and slow modes for spinning disks that are located off the 

projectile axis of symmetry in the buttline direction and the waterline direction. It is important to 

make this distinction because rotating the disk about the KB

r axis from the projectile axis of 

symmetry produces two unique angles relative to the projectile in each direction. Analysis of the 

fast and slow modes for a system with the disk located off the projectile axis of symmetry in the 

buttline direction and varying its position in the direction of the stationline produce the same 

mode path as that for a disk located on the projectile axis of symmetry and varying its stationline 

location. The path is the same, however, the locations of the modes along the path do vary 

slightly with stationline position. This slight variation is undetectable on the scale used for 

graphs shown in this report. Therefore, Figure 9 can be referred to for a root-locus of the 

previously mentioned system for a disk located anywhere along the stationline and anywhere off 

the axis in the buttline direction. The stationline location does, however, change the fast and 

slow modes of the system for a disk that is located off-axis in the waterline direction as shown in 

Figure 10. 
19

-5 -4 -3 -2 -1 0 1 2

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

Re(s) x 1000

Im(s)

SL = 0.0627

SL = 0.5627

SL = 1.0627

Figure 10. Root-locus for a disk located 3/4 in off-axis of symmetry in the waterline direction (spin ratio = 10, 

mass ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk angle; □ = 180o disk angle; x = rigid projectile). 

In Figure 10, the fast and slow modes are plotted in increments of disk angle, and to add clarity 

to the graphs only the eigenvalues containing positive imaginary parts are shown. The root-locus 

for the stationline of the disk center of gravity that coincides with the stationline of the projectile 

center of gravity is the same as that for disks either located on the axis of symmetry or off-axis in 

the buttline direction. The plots are only shown for stationline positions varied from the aft of 

the projectile to its center of gravity. Disks that are located forward of the projectile center of 

gravity, at distances from the center of gravity equal to those that are aft of it, produce the same 

modes for the respective distances. Notice that the fast modes for SL = 0.0627 ft and 

SL = 0.5627 ft become slow modes at a disk angle of 180° (designated by the square) and the 

slow modes change to fast modes. It is also interesting that for these two curves the directions of 

the root-locus path change for disk angles slightly <180° and slightly >90° (designated by the 

circle). The previously mentioned studies show that for a disk located at the projectile center of 

gravity stationline, the epicyclic dynamics of the system are not affected by the off-axis disk 

location. They also show that a spinning disk mounted at an angle of 90° with a stationline 

position equal to the stationline of the center of gravity of the projectile is the same system as a 
20

rigid projectile without a disk. The slow mode of the system with a disk located off-axis in the 

buttline direction is predominantly on the left side of the s-plane, whereas for a disk located offaxis in the direction of the waterline the slow mode is divided between the left and right sides of 

the s-plane. 

Figures 11 and 12 show the root-locus for various angles of a disk located off-axis in various 

directions lying between the waterline and buttline. The fast mode shown in Figure 11 and the 

slow mode in Figure 12 are for a disk located 1/4 in off the projectile axis of symmetry with a 

spin ratio of 10 and mass ratio of 1/10 mounted at 0.0627 ft along the stationline. The lines 

running vertically represent the epicyclic mode locations for different directions of offset defined 

in the figure legend. The center of gravity location is the same regardless of the disk angle. The 

vertical line running roughly down the center of the plot represents the epicyclic modes for a disk 

located solely in the positive buttline direction. Each horizontal line in the graph represents the 

epicyclic mode locations for a specific disk angle designated by the values shown. The 

horizontal lines located on the left-hand side of the center line represent angles from 0° to 180°, 

and those on the right-hand side represent angles from 180° to 360°. If the fast mode of a disk 

spinning at an angle of 120° and located 45° from the positive buttline direction were desired to 

be known, one would follow the dash-dot line until it crossed the bottom half of the 120° curve 

and arrive at a value of –0.53E-3

+0.017i. These plots definitively show that the epicyclic modes 

of a projectile containing a rotating internal part are significantly affected by the location and 

angle of the rotating internal part. Figure 13 shows a system with a disk mounted at 0.0627 ft 

along the stationline and at various distances off-axis in a direction equidistant from the buttline 

and waterline. Increasing the distance that the disk is located off-axis increases the effects that 

the spinning disk has on the dynamics of the system. Greater distances off-axis cause the fast 

and slow modes to swap at disk angles near 180°, whereas smaller distances, such as 1/4 in, 

shown in Figures 11 and 12, do not. Figure 14 shows the root-locus for a disk located 3/4 in off 

the projectile axis of symmetry in a direction equidistant from the waterline and buttline 

directions for various spin ratios. Like disks located on the axis of symmetry, increasing the spin 

ratio increases the range of the epicyclic modes. For disk angles between 0° and 180° the fast 

mode is predominantly on the left-hand side of the imaginary axis, and the slow mode is 

predominantly on the right-hand side. The opposite is true for disk angles between 

180° and 360°. Increasing the spin ratio moves the modes farther from the imaginary axis for 

disk angles other than 0° and 360° and expands their range along the imaginary axis.
21

Figure 11. Fast mode root-locus for a disk located 1/4 in off-axis of symmetry in the waterline 

and buttline directions (spin ratio = 10; mass ratio = 1/10). 
22

Figure 12. Slow mode root-locus for a disk located 1/4 in off-axis of symmetry in the waterline 

and buttline directions (spin ratio = 10; mass ratio = 1/10). 
23

-3 -2 -1 0 1 2 3

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

Re(s) x 1000

Im(s)

Off-Axis = 1/4"

Off-Axis = 3/4"

Off-Axis = 1-1/4"

Figure 13. Root-locus for a disk located at various distances off-axis of symmetry 

in a direction equidistant from the waterline and buttline directions (spin 

ratio = 10; mass ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk 

angle; □ = 180o disk angle; x = rigid projectile). 

-8 -6 -4 -2 0 2 4 6 8

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045

Re(s) x 1000

Im(s)

Spin Ratio = 5

Spin Ratio = 10

Spin Ratio = 20

Figure 14. Root-locus for a disk located 3/4 in off-axis of symmetry in a direction 

equidistant from the waterline and buttline directions for various spin ratios 

(mass ratio = 1/10) (◊ = 0o, 360o disk angle; ○ = 90o disk angle; □ = 180o 

disk angle; x = rigid projectile).
24

6. Conclusions 

The equations of motion for a projectile containing an axisymmetric rotating internal disk that 

spins at a constant rate have been developed. The model allows for the disk to be located off the 

axis of symmetry of the projectile and oriented at an arbitrary angle relative to the projectile axis 

of symmetry. Projectile linear theory has been modified to accommodate projectile 

configurations that contain an internal rotating disk. The addition of an internal rotating disk 

alters several of the coefficients in the epicyclic dynamic equations leading to modified fast and 

slow epicyclic modes.Using modified projectile linear theory, the effect of disk mass, speed, and 

location is systematically studied. If the mass ratio times the spin ratio is held constant the same 

epicyclic modes are produced for a disk located at any position along the projectile axis of 

symmetry. For a specified mass and spin ratio the stationline position of the disk has no effect 

on the dynamics of the system for a disk located either on the projectile axis of symmetry. 

However, the stationline position does affect the dynamics if the disk is off-axis. In these 

instances, the farther the disk is located away from the center of gravity of the projectile the more 

influence its mass, spin rate, and angle relative to the projectile axis of symmetry have on the 

epicyclic dynamics of the system. 
25

7. References 

1. Murphy, C. Symmetric Missile Dynamic Instabilities. Journal of Guidance and Control

1981, 4 (5), 464–471. 

2. Soper, W. Projectile Instability Produced by Internal Friction. AIAA Journal 1978, 16 (1), 

8–11. 

3. Murphy, C. Influence of Moving Internal Parts on Angular Motion of Spinning Projectiles. 

Journal of Guidance and Control 1978, 1 (2), 117–122. 

4. D’Amico, W. Comparison of Theory and Experiment for Moments Induced by Loose 

Internal Parts. Journal of Guidance and Control 1987, 10 (1), 14–19. 

5. Hodapp, A. Passive Means for Stabilizing Projectiles With Partially Restrained Internal 

Members. Journal of Guidance and Control 1989, 12 (2), 135–139. 

6. Goddard, R. Apparatus for Steering Aircraft. U.S. Patent 2,594,766, April 1952. 

7. Barrett, R.; Stutts, J. Modeling, Design, and Testing of a Barrel-Launched Adaptive 

Munition. Proceedings of the 4th Annual Society of Photo-Optical Engineers Symposium on 

Smart Structures, San Diego, CA; Society of Photo-Optical Engineers: New York, 

March 1997. 

8. Schmidt, E.; Donovan, V. Technique to Reduce Yaw and Jump of Fin-Stabilized Projectiles. 

Journal of Spacecraft and Rockets 1998, 35 (1), 110–111. 

9. Costello, M.; Agarwalla, R. Improved Dispersion of a Fin Stabilized Projectile Using a 

Passive Moveable Nose. Journal of Guidance, Control, and Dynamics 2000, 23 (5), 

900–903. 

10. Smith, J.; Smith, K.; Topliffe, R. Feasibility Study for Application of Modular Guidance 

and Control Units to Existing ICM Projectiles; ARLCD-CR-79001; U.S. Army Armament 

Research and Development Command: Picatinny Arsenal, Dover, NJ, 1978. 

11. Costello, M.; Peterson, A. Linear Theory of a Dual Spin Projectile in Atmospheric Flight. 

Journal of Guidance, Control, and Dynamics 2000, 23 (4), 789–797. 

12. Burchett, B.; Peterson, A.; Costello, M. Prediction of Swerving Motion of a Dual-Spin 

Projectile With Lateral Pulse Jets in Atmospheric Flight. Mathematical and Computer 

Modeling 2002, 35 (7–8), 821–834. 
26

13. Von Mises, R. Theory of Flight. Dover: New York, 1959; Chapter 1. 

14. Murphy, C. H. Free Flight Motion of Symmetric Missiles; Report 1216; U.S. Army Ballistic 

Research Laboratory: Aberdeen Proving Ground, MD, 1963. 
27

Appendix  A. Coefficients of the Epicyclic Dynamic Equations in C-Form 

 

 This appendix appears in its original form, without editorial change. 
28

The coefficients of the characteristic equation in C-programming language form are given in 

equations A1–A6. The equations were derived using Wolfram Research Inc.’s Mathematica and 

the following variable substitutions need to be made in order for their appearance to be 

consistent with the text. 

Variable Substitutions 

dx,dy, dz PD PD PD ⇔ x , y , z dax, day, daz PA PA PA ⇔ x , y , z dmx, dmy, dmz 

⇔ xPM , y PM , z PM

phi, theta, psi ⇔ φ,θ ,ψ md ⇔ mD mp ⇔ mP

CYPA ⇔ CYPA CNA ⇔ CNA CMQ ⇔ CMQ

CLP ⇔ CLP CDD ⇔ C DD CX0, CY0, CZ0 

0 0 0 , , ⇔ CX CY CZ

CX2 ⇔ CX 2 CYB1, CZB1 1 1 , ⇔ CYB CZB rho ⇔ ρ

omega ⇔ Ω Cos(disk), Sin(disk), 

0 n X nY nZ ⇔ , ,

JPXX, JPYY, JPZZ, JPXY, JPXZ, JPYZ PXX PYY PZZ PXY PXZ PYZ ⇔ I , I , I , I , I , I

JDXX, JDYY, JDZZ, JDXY, JDXZ, JDYZ DXX DYY DZZ DXY DXZ DYZ ⇔ I , I , I , I , I , I

A1) A = -(CNA*D*rho*S)/(2.*mc) 

A2) B = (D*(-((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*(-(CYPA*D*dy*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Cos(phi))/(2.*(md + mp)) + & 

(CYPA*D*dz*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dy*md*rho*S*V*Sin(phi))/(2.*(md + mp)))) + (JPXX + ((Power(dy,2) + 

& 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CYPA*D*dmx*pt*rho*S)/4. + & 

(CYPA*D*dx*md*pt*rho*S*Power(Cos(phi),2))/(4.*(md + mp)) + & 

(CYPA*D*dx*md*pt*rho*S*Power(Sin(phi),2))/(4.*(md + & 

mp))))/(-Power(Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + mp))*Sin(phi),2) 

+ & 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 
29

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))) + & 

((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((-((Cos(phi)*(JPXY - 

(dx*dy*md*mp)/(md & 

+ mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md 

+ & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md 

+ & 

mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + 

& 

mp))*Sin(phi))*(-(CYPA*D*dy*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Cos(phi))/(2.*(md + mp)) + & 

(CYPA*D*dz*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dy*md*rho*S*V*Sin(phi))/(2.*(md + mp)))) + (JPXX + ((Power(dy,2) + 

& 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CYPA*D*dmx*pt*rho*S)/4. + & 

(CYPA*D*dx*md*pt*rho*S*Power(Cos(phi),2))/(4.*(md + mp)) + & 

(CYPA*D*dx*md*pt*rho*S*Power(Sin(phi),2))/(4.*(md + mp))))) + & 

(-((-(CYPA*D*dy*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Cos(phi))/(2.*(md + mp)) + & 

(CYPA*D*dz*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dy*md*rho*S*V*Sin(phi))/(2.*(md + mp)))*((JPXZ - 

(dx*dz*md*mp)/(md + & 

mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 
30

JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*((CNA*dax*rho*S*V)/2. - & 

(CNA*dx*md*rho*S*V*Power(Cos(phi),2))/(2.*(md + mp)) - & 

(CNA*dx*md*rho*S*V*Power(Sin(phi),2))/(2.*(md + mp))))*(-

Power(Cos(phi)*(JPXY & 

- (dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))))/((-Power(Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + mp))*Sin(phi),2) 

+ & 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))*(-Power(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp),2) + (-Power(Cos(phi)*(JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 
31

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))*(-Power((JPXZ - & 

(dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + 

(JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md 

+ mp) & 

+ JPYY*(md + mp))*Power(Cos(phi),2) - 2*dy*dz*md*mp*Cos(phi)*Sin(phi) + 

& 

(JPYY*md + JPYY*mp + Power(dx,2)*md*mp + Power(dz,2)*md*mp + 

JDYY*(md + & 

mp)*Power(Cos(disk),2) + JDXX*(md + 

mp)*Power(Sin(disk),2))*Power(Sin(phi),2) & 

+ JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))))))/V 

A3) C = (D*(-((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*(-(CYPA*D*dz*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) - & 

(CNA*dy*md*rho*S*V*Cos(phi))/(2.*(md + mp)) - & 

(CYPA*D*dy*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Sin(phi))/(2.*(md + mp)))) + (JPXX + ((Power(dy,2) + 

& 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CNA*dax*rho*S*V)/2. + & 

(CNA*dx*md*rho*S*V*Power(Cos(phi),2))/(2.*(md + mp)) + & 

(CNA*dx*md*rho*S*V*Power(Sin(phi),2))/(2.*(md + mp))))/(-

Power(Cos(phi)*(JPXY & 

- (dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))) + & 

((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 
32

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((-((Cos(phi)*(JPXY - 

(dx*dy*md*mp)/(md & 

+ mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md 

+ & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md 

+ & 

mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + 

& 

mp))*Sin(phi))*(-(CYPA*D*dz*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) - & 

(CNA*dy*md*rho*S*V*Cos(phi))/(2.*(md + mp)) - & 

(CYPA*D*dy*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Sin(phi))/(2.*(md + mp)))) + (JPXX + ((Power(dy,2) + 

& 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CNA*dax*rho*S*V)/2. + & 

(CNA*dx*md*rho*S*V*Power(Cos(phi),2))/(2.*(md + mp)) + & 

(CNA*dx*md*rho*S*V*Power(Sin(phi),2))/(2.*(md + mp))))) + & 

(-((-(CYPA*D*dz*md*pt*rho*S*Cos(phi))/(4.*(md + mp)) - & 

(CNA*dy*md*rho*S*V*Cos(phi))/(2.*(md + mp)) - & 

(CYPA*D*dy*md*pt*rho*S*Sin(phi))/(4.*(md + mp)) + & 

(CNA*dz*md*rho*S*V*Sin(phi))/(2.*(md + mp)))*((JPXZ - 

(dx*dz*md*mp)/(md + & 

mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 
33

JDYY*Power(Sin(disk),2))*(-(CYPA*D*dmx*pt*rho*S)/4. + & 

(CYPA*D*dx*md*pt*rho*S*Power(Cos(phi),2))/(4.*(md + mp)) + & 

(CYPA*D*dx*md*pt*rho*S*Power(Sin(phi),2))/(4.*(md + & 

mp))))*(-Power(Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + mp))*Sin(phi),2) 

+ & 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))))/((-Power(Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + mp))*Sin(phi),2) 

+ & 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))*(-Power(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp),2) + (-Power(Cos(phi)*(JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 
34

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))*(-Power((JPXZ - & 

(dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + 

(JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md 

+ mp) & 

+ JPYY*(md + mp))*Power(Cos(phi),2) - 2*dy*dz*md*mp*Cos(phi)*Sin(phi) + 

& 

(JPYY*md + JPYY*mp + Power(dx,2)*md*mp + Power(dz,2)*md*mp + 

JDYY*(md + & 

mp)*Power(Cos(disk),2) + JDXX*(md + 

mp)*Power(Sin(disk),2))*Power(Sin(phi),2) & 

+ JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))))))/V 

A4) D = D 

A5) E = (D*(-((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*(JDXX*omega*Power(Cos(disk),2)*Sin(disk)*Sin(phi) + & 

JDXX*omega*Power(Sin(disk),3)*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CMQ*Power(D,2)*rho*S*V)/4. - & 

2*JPYZ*pt*Power(Cos(phi),2) + (2*dy*dz*md*mp*pt*Power(Cos(phi),2))/(md + 

mp) & 

+ 2*JDYY*pt*Cos(phi)*Sin(phi) + & 

(2*Power(dy,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

(2*Power(dz,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

2*JDYY*pt*Power(Cos(disk),2)*Cos(phi)*Sin(phi) - & 

2*JDXX*pt*Cos(phi)*Power(Sin(disk),2)*Sin(phi) + 

2*JPYZ*pt*Power(Sin(phi),2) & 

- (2*dy*dz*md*mp*pt*Power(Sin(phi),2))/(md + mp)))/(-Power(Cos(phi)*(JPXY 

- & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 
35

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))) + & 

((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((-((Cos(phi)*(JPXY - 

(dx*dy*md*mp)/(md & 

+ mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md 

+ & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md 

+ & 

mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + 

& 

mp))*Sin(phi))*(JDXX*omega*Power(Cos(disk),2)*Sin(disk)*Sin(phi) + & 

JDXX*omega*Power(Sin(disk),3)*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CMQ*Power(D,2)*rho*S*V)/4. - & 

2*JPYZ*pt*Power(Cos(phi),2) + (2*dy*dz*md*mp*pt*Power(Cos(phi),2))/(md + 

mp) & 

+ 2*JDYY*pt*Cos(phi)*Sin(phi) + & 

(2*Power(dy,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

(2*Power(dz,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

2*JDYY*pt*Power(Cos(disk),2)*Cos(phi)*Sin(phi) - & 

2*JDXX*pt*Cos(phi)*Power(Sin(disk),2)*Sin(phi) + 

2*JPYZ*pt*Power(Sin(phi),2) & 
36

- (2*dy*dz*md*mp*pt*Power(Sin(phi),2))/(md + mp)))) + & 

(-((JDXX*omega*Power(Cos(disk),2)*Sin(disk)*Sin(phi) + & 

JDXX*omega*Power(Sin(disk),3)*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + & 

mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(JDYY*pt*Power(Cos(phi),2)) - & 

JPXX*pt*Power(Cos(phi),2) - 

(2*Power(dy,2)*md*mp*pt*Power(Cos(phi),2))/(md + & 

mp) - JDXX*pt*Power(Cos(disk),2)*Power(Cos(phi),2) + & 

JDYY*pt*Power(Cos(disk),2)*Power(Cos(phi),2) - & 

JDXX*omega*Power(Cos(disk),3)*Power(Cos(phi),2) + & 

JDXX*pt*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

JDYY*pt*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

JDXX*omega*Cos(disk)*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

4*JPYZ*pt*Cos(phi)*Sin(phi) + (4*dy*dz*md*mp*pt*Cos(phi)*Sin(phi))/(md + 

mp) & 

+ JDYY*pt*Power(Sin(phi),2) - JPXX*pt*Power(Sin(phi),2) - & 

(2*Power(dz,2)*md*mp*pt*Power(Sin(phi),2))/(md + mp) - & 

JDXX*pt*Power(Cos(disk),2)*Power(Sin(phi),2) - & 

JDYY*pt*Power(Cos(disk),2)*Power(Sin(phi),2) - & 

JDXX*omega*Power(Cos(disk),3)*Power(Sin(phi),2) - & 

JDXX*pt*Power(Sin(disk),2)*Power(Sin(phi),2) - & 

JDYY*pt*Power(Sin(disk),2)*Power(Sin(phi),2) - & 

JDXX*omega*Cos(disk)*Power(Sin(disk),2)*Power(Sin(phi),2)))*(-

Power(Cos(phi)*( & 

JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-

JPXZ & 

+ (dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))))/((-Power(Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + mp))*Sin(phi),2) 

+ & 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 
37

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))*(-Power(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp),2) + (-Power(Cos(phi)*(JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))*(-Power((JPXZ - & 

(dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + 

(JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md 

+ mp) & 

+ JPYY*(md + mp))*Power(Cos(phi),2) - 2*dy*dz*md*mp*Cos(phi)*Sin(phi) + 

& 

(JPYY*md + JPYY*mp + Power(dx,2)*md*mp + Power(dz,2)*md*mp + 

JDYY*(md + & 

mp)*Power(Cos(disk),2) + JDXX*(md + 

mp)*Power(Sin(disk),2))*Power(Sin(phi),2) & 

+ JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))))))/V 

A6) F = -((D*(-((-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 
38

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp))*(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md 

+ & 

mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + 

& 

mp))*Sin(phi))*(JDXX*omega*Power(Cos(disk),2)*Sin(disk)*Sin(phi) + & 

JDXX*omega*Power(Sin(disk),3)*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(CMQ*Power(D,2)*rho*S*V)/4. - & 

2*JPYZ*pt*Power(Cos(phi),2) + (2*dy*dz*md*mp*pt*Power(Cos(phi),2))/(md + 

mp) & 

+ 2*JDYY*pt*Cos(phi)*Sin(phi) + & 

(2*Power(dy,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

(2*Power(dz,2)*md*mp*pt*Cos(phi)*Sin(phi))/(md + mp) - & 

2*JDYY*pt*Power(Cos(disk),2)*Cos(phi)*Sin(phi) - & 

2*JDXX*pt*Cos(phi)*Power(Sin(disk),2)*Sin(phi) + 

2*JPYZ*pt*Power(Sin(phi),2) & 

- (2*dy*dz*md*mp*pt*Power(Sin(phi),2))/(md + mp)))) + & 

(-((JDXX*omega*Power(Cos(disk),2)*Sin(disk)*Sin(phi) + & 

JDXX*omega*Power(Sin(disk),3)*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + & 

mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + (JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(-(JDYY*pt*Power(Cos(phi),2)) - & 

JPXX*pt*Power(Cos(phi),2) - 

(2*Power(dy,2)*md*mp*pt*Power(Cos(phi),2))/(md + & 

mp) - JDXX*pt*Power(Cos(disk),2)*Power(Cos(phi),2) + & 

JDYY*pt*Power(Cos(disk),2)*Power(Cos(phi),2) - & 

JDXX*omega*Power(Cos(disk),3)*Power(Cos(phi),2) + & 

JDXX*pt*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

JDYY*pt*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

JDXX*omega*Cos(disk)*Power(Cos(phi),2)*Power(Sin(disk),2) - & 

4*JPYZ*pt*Cos(phi)*Sin(phi) + (4*dy*dz*md*mp*pt*Cos(phi)*Sin(phi))/(md + 

mp) & 

+ JDYY*pt*Power(Sin(phi),2) - JPXX*pt*Power(Sin(phi),2) - & 

(2*Power(dz,2)*md*mp*pt*Power(Sin(phi),2))/(md + mp) - & 

JDXX*pt*Power(Cos(disk),2)*Power(Sin(phi),2) - & 

JDYY*pt*Power(Cos(disk),2)*Power(Sin(phi),2) - & 

JDXX*omega*Power(Cos(disk),3)*Power(Sin(phi),2) - & 
39

JDXX*pt*Power(Sin(disk),2)*Power(Sin(phi),2) - & 

JDYY*pt*Power(Sin(disk),2)*Power(Sin(phi),2) - & 

JDXX*omega*Cos(disk)*Power(Sin(disk),2)*Power(Sin(phi),2)))*(-

Power(Cos(phi)*( & 

JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-

JPXZ & 

+ (dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + & 

mp))))/(V*(-Power(-((Cos(phi)*(JPXY - (dx*dy*md*mp)/(md + mp) + (JDXX - 

& 

JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + (dx*dz*md*mp)/(md + & 

mp))*Sin(phi))*((JPXZ - (dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk))*Sin(phi))) + 

& 

((JPXX + ((Power(dy,2) + Power(dz,2))*md*mp)/(md + mp) + & 

JDXX*Power(Cos(disk),2) + JDYY*Power(Sin(disk),2))*((-(dy*dz*md*mp) + & 

JPYZ*(md + mp))*Power(Cos(phi),2) + Cos(phi)*(-(JDYY*md) - JDYY*mp - & 

Power(dy,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2))*Sin(phi) - (-(dy*dz*md*mp) + JPYZ*(md 

+ & 

mp))*Power(Sin(phi),2)))/(md + mp),2) + (-Power(Cos(phi)*(JPXY - & 

(dx*dy*md*mp)/(md + mp) + (JDXX - JDYY)*Cos(disk)*Sin(disk)) + (-JPXZ + 

& 

(dx*dz*md*mp)/(md + mp))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(Power(Cos(phi),2)*(JPYY*md + JPYY*mp + & 

Power(dx,2)*md*mp + Power(dz,2)*md*mp + JDYY*(md + 

mp)*Power(Cos(disk),2) + & 

JDXX*(md + mp)*Power(Sin(disk),2)) + Sin(phi)*(2*dy*dz*md*mp*Cos(phi) + 

& 

((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md + mp) + JPYY*(md + & 

mp))*Sin(phi)) - JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))*(-Power((JPXZ - & 

(dx*dz*md*mp)/(md + mp))*Cos(phi) + (JPXY - (dx*dy*md*mp)/(md + mp) + 

(JDXX - & 

JDYY)*Cos(disk)*Sin(disk))*Sin(phi),2) + ((JPXX + ((Power(dy,2) + & 

Power(dz,2))*md*mp)/(md + mp) + JDXX*Power(Cos(disk),2) + & 

JDYY*Power(Sin(disk),2))*(((Power(dx,2) + Power(dy,2))*md*mp + JDYY*(md 

+ mp) & 
40

+ JPYY*(md + mp))*Power(Cos(phi),2) - 2*dy*dz*md*mp*Cos(phi)*Sin(phi) + 

& 

(JPYY*md + JPYY*mp + Power(dx,2)*md*mp + Power(dz,2)*md*mp + 

JDYY*(md + & 

mp)*Power(Cos(disk),2) + JDXX*(md + 

mp)*Power(Sin(disk),2))*Power(Sin(phi),2) & 

+ JPYZ*(md + mp)*Sin(2*phi)))/(md + mp))))) 
41

Appendix B. Coefficients of the Epicyclic Dynamic Equations in SymbolicForm for a Special Case 

Application of the linear theory assumptions to the simplified system of a projectile containing a 

rotating internal part located on the projectile axis of symmetry at a disk angle, Ω, equal to 0o 

yields the coefficients of the epicyclic equations shown in equations B-1 through B-6: 

C

NA

m

SDC A

2

ρ = − ; (B-1) 

( )















 −































 = V

x x C p

I

m D

m

SD B m YPA

m

PM PD

CYY

C

C

C

D

2

~

2

ρ ; (B-2) 

( )















 −































 = D

x x C

I

m D

m

SD C m NA

m

PA PD

CYY

C

C

C

D

2 2

ρ ; (B-3) 

D = D ; (B-4) 

2 2

2

MQ

CYY

C

C

C

I

m D

m

SD E 





























 = ρ ; (B-5) 

( )

CYY

PXX DXX

I

I p I p

V

D F + + Ω 



 



 =

~ ~

. (B-6) 

Magnus force is assumed to be small in comparison to other aerodynamic forces and is dropped 

from equations B-1 through B-6. However, due to the magnitude amplification resulting from 

the cross product between Magnus force and its respective moment arm, the Magnus moment is 

retained in equation B-2. The equations for the simplified system previously presented are 

essentially the same as those derived for a dual-spin projectile by Costello and Peterson.1 They 

are identical if the aerodynamic coefficients applied on the aft section of a dual-spin projectile 

are neglected.

 1

 Costello, M.; Peterson, A. Linear Theory of a Dual Spin Projectile in Atmospheric Flight. Journal of Guidance, Control, 

and Dynamics 2000, 23 (4), 789–797. 
42

INTENTIONALLY LEFT BLANK. 
43

List of Abbreviations and Symbols 

I I KI I J v v v , = inertial frame unit vectors ,

B B KB I J v v v , = body frame unit vectors ,

N N KN I J v v v , = no-roll frame unit vectors ,

PD PD PD x , y ,z = body frame components of distance vector from projectile center of mass 

to disk center of mass 

P P P x , y ,z = body frame components of distance vector from projectile base to 

projectile center of mass 

D D D x , y , z = body frame components of distance vector from projectile base to disk 

center of mass 

E E E x , y , z = nominal body frame components of distance vector from projectile base to 

projectile center of mass 

PA PA PA x , y ,z = no-roll frame components of distance vector from projectile center of mass 

to center of pressure 

PM PM PM x , y ,z = no-roll frame components of distance vector from projectile center of mass 

to Magnus center of pressure 

u, v,w = translational velocity components of the two-body system center of mass 

resolved in the body frame 

u v w

~, ~, ~ = translational velocity components of the two-body system center of mass 

resolved in the no-roll frame 

p, q,r = roll, pitch, and yaw components of the angular velocity vector of projectile 

body expressed in the body frame 

p q r

~, ~, ~ = roll, pitch, and yaw components of the angular velocity vector of projectile 

body expressed in the no-roll frame 

φ,θ ,ψ = Euler roll, pitch, and yaw angles 

x, y,z = position vector components of the two-body system center of mass 

expressed in the inertial frame 

ω N / I

v = angular velocity vector of no-roll reference frame with respect to the 

inertial frame
44

ω D / I

v = angular velocity vector of disk body with respect to the inertial frame 

ω D / P

v = angular velocity vector of disk body with respect to the projectile body 

frame 

ωB / I

v = angular velocity vector of projectile body with respect to the inertial frame 

C I v /

v = translational velocity vector of two-body system with respect to the 

inertial frame 

aC / I

v = translational acceleration vector of two-body system with respect to the 

inertial frame 

mD = disk body mass 

mP = projectile body mass 

mC = two-body system mass 

mE = nominal projectile mass 

aD / I

v = translational acceleration vector of disk mass center with respect to the 

inertial frame 

aP / I

v = translational acceleration vector of projectile body mass center with 

respect to the inertial frame 

WC

v = weight vector of two-body system 

WP

v = weight vector of projectile body 

WD

v = weight vector of disk body 

FR

v = reaction force vector 

α P / I

v = angular acceleration vector of projectile body with respect to the inertial 

frame 

P D r →

v = distance vector from projectile center of mass to disk center of mass 

FA

v = aerodynamic forces vector 

X A YA Z A , , = aerodynamic force vector components expressed in the body frame 

X A YA Z A

~

, ~, ~ = aerodynamic force vector components expressed in the no-roll frame 

XW YW ZW , = projectile weight vector co , mponents expressed in the body frame
45

XW YW ZW

~

, ~, ~ = projectile weight vector components expressed in the no-roll frame 

X Y Z

~

,

~

, ~ = total external force components on the projectile body expressed in the 

no-roll frame 

TN = transformation matrix from the body frame to the no-roll frame 

LR MR NR , , = reaction moment components due to cross product of distance vector from 

projectile center of mass to disk center of mass with reaction force vector 

expressed in the body frame 

LR M R NR

~, ~ , ~ = reaction moment components due to cross product of distance vector from 

projectile center of mass to disk center of mass with reaction force vector 

expressed in the fixed frame 

HP / I

v = angular momentum vector of the projectile body with respect to the 

inertial frame 

HD / I

v = angular momentum vector of the disk body with respect to the inertial 

frame 

M R

v = reaction moment vector 

M A

v = total aerodynamic moment vector 

HPx HPy HPz , = angular momentum derivative vector components of the projectile body ,

expressed in the body frame 

HPx HPy HPz

~, ~, ~ = angular momentum derivative vector components of the projectile body 

expressed in the no-roll frame 

HDx HDy HDz , , = angular momentum derivative vector components of the projectile body 

expressed in the body frame 

HDx HDy HDz

~, ~, ~ = angular momentum derivative vector components of the projectile body 

expressed in the no-roll frame 

LA M A NA , , = total aerodynamic moment vector components expressed in the body 

frame 

LA M A NA

~

, ~ , ~ = total aerodynamic moment vector components expressed in the no-roll 

frame 

PI = mass moment of inertia of the projectile body about its mass center with 

respect to the body frame
46

DI = mass moment of inertia of the disk body about its mass center with respect 

to the disk frame 

EI = nominal mass moment of inertia of the projectile body about its mass 

center with respect to the body frame 

L M N

~, ~ ,~ = total aerodynamic moment vector components expressed in the no-roll 

frame 

CX 0 = zero yaw axial force aerodynamic coefficient for the projectile body 

CX 0 = zero yaw axial force aerodynamic coefficient for the projectile body 

CX 2 = yaw drag axial force aerodynamic coefficient for the projectile body 

CY 0 = trim side force aerodynamic coefficient for the projectile body 

CYB1 = normal force aerodynamic coefficient along N J

v axis for the projectile 

body 

CZ 0 = trim vertical force aerodynamic coefficient for the projectile body 

CZB1 = normal force aerodynamic coefficient along KN

v axis for the projectile 

body 

CYPA = Magnus force aerodynamic coefficient for the projectile body 

CDD = roll moment aerodynamic coefficient due to fin cant for the projectile 

body 

CLP = roll damping moment aerodynamic coefficient due to fin cant for the 

projectile body 

CMQ = pitch damping moment aerodynamic coefficient due to fin cant for the 

projectile body 

V = velocity magnitude 

α = longitudinal aerodynamic angle of attack 

β = lateral aerodynamic angle of attack 

ρ = density of air 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

1

 2 DEFENSE TECHNICAL 

 INFORMATION CENTER 

 DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FT BELVOIR VA 22060-6218 

 1 COMMANDING GENERAL 

 US ARMY MATERIEL CMD 

 AMCRDA TF 

 5001 EISENHOWER AVE 

 ALEXANDRIA VA 22333-0001 

 1 INST FOR ADVNCD TCHNLGY 

 THE UNIV OF TEXAS AT AUSTIN 

 3925 W BRAKER LN STE 400 

 AUSTIN TX 78759-5316 

 1 US MILITARY ACADEMY 

 MATH SCI CTR EXCELLENCE 

 MADN MATH 

 THAYER HALL 

 WEST POINT NY 10996-1786 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRL D 

 DR D SMITH 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRL CS IS R 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRL CI OK TL 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRL CS IS T 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 2 DIR USARL 

 AMSRL CI LP (BLDG 305) 

 AMSRL CI OK TP (BLDG 4600)