Question on Aircraft Drone Modeling - Simulation - Data Selection and Entry

Okay, I'll be blunt. I'm trying to learn about aircraft aerodynamics by working through a textbook on drone design. I have no instructor, and I'm 30 years distant from my own college degree, so I've forgotten a step or two. That means when matlab has a segment of code that say that the term that is added is m[J K L] (vertically aligned) I'm at something of a loss regarding what I'm adding. Do I add m(J) to the topmost line, then add m(K) to the second line equation, and m(L) to the bottom line?

The first part assignment (Chapter 2) was to produce a model of the an aircraft in Matlab and get it to move on the screen. I did this. The code is below:

The pn, pe, and pd variables are apparently location variables. The dpsi is a change in the horizontal orientation of the aircraft (it turns right or left). The pd = pd - 1 threw me for a moment, until I realized what I'd done was to increment the direciton of motion with a negative one to make the plane move in the direction it is pointing. I may need to correct this. The following segment of code has no tie to any real estimate of forces acting on the drone. It is purely my data producing an aircraft that moves on the scree.

pn=0;

pe=0;

pd=0;

phi=0;

theta=0;

psi=0;

handle=[];

%Simulation parameters

dpn=0.5; %Change in position.

dpsi=0.01; %Change in yaw angle.

simlength=125;%Number of Sim Steps

%Draw and update the plane's position.

for k = 1:simlength

pn = 75+dpn;

pe = pe+dpn;

%THE CODE THAT PRODUCES A PROBLEM IS BELOW

pd = pd -1;

psi = psi+dpsi;

%Draw or update the plane.

handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle);

pause(0.1); %Pause to visualize the movement

end

function handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle)

%define points on plane in local NED coordintates

NED = airplanepoints;

%rotate plane by (phi; theta; psi)

NED = rotate(NED, phi, theta, psi);

%translate plane to [pn; pe; pd]

NED = translate(NED, pn, pe, pd);

%transform vertices from NED to XYZ

R = [...

0, 1, 0;...

1, 0 , 0;...

0, 0, -1;...

];

XYZ = R* NED';

%plot plane

if isempty(handle)

handle = plot3(XYZ(1,:),XYZ(2,:),XYZ(3,:),'b');

xlim([-100 150]);

ylim([-100 150]);

zlim([0 150]);

grid on;%show grid

xlabel('X');ylabel('Y');zlabel('Z');%label axes

else

set(handle, 'XData', XYZ(1,:), 'YData',XYZ(2,:), 'ZData',XYZ(3,:));

drawnow

end

end

function XYZ = airplanepoints

%define points on the aircraft in local NED

XYZ = [...

0 0 0;%point1

-2 1 1;%point2

-2 1 -1;%point3

0 0 0;%point1

-2 -1 1;%point4

-2 -1 -1;%point5

0 0 0;%point1

-2 1 1;%point2

-2 -1 1;%point4

0 0 0;%point1

-2 1 -1;%point3

-2 -1 -1;%point5

0 0 0;%point1

-2 1 1;%point2

-18 0 0;%point6

-2 1 -1;%point3

-18 0 0;%point6

-2 -1 -1;%point5

-18 0 0;%point6

-2 -1 1;%point4

-18 0 0;%point6

-2 1 1;%point2

0 0 0;%point1

-5 0 0;%point7

-5 -10 0;%point8

-8 -10 0;%point9

-8 10 0;%point10

-5 10 0;%point11

-5 0 0;%point7

-15.5 0 0;%point12

-15.5 2 0;%point13

-17.5 2 0;%point14

-17.5 -2 0;%point15

-15.5 -2 0;%point16

-15.5 0 0;%point12

-18 0 0;%point6

-18 0 -3;%point17

-15.5 0 0;%point15

-18 0 0;%point16

];

end

function XYZ=rotate(XYZ,phi,theta,psi)

%define rotation matrix

R_roll = [

1, 0, 0;

0, cos(phi), -sin(phi);

0, sin(phi), cos(phi)];

R_pitch = [

cos(theta), 0, sin(theta);

0, 1, 0;

-sin(theta), 0, cos(theta)];

R_yaw = [

cos(psi), -sin(psi), 0;

sin(psi), cos(psi), 0;

0, 0, 1];

R = R_roll*R_pitch*R_yaw;

%rotate vertices

XYZ =R*XYZ';

XYZ = XYZ';

end

function XYZ = translate(XYZ, pn, pe, pd)

XYZ = XYZ' + repmat([pn;pe;pd],1,size(XYZ,1));

XYZ=XYZ';

end

Alright, this is the end of the purely "draw a drone" segment ofthe code. The book then says to produce code that will permit the aircraft to respond to real world forces. That code follows. I haven't even tried to run it. I'm not sure that I even have the right variables transferred to the functions. The equations generally serve to translate forces into a reference frame that is relevant to the aircraft. I haven't check which reference frame is the starting one and which is the final one. I leave the movement code alone and maintain the same size box as in the prior simulation. I'm hoping someone wtho knows something about this book, Small, Unmanned Aircraft Theory and Practice will find this and assist. I need some instructor level guidance. Maybe theres an aircraft designer who understands this subject. I'll post this and hope. (The G1, G2, G3, etc. in the following code is a capital gamma with those numbers afterward in the text. There's also a Jy at the end of the line of data that is transferred. I don't know what that is. Does there need to be a Jx and a Jz transferred as well?)

pn=0;

pe=0;

pd=0;

phi=0;

theta=0;

psi=0;

handle=[];

%Simulation parameters

dpn=0.5; %Change in position.

dpsi=0.01; %Change in yaw angle.

%SIMULATION

u = 1;

v = 1;

w = 1;

theta = 3.1415/8;

psi = 3.1415/9;

phi = 3.1415/10;

l = 2;

m =3;

n=4;

p = 5;

q = 6;

r = 7;

G1 = 1;

G2 = 1.1;

G3 = 1.2;

G4 = 1.3;

G5 = 1.4;

G6 = 1.5;

G7 = 1.6;

G8 = 1.7;

Jy = 2;

Mavequations(u,v,w,theta,psi,phi,l,m,n,p,q,r,G1,G2,G3,G4,G5,G6,G7,G8,Jy);

%simlength=125;%Number of Sim Steps

%Draw and update the plane's position.

%for k = 1:simlength

% pn = 75+dpn;

% pe = pe+dpn;

%THE CODE THAT PRODUCES A PROBLEM IS BELOW

% pd = pd -1;

% psi = psi+dpsi;

% %Draw or update the plane.

% handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle);

% pause(0.1); %Pause to visualize the movement

%

function MAVEquations(u,v,w,theta,psi,phi,l,m,n,p,q,r,G1,G2,G3,G4,G5,G6,G7,G8,Jy)

pn = cos(theta)*cos(psi)*sin(phi)*sin(theta)*cos(psi)*u-cos(phi)*sin(psi)*cos(phi)*sin(theta)*cos(psi)*v+sin(phi)*sin(psi)*w;

pe = cos(theta)*sin(psi)*sin(phi)*sin(theta)*sin(psi)*u +cos(phi)*cos(psi)*cos(phi)*sin(theta)*sin(psi)*v+sin(phi)*cos(psi)*w;

pd = -sin(theta)*u + sin(phi)*cos(theta)*v + cos(phi)*cos(theta)*w;

ud = r*v-q*w +(fx+fy+fz)/m;

vd = p*w-r*u +(fx+fy+fz)/m;

wd = q*u-p*v +(fx+fy+fz)/m;

phid =1*p + sin(phi)*tan(theta)*q + cos(phi)*cos(theta)*r;

thetad = 0*p + cos(phi)*q - sin(phi);

psid = 0 * p + sin(phi)/cos(theta) *q + cos(phi)/cos(theta)*r;

pd = G1*p*q-G2*q*r + G3*l+G4*n;

qd = G5*p*r-G6(p^2-r^2)+m/Jy;

rd = G7*p*q-G1*q*r+G4*l+G8*n;

handle = drawPlanBody(pn,pe,pd,phi,theta,psi,handle);

pause(0.1);

end

function handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle)

%define points on plane in local NED coordintates

NED = airplanepoints;

%rotate plane by (phi; theta; psi)

NED = rotate(NED, phi, theta, psi);

%translate plane to [pn; pe; pd]

NED = translate(NED, pn, pe, pd);

%transform vertices from NED to XYZ

R = [...

0, 1, 0;...

1, 0 , 0;...

0, 0, -1;...

];

XYZ = R* NED';

%plot plane

if isempty(handle)

handle = plot3(XYZ(1,:),XYZ(2,:),XYZ(3,:),'b');

xlim([-100 150]);

ylim([-100 150]);

zlim([0 150]);

grid on;%show grid

xlabel('X');ylabel('Y');zlabel('Z');%label axes

else

set(handle, 'XData', XYZ(1,:), 'YData',XYZ(2,:), 'ZData',XYZ(3,:));

drawnow

end

end

function XYZ = airplanepoints

%define points on the aircraft in local NED

XYZ = [...

0 0 0;%point1

-2 1 1;%point2

-2 1 -1;%point3

0 0 0;%point1

-2 -1 1;%point4

-2 -1 -1;%point5

0 0 0;%point1

-2 1 1;%point2

-2 -1 1;%point4

0 0 0;%point1

-2 1 -1;%point3

-2 -1 -1;%point5

0 0 0;%point1

-2 1 1;%point2

-18 0 0;%point6

-2 1 -1;%point3

-18 0 0;%point6

-2 -1 -1;%point5

-18 0 0;%point6

-2 -1 1;%point4

-18 0 0;%point6

-2 1 1;%point2

0 0 0;%point1

-5 0 0;%point7

-5 -10 0;%point8

-8 -10 0;%point9

-8 10 0;%point10

-5 10 0;%point11

-5 0 0;%point7

-15.5 0 0;%point12

-15.5 2 0;%point13

-17.5 2 0;%point14

-17.5 -2 0;%point15

-15.5 -2 0;%point16

-15.5 0 0;%point12

-18 0 0;%point6

-18 0 -3;%point17

-15.5 0 0;%point15

-18 0 0;%point16

];

end

function XYZ=rotate(XYZ,phi,theta,psi)

%define rotation matrix

R_roll = [

1, 0, 0;

0, cos(phi), -sin(phi);

0, sin(phi), cos(phi)];

R_pitch = [

cos(theta), 0, sin(theta);

0, 1, 0;

-sin(theta), 0, cos(theta)];

R_yaw = [

cos(psi), -sin(psi), 0;

sin(psi), cos(psi), 0;

0, 0, 1];

R = R_roll*R_pitch*R_yaw;

%rotate vertices

XYZ =R*XYZ';

XYZ = XYZ';

end

function XYZ = translate(XYZ, pn, pe, pd)

XYZ = XYZ' + repmat([pn;pe;pd],1,size(XYZ,1));

XYZ=XYZ';

end

Yes, this is a Hail Mary appeal. I hope that the right person will find it and help me to develop an answer. Thank you very much.