Finding Attitude of a Spin Axis From Roll Angles 

by Andrew A. Thompson 

ARL-TR-5469 March 2011 

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-TR-5469 March 2011 

Finding Attitude of a Spin Axis From Roll Angles 

Andrew A. Thompson 

Weapons and Materials Research Directorate, ARL

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)

March 2011 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

1 January 2009–1 January 2010 

4. TITLE AND SUBTITLE 

Finding Attitude of a Spin Axis From Roll Angles 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Andrew A. Thompson 

5d. PROJECT NUMBER 

AH80 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

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

U.S. Army Research Laboratory 

ATTN: RDRL-WML-A 

Aberdeen Proving Ground, MD 21005-5066 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-5469 

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 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

This report describes a method for determining the attitude of a spin axis from rotational measurements between known 

locations. The described method goes beyond “finding down” to solve for the attitude. This solution was conceived as a 

method to enhance a GPS receiver on a spinning body. The method was made more general by only dealing with locations and 

directions and is potentially useful for a wider class of applications. The generalized problem has been solved and simulations 

have been run to verify the implementation of the solution. Pressure waves and electromagnetic energy could be used; systems 

could combine magnetic, visual, radar, IR, and acoustic information to determine attitude. The method described has received 

U.S. Patent 7,388,538 B1. 

15. SUBJECT TERMS

GPS, orientation, finding down, attitude estimation, roll angles 

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION 

OF ABSTRACT 

UU 

18. NUMBER 

OF PAGES 

32 

19a. NAME OF RESPONSIBLE PERSON

Andrew A. Thompson 

a. REPORT 

Unclassified 

b. ABSTRACT 

Unclassified 

c. THIS PAGE 

Unclassified 

19b. TELEPHONE NUMBER (Include area code) 

410-278-6805 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
iii

Contents 

List of Figures iv

Acknowledgments v

1. Introduction 1

2. Coordinate Estimation 1

3. Development Background 2

4. Nonlinear Least Squares 4

5. Roll Angle Measurements 5

6. Simulation 7

7. Precise Derivative 7

8. Closing Remarks 10

Appendix. Matlab Code 11

Distribution List 21 
iv

List of Figures 

Figure 1. Roll angle.........................................................................................................................3

Figure 2. Agreement for randomly chosen directions and spin axis...............................................9

Figure 3. Inner products for perpendicular directions. ...................................................................9

Figure 4. Inner products for 120° separation of directions. ..........................................................10
v

Acknowledgments 

The author would like to thank Tom Harkins for his diligent review and the suggestions he made, 

which resulted in an improvement in the quality of this report. 
vi

INTENTIONALLY LEFT BLANK. 
1

1. Introduction 

This report describes a method for determining the attitude of a spin axis from rotational 

measurements between known locations. This solution was conceived as a method to enhance a 

GPS receiver on a spinning body. The method was made more general by only dealing with 

locations and directions and is potentially useful for a wider class of applications. The 

generalized problem has been solved and simulations have been run to verify the implementation 

of the solution. 

The ability of a single sensor to determine both position and attitude is appealing from a tactical 

perspective. Software modification and the addition of appropriate antennas to a GPS receiver 

would realize a measurement that includes position and attitude. This package would provide all 

the information needed for input to control mechanisms for improving the precision of 

projectiles, indirectly increasing the lethality. The increase in lethality is a direct result of 

control algorithms reducing the delivery error. The accuracy of these algorithms is a function of 

the position and attitude accuracy. Obtaining an acceptable estimate of the spin angle has been a 

stumbling block in this area and is often referred to as the “finding down” problem. The 

described method goes beyond “finding down” to solve for the attitude. 

Another example of a possible use of this method would be to find the attitude of a spacecraft. 

Using stars as known sources of energy and rotational information, the spin axis of the spacecraft 

could be found. Similarly, this method could be applied to estimate the earth’s or any other 

planet’s spin axis. It is also possible to imagine systems using information from many spectral 

bands to obtain the needed roll information. Both pressure waves and electromagnetic energy 

could be used; systems could combine magnetic, visual, radar, IR, and acoustic information to 

determine attitude. The method described has received U.S. Patent 7,388,538 B1. 

2. Coordinate Estimation 

This section discusses some of the issues with GPS and coordinate adjustment. Coordinate 

adjustment is used within GPS receivers to process pseudorange measurements. 

Coordinate adjustment is central to precise estimation of location. Surveying and geodesy both 

have contributed to this body of literature. Wolf and Ghilani1 describe the way coordinate 

estimation is used in surveying. Strang and Borre2 discuss the way coordinate adjustment is used 

 1

Wolf, P. R.; Ghilani, C. D. Adjustment Computations: Statistics and Least Squares in Surveying and GIS; John Wiley and 

Sons: New York, NY, 1997. 

2

Strang, G.; Borre, K. Linear Algebra, Geodesy, and GPS; Wellesly-Cambridge Press: Wellesly, MA, 1997. 
2

in geodesy. Thompson3 discusses how angles, recast as inner product measurements, can be 

used within the umbrella of coordinate adjustment. Finding equations for the measurements in 

terms of the unknown location establishes the necessary function. The gradient of each of these 

functions is used as a row in a least-squares observation matrix. Using an iterative process, the 

location can be estimated. 

3. Development Background 

GPS receivers have advanced in complexity and features. Of particular interest are the receivers 

containing hundreds of thousands of correlators. These correlators are used to search the 

frequency-time delay space and find the GPS signal associated with a given satellite. After 

satellites have been acquired, many of these correlators are idle. 

A spinning round will create a Doppler effect on the received GPS signal for any radially 

oriented, non-circumferential antenna. Spinning rounds typically use ring antennas to bypass 

this situation. Consider the situation with a spinning round using both a ring antenna and a patch 

(or spin-axis-aligned antenna.) The GPS receiver will use the ring antenna to acquire and track 

the GPS signals; thus, the signal tracking information is available within the receiver. The patch 

antenna will have a frequency shift associated with its received signal for its rotation cycle, 

except when it is moving perpendicular to the satellite direction. By using the tracking 

information from the ring antenna and interrogating the patch antenna, the point of alignment can 

be found. The idle correlators can be used for this task. Using this information and timing 

information, roll angles between satellites can be found. With knowledge of position and 

direction to each satellite and three such measurements of roll angle, it is possible to find the 

orientation of the spin axis. Also note that a directional antenna could be used, removing the 

need for a Doppler shift. Signal parameters like signal-to-noise ratio could also be used to find 

roll angles between satellites or beacons. 

The previous discussion gives a possible implementation of the method. The requirements are to 

get roll angles between known directions. To visualize roll angles, consider a half plane attached 

to a spin axis. As this plane rotates, points in the containing space enter and leave this plane; the 

roll angle between two points is the amount of roll between one point leaving this rotating plane 

and the other point entering. Figure 1 is a graphic representation of this idea. The green 

rectangle represents a plane of observation for the sensor. An object in this plane is observed. 

As the observing plane spins to the position of the yellow rectangle, another object is observed. 

The roll angle between these observations is recorded as the measurement. Due to the geometry, 

 3

Thompson, A. Advances in Applied and Computational Mathematics; Nova Science Publishers: New York, NY, 2006; 

Chapter 18. 
3

-1 -0.5 0 0.5 1 0

0.5

1

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

Figure 1. Roll angle. 

the constraints imposed by the roll angles are sufficient to define the orientation or attitude of the 

spin axis. The following extreme situations are discussed in an attempt to clarify the nature of 

the information being used. Consider two directions or vectors and the plane formed by them. If 

the spin-axis vector is orthogonal to this plane, then the roll angle between the two directions is 

the same as the angle between the two directions. If the spin axis is in the same plane between 

the two directions, then the roll angle will be 180°, shown in figure 1 as the green and cyan 

rectangles. However, if the spin axis is located in the same plane so that both directions are to 

one side, the roll angle between the two will be 0°, as they will both be in the same half plane. 

The roll angle between two directions depends on the orientation of the spin axis. Three roll 

measurements will define the attitude of the spin axis; thus, four known directions are required. 

For GPS systems, four satellites are always in view. 

It is a small extension to consider range systems with beacons at known locations being used as 

known directions. Lasers, infrared, and other portions of the electromagnetic spectrum can be 

used as energy sources; obviously, stars could make reliable energy sources. Sound waves and 

pressure waves of various types can be used as long as the direction to the energy source is 

known. The ability to monitor the precise attitude of a spinning shaft has uses in commercial 

enterprise as well as military and space applications. 
4

4. Nonlinear Least Squares 

A brief review of least squares will be given first. Start with the over-determined equation 

FX  Y , where X is the vector of unknown coordinates, Y represents measurements, and F is a 

matrix model of a known linear relationship. Then project both sides onto the column space of 

the matrix F. Multiply each side by the transpose of F and then, assuming X is multiplied by a 

square matrix of full rank, multiply each side by the inverse. The result will be 

X F F F Y t 1 t ( )

  . (1) 

This result, discussed in many textbooks, will be used as the core of an iterative procedure for 

nonlinear least squares. Both F and Y are known values, while X represents the unknown 

parameters. 

Suppose the functional relationship, f , between the measurements and the coordinates ( ) X is 

known. This function is typically nonlinear. To start, the value of X needs to be guessed or

approximated. Denote this value as X0 . The gradient matrix for the measurements with respect 

to the coordinates will now be denoted by F. If X0 is reasonably close to X, then the following 

linearized approximation is valid: 

0 0 0 ( ) ( ) ( ). X f X fX F X X    (2) 

The left side is the measurement vector, while the first term on the right is the evaluation of the 

functional relationship at the point X0 . Moving this term to the left side allows the interpretation 

of the left side as the deviation of the measurements from the approximated coordinates. The

gradient evaluated at the point is multiplied by the unknown deviation from the approximation of

the location. The previous equation can be written in the form Y = FX as follows: 

0 0 0 ( ) ( ) ( ). X f X fX F X X    (3) 

Linear least squares can be used to solve for the location deviation. After the solution is added to 

the current approximation (and becomes the new current approximation), the process can be 

repeated until the left side gets close to 0. If this is understood, what remains is to discuss 

realizations of f and F for different measurements in terms of coordinates. Generally, within 

the F matrix, different types of measurements can be used; however, only roll angles will be

considered. Each measurement will define a row of the matrix.
5

5. Roll Angle Measurements 

First note that all the quantities are vectors on the unit sphere. All directions are normalized. 

The unknown quantity is the vector representing the spin axis. Consider the ideal case where the 

sensor only observes energy in a half plane attached to the spin axis. Recall that the 

Gram-Schmidt process is a way to set up an orthogonal basis for a set of vectors. For two 

vectors, we can use this to create two orthogonal vectors. Consider two known directions and a 

spin axis. The roll angle is the angle between the orthogonal components of the known 

directions to the spin angle. First, each of these vectors will be described, and then their inner 

product will give the angle between them. Let x represent the unknown coordinates of the spin 

axis; let d d 1 2 and be two known unit directions. The orthogonal or radial components of each 

direction are 

1 1 1 , , r d   d x x (4) 

and 

2 2 2 , , r d   d x x (5) 

where the brackets indicate the inner product function. Assume these new vectors are 

normalized using normalization factors n1 and n2. The inner product between these will yield the 

cosine of the angle between them. 

11 2 2 11 2 2

1 2 12

, , cos( ) , , ,dd dd xx xx dp dp x x

n n nn

       (6) 

where the inner products are defined by 1 p and 2 p . In the sequel, the additional subscript 

indicates the component of the represented quantity. This can be written out as 

           11 1 1 21 2 1 12 1 2 22 2 2 13 1 3 23 2 3

1 2

1 d p x d p x d p x d p x d p x d p x

n n          . (7) 

Next, this expression will be expanded and then rearranged into a more concise group of symbols. 

Temporarily ignoring the normalization factors, multiply out each of the terms as follows: 

2 2

11 21 1 21 1 2 11 1 1 2 1 12 22 1 22 2 2 12 2 1 2 2

2

13 23 1 23 3 2 13 3 1 2 3 .

d d pd x pd x p p x d d pd x pd x p p x

d d pd x pd x p p x

     

  

 (8) 

Now the regrouping of terms results in 

 

   

11 21 12 22 13 23 1 21 1 22 2 23 3

222

2 11 1 12 2 13 3 1 2 1 2 3 .

dd dd dd p dx d x d x

p d x d x d x pp x x x

    

    

 (9) 
6

This result can be written in terms of inner products as follows: 

2

1 2 1 2 2 1 12    d d p d p d pp ,,, . x xx (10) 

Recall all the original vectors, including the sensor spin-axis vector, were considered to be 

normalized. The preceding equation can be recast as 

1 2 12 21 12   d d pp pp pp , . (11) 

Reintroducing the normalization factor, 

1 2 12 12 1 2

12 12

, , ,, . d d pp dd d d x x

nn nn

    (12) 

Now, the relationship between the cosine of the roll angle and the coordinates of the normalized 

spin axis has been established. Using this functional relationship, the gradient of the cosine of 

the roll angle with respect to the coordinates of the unknown spin axis can be determined. Next, 

an expression for the partial with respect to the ith component of the spin axis is needed as 

follows: 

12 1 2

1 2 12 21

1 2

, ,,

, , . i i

dd d d x x

n n dd dd x x

x n n i

     

     

(13)

This assumes the normalization factors make negligible contributions to the partial derivative. 

Assuming the normalization factors are constant simplifies each least-squares iteration, since the 

partials associated with the normalization factors are ignored. Since the exact partial is not being 

used, there will be a suboptimal convergence. The tradeoff between the number of steps saved 

per iteration and the additional number of iterations needed can be used to analyze this situation. 

The numerator of the roll angle expression contains the necessary information; thus, it seems the 

normalization factors do not contain information pertinent to the roll angle. In some situations, 

the extra complexity due to the partials of the normalization factors may be tolerated; however, 

as the simulations have shown, the normalization factors can be treated as constants in the 

partial. Although this simplification reduces implementation complexity, it results in a 

suboptimal convergence. At this point, it is possible to form a row in the previously described F

matrix that corresponds to the cosine of the observed roll angle. To estimate the three 

coordinates of the spin axis, a minimum of three rows is required. By dealing directly with the 

cosine of the angle, the partial is simplified in comparison to the partial of the angle with respect 

to the coordinates. This perspective is no problem mathematically. 
7

6. Simulation 

A Matlab toolbox was designed to test these ideas and verify that the method would produce 

accurate results. After completing this validation, the set of GPS toolbox functions from Kai 

Borre was used to test the method under realistic conditions. In these simulations, the orientation 

of the spin axis was found. The method works; what remains is to work out details with specific 

applications. Matlab code for the simulation is included as the appendix. 

7. Precise Derivative 

In a previous discussion of using roll angles to determine attitude, it was stated that a reduced 

form of the derivative could be used in the linearization step of the iterative process. The 

simplification was justified through the use of a simulation. This simplification is investigated to 

give a more quantitative perspective to the former qualitative statement. The expression under 

investigation follows and represents a roll angle about a spin axis (x) between two known 

directions as follows: 

12 1 2

1 2

, ,, . dd dx d x

n n

 (14) 

Here the d terms are known directions of unit magnitude, and x is an unknown direction of unit 

magnitude. The terms in the denominator are both normalization factors and are as follows: 

. . iii n d dxx   (15) 

Adding an additional subscript to identify the component, this can be written as 

 

  

 

1

2 2 2

1 11 2 2 3 3 1 2 11 2 2 3 3 2

2

3 11 2 2 3 3 3

( )( ) ,

( )

ii i i i i i i

i

ii i i

dx d d d d d d x

n

d x xx x x x

d dx dx dx x

    

    

   

 (16) 

or more concisely, 

 

1

3 2 2

11 2 2 33

1

( ). i ij i i i j

j

n d dx dx dx x 

 

        (17) 
8

The partial of equation 17 is described by the following equation: 

  11 2 2

11 2 2 3 3

3 3

( 2( )2 2

1 ) . 2

ik i i

ij i i i j ij j ik k

kj kj i ij k i

j i

d dx dx

d dx d x dx x dx dx

n dx dx

x n

 

     

                

 

 (18) 

The numerator of the first expression has the following derivative: 

12 1 2

12 21

, ,, , ,. i i

i

dd dx d x

d dx d dx

x

 

 

 (19) 

Now that the individual pieces have been defined, they will be assembled using the following 

equation. First, let p represent the numerator of the first expression. The expression of interest 

can be written as 

1 2

. p

n n

 (20) 

The partial derivative of it can be found using the product rule as follows: 

1 1

12 1 2 11 12 12

1 2 21 1 2 .

ii i i

pnn n n p n n pn n pn n

x x xx

             (21) 

This can now be evaluated by substituting the former expressions in the previous equation. If the 

directions and spin axis are known, then each of the terms will indicate the magnitude of the 

particular partial’s influence in relation to the total. If the second and third terms are relatively 

large, this can differ significantly from the abbreviated partial used in the simulation. The 

appendix includes code for evaluating the partial. 

A simulation was designed to compare the partial of p with the partial of 

n1n2

p . The measure of 

agreement was the inner product of the normalized partials. If they both point in the same 

direction, the inner product is 1; if the inner product is 0, the partials are perpendicular, and 

negative values will result in divergence in the estimation process. Three cases were 

investigated. In the first case, the two known directions and the spin axis were chosen randomly. 

The results are shown in figure 2. 

Although the majority of the cases are in the correct direction and thus lead to convergence, there 

are instances approaching complete divergence. From figure 2, it can be conjectured that using 

the partial of p will result in suboptimal convergence that will take a step in the wrong direction 

on ~10% of the iterations and move in the appropriate direction on ~70% of the steps, as 

indicated by positive inner products. 
9

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0

100

200

300

400

500

600

Random Case

Inner Product

Count out of 1000

Figure 2. Agreement for randomly chosen directions and spin axis. 

A second simulation was performed with the two known directions perpendicular to each other. 

Figure 3 illustrates the results of the simulation. Comparing this to figure 2, it can be seen that 

the number and magnitude of divergent instances have decreased. In this case, ignoring the 

partials of the normalization factors has a smaller effect on the convergence of the estimator, and 

less than 2% of the cases are divergent. Figure 4 shows that if the two known directions are 120° 

apart, 10% of the cases diverge. 

-0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0

100

200

300

400

500

Perpendicular Directions

Inner Product

Count out of 1000

Figure 3. Inner products for perpendicular directions. 
10

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0

100

200

300

400

500

120 Degree Directions

Inner Product

Count out of 1000

Figure 4. Inner products for 120° separation of directions. 

Although the estimator converged while ignoring the partials of the normalization factors, the 

results of this section indicate it may be at the cost of extra iterations. An implementation of the 

method that does not include the partials of the normalization factors should present a rationale 

for their omission. 

8. Closing Remarks 

This method can have widespread application. Attitude estimation is a problem in many fields, 

from detecting small shifts in snow to predicting avalanches to aerospace applications. Many 

systems’ sole purpose is to estimate the attitude of an object. A method that estimates attitude 

from known directions and roll angle measurements has been developed, and simulations show 

that if the requirements are met, accurate estimates of the axis orientation can be made. For 

artillery systems, this is a solution to the “find down” problem. For systems already equipped 

with GPS for position information, this will allow the receivers to be extended and also deliver 

attitude information. This method can be used on a spinning sensor where the directions to 

known sources of energy are known. It would also be possible to simulate a spin axis with a 

steerable antenna. A single sensor providing both location and attitude information would solve 

many tactical and test range problems. 
11

Appendix. Matlab Code 

 

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

function [dip, num]=ipiter(a,meas) 

%for i=1:n a(:,i)=a(:,i)/norm(a(:,i)); end 

aip=a'*a; %inner products 

[m n]=size(a); 

%sweep the spin axis component from the other vectors 

sa=a; 

for i=2:n 

 sa(:,i)=a(:,i)-aip(1,i)*a(:,1); 

end 

anorm=sqrt(sum(sa.^2)); %norm of the vectors orth to spin axis 

%get the measurements c notes 

num=anorm; 

for i=2:n-1 

 num(i)=(aip(i,i+1)-aip(1,i)*aip(1,i+1))/(anorm(i)*anorm(i+1)); 

end 

%c the angular situation 

%acos(num)*180/pi; 

% find differential 

dip=zeros(n-2,3); %initialize 

for j=2:n-1 

 denom=anorm(j)*anorm(j+1); 

 for i=1:3 

 dip(j-1,i)=-(a(i,j)*aip(1,j+1)+a(i,j+1)*aip(1,j))/denom; 

 end 

end 

routine ipconverge 

tst=astart(:,1)'; 

%a(:,1)=[-.2;-.9;.0001]; 

it=0; 

[r n_ipc]=size(meas); 

n_ipc=n_ipc-1; 

while (1-abs(tst*a(:,1))>.000001) 

%for i=1:20 

 it=it+1; 

 [dip num]=ipiter(a,meas); 

 y=(meas(2:n_ipc)-num(2:n_ipc))'; 

 xtx=inv(dip'*dip); 

 yadj=xtx*dip'*y; 

 a(:,1)=a(:,1)+yadj; 

 a(:,1)=a(:,1)/norm(a(:,1)); 

 if it>100 break; end 

end 
13

%function y=gpsipsim() 

close all; 

clear all; 

% ------------------------------------------------------------- 

% --- Read ephemeris (RINEX nav-file or YUMA almanac) --- 

% --- and delete double occurences of SVN's --- 

% ------------------------------------------------------------- 

eph = sortrows(rdyuma('yuma79.txt')); 

% ----------------------------------------- 

% --- Station coordinates (WGS84 / XYZ) --- 

% ----------------------------------------- 

station(1) = 3923625; % Delft (Geodesy-building) 

station(2) = 298462; 

station(3) = 5002803; 

cutoff = 10; % Cutoff elevation in degrees 

% ----------------- 

% --- Time/span --- 

% ----------------- 

tsat = mktsat ('28-FEB-2001 12:00:00','28-FEB-2001 14:00:00',300); 

% ------------------------------------------------------- 

% --- Compute satellite positions (XYZ WGS84 & AZ/EL) --- 

% ------------------------------------------------------- 

[xsat,ysat,zsat,azim,elev,rotmat] = cpaziele (tsat,eph,station); 

%get local directions form the rotation matrix 

local_vert=-rotmat(3,:)'; 

local_north=rotmat(1,:)'; 

%pick a random angle to rotate about the local vert 

%make a quaternion rotation about the vert 

%find the new axis by rotating local north about the vert 

d_r=2*pi*rand; 

q_d=q_aa2quat(local_vert,d_r); 

axis=q_rotate(q_d,local_north); 

%follow above procedure choose rotation 

%make rotation quat to rotate about previous horizantal axis 

%calculate new axis by rotaing the verticle axis 

d_r=45*pi/180; 
14

q_d=q_aa2quat(axis,d_r); 

s_axis=q_rotate(q_d,local_vert); 

a=zeros(3,29); 

[r n_epochs]=size(tsat); 

con_its=[]; 

est_err=[]; 

dxtx=[]; 

itlist=[]; 

i_noncon=0; 

for i=1:n_epochs 

 vis_sat = find (elev(:,i) > cutoff); 

 n_direct=size(vis_sat); 

 for j=1:size(vis_sat) 

 tmp=[xsat(vis_sat(j),i);ysat(vis_sat(j),i);zsat(vis_sat(j),i)]; 

 tmp=tmp-station'; 

 tmp=tmp/norm(tmp); 

 a(:,j)=tmp; 

 end 

 a=[s_axis a]; 

 n=n_direct+1; 

 a=a(:,1:n); 

 aip=a'*a; 

 sa=a; 

 for i=2:n 

 sa(:,i)=a(:,i)-aip(1,i)*a(:,1); 

 end 

 anorm=sqrt(sum(sa.^2)); 

 num=a(1,:); 

 aip=a'*a; 

 for i=2:n_direct 

 num(i)=(aip(i,i+1)-aip(1,i)*aip(1,i+1))/(anorm(i)*anorm(i+1)); 

 end 

 meas=num; 

 astart=a; 

 d_r=2*pi*rand; 

 q_d=q_aa2quat(s_axis,d_r); 

 axis3=q_rotate(q_d,axis); 

 d_r=10/180*pi; 

 q_d=q_aa2quat(axis3,d_r); 

 axis_p=q_rotate(q_d,s_axis); 

 a(:,1)=axis_p; 

 ipconverge 

 con_its=[con_its it]; 

 err=astart(:,1)-a(:,1); 

 est_err=[est_err err]; 

 dxtx=[dxtx det(xtx)]; 
15

 if it>100 

 i_noncon=i_noncon+1; 

 noncon(i_noncon).a=astart; 

 noncon(i_noncon).sa=axis_p; 

 noncon(i_noncon).meas=meas; 

 i 

 xtx 

 anorm 

 aip 

 a 

 end 

end 

 

example of use 

gpsipsim 

ipinvest 

a=noncon.a 

ipinvest 

plot_hs 

a=noncon.a 

r_set 

plot_hs 

con_its 

gpsipsim 

con_its 

help find 

find(con_its>100) 

noncon 

size(noncon) 

nonconinv 

plot_hs 

pplot_hs 

plot_hs 

noncon 

r_noncon=noncon 

con_its 

rotatenocon 

rotatnocon 

rotatnoncon 

noncontemp=noncon 

noncon=r_noncon; 

nonconinv 

plot_hs 
16

%function y=ipadjust() 

n=5; %n must b more than 3 

a=randn(3,n); 

for i=1:n a(:,i)=a(:,i)/norm(a(:,i)); end 

aip=a'*a %inner products 

%sweep the spin axis component from the other vectors 

sa=a; 

for i=2:n 

 sa(:,i)=a(:,i)-aip(1,i)*a(:,1); 

end 

anorm=sqrt(sum(sa.^2)); %norm of the vectors orth to spin axis 

%get the measurements c notes 

num=anorm; 

for i=2:n-1 

 num(i)=(aip(i,i+1)-aip(1,i)*aip(1,i+1))/(anorm(i)*anorm(i+1)); 

end 

%c the angular situation 

acos(num)*180/pi 

% find differential 

dip=zeros(n-2,3); %initialize 

for j=2:n-1 

 denom=anorm(j)*anorm(j+1); 

 for i=1:3 

 dip(j-1,i)=-(a(i,j)*aip(1,j+1)+a(i,j+1)*aip(1,j))/denom; 

 end 

end 

%function y=ipconverge_d(d,rotaxis) 

d_r=d/180*pi 

q_d=q_aa2quat(rotaxis,d_r) 

axis=q_rotate(q_d,v) 

a(:,1)=axis 

ipconverge 

y=it 

%function y=ipinvest(a) 

aip=a'*a %inner products 
17

%sweep the spin axis component from the other vectors 

sa=a; 

[r n]=size(a); 

for i=2:n 

 sa(:,i)=a(:,i)-aip(1,i)*a(:,1); 

end 

anorm=sqrt(sum(sa.^2)); %norm of the vectors orth to spin axis 

%get the measurements c notes 

num=anorm; 

for i=2:n-1 

 num(i)=(aip(i,i+1)-aip(1,i)*aip(1,i+1))/(anorm(i)*anorm(i+1)); 

end 

astart=a; 

meas=num; 

tst=astart(:,1)'; 

a(:,1)=starta 

% 

it=0; 

[r n_ipc]=size(meas); 

n_ipc=n_ipc-1; 

seq=[]; 

while (1-abs(tst*a(:,1))>.000001) 

%for i=1:20 

 it=it+1; 

 [dip num]=ipiter(a,meas); 

 y=(meas(2:n_ipc)-num(2:n_ipc))'; 

 xtx=inv(dip'*dip); 

 yadj=xtx*dip'*y; 

 %af=yadj-yadj'*a(:,1)*a(:,1); this method oscilated toom uch 

 %yadj_1=norm(yadj)*af/norm(af); 

 temp=a(:,1)+yadj; %method2 

 temp=temp/norm(temp); %m2 

 t_2=temp-a(:,1); %m2 

 t_2=norm(yadj)/norm(t_2)*t_2; %m2 

 a(:,1)=a(:,1)+1*t_2; %m2 

 a(:,1)=a(:,1)/norm(a(:,1)); %m2 

 %a(:,1)=temp/norm(temp); %method 1 

 seq=[seq a(:,1)]; 

 if it>100 break; end 

end 

figure;plot3(seq(1,:),seq(2,:),seq(3,:),'-+') 

grid on; 

hold on;plot3(astart(1,1),astart(2,1),astart(3,1),'ro') 

plot3(seq(1,1),seq(2,1),seq(3,1),'go') 
18

%this was to plot eigenvalues from [d v]=eig(xtx) 

%x=astart(:,1);for iz=1:3 plot3([x(1) x(1)+.02*v(1,iz)],[x(2) x(2)+.02*v(2,iz)], [x(3) 

x(3)+.02*v(3,iz)]); end 

%function y=ipstart() 

ipadjust 

astart=a; 

meas=num 

v=astart(:,1) 

vperp=[v(2);-v(1);0] 

deg=2*pi*rand 

pertaxis=q_aa2quat(v,deg) 

rotaxis=q_rotate(pertaxis,vperp) 

d5=5/180*pi 

q5=q_aa2quat(rotaxis,d5) 

axis_5=q_rotate(q5,v) 

a(:,1)=axis_5 

ipconverge 

d10=10/180*pi 

q10=q_aa2quat(rotaxis,d10) 

axis_10=q_rotate(q10,v) 

a(:,1)=axis_10 

ipconverge 

%function y=iptest() 

n=3; %n must b more than 3 

a=randn(3,n); 

for i=1:n a(:,i)=a(:,i)/norm(a(:,i)); end 

a; 

a2=a.^2; 

sum(a2) 

aip=a'*a 

ip1=aip(1,3); 

ip2=aip(2,3); 

sa1=a(:,1)-ip1*a(:,3); 

sa1'*a(:,3) 

sa2=a(:,2)-ip2*a(:,3); 

s_ip=sa1'*sa2/(norm(sa1)*norm(sa2)); 
19

acos(s_ip)*180/pi 

acos(aip(1,2))*180/pi 

sa2'*a(:,3) 

ns1=norm(sa1); 

ns2=norm(sa2); 

num=aip(1,2)-aip(1,3)*aip(2,3); 

acos(num/(ns1*ns2))*180/pi 

% find differential 

dip=sa1; 

for i=1:3 

 dip(i)=-(a(i,1)*aip(2,3)+a(i,2)*aip(1,3)); 

end 

dip=dip/(ns1*ns2) 

%function y=nonconinv() 

[r c]=size(noncon); 

for i=1:c 

 a=noncon(i).a 

 starta=noncon(i).sa 

 ipinvest 

end 

%function y=plot_hs() 

n=20 

theta = pi*(-n:2:n)/n; 

phi = (pi/2)*(0:2:n)'/n; 

X = cos(phi)*cos(theta); 

Y = cos(phi)*sin(theta); 

Z = sin(phi)*ones(size(theta)); 

%figure 

plot3(X,Y,Z,'y') 

[r,c]=size(Z); 

hold on 

for i=1:r 

 plot3(X(i,:),Y(i,:),Z(i,:),'y') 

end 

%function y=r_set() 

q1=q_aa2quat([1;1;0],pi/4) 
20

ar=a;for j=1:9 ar(:,j)=q_rotate(q1,a(:,j)); end 

a=ar 

ipinvest 

%function y=rotatnoncon() 

r_noncon=noncon; 

q1=q_aa2quat([1;1;0],pi/4); %rotational axis and angular rotation 

[r c]=size(noncon); 

for i=1:c 

 ar=noncon(i).a;for j=1:9 ar(:,j)=q_rotate(q1,a(:,j)); end 

 r_noncon(i).a=ar; 

 r_noncon(i).sa=q_rotate(q1,noncon(i).sa); 

end 
NO. OF 

COPIES ORGANIZATION 

21

 1 DEFENSE TECHNICAL 

 (PDF INFORMATION CTR 

 only) DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FORT BELVOIR VA 22060-6218 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 IMNE ALC HRR 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL CIM L 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL CIM P 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL D 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 1 DIR USARL 

 RDRL CIM G (BLDG 4600) 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

22

 2 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL SER L 

 M DUBEY 

 B PIEKARSKI 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL SEE E 

 Z SZTANKAY 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL SES 

 J EICKE 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL SES A 

 J PRICE 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL SES S 

 A EDELSTEIN 

 D FLIPPEN 

 A LADAS 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP F A 

 B CHRISTOPHERSON 

 W KONICK 

 D WARD 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL WMM B 

 A FRYDMAN 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIR US ARMY CECOM RDEC 

 AMSEL RD C2 CS CS 

 J VIG 

 FORT MONMOUTH NJ 07703-5601 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR QEM E 

 M BOMUS 

 BLDG 65S 

 PICATINNY ARSENAL NJ 07806-5000 

 2 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEM A 

 S CHUNG 

 T RECCHIA 

 BLDG 95 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEM C 

 J POTUCEK 

 BLDG 61S 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEM C 

 A MOCK 

 BLDG 171A 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP 

 S PEARCY 

 BLDG 94 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP 

 M CILLI 

 BLDG 382 

 PICATINNY ARSENAL NJ 07806-5000 

 3 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP E 

 D CARLUCCI 

 M HOLLIS 

 J KALINOWSKI 

 BLDG 94 

 PICATINNY ARSENAL NJ 07806-5000 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

23

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP F 

 H RAND 

 BLDG 61S 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CMDR US ARMY TACOM ARDEC 

 AMSRD AR AEP I 

 S LONGO 

 BLDG 95N 

 PICATINNY ARSENAL NJ 07806-5000 

 2 CDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP S 

 Q HUYNH 

 T ZAPATA 

 BLDG 94 

 PICATINNY ARSENAL NJ 07806-5000 

 2 CDR US ARMY TACOM ARDEC 

 AMSRD AAR AEP S 

 C PEREIRA 

 BLDG 407 

 PICATINNY ARSENAL NJ 07806-5000 

 2 CDR US ARMY TACOM ARDEC 

 AMSRD AAR AEM L 

 G KOLASA 

 A MOLINA 

 BLDG 65S 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CDR US ARMY TACOM ARDEC 

 AMSRD AAR AEM M 

 M PALATHINGAL 

 BLDG 65S 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CDR US ARMY TACOM ARDEC 

 SFAE SDR SW IW B 

 D AHMAD 

 BLDG 151 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CDR US ARMY TACOM ARDEC 

 SFAE AMO MAS 

 C GRASSANO 

 BLDG 354 

 PICATINNY ARSENAL NJ 07806-5000 

 1 DIR US ARMY RTTC 

 STERT TE F TD 

 R EPPS 

 REDSTONE ARSENAL AL 35898-8052

 1 CDR US ARMY TACOM ARDEC 

 SFAE ASM TMA MS 

 R KOWALSKI 

 BLDG 171A 

 PICATINNY ARSENAL NJ 07806-5000 

 2 PRODUCT MANAGER FOR MORTARS 

 SFAE AMO CAS MS 

 G BISHER 

 D SUPER 

 BLDG 162 S 

 PICATINNY ARSENAL NJ 07806-5000 

 2 CDR US AMRY TACOM ARDEC 

 SFAE AMO CAS 

 A HERRERA 

 M MORATZ 

 BLDG 171A 

 PICATINNY ARSENAL NJ 07806-5000 

 1 CDR OFC OF NAVAL RSCH 

 CODE 333 

 P MORRISSON 

 800 N QUINCY ST RM 507 

 ARLINGTON VA 22217-5660 

 1 DIR NAVAL AIR SYSTEMS CMD 

 TEST ARTICLE PREP DEP 

 CODE 5 4 

 R FAULSTICH 

 BLDG 1492 UNIT 1 

 47758 RANCH RD 

 PATUXENT RIVER MD 20670-1456 

 1 CDR NAWC WEAPONS DIV 

 CODE 543200E 

 G BORGEN 

 BLDG 311 

 POINT MUGU CA 93042-5000 

 2 PROGRAM MANAGER ITTS 

 PEO STRI 

 AMSTI EL 

 C GOODWIN 

 D SCHNEIDER 

 12350 RESEARCH PKWY 

 ORLANDO FL 32826-3276 

 1 CDR US ARMY RDEC 

 AMSRD AMR SG SD 

 P JENKINS 

 BLDG 5400 

 REDSTONE ARSENAL AL 35898-5247 
24

 1 CDR US ARMY AVN MIS CMND 

 AMSRD AMR SG SP 

 P RUFFIN 

 REDSTONE ARSENAL AL 35898 

 2 ARROW TECH ASSOCIATES 

 W HATHAWAY 

 M SOSTRIN 

 1233 SHELBURNE RD STE 8 

 SOUTH BURLINGTON VT 05403 

 1 SAIC 

 D HALL 

 1150 FIRST AVE STE 400 

 KING OF PRUSSIA PA 19406 

 1 AAI CORPORATION 

 M/S 113/141 

 C BEVARD 

 124 INDUSTRY LN 

 HUNT VALLEY MD 21030 

 1 DEPT OF MECHL ENGRG 

 DREXEL UNIV 

 B C CHANG 

 3141 CHESTNUT ST 

 PHILADELPHIA PA 19104 

 3 DEPT OF ELECTRICAL AND 

 COMPUTER ENGRG 

 V ANDERSON 

 G ARCE 

 X MA 

 140 EVANS HALL 

 NEWARK DE 19716 

 1 GEORGIA TECH RESEARCH INST 

 GTRI/ATAS 

 A LOVAS 

 SMYRNA GA 30080 

 1 JOHNS HOPKINS UNIV 

 APPLIED PHYSICS LAB 

 W D’AMICO 

 1110 JOHNS HOPKINS RD 

 LAUREL MD 20723-6099 

 4 CHLS STARK DRAPER LAB 

 J CONNELLY 

 T EASTERLY 

 A KOUREPENIS 

 J SITOMER 

 555 TECHNOLOGY SQ 

 CAMBRIDGE MA 02139-3563 

ABERDEEN PROVING GROUND

 3 CDR US ARMY TACOM ARDEC 

 AMSRD AAR AEF T 

 R LIESKE 

 J MATTS 

 J WHITESIDE 

 BLDG 120 

 1 CDR USAEC 

 CSTE AEC SVE B 

 D SCOTT 

 BLDG 4120 

 33 DIR USARL 

 RDRL CII C 

 B BODT 

 RDRL SLB D 

 J COLLINS 

 L MOSS 

 RDRL VT 

 S WILKERSON 

 RDRL VTA 

 M FIELDS 

 RDRL WML A 

 B FLANDERS 

 A THOMPSON (4 CPS) 

 D WEBB 

 R PEARSON 

 P WYANT 

 R YAGER 

 B OBERLE 

 M ARTHUR 

 RDRL WML F 

 E BUKOWSKI 

 J CONDON 

 B DAVIS 

 T HARKINS 

 D HEPNER 

 K HUBBARD 

 M ILG 

 D LYON 

 D MCGEE 

 P MULLER 

 P PEREGINO 

 RDRL WML G 

 T BROWN 

 RDRL WMS 

 T ROSENBERGER 

 RDRL WML E 

 F FRESCONI 

 B GUIDOS 

 P WEINACHT 

 RDRL WMB D 

 M NUSCA