Model Motion of Circling Airplane

Start with an airplane moving at 150 kmph in a circle of radius 10 km and descending at the same time at a rate of 20 m/sec. Compute the motion of the airplane from its instantaneous acceleration as an argument to the step method. Set the initial orientation of the platform to the identity, coinciding with the global coordinate system.

Set up the scenario

Specify the initial position and velocity of the airplane. The airplane has a ground range of 10 km and an altitude of 20 km.

range = 10000;
alt = 20000;
initPos = [cosd(60)*range;sind(60)*range;alt];
originPos = [1000,1000,0]';
originVel = [0,0,0]';
vs = 150.0;
phi = atan2d(initPos(2)-originPos(2),initPos(1)-originPos(1));
phi1 = phi + 90;
vx = vs*cosd(phi1);
vy = vs*sind(phi1);
initVel = [vx,vy,-20]';
platform = phased.Platform('MotionModel','Acceleration',...
    'AccelerationSource','Input port','InitialPosition',initPos,...
    'InitialVelocity',initVel,'OrientationAxesOutputPort',true,...
    'InitialOrientationAxes',eye(3));
relPos = initPos - originPos;
relVel = initVel - originVel;
rel2Pos = [relPos(1),relPos(2),0]';
rel2Vel = [relVel(1),relVel(2),0]';
r = sqrt(rel2Pos'*rel2Pos);
accelmag = vs^2/r;
unitvec = rel2Pos/r;
accel = -accelmag*unitvec;
T = 0.5;
N = 1000;

Compute the trajectory

Specify the acceleration of an object moving in a circle in the x-y plane. The acceleration is v^2/r towards the origin.

posmat = zeros(3,N);
r1 = zeros(N);
v = zeros(N);
for n = 1:N
    [pos,vel,oax] = platform(T,accel);
    posmat(:,n) = pos;
    vel2 = vel(1)^2 + vel(2)^2;
    v(n) = sqrt(vel2);
    relPos = pos - originPos;
    rel2Pos = [relPos(1),relPos(2),0]';
    r = sqrt(rel2Pos'*rel2Pos);
    r1(n) = r;
    accelmag = vel2/r;
    accelmag = vs^2/r;
    unitvec = rel2Pos/r;
    accel = -accelmag*unitvec;
end

Display the final orientation of the local coordinate system.

disp(oax)
   -0.3658   -0.9307   -0.0001
    0.9307   -0.3658   -0.0010
    0.0009   -0.0005    1.0000

Plot the trajectory and the origin position

posmat = posmat/1000;
figure(1)
plot3(posmat(1,:),posmat(2,:),posmat(3,:),'b.')
hold on
plot3(originPos(1)/1000,originPos(2)/1000,originPos(3)/1000,'ro')
xlabel('X (km)')
ylabel('Y (km)')
zlabel('Z (km)')
grid
hold off