Engineering Lab report - Matlab

matlab_files/Euler_Digitizer.m

% BIOEN 4250 FALL 2018 LAB 2 - Rigid Body Kinematics % PROF. Lucas Timmins % TAs: Allen Lin and Kelly Smith % Code Written by E. David Bell, modified from code developed by TREVOR LUJAN % ########### RIGID BODY KINEMATICS ############### % STUDENT NAME : % STUDENT ID : % STUDENT GROUP : % Clear all stored variables and set output format to short, % which is a 5 digit format clear all; format short; % Load the file containing the digitized points for Position 1. This 3 % column matrix will be used to determine the transformations between the % marker clusters and embedded systems, to be inputted in the overall % transformation matrix. D1 = xlsread('Position 1.xlsx'); % Load the file containing the digitized points for Position 2. This 3 % column matrix will be used to determine the transformations between the % marker clusters specific to Position 2. D2 = xlsread('Position 2.xlsx'); % Load the file containing the digitized points for Position 3. This 3 % column matrix will be used to determine the transformations between the % marker clusters specific to Position 3. D3 = xlsread('Position 3.xlsx'); % FORMULATE THE COORDINATE SYSTEM FOR THE KINEMATIC MARKERS % "KINEMATIC MARKERS" ARE THE DIVITS ON THE L SHAPED PLATES THAT ARE TRACKED % AT EACH POSITION - ONE SET FOR FEMUR (upper section)AND ONE SET FOR TIBIA % (lower section). THEY ARE REFFERED TO AS "MARKERS" IN ORDER TO % DISTINGUISH THEM FROM THE POINTS USED TO MAKE THE EMBEDED COORDINATE % SYSTEMS (THOSE DIVITS ON THE MAIN PIECES OF THE FIXTURE)AND BECAUSE THE % SAME ANALYSIS PROCEDURE WOULD BE DONE HAD THIS DATA BEEN AQUIRED FROM % VISUAL TRACKING SOFTWARE. % *********** ESTABLISH ORTHONORMAL COORDINATE SYSTEM ON FEMUR USING KINEMATIC MARKERS **************** %%%%%%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % dft: digitized femur top : REFERS TO THE DIVIT AT THE TOP OF THE L PIECE INITIALLY % dfb: digitized femur bottom : REFERS TO THE MIDDLE DIVIT ON THE L PIECE % dff: digitized femur front : REFERS TO THE DIVIT ON THE L PIECE CLOSEST TO THE JOINT % e_dfx: base vector along x digitized femur % (2 PTS) C: (What does the following 3 lines of code do?) dft = D1(1,1:3); dfb = D1(2,1:3); dff = D1(3,1:3); % (4 PTS) C: (What does the following do (2 pts) and specifically, why are %the cross products taken? (2 pts)) dfx = dff - dfb; tempdfy = dft - dfb; dfz = cross(dfx,tempdfy); dfy = cross(dfz,dfx); % (2 PTS) C: (What does the following code do?) e_dfx = dfx/mag(dfx); e_dfy = dfy/mag(dfy); e_dfz = dfz/mag(dfz); % *********** ESTABLISH ORTHONORMAL COORDINATE SYSTEM ON TIBIA USING KINEMATIC MARKERS **************** %load tibia kinematic marker coordinates %%%%%%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % dtt: digitized tibia top : REFERS TO THE DIVIT AT THE TOP OF THE L PIECE INITIALLY % dtb: digitized tibia bottom : REFERS TO THE MIDDLE DIVIT ON THE L PIECE % dtf: digitized tibia front : REFERS TO THE DIVIT ON THE L PIECE CLOSEST TO THE JOINT % e_dtx: base vector along x digitized tibia % (2 PTS) C: (How is this section of the code different from the similar % section above?) dtt = D1(4,1:3); dtb = D1(5,1:3); dtf = D1(6,1:3); dtx = dtf - dtb; tempdty = dtt - dtb; dtz = cross(dtx,tempdty); dty = cross(dtz,dtx); e_dtx = dtx/mag(dtx); e_dty = dty/mag(dty); e_dtz = dtz/mag(dtz); % *********** ESTABLISH ORTHONORMAL REFERENCE COORDINATE SYSTEM FOR TROUBLESHOOTING**************** %This section of the code is commented out because it is not normally %needed, but please still address the one comment requested below. %%%%%%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % drt - digitized reference top: Top back corner of the mounting frame adjacent to the "femur" % drb - digitized reference bottom: Bottom back corner of the mounting frame % drf - digitized reference front: Bottom front corner of the mounting frame % e_drx: base vector along x digitized reference % C: (6 pts) (Explain where the reference system was taken (2 pts) and how it would % be usefull for trouble shooting (4 pts). ) % drt = D1(7,1:3); % drb = D1(8,1:3); % drf = D1(9,1:3); % % drz = drt - drb; % tempdry = drf - drb; % drx = cross(tempdry,drz); % dry = cross(drz,drx); % % e_drx = drx/mag(drx); % e_dry = dry/mag(dry); % e_drz = drz/mag(drz); % ****** EMBEDDED COORDINATE SYSTEMS IN FEMUR AND TIBIA [DIGITIZED] ******* % NOTE: The mechanical surrogate "knee" is considered to be a left knee in % relation to which side is medial vs lateral. Therefore, the side with the % mounting bracket attached is the lateral side, and the side with the % "markers" attached is the medial side. %%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %e: embedded, m: medial, l: lateral, o: origin, d: distal %a: anterior, p: posterior f: femur, t: tibia, Prox: proximal %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % *********** EMBEDDED FEMUR (femur = upper segment of fixture) *********** % Load embedded femur, lateral coordinates efl = D1(11,1:3); % Load embedded femur, medial efm = D1(10,1:3); % Create a vector in the femoral embedded medial direction following the % coordinate axis conventions in the report instructions efx = efm-efl; % (2 PT) C: (Which points are loaded in the following two lines of code?) efda = D1(13,1:3); efdp = D1(12,1:3); % (4 PTS) C: (Explain why a mean is taken) A = [efdp;efda]; efd = mean(A); % (2 PT) C: (What data points are loaded in the following 2 lines of code?) efaProx = D1(15,1:3); efpProx = D1(14,1:3); % (4 PTS) C:(Explain why a mean is taken) A = [efaProx;efpProx]; efProx = mean(A); % (4 PTS) C: (What does "temp_efz" represent, AND why was it not %calculated directly from digitized points like has been done in other %parts of the code?) temp_efz = efProx-efd; efy = cross(efx, temp_efz); efz = cross(efx,efy); % (2 PTS) C: (What does the following code do?) e_efx = efx/mag(efx); e_efy = efy/mag(efy); e_efz = efz/mag(efz); % *********** EMBEDDED TIBIA (tibia = lower segment of fixture)************ % SAME NAMING CONVENTIONS AS EMBEDDED FEMUR % etl = D1(17,1:3); etm = D1(16,1:3); etx = etm-etl; etaProx = D1(19,1:3); etpProx = D1(18,1:3); A = [etaProx;etpProx]; etProx = mean(A); etda = D1(21,1:3); etdp = D1(20,1:3); A = [etda;etdp]; etd = mean(A); temp_etz = etProx-etd; ety = cross(etx, temp_etz); etz = cross(etx,ety); e_etx = etx/mag(etx); e_ety = ety/mag(ety); e_etz = etz/mag(etz); % *********** MARKER TO EMBEDDED TRANSFORMATION MATRICES **************** % (4 PTS) C:(What does the following function produce, {i.e. what is %"T_fe_fm")(2 pts), and what it will be used for? (2 pts)) T_fe_fm = transform(e_efx, e_efy, e_efz,e_dfx, e_dfy, e_dfz,efd,dfb); % (4 PTS) C: (What does the following function produce, {i.e. what is %"T_tm_te")(2 pts), and what it will be used for? (2 pts)) T_tm_te = transform(e_dtx, e_dty, e_dtz, e_etx, e_ety, e_etz, dtb, etProx); %%% PROCESS "MARKER" MOTION TRACKING DATA AND CALCULATE THE EULER ANGLES%%% %%%%%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%%%% %mff: marker femur front mtf: marker tibia front %mft: marker femur top mtt: marker tibia top %mfb: marker femur bottom mtb: marker tibia bottom %e_mfx: base vector along x marker femur %e_mtz: base vector along z marker tibia %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (4 PTS) C: (In general, what does each cycle of the following loop do?) i=1; for i=1:1:3 if i==1 M = D1; end if i==2 M = D2; end if i==3 M = D3; end % (2 PTS) C: (What is being loaded in the following 3 lines of code?) mft = M(1,:); mfb = M(2,:); mff = M(3,:); % (2 PTS) C: (What do the following lines of code do, and why is the %cross product taken?) mfx = mff - mfb; temp_mfy = mft - mfb; mfz = cross(mfx,temp_mfy); mfy = cross(mfz,mfx); % (2 PTS) C: (What does the following code do?) e_mfx = mfx/mag(mfx); e_mfy = mfy/mag(mfy); e_mfz = mfz/mag(mfz); % (4 PTS) C:(In general, what does the following block of code do? {up %the next comment} mtt = M(4,:); mtb = M(5,:); mtf = M(6,:); mtx = mtf - mtb; temp_mty = mtt - mtb; mtz = cross(mtx,temp_mty); mty = cross(mtz,mtx); e_mtx = mtx/mag(mtx); e_mty = mty/mag(mty); e_mtz = mtz/mag(mtz); % (4 PTS) C: (What does the following function produce, {i.e. what is %"T_fm_tm")(2 pts), and what it will be used for? (2 pts)) T_fm_tm = transform(e_mfx, e_mfy, e_mfz, e_mtx, e_mty, e_mtz, mfb, mtb); % (4 PTS) C: (What does the following function produce, {i.e. what is %"T_fe_te")(2 pts), and what it will be used for? (2 pts)) T_fe_te = [T_tm_te]*[T_fm_tm]*[T_fe_fm]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following set of "if" statements simply saves the values of % T_fm_tm and T_fe_te as seperate variables so their values from each % loop iteration will be available to the user via the command window % after the program finishes. This variables should be useful for the % manual calculation of the contents of Table 1 in the lab report. if i==1 T_fm_tm1 = T_fm_tm; T_fe_te1 = T_fe_te; end if i==2 T_fm_tm2 = T_fm_tm; T_fe_te2 = T_fe_te; end if i==3 T_fm_tm3 = T_fm_tm; T_fe_te3 = T_fe_te; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% CALCULATE EULER ANGLES FROM EMBEDED TRANSFORMATIONS %%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% Naming Conventions %%%%%%%%%%%%%%%%%%%%%%%%%% %gamma#= rotation about x-axis(rad) gamma#Deg= rotation about x-axis(deg) %theta#= rotation about y-axis(rad) theta#Deg= rotation about y-axis(deg) %phi# = rotation about z-axis(rad) phi#Deg= rotation about z-axis(deg) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NOTE: Remember, each transformation has 2 complete sets of Euler % angles (thus gamma1 is from the first set, and gamma2 is from the % second set). % (4 PTS) C: (In general, what does the code from here to the end of % the loop do?) T=T_fe_te; %(2 PTS) C: (What does this following 3 lines calculate?) theta1=-asin(T(3,1)); gamma1=atan2((T(3,2)/cos(theta1)),(T(3,3)/cos(theta1))); phi1=atan2((T(2,1)/cos(theta1)),(T(1,1)/cos(theta1))); %(2 PTS) C: (Why are a second set of calculations performed in the %following 3 lines that is so similar to that in the 3 lines above?) theta2=pi-theta1; gamma2=atan2((T(3,2)/cos(theta2)),(T(3,3)/cos(theta2))); phi2=atan2((T(2,1)/cos(theta2)),(T(1,1)/cos(theta2))); %(2 PTS) C: (What are the two things accomplished by the following %6 lines of code?) theta1Deg(i) = theta1*(180/pi); gamma1Deg(i) = gamma1*(180/pi)'; phi1Deg(i) = phi1*(180/pi); theta2Deg(i) = theta2*(180/pi); gamma2Deg(i) = gamma2*(180/pi); phi2Deg(i) = phi2*(180/pi); end % (1 pt) C: (What is the following variable use for?) Position = [1, 2, 3]; % (5 PTS) C: (Provide an explanation of the subplots) figure(1); title('Embedded Euler Angle Plots'); orient landscape; subplot(2,3,1) plot(Position,gamma1Deg) title ('Gamma-1 (rotation about x-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); subplot(2,3,2) plot(Position,theta1Deg) title ('Theta-1 (rotation about y-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); subplot(2,3,3) plot(Position,phi1Deg) title ('Phi-1 (rotation about z-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); subplot(2,3,4) plot(Position,gamma2Deg) title ('Gamma-2 (rotation about x-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); subplot(2,3,5) plot(Position,theta2Deg) title ('Theta-2 (rotation about y-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); subplot(2,3,6) plot(Position,phi2Deg) title ('Phi-2 (rotation about z-axis)'),ylabel('Angle(Deg)'),xlabel('Position Number'); %For clarification, in the following 2 lines, "F1" is the handle the %program uses to refer to figure(1), and the set function, is defining %the property "Position"(which is how the program defines where, and how %large the figure will be on the screen). Feel free to adjust this if you %wish to get the figure to be the size you want. Dave (your TA) will be %happy to explain how to do this, during office hours, if you want %(along with any other coding questions you may have for this lab). F1 = gcf; set(F1,'Position', [19 187 1284 477]) %"SUPTITLE" Adds a master title to the entire figure above all the other %subplots. Dave got this subfunction on a webpage where people can submit %specialized Matlab codes they wrote themselves. %YOU NEED TO USE THIS SUBFUNCTION TO GET THIS CODE TO WORK, BUT YOU DO NOT %NEED TO INCLUDE IT IN YOUR APPENDIX %Please remember to update this title to designate the figure as being from %your specific group data, or the TA provided data. suptitle('Embedded Euler Angle Plots') message='NORMAL TERMINATION!'