Help with a rotation matrix

I have a code which requires 4 different rotation matrices, and since all my variables are in the form of arrays I'm trying to make these matrices also in the form of arrays. So far, the code I have below is what I believe should work. The parameters at the top are parameters for each of the 4 orbits, laid out in 1x4 array. R1, R2 and R3 are not giving me the answers I expect in all the cells (R1,2 and 3 will be combined to make R the final matrix). The matrices should have more zeros in them but I am getting a very very small number (10^-17) instead, I don't know where the maths is going wrong? Thanks! (.png attached at top).

clear
close

i = [90 90 90 90];                   %[deg] Inclinations for each orbit
o = [0 0 0 0];                       %[deg] Small omega of each orbit
O = [0 45 90 135];                   %[deg] RAANs for each orbit

R1 = [cos(o).*cos(O)-sin(o).*cos(i).*sin(O);         % Column 1 of rotation matrix for the 4 orbits
      cos(o).*sin(O)+sin(o).*cos(i).*cos(O);
      sin(o).*sin(i)];
R2 = [-sin(o).*cos(O)-cos(o).*cos(i).*sin(O);        % Column 2 of rotation matrix for the 4 orbits
      -sin(o).*sin(O)+cos(o).*cos(i).*cos(O);
      cos(o).*sin(i)];
R3 = [sin(i).*sin(O);                                % Column 3 of rotation matrix for the 4 orbits
      -sin(i).*cos(O);
      cos(i)];

[Image Analyst]

From Wikipedia this is a rotation matrix:

I'm not quite sure why your rotation matrices (R1, R2, and R4) have so many terms.

Plus I don't know how the i, o, and O are related to points in the plane that might be rotated through some angle. Why are they not Nx2 arrays of (x,y) coordinates?

A diagram would be helpful.

[Scott Hunter]

I have attached the formula for the rotation matrix at the top.

i is the angle of inclination of the plane with respect to the x-y plane, o and O are angles to do with the rotation of the plane about the z axis, however o is always zero in my case.

[Bruno Luong]

you seem to apply sin/cos to degree angles which is wrong, since argument must be in radian.

[Scott Hunter]

I initially had the argument in radians (just did deg2rad) and I was still getting incorrect results. I then read somewhere MATLAB uses degrees for trig. functions, is that incorrect?

[Scott Hunter]

Could it possibly be to do with the fact I have element-wise multiplication?

[Bruno Luong]

"I then read somewhere MATLAB uses degrees for trig. functions, is that incorrect?"

Yes incorrect.

"Could it possibly be to do with the fact I have element-wise multiplication?"

No, element wise is correct.

But take a step back : Why you think R2 and R3 are wrong?

[Scott Hunter]

Okay, I've reverted it to radians thank you.

I think R2 and R3 are incorrect because looking at the angles and the formula in the .png attached to my question, if we take R2(2,1) for example, using i, o and O from the parameters the answer should be 0, however when I view the variable I get 6.12323399573677e-17. granted this is almost 0, but not 0.

Edit: I've acutally just spotted R1 is also incorrect, for R1(1,3) I should get 0 but also get the same answer as above.

[Bruno Luong]

You should expect such small error when working with finite precision numbers

>> cos(pi/2)

ans =

     6.1232e-17

>>

[Scott Hunter]

wow okay, so really my code is fine? Thank you!

[James Tursa]

@IA: "... I'm not quite sure why your rotation matrices (R1, R2, and R4) have so many terms ..."

Because the formulas he is using are the result of three 3D rotation matrices multiplied together, the result of which is still a rotation matrix.

[James Tursa]

@Scott: For your downstream code, I would strongly advise arranging all of your columns into a 3D matrix that is size 3x3xN, where N is the number of rotation matrices that you are working with (in this case 4) and the 3x3 pages are the individual rotation matrices. That way the first two dimensions form pages of your rotation matrices, which makes it compatible with a lot of code (e.g. in the FEX) that has been written to deal with pages that are stored this way.

[Scott Hunter]

@James: By columns do you mean my R1, 2 and 3? I am not sure I understand how I would go about implementing that into the code I have to follow the rotation matrices. I only just learned how to use arrays for variables, I had my whole thing written out before (oops).

[James Tursa]

R1 contains four 1st columns, R2 contains four 2nd columns, and R3 contains four 3rd columns. I am suggesting that the first page of your 3D R array would be:

[R1(:,1) R2(:,1) R3(:,1)]

the second page of R would be

[R1(:,2) R2(:,2) R3(:,2)]

etc.

Rather that constructing R manually this way, one could pre-allocate R and then just fill it in appropriately. E.g.,

N = numel(i); % the number of rotation matrices
R = zeros(3,3,N); % pre-allocate the array of rotation matrices
R(:,1,:) = R1; % fill in the 1st columns of R
R(:,2,:) = R2; % fill in the 2nd columns of R
R(:,3,:) = R3; % fill in the 3rd columns of R

Then R(:,:,1) will be the first rotation matrix, R(:,:,2) will be the second rotation matrix, etc. In general, R(:,:,k) will be the k'th rotation matrix. This arrangement will pay dividends in your downstream code.

[Scott Hunter]

Oh right I see thank you. Why is this a benefit?

[James Tursa]

Many FEX submissions assume this layout so you could use them with your R array. E.g.,

https://www.mathworks.com/matlabcentral/fileexchange/27680-multiple-eigen-values-for-2x2-and-3x3-matrices

https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver

https://www.mathworks.com/matlabcentral/fileexchange/68976-multipleqr

https://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver

https://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

https://www.mathworks.com/matlabcentral/fileexchange/37515-mmx

https://www.mathworks.com/matlabcentral/fileexchange/8773-multiple-matrix-multiplications-with-array-expansion-enabled

And several MATLAB functions also assume this layout. E.g.,

https://www.mathworks.com/help/aerotbx/ug/dcm2quat.html

https://www.mathworks.com/help/aerotbx/ug/quat2dcm.html

https://www.mathworks.com/help/aerotbx/ug/dcm2angle.html

https://www.mathworks.com/help/aerotbx/ug/angle2dcm.html

Answer 1

[Bruno Luong]

You seem to apply sin/cos to degree angles which is wrong, since argument must be in radian.

Also the code is correct and might returns small error when expected result is 0.

   >> cos(pi/2)

    ans =

       6.1232e-17

PS: I put my comment here for you to accept if it solves your pb

[Scott Hunter]

I will do! Maybe add in that my confusion was due to the error in accuracy?

[Image Analyst]

You can use degrees. Just switch to the "d" versions of the trig functions: sind() and cosd().

Answer 2

[sanjay]

cos(pi/2)

ans =

6.1232e-17

[Steven Lord]

That's correct. If you want to compute the cosine or sine of an angle represented as a multiple of pi radians, you can use cospi and sinpi.

cos(pi/2)

ans = 6.1232e-17

cospi(1/2)

ans = 0

sin(pi)

ans = 1.2246e-16

sinpi(1)

ans = 0

[Stephen23]

Or use degrees:

cosd(90)

ans = 0

sind(180)

ans = 0