Tip:
Highlight text to annotate it
X
continued... % Original Scripts ©BiomEng 2006-2013 [Embed as a video - Do not use my script without a reference]
% Biomechanics - Objects Moving Through Air (Ball bearing) % ACCELERATION X % Ax [N3 N4]=size(A_BB1(:,1)); subplot(1,1,1);plot(T2,A_BB1(:,1),'-r',T2,FA_BB1(1:N3,1),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') legend('Ax-Raw','Ax-Filtered',2,'Location','Best'); %subplot(2,1,2);plot(T2,FA_BB1(1:N3,1),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') %legend('Ax-Filtered',1,'Location','Best');
figure for dj1=(1:N3)'; EAx1=max(abs(A_BB1(dj1,1)-FA_BB1(dj1,1))); EAx2a=sum(abs(A_BB1(dj1,1)-FA_BB1(dj1,1))); EAx3a=sqrt(sum(abs(A_BB1(dj1,1)-FA_BB1(dj1,1)).^2)); EAx2=EAx2a/N3; EAx3=EAx3a/N3; end
%ACCELERATION Y % Ay [N3 N4]=size(A_BB1(:,1)); subplot(1,1,1);plot(T2,-A_BB1(:,2),'-r',T2,-FA_BB1(1:N3,2),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') legend('Ay-Raw','Ay-Filtered',2,'Location','Best'); %subplot(2,1,2);plot(T2,-FA_BB1(1:N3,2),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') %legend('Ay-Filtered',1,'Location','Best');
figure for dj1=(1:N3)'; EAy1=max(abs(A_BB1(dj1,2)-FA_BB1(dj1,2))); EAy2a=sum(abs(A_BB1(dj1,2)-FA_BB1(dj1,2))); EAy3a=sqrt(sum(abs(A_BB1(dj1,2)-FA_BB1(dj1,2)).^2)); EAy2=EAy2a/N3; EAy3=EAy3a/N3; end
% ACCELERATION Z % Az [N3 N4]=size(A_BB1(:,1)); subplot(1,1,1);plot(T2,A_BB1(:,3),'-r',T2,FA_BB1(1:N3,3),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') legend('Az-Raw','Az-Filtered',2,'Location','Best'); %subplot(2,1,2);plot(T2,FA_BB1(1:N3,3),'-b','Linewidth',2); xlabel('Time (msec)'); ylabel('Acceleration (m/sec^2)') %legend('Az-Filtered',1,'Location','Best');
figure for dj1=(1:N3)'; EAz1=max(abs(A_BB1(dj1,3)-FA_BB1(dj1,3))); EAz2a=sum(abs(A_BB1(dj1,3)-FA_BB1(dj1,3)))/N3; EAz3a=sqrt(sum(abs(A_BB1(dj1,3)-FA_BB1(dj1,3)).^2)/N3); EAz2=EAz2a/N3; EAz3=EAz3a/N3; end
% Times tqs2=(t(T1(1,1)-1:T1(1,2)+1,1)); tqs1=tqs2-tqs2(1,1);
% QUINTIC SPLINES % Velocity X title('DX-VX-ERROR-Q SPLINE','FontWeight','bold'); splinetool(tqs1,BBXCD) subplot(3,1,1); DX=plot(tqs1,BBXCD,'-b','Linewidth',2); spDX=spaps(tqs1,BBXCD,1.9304e-007,3); legend('QSPL-dX',1,'Location','Best'); xlabel('Time (sec)'); ylabel('Displacement (m)')
subplot(3,1,2); spVX=fnder(spDX); fnplt(spVX,2); legend('QSPL-VX',1,'Location','Best'); xlabel('Time (sec)') ylabel('Velocity (m/sec)') subplot(3,1,3); plot(tqs1,BBXCD-fnval(spDX,tqs1),'Linewidth',2); legend('error-QSpl for Vx',1,'Location','Southeast') xlabel('Time (sec)') ylabel('Error in QSpl for Vx') figure
% QUINTIC SPLINES % Velocity Y subplot(3,1,1); title('Dy-Vy-ERROR-Q SPLINE','FontWeight','bold'); subplot(3,1,1); DY=plot(tqs1,-BBYCD,'-b','Linewidth',2); spDY=spaps(tqs1,-BBYCD,1.9304e-007,3); legend('QSPL-dY',1,'Location','Best'); xlabel('Time (sec)'); ylabel('Displacement (m)')
subplot(3,1,2); spVY=fnder(spDY); fnplt(spVY,2); legend('QSPL-VY',1,'Location','Best'); xlabel('Time (sec)');ylabel('Velocity (m/sec)') subplot(3,1,3); plot(tqs1,-BBYCD-fnval(spDY,tqs1),'Linewidth',2); legend('error-QSpl for Vy',1,'Location','Southeast') xlabel('Time (sec)');ylabel('Error in QSpl for Vy') figure
% QSPLINE ACCELERATION+VELOCITY % X AXIS title('Ax-QSplines','FontWeight','bold'); subplot(2,1,1); fnplt(spVX); legend('Vx-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Velocity (m/sec)') subplot(2,1,2); spAX=(fnder(spVX)); fnplt(spAX); legend('Ax-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Acceleration (m/sec^2)')
% QSPLINE ACCELERATION+VELOCITY % Y AXIS title('Ay-QSplines','FontWeight','bold'); subplot(2,1,1); fnplt(spVY); legend('Vy-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Velocity (m/sec)') subplot(2,1,2); spAY=(fnder(spVY)); fnplt(spAY); legend('Ay-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Acceleration (m/sec^2)') 0:02:10.000,0:02:20.000 % Velocity Z title('DZ-VZ-ERROR-Q SPLINE','FontWeight','bold'); subplot(3,1,1); DZ=plot(tqs1,BBZCD+cor/1000,'-b','Linewidth',2); spDZ=spaps(tqs1,BBZCD,1.9304e-007,3); legend('QSPL-dZ',1,'Location','Best'); xlabel('Time (sec)'); ylabel('Displacement (m)') subplot(3,1,2); spVZ=fnder(spDZ); fnplt(spVZ,2);
legend('QSPL-VZ',1,'Location','Best'); xlabel('Time (sec)') ylabel('Velocity (m/sec)') subplot(3,1,3); plot(tqs1,BBZCD-fnval(spDZ,tqs1),'Linewidth',2); legend('error-QSpl for Vz',1,'Location','Southeast') xlabel('Time (sec)') ylabel('Error in QSpl for Vz') figure
% QSPLINE ACCELERATION+VELOCITY % Z AXIS title('Az-QSplines','FontWeight','bold'); subplot(2,1,1); fnplt(spVZ); legend('Vz-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Velocity (m/sec)') subplot(2,1,2); spAZ=(fnder(spVZ)); fnplt(spAZ); legend('Az-QS',1,'Location','Best'); xlabel('Time (sec)') ylabel('Acceleration (m/sec^2)') 0:02:40.000,0:03:00.000 % FIND INSTANT TIME OF RELEASE FROM % RATE OF CHANGE OF LENGTH BB-THUMB LBTx=BBXCD-THUXCD; LBTy=BBYCD-THUYCD; LBTz=BBZCD-THUZCD; LBWx=BBXCD-WRJCX; LBWy=BBYCD-WRJCY; LBWz=BBZCD-WRJCZ; h1=0.00104167; [LENGTH1 LENGTH2]=size(LBTx); for n1=(2:2:LENGTH1-1)'; % OR BBSIDE-1 D1_BT=(LBTx(n1+1,1)-LBTx(n1-1,1))/(2*h1); D1_BT(:,2)=(LBTy(n1+1,1)-LBTy(n1-1,1))/(2*h1); D1_BT(:,3)=(LBTz(n1+1,1)-LBTz(n1-1,1))/(2*h1); D2_BT=(LBTx(n1+1,1)-2*LBTx(n1,1)+LBTx(n1-1,1))/(h1^2); D2_BT(:,2)=(LBTy(n1+1,1)-2*LBTy(n1,1)+LBTy(n1-1,1))/(h1^2); D2_BT(:,3)=(LBTz(n1+1,1)-2*LBTz(n1,1)+LBTz(n1-1,1))/(h1^2); end 0:03:00.000,0:03:10.000 [length1 length2]=size(tqs1); tqs4=tqs1(2:2:length1-1,1); tqs3=tqs4-tqs4(1,1); subplot(1,1,1); plot(tqs1 ,-LBTy,'Linewidth',2); xlabel('Time (sec)') ylabel('Displacement (m)') figure
subplot(2,1,1); plot(tqs3,-D1_BT(:,1),'Linewidth',2); xlabel('Time (sec)') ylabel('1st derivative (m)') subplot(2,1,2); plot(tqs3,-D1_BT(:,2),'Linewidth',2); xlabel('Time (sec)') ylabel('1st derivative (m)') figure
CHECK=[Field Time]; CHECK1=[tqs1 tqs2]; CHECK2=[tqs3 tqs4]; BBX(T1(1,1)); % 866 BBX(T1(1,2)); % 1341 BBZ(T1(1,1)); % 866 BBZ(T1(1,2)); % 1341 BBY(T1(1,1)); % 866 BBY(T1(1,2)); % 1341
% INSTANT TIME OF RELEASE+EQUATIONS OF MOTION UNDER CONST ACCELERATION g=9.81; mass=0.0165; FINIT=T1(1,2); tst=t(INIT,1); % STO 866 tfin=t(T1(1,2)); % STO 1341 tp(1)=tst; VXY=sqrt(VXINIT^2+VYINIT^2); VYZ=sqrt(VYINIT^2+VZINIT^2); VXYZ=sqrt(VXINIT^2+VYINIT^2+VZINIT^2); Vp(1,:)=[VXINIT,-VYINIT,VZINIT]; Vp1(1,:)=Vp(1,:); Xp(1,:)=[BBX(INIT,1)/1000,BBY(INIT,1)/1000,(BBZ(INIT,1)+cor)/1000]; Xp1(1,:)=Xp(1,:); nsteps=T1(1,2)-INIT; % NO DRAG F(1:nsteps,1)=0; F(:,2)=0; F(:,3)=mass*(-g); dt=h1;
% INSTANT TIME OF RELEASE+EQUATIONS OF MOTION DRAG+LIFT+MAGNUS r0=1.2; %kgr/m^3 AIR DENSITY(AREA) (1000mbar-15oC) %? = 105 / (287 × 313) = 1.113 kg/m3 %? = 1.8 10-5 kg/ms (Ns/m2) (15o C) (P) mi=1.80*10^-5; % kgr/m*sec= gr/cm.sec=1 poise=1 Pa.sec1 Pa s = 1 N s/m2 = 1 kg/m s % AIR VISCOCITY diam=0.016; % DIAMETER IN m A=3.14/4*(diam)^2; % CROSS SECTION FOR A SPHERE(A CIRCLE) 0:04:10.000,0:04:20.000 %REYNOLDS NUMBER Re=diam*VXYZ*r0/mi;%diam*VXYZ*r0/mi; SOS=340; %m/sec Mach=VXYZ/SOS; AERODYN=[Re Mach]; Cd=0.78; % estimated drag coeff. (3d) Cl=0.01; % estimated lift(negative) coeff. (3d) 0:04:20.000,0:04:30.000 % Lift F1(1:nsteps,1)=-Cd*0.5*r0*A*VXINIT^2; F1(:,2)= Cd*0.5*r0*A*VYINIT^2; F1(:,3)= -mass*g-Cd*0.5*r0*A*(VZINIT)^2-Cl*0.5*r0*A*(VYINIT)^2; % IF ROTATION CLOCKWISE> Magnus UP - LIFT FORCE UPWARDS % IF ROTATION ANTI-CLOCKWISE> Magnus DOWN - LIFT FORCE DOWNWARDS for i=1:nsteps; t(i+1,:)=t(i)+dt; Vp(i+1,:)=Vp(i,:)+F(i,:)./mass*dt; % using equations of motion under const accel. Xp(i+1,:)=Xp(i,:)+Vp(i,:)*dt; % using equations of motion under const accel. Vp1(i+1,:)=Vp1(i,:)+F1(i,:)./mass*dt; % using equations of motion +drag lift+magnus Xp1(i+1,:)=Xp1(i,:)+Vp1(i,:)*dt; % using equations of motion +drag lift +magnus end
hold all plot3(BBX(INIT:T1(1,2),1)/1000,BBY(INIT:T1(1,2),1)/1000,(BBZ(INIT:T1(1,2),1)+cor)/1000,'o') % real trajectory plot3(Xp(:,1),Xp(:,2),Xp(:,3),'b','Linewidth',2) % predicted trajectory CONST ACCEL plot3(Xp1(:,1),Xp1(:,2),Xp1(:,3),'r','Linewidth',2) % predicted trajectory +DRAG+ LIFT hold off legend('VICON','Constant acceleration','Drag-Lift',3,'Location','Best'); xlabel('Displacement-x (m)') ylabel('Displacement-y (m)') zlabel('Displacement-z (m)') figure
% TARGET HIT subplot(1,1,1); dbrd=plot(DB1Xa(1,1)/1000,(cor+DB1Z(1,1))/1000,'+b',DB2Xa(1,1)/1000,(cor+DB2Za(1,1))/1000,'+b',DB4X(1,1)/1000,(cor+DB4Z(1,1))/1000,'+b',DB3X(1,1)/1000,(cor+DB3Za(1,1))/1000,'+b',C(1,1)/1000,(C(1,3)+cor)/1000,'+b',(BBXCD(BBSIZE,1)),(BBZCD(BBSIZE,1)+cor/1000),'or',Xp(nsteps+1,1),Xp(nsteps+1,3),'bs',Xp1(nsteps+1,1),Xp1(nsteps+1,3),'rs'); legend('Board-Marker 1','Board-Marker 2','Board-Marker 3','Board-Marker 4','Center of Board','Camera','Constant acceleration','Drag & Lift',8,'Location','Best') xlabel('x (m)') ylabel('z (m)') rotate(dbrd,[0 1 0],180); figure 0:04:50.000,0:05:10.000 VQS=[spVX.coefs',-spVY.coefs',spVZ.coefs']; [N11 N22]=size(VQS(:,1)); [X11 X22]=size(Xp(:,1)); for djj=(2:2:N11-4)'; for dj=(1:N1)'; EQSVx1=max(abs(V_BB1(dj,1)-VQS(djj,1))); EQSVx2a=(sum(abs(V_BB1(dj,1)-VQS(djj,1))));EQSVx2=EQSVx2a/N11; EQSVx3a=sqrt(sum(abs(V_BB1(dj,1)-VQS(djj,1)).^2));EQSVx3=EQSVx3a/N11; EQSVy1=max(abs(V_BB1(dj,2)-VQS(djj,2))); EQSVy2a=(sum(abs(V_BB1(dj,2)-VQS(djj,2))));EQSVy2=EQSVy2a/N11; EQSVy3a=sqrt(sum(abs(V_BB1(dj,2)-VQS(djj,2)).^2));EQSVy3=EQSVy3a/N11; EQSVz1=max(abs(V_BB1(dj,3)-VQS(djj,3))); EQSVz2a=(sum(abs(V_BB1(dj,3)-VQS(djj,3))));EQSVz2=EQSVz2a/N11; EQSVz3a=sqrt(sum(abs(V_BB1(dj,3)-VQS(djj,3)).^2));EQSVz3=EQSVz3a/N11; end end
MATRIX3=[T1(1,1), INIT, T1(1,2);fcx, fcy, fcz;httx,htty, httz;px, py, pz]; MATRIX4=[BBX(T1(1,2),1)/1000 BBY(T1(1,2),1)/1000 (BBZ(T1(1,2),1)+cor)/1000;Xp(nsteps+1,1), Xp(nsteps+1,2), Xp(nsteps+1,3);Xp1(nsteps+1,1), Xp1(nsteps+1,2), Xp1(nsteps+1,3);EVx1, EVx2, EVx3;EQSVx1, EQSVx2, EQSVx3;EVy1, EVy2, EVy3;EQSVy1, EQSVy2, EQSVy3;EVz1, EVz2, EVz3;EQSVz1, EQSVz2, EQSVz3];