3 - Purdue University



3.4 Attitude Dynamics

Chris Heidelberger

Nomenclature

a = semi-major axis of orbit, km

c.m. = center of mass

e = eccentricity of orbit

ej = j’th quaternion element, rad

E = Eccentric anomaly, rad

H = angular momentum vector, km2/s

kj = j’th non-dimensional shape factor

M = Mean anomaly, 1/s

R = radius of orbit, km

t = current time, s

tp = time of perigee, s

wj = j’th angular velocity element (3), rad/s

W = angular velocity vector, rad/s

V = tangential velocity vector, km/s

3.4.1 Introduction

With the trajectory in place for all three vehicles and a mission requirement met by spinning with a tether, we must now know how the vehicles are going to accomplish other required tasks while in transfer. Attitude dynamics gives insight into how the vehicles are moving in space and the stability of that motion. The two main mission requirements affected by the motion of the vehicle are communications pointing and solar array pointing. Other issues are also affected by the motion; included in this discussion are heating due to solar radiation and end-of-mission translational maneuvers. Finding a solution to these problems involves first determining the orientation and rotation information for the vehicle’s motion in orbit (attitude dynamics) and then determining how this motion and orientation affect the issues. The final solution to these problems involves finding the best motion for the vehicle to follow while in transfer and the best initial orientation for this motion.

3.4.3 Development of Attitude Equations of Motion

To determine the position and attitude information throughout the transfer orbit, we numerically integrate equations of motion (EOMs). These EOMs are for an asymmetric, rigid body in an elliptical orbit; also note that these EOMs are general equations that apply to any rigid body in any closed orbit. There are eight EOMs for this system: three equations for the angular velocities, four equations for the attitude information, and one equation for the radius of the orbit. There are two main systems of numbers used to describe the orientation of a body: Euler angles and quaternions. Euler angles are the three direct angles of rotation from the axes of the reference frame. These are very easy to visualize and to integrate, but their three differential equations involve singularities at specific orientations. Quaternions take the three rotations from the Euler sequence and transform them into one rotation about a single axis. The quaternion system consists of a vector of three scalars that are similar to the angles of rotation from the Euler sequence, and a scalar value that parameterizes the total rotation. Integrating quaternions, however, requires four differential equations instead of three and is not visually intuitive. We use quaternions for orientation information, though, because there are no singularities involved, and numerically integrating the four equations isn’t much more difficult than integrating three equations.

Assuming orientation doesn’t affect translation but translation does affect orientation, the seven EOMs for the body are dependent on the EOM for radius; the radius EOM, though, is independent of the body EOMs. Real applications in industry commonly make this assumption because the effect on translation by orientation is very minute for most systems. With this assumption we can derive the EOMs for angular velocities. Angular velocity changes with angular momentum and the torques on the body. There are many torques acting on a system in orbit in space; these include, among others, gravity gradient torque, solar radiation pressure, magnetic fields, solar wind, and Newtonian drag. Solar radiation pressure is the largest of these torques and must be taken into account in both the trajectory and attitude dynamics of a mission.1 We model it as a constant force acting only on the solar arrays (they are much larger in area and are constantly perpendicular to the force because they are tracking the sun). While gravity gradient is among the smallest of these, it is the most significant torque on the rotation of the system and is dependent only on the position of a body in orbit and the inertia properties of the body.

In order to model the gravity gradient torque on a vehicle in orbit we need to find the radius from the sun to the c.m. of the vehicle, R. Since the vehicle is in an elliptic orbit, R isn’t constant or linear with respect to time; therefore we need a differential equation for R with respect to the independent variable of integration. Time isn’t necessarily the best independent variable to use for a couple reasons: since the vehicle will be transferring over a period of millions of seconds this would involve numerically integrating equations containing very large numbers, which requires a lot of time and computing power. Radius also isn’t easily determined as a function of time. An independent variable that is small in magnitude over the entire integration period and still easily related to radius can be either true anomaly or eccentric anomaly, E. True anomaly is the angle centered on the attracting focus of the orbit and measured from the periapsis of the orbit to the vehicle’s radius vector. It is easier to use, but is more difficult to relate to time. Eccentric anomaly is the angle centered at the center of the ellipse (between the two foci) and measured the same way as true anomaly. It is easily related to time and radius so we use it as the independent variable in integration. Differentiating Kepler’s equation

M*( t – tp ) = E - e*sin (E) (3.4.1)

with respect to time gives an equation for the eccentric anomaly rate with respect to time. Using this rate we can find differential equations for radius, angular velocities, and orientation with respect to eccentric anomaly. We can also non-dimensionalize the equations by dividing by the mean anomaly. Equations 3.4.2-8 are the EOMs that result from this analysis and eq. 3.4.9 is the differential equation for the radius of the orbit.

w’1=[1-e*cos(E)]*k1*[w2*w3 - 12*a3/R3*(e1*e2-e3*e4)*(e3*e1+e2*e4)] (3.4.2)

w’2=[1-e*cos(E)]*k2*[w3*w1 - 6*a3/R3*(1-2*e22-2*e32)*(e3*e1+e2*e4)] (3.4.3)

w’3=[1-e*cos(E)]*k3*[w1*w2 - 6*a3/R3*(1-2*e22-2*e32)*(e1*e2-e3*e4)] (3.4.4)

e’1=0.5*[1-e*cos(E)]*{w1*e4 - w2*e13 + [w3 + sqrt(a3/R3)]*e2} (3.4.5)

e’2=0.5*[1-e*cos(E)]*{w1*e3 + w2*e4 – [w3 + sqrt(a3/R3)]*e1} (3.4.6)

e’3=0.5*[1-e*cos(E)]*{-w1*e2 + w2*e1 + [w3 - sqrt(a3/R3)]*e4} (3.4.7)

e’4=-0.5*[1-e*cos(E)]*{w1*e1 + w2*e2 + [w3 - sqrt(a3/R3)]*e3} (3.4.8)

R’=a*e*sin(E) (3.4.9)

3.4.4 Stability Analysis

Now that we have the EOMs for any system in a closed orbit we simulate the motion of our vehicles. Running a simulation requires the initial conditions of the orientation and rotation of the body, the orbital elements for the transfer, and the mass properties of the body. For the stability analyses, we run test cases for two nominal motions and then cases with initial perturbations to those motions. To determine stability we look at four properties of the motion: nutation or tilt of the spin plane, precession or rotation of the spin plane, spin of the body in the near-axial direction, and the spin of the body about the c.m. for artificial gravity. All of these parameters are presented in the pseudo-inertial reference frame of the solar system. The parameters are integrated with respect to the orbital reference frame and then converted to the inertial frame.

Habitation module stability

Using the mass properties of the Habitation module (Hab) and the orbital elements of the free return trajectory, we can determine its stability. The first nominal motion (nominal motion 1) is one where the body is spinning in the orbital plane with its angular velocity vector perpendicular to the orbit plane (fig. 3.4.1). Figure 3.4.2 shows the results of this simulation.

[pic]

Fig. 3.4.1 Nominal motion 1, spinning in the orbital plane.

[pic]

Fig. 3.4.2 Stability parameters of Hab module for nominal motion 1.

It can be seen that tilt, precession and axial spin all remain constant through the entire transfer. This is expected because we have a large body spinning at a relatively fast rate so H is very large. We are very far from the sun so the gravity gradient torque is very small, and both solar pressure and gravity moments are counteracted on each half in a single rotation. This nominal motion, therefore, is stable.

Next, we perturb this motion by an initial offset of 0.5 and 1.0 degrees. Misaligned thrusters during spin up, solar pressure on the body while spinning, solar flares, and micrometeorite hits contribute to perturbations. These perturbations can affect the motion in two ways. The first is the tilt of the body’s spin plan as in a pitch maneuver while the second would act like a roll maneuver of the spin plane. The simulations for the first type are shown in fig. 3.4.3.

[pic]

Fig. 3.4.3 Perturbations of nominal motion 1 of 0.5 and 1.0 degree.

Though the body does change its orientation over the orbit, it is only on the order of 10-5 degrees. The body now spins on a plane that is inclined by 0.5 degrees. When the body is perturbed by twice this amount, 1.0 degrees, the stability results are almost exactly the same except the plane is tilted by 1.0 degrees. This shows that this motion and initial orientation are stable because perturbations stay in the same magnitude and are stable. The perturbed solutions don’t converge to the nominal motion, though, so the system isn’t asymptotically stable. Analysis of the second type of perturbation yields identical results for tilt and body spin, but the negative response for precession and axial spin. This means the motion is stable for this type of perturbation and is therefore a stable motion.

The second nominal motion is one where the body spins perpendicular to the orbital plane and W lies in the orbital plane (referred to as nominal motion 2 or spinning out of the plane). There are an infinite number of variations on this motion, depending on where w is pointing. The vehicle can start out spinning in the direction of V, where W is perpendicular to V; or the vehicle can spin perpendicular to V so that W is aligned with V. Analysis shows that all of these motions are similar in their dynamic response and stability so we just present the case where the angular velocity vector is aligned with the tangential velocity vector (fig. 3.4.4).

[pic]

Fig. 3.4.4 Stability of Hab module in nominal motion 2.

Although this nominal motion does change through the transfer it is only on the order of 10-6 degrees. Rotating the spin plane does not change the dynamics of the body because the spin plane remains inertially fixed, therefore the spin plane can be changed to optimize the communication, solar array, and thermal pointing.

The same types of disturbances from the first nominal motion can cause initial perturbations in this nominal motion, too. There are also two different types perturbations that are created. The first one is just a rotation of the spin plane and is just a special case of the nominal motion. The second one is a tilt in the spin plane making the angular velocity vector to be no longer aligned with the tangential velocity vector. Figure 3.4.5 shows this analysis.

[pic]

Fig. 3.4.5 Stability of Hab with perturbations to nominal motion 2.

The perturbations stay in the same magnitude and are stable; therefore, this nominal motion is stable. Because the solutions don’t converge to the nominal solution, the system is not asymptotically stable. An important thing to note is that the spin plane in this nominal motion is inertially fixed (fig. 3.4.6). This means it rotates with respect to the orbital reference frame, which causes problems when the vehicles have to track the sun or Earth.

[pic]

Fig.3.4.6 Nominal motion 2, inertially fixed spin plane and pointing issues.

Earth Return Assembly stability

The Earth Return Assembly (ERA) on the initial transit will not need to be tethered for artificial gravity so it can use one of three methods to maintain stability throughout the orbit. The first and simplest method is spin stabilization; the vehicle spins about its symmetric axis creating a large angular momentum vector that resists the external torques making the vehicle stable. The second method is three-axis stabilization. This involves using the vehicles mass properties to balance the gravity torques making the vehicle stable. The final method uses active and passive control systems to perform maneuvers and counteract adverse torques. This method is more complex and would weigh more, but more reliably creates stability. Three-axis stabilization is the least propulsion intensive because the vehicle does not have to spin up; it also does not have the intensive communications and solar array pointing problems because it isn’t spinning. Figure 3.4.7 shows that the ERA is stable just by being in orbit in a certain orientation.

[pic]

Fig. 3.4.7 ERA with three-axis stabilization; perturbations of 0.5o and 1.0o

Since it is stable with any control, three-axis stabilization does not need extra control system mass, which is expensive. This orientation, however, is contained in a special set orientations. This set corresponds to the cases where the axis of symmetry is aligned with the radius vector (asymptotically stable) and the marginally stable case where the vehicle’s axis of symmetry is perpendicular to the radius vector and still in the orbit plane (fig. 3.4.7 and 3.4.8). If the axis of symmetry is perpendicular to the orbit plane and R, it is not a stable orientation. After the transfer orbit injection from the propulsion system, the ERA will be close to the second of the two stable orientations (it will be off by the magnitude of the flight path angle) so we choose this orientation as our initial condition if we use three-axis stabilization.

[pic]

Fig. 3.4.8 Nominal motion for ERA with three-axis stabilization

Earth Return Vehicle stability

The Earth Return Vehicle (ERV) on the return trip has similar mass properties and is dynamically similar to the Hab, which gives similar stability results. The main difference between the Hab and the ERV is that the ERV starts spinning at 1.24 rpm’s to give 0.38g and then spins up 0.2 rpm’s per day until it is spinning at 2 rpm’s for 1g. This spin up only makes the vehicle more stable.

3.4.5 Initial Conditions for Pointing

Solar array pointing

In order to draw power from the solar arrays they need to be pointed at the sun. If the vehicle is spinning, this creates quite a problem. The ERA does not need to have its arrays deployed during transfer so solar array pointing will not factor in the design of the ERA. The spinning vehicles, Hab and ERV, both have solar arrays deployed for either primary or secondary power systems. If the vehicle spins out of the orbit plane, the arrays have to track the sun while the vehicle moves around it, compensate for the spin of the vehicle, and account for the rotation of the spin plane in the orbit plane. There are times in this nominal motion, however, when the arrays can only point towards the sun during half of the rotation; during the other half of the rotation, the arrays are pointed away from the sun (fig. 3.4.6). This requires a propulsive maneuver to compensate for the rotation of the spin plane with respect to the orbit frame. This is not an attractive option since propulsive maneuvers are very expensive. The first nominal motion, when the vehicle spins in the plane, only requires the arrays to track the movement away from the sun and compensate for the spin of the vehicle. The best way to solve this is to place the arrays on a de-spun section to overcome the vehicle spin and have small motors to track the sun. Spinning in the plane is therefore the better solution for the Hab and ERV since it is easier to track the sun and requires no propulsive maneuvers.

Communications Pointing

In any manned mission, constant communication is major requirement. This does not necessarily mean constant high-bandwidth communications, though. Over the six-month transfer orbit, there are periods when only intermittent voice or data commands are needed. For low speed telemetry and voice transmission, an omni-directional antenna is acceptable and has no pointing requirements. Also due to the long-term transfer, high-speed voice and video transmissions are a requirement. This means the vehicles need a high-gain communications dish. These dishes require a direct line of sight to receive and transmit, though, so there are very stringent pointing requirements. Since we do not know when our vehicle needs full bandwidth and when it doesn’t, we have to design our pointing system to allow for constant line of sight pointing to the receiving dishes on Earth. For the ERA, which does not have a spinning requirement, the vehicle can be three-axis stabilized and it is easy to have a small motor that tracks Earth in its orbit since it is not moving very fast relative to the ERA. This is also the case if active control systems are used on the ERA. If the ERA were spin-stabilized, there would be large pointing problems that would have to be overcome with complicated motors and controllers. With spinning vehicles such as the Hab or ERV, the dish would have to be placed on a de-spun section and would also have to track the earth. When the vehicles spin out of the orbit plane, the dish must track in three dimensions and has to change its pointing angle by large amounts in one spin during certain periods in a transfer (fig. 3.4.6). These two issues require complicated motor and control systems. Vehicles spinning in the orbital plane would only have to track in two dimensions (with small three dimensional variations when the transfer orbit is a little off of the Ecliptic plane). If we just place the dish on the same de-spun section as the solar arrays, the dish only needs to track the earth in its orbit. Spinning in the plane, then, is the best solution for communications pointing on the Hab and ERV.

Thermal considerations

While in orbit, the vehicles experience heating from solar radiation. The best way to keep the vehicle inside a temperature threshold is to have it rotate such that the entire vehicle gets an even dosage of radiation. This method is referred to as the even-bake or rotisserie method. The other method is to maintain an orientation where the radiation is only heating a very well insulated and non-volatile area of the vehicle. With the ERA, which does not need to spin and does not have fragile humans, the entry heat shield can block most of the heat and some of the radiation from the volatile systems and the spent propulsion system is non-volatile and can block both the radiation and heat. The Hab and ERV are both spinning so it is easiest to use the even-bake method. Spinning out of the plane creates periods during the transfer where one side of the vehicle is getting the bulk of the radiation (fig. 3.4.6). Spinning in the plane allows for even distribution throughout the entire transfer.

Translational maneuvers

At the end of the missions for the Hab and ERV, the counterweight needs to be detached. The vehicles also need trajectory correctional translation maneuvers that are in the orbit plane. The method for achieving both of these goals is cutting the tether at a certain point in the spin that accomplishes the required maneuvers. This cutting, though, requires the vehicle to spin in the orbital plane. If the vehicle is out of the plane, it needs a propulsive maneuver to put it into this orientation. If the vehicle is spinning in the plane, no maneuvers are required.

3.4.6 Conclusions

We have shown that, due the large size of the vehicles and the low amount of torque on them, the three vehicles are stable in their nominal motions. We have also shown that, for the Hab and ERV, spinning in the orbit plane (fig. 3.4.1) is the best solution for communications and solar array pointing, thermal issues, and translational maneuvers when cutting the tether. Three-axis stabilization of the ERA is the best method for controlling the vehicle during its transfer with the initial orientation such that the axis of symmetry is perpendicular to the radius vector (fig. 3.4.8).

3.4.7 Acknowledgements

The basis for all of the equations and foundations for the work in this section came from the work and notes by Professor Howell for the course AAE490c. I would also like to personally thank Professor Howell for the extra assistance, guidance, and time given to me to finish this section.

References

1) Longuski, James M., Todd, Richard E., Kőnig, Wolfgang W., “Survey of Non-gravitational Forces and Space Environmental Torques: Applied to the Galileo,” Journal of Guidance, Control, and Dynamics, Vol.15, No.3, May-June 1992, pp. 545-553

3.4 Appendix

Chris Heidelberger

3.4.1 Matlab files

Run file

%run script for 450 by Chris Heidelberger

clear all

global k1;

global k2;

global k3;

global om;

global a1;

global e1;

global mew;

%-----------------------Hab

I1=2.1645e9;

I2=2.721e5;

I3=2.1645e9;

%-----------------------ERA

%I1=1.5e7;

%I2=4.12e6;

%I3=1.5e7;

%-----------------------ERV

%I1=5.467e9;

%I2=3.36e5;

%I3=5.467e9;

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

k1=(I2-I3)/I1;

k2=(I3-I1)/I2;

k3=(I1-I2)/I3;

mew=1.327e11;

re=1.49596e8;

ae=re;

ee=0;

ome=sqrt(mew/(ae^3));

%-----------------------Hab

r0=re;

a1=1.2992638*re;

e1=0.24863708;

om=sqrt(mew/(a1^3));

thst0=15.74*pi/180; %initial/final orbit locations (deg->rad)

thstf=143.36*pi/180;

Ea0=2*atan(sqrt((1-e1)/(1+e1))*tan(thst0/2));

Eaf=2*atan(sqrt((1-e1)/(1+e1))*tan(thstf/2));

w10=0*2*pi/60*om; %initial conditions (rpm->rad/s)

w20=0*2*pi/60*om; %non-dim. by omega=sqrt(mew/a3)

w30=2*2*pi/60*om;

psi0=0.5*pi/180;

lambda=[1,0,0];

e0=sin(psi0/2).*lambda;

e40=cos(psi0/2);

Ea1=Ea0:.01:Eaf;

coordx=(Ea1-Ea0)/(Eaf-Ea0).*100;

%-----------------------ERA

%r0=re;

%a1=1.2524363*re;

%e1=.21018772;

%om=sqrt(mew/(a1^3));

%w10=0*2*pi/60;

%w20=0*2*pi/60;

%w30=0*2*pi/60;

%thst0=352.71835*pi/180;

%thstf=201.00936*pi/180;

%Ea0=2*atan(sqrt((1-e1)/(1+e1))*tan(thst0/2));

%Eaf=2*atan(sqrt((1-e1)/(1+e1))*tan(thstf/2));

%psi0=0.5*pi/180;

%lambda=[1,0,0];

%e0=sin(psi0/2).*lambda;

%e40=cos(psi0/2);

%Ea1=Ea0:.01:Eaf+2*pi;

%coordx=(Ea1-Ea0)/(Eaf-Ea0+2*pi).*100;

%initial and final conditions for ode45

w01=[w10 w20 w30 e0(1) e0(2) e0(3) e40 r0 e0(1) e0(2) e0(3) e40];

options=odeset('RelTol',1e-8,'AbsTol',1e-8);

[EA,ww1]=ode45('elip_non_inert_eom',Ea1,w01,options);

w11=ww1(:,1);

w21=ww1(:,2);

w31=ww1(:,3);

e11=ww1(:,4);

e21=ww1(:,5);

e31=ww1(:,6);

e41=ww1(:,7);

r11=ww1(:,8);

e12=ww1(:,9);

e22=ww1(:,10);

e32=ww1(:,11);

e42=ww1(:,12);

w11r=w11.*60./2./pi.*om;

w21r=w21.*60./2./pi.*om;

w31r=w31.*60./2./pi.*om;

c111=1 - 2.*e31.^2 - 2.*e11.^2;

c131=2.*(e31.*e11+e21.*e41); %essential dir.cos elements

c231=2.*(e21.*e31-e11.*e41);

c311=2.*(e31.*e11-e21.*e41);

c321=2.*(e21.*e31+e11.*e41);

c331=1 - 2.*e11.^2 - 2.*e21.^2;

c112=1 - 2.*e32.^2 - 2.*e12.^2;

c132=2.*(e32.*e12+e22.*e42); %essential dir.cos elements

c232=2.*(e22.*e32-e12.*e42);

c312=2.*(e32.*e12-e22.*e42);

c322=2.*(e22.*e32+e12.*e42);

c332=1 - 2.*e12.^2 - 2.*e22.^2;

nut1=acos(c331);

prec1=atan2(c231,c131); %precession and spin angle from dircos

spin1=atan2(c321,-c311); %using atan2 fn. which does quadrant check

nutd1=nut1.*(180/pi); %convert to degrees

precd1=prec1.*(180/pi);

spind1=spin1.*(180/pi);

nut2=acos(c332);

prec2=atan2(c232,c132); %precession and spin angle from dircos

spin2=atan2(c322,-c312); %using atan2 fn. which does quadrant check

nutd2=nut2.*(180/pi); %convert to degrees

precd2=prec2.*(180/pi);

spind2=spin2.*(180/pi);

%r1=r11;

%len=85/1000; %length of tether to c.m. (m->km)

%thst1=2*atan(sqrt((1+e1)/(1-e1))*tan(EA./2));

%t0tp1=sqrt(a1^3/mew)*(Ea0-e1*sin(Ea0));

%t1tp1=sqrt(a1^3/mew).*(Ea1-e1.*sin(Ea1))-t0tp1;

%t0tpe=sqrt(ae^3/mew)*(thst0-ee*sin(thst0));

%thste=(ome.*(t1tp1+t0tpe))';

%beta=(thste-thst1)+0.0001;

%rec=(re^2+r1.^2-2*re.*r1.*cos(beta)).^(0.5);

%alpha=pi/2-acos((re^2-r1.^2-rec.^2)./(-2.*r1.*rec));

%epsi=-(prec1)+alpha;

%epsi2=pi-epsi;

%res=(len^2+rec.^2-2*len.*rec.*cos(epsi)).^(0.5);

%res2=(len^2+rec.^2-2*len.*rec.*cos(epsi2)).^(0.5);

%kappa=acos((rec.^2-len^2-res.^2)./(-2*len.*res)).*(180/pi);

%kappa2=acos((rec.^2-len^2-res2.^2)./(-2*len.*res2)).*(180/pi);

figure(1);

subplot(2,2,1);

hold on;

plot(coordx,nutd2);

ylabel('Tilt (deg)');

subplot(2,2,2);

hold on;

plot(coordx,precd2+90);

ylabel('Precession (deg)');

subplot(2,2,3);

hold on;

plot(coordx,spind2-90);

xlabel('% of transfer');

ylabel('Axial Spin (deg)');

subplot(2,2,4);

hold on;

plot(coordx,w31r,':',coordx,w21r,coordx,w11r,'--');

xlabel('% of transfer');

ylabel('Body Spin');

axis([0 100 -1 3]);

%figure(2);

%plot(coordx,kappa,coordx,kappa2,':');

%xlabel('% of transfer');

%ylabel('comm. angle');

%figure(3);

%hold on;

%plot(t1tp1,kappa-kappa2);

%xlabel('time');

%ylabel('difference in angle');

EOM file

%script to set eom into state variable form

function ydot = elip_non_inert_eom(Ea,y)

global k1;

global k2;

global k3;

global a1;

global e1;

global om;

global mew;

ydot(1)=(1-e1*cos(Ea))*k1*(y(2)*y(3) - 12*a1^3/y(8)^3*(y(4)*y(5)-y(6)*y(7))*(y(6)*y(4)+y(5)*y(7)));

ydot(2)=(1-e1*cos(Ea))*k2*(y(3)*y(1) - 6*a1^3/y(8)^3*(1-2*y(5)^2-2*y(6)^2)*(y(6)*y(4)+y(5)*y(7)));

ydot(3)=(1-e1*cos(Ea))*k3*(y(1)*y(2) - 6*a1^3/y(8)^3*(1-2*y(5)^2-2*y(6)^2)*(y(4)*y(5)-y(6)*y(7)));

ydot(4)=0.5*(1-e1*cos(Ea))*(y(1)*y(7) - y(2)*y(6) + (y(3)+sqrt(a1^3/y(8)^3))*y(5));

ydot(5)=0.5*(1-e1*cos(Ea))*(y(1)*y(6) + y(2)*y(7) - (y(3)+sqrt(a1^3/y(8)^3))*y(4));

ydot(6)=0.5*(1-e1*cos(Ea))*(-y(1)*y(5) + y(2)*y(4) + (y(3)-sqrt(a1^3/y(8)^3))*y(7));

ydot(7)=-0.5*(1-e1*cos(Ea))*(y(1)*y(4) + y(2)*y(5) + (y(3)-sqrt(a1^3/y(8)^3))*y(6));

ydot(8)=abs(a1*e1*sin(Ea));

ydot(9)=0.5*(1-e1*cos(Ea))*(y(1)*y(12) - y(2)*y(11) + (y(3))*y(10));

ydot(10)=0.5*(1-e1*cos(Ea))*(y(1)*y(11) + y(2)*y(12) - (y(3))*y(9));

ydot(11)=0.5*(1-e1*cos(Ea))*(-y(1)*y(10) + y(2)*y(9) + (y(3))*y(12));

ydot(12)=-0.5*(1-e1*cos(Ea))*(y(1)*y(9) + y(2)*y(10) + (y(3))*y(11));

ydot=ydot';

return

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download