Appendix B: nemic.c
/* * File: * nemic.c * * Description: * This program mathematically analyzes and simulates * the rattleback by taking various input parameters and * using them in the calculation. The programs uses a * series of equations developed by Professor Mitiguy and * Kane. The series of equations can be split into two * parts: the kinematic & kinetic analysis. At the end * of the program, Mathmatica is called for the final * calculations. * * Compilation: * [me@localhost](~)> gcc -o nemic nemic.c -lm * * Syntax: * [me@localhost](~)> ./nemic * * $ SCC: nemic.c, v 2.3.3 03/12/2000 01:14:35 tobkin Exp $ * * */ #includeAppendix C: nemicsolver.m#include #include float inputparams(char *descr) { float var; printf("Input the %s: ", descr); scanf("%f", &var); return var; } int main(void) { int i; float tmax, a, b, c, h; float A, B, C, D, M; float g, sigma, alpha; float beta, gamma, omega1; float omega2, omega3, u1; float u2, u3, epsilon, epsilont; float x1, x2, x3, u1t, u2t, u3t; float x1t, x2t, x3t, v1, v2; float v3, foo1, foo2, foo3; float F1, F2, F3, R1, R2, R3; float I11, I22, I33, I12, I23; float I13, I31, I32, I21, S1; float S2, S3, Q1, Q2, Q3; float g1, g2, g3, g0; float pi=3.141592; char words[80]; for (i=0;i<30;i++) { printf("\n"); } printf(" _ \n"); printf(" (_) \n"); printf(" _ __ ___ _ __ ___ _ ___ \n"); printf(" | '_ \\ / _ \\ '_ ` _ \\| |/ __| \n"); printf(" | | | | __/ | | | | | | (__ \n"); printf(" |_| |_|\\___|_| |_| |_|_|\\___| \n\n"); printf(" Version 2.3.3 \n\n"); printf("Welcome to the rattleback analyzer program. \n\n"); printf("Instructions: \n"); printf("1. Gather information about the rattleback to being \n"); printf(" calculated and simulated. \n"); printf("2. Input the various values when prompted for in the \n"); printf(" program. \n\n"); strcpy(words, "time the simulation should run (20)"); tmax=inputparams(words); strcpy(words, "ellipsoid shape (x-axis, .2)"); a=inputparams(words); strcpy(words, "ellipsoid shape (y-axis, .03)"); b=inputparams(words); strcpy(words, "ellipsoid shape (z-axis, .02)"); c=inputparams(words); strcpy(words, "height of coordinate system above the center of gravity (.01)"); h=inputparams(words); strcpy(words, "mass moment of inertia A (.0002)"); A=inputparams(words); strcpy(words, "mass moment of inertia B (.0016)"); B=inputparams(words); strcpy(words, "mass moment of inertia C (.0017)"); C=inputparams(words); strcpy(words, "mass moment of inertia D (-.000020)"); D=inputparams(words); strcpy(words, "mass of the ellipsoid (1)"); M=inputparams(words); strcpy(words, "gravitational pull (9.81)"); g=inputparams(words); strcpy(words, "air resistance coefficient (0)"); sigma=inputparams(words); alpha = ((.5 * pi) / 180); beta = ((.5 * pi) / 180); gamma = 0; omega1 = 0; omega2 = 0; omega3 = 0; u1 = (-cos(alpha) * sin(beta)); u2 = (sin(alpha)); u3 = (cos(alpha) * cos(beta)); epsilon = (sqrt(pow((a * u1),2) + pow((b * u2),2) + pow((c * u3),2))); x1 = ((pow(a,2) * u1) / epsilon); x2 = ((pow(b,2) * u2) / epsilon); x3 = ((pow(c,2) * u3) / epsilon); u1t = ((omega2 * u2) - (omega2 * u3)); u2t = ((omega3 * u3) - (omega3 * u1)); u3t = ((omega1 * u1) - (omega1 * u2)); epsilont = ((pow(a,2) * u1 * u1t) + (pow(b,2) * u2 * u2t) + (pow(c,2) * u3 * u3t)); x1t = ((pow(a,2) * ((epsilon * u1t) - (epsilont * u1))) / pow(epsilon,2)); x2t = ((pow(b,2) * ((epsilon * u2t) - (epsilont * u2))) / pow(epsilon,2)); x3t = ((pow(c,2) * ((epsilon * u3t) - (epsilont * u3))) / pow(epsilon,2)); v1 = ((omega2 * (h - x3)) + (omega3 * x2)); v2 = ((-omega3 * x1) - (omega1 * (h - x3))); v3 = ((-omega1 * x2) + (omega2 * x1)); foo1 = ((omega2 * (v3 - x3t)) - (omega3 * (v2 - x2t))); foo2 = ((omega3 * (v1 - x1t)) - (omega1 * (v3 - x3t))); foo3 = ((omega1 * (v2 - x2t)) - (omega2 * (v1 - x1t))); F1 = ((-sigma * omega1) + ((M * g) * (((x3 - h) * u2) - (x2 * u3)))); F2 = ((-sigma * omega2) + ((M * g) * (((h - x3) * u1) - (x1 * u3)))); F3 = ((-sigma * omega3) + ((M * g) * ((x2 * u1) - (x1 * u2)))); R1 = (((D * omega1) + ((B - C) * omega2)) * omega3); R2 = ((((C - A) * omega1) - (D * omega2)) * omega3); R3 = (D * (pow(omega2,2) - pow(omega1,2)) + ((A - B) * omega1 * omega2)); I11 = (A + (M * (pow(x2,2) + pow((h - x3),2)))); I22 = (B + (M * (pow(x1,2) + pow((h - x3),2)))); I33 = (C + (M * (pow(x1,2) + pow(x2,2)))); I12 = (D - (M * x1 * x2)); I23 = (M * (h - x3) * x2); I13 = (M * (h - x3) * x1); I31 = I13; I32 = I23; I21 = I12; S1 = (M * (((h - x3) * foo2)) + (x2 * foo3)); S2 = (M * (((x3 - h) * foo1)) - (x1 * foo3)); S3 = (M * ((x1 * foo2) - (x2 * foo1))); Q1 = F1 + R1 + S1; Q2 = F2 + R2 + S2; Q3 = F3 + R3 + S3; g1 = (((Q1 * I22 * I33) + (I12 * I23 * Q3) + (I13 * Q2 * I32)) - ((Q3 * I22 * I13) + (I32 * I23 * Q1) * (I33 * Q2 * I12))); g2 = (((I11 * Q2 * I33) + (Q1 * I23 * I31) + (I13 * I21 * Q3)) - ((I31 * Q2 * I13) + (Q3 * I23 * I11) * (I33 * I21 * Q1))); g3 = (((I11 * I22 * Q3) + (I12 * Q2 * I31) + (Q1 * I21 * I32)) - ((I31 * I22 * Q1) + (I32 * Q2 * I11) * (Q3 * I21 * I12))); g0 = (((I11 * I22 * I33) + (I12 * I23 * I31) + (I13 * I21 * I32)) - ((I31 * I22 * I13) + (I32 * I23 * I11) + (I33 * I21 * I12))); printf("\n\nQ1 = %f\n", Q1); printf("Q2 = %f\n", Q2); printf("Q3 = %f\n", Q3); printf("g1 = %f\n", g1); printf("g2 = %f\n", g2); printf("g3 = %f\n", g3); printf("g0 = %f\n", g0); return 0; } /* EOF */
tspan = [0,20]
y10 = .5*pi/180
y20 = .5*pi/180
y30 = 0
y40 = 0
y50 = 0
y60 = -5
y0 = [y10;y20;y30;y40;y50;y60]
[t,y] = ode45('nemic',tspan,y0);
figure
plot(t,y(:,1))
grid on
xlabel('t (sec)')
ylabel('alpha (rad)')
title('Roll angle ')
figure
plot(t,y(:,2))
grid on
xlabel('t (sec)')
ylabel('beta (rad)')
title('Pitch angle')
figure
plot(t,y(:,3))
grid on
xlabel('t (sec)')
ylabel('gamma (rad)')
title('Yaw angle')
figure
plot(t,y(:,4))
grid on
xlabel('t (sec)')
ylabel('Omega1 (rad/sec)')
title('Angular velocity of roll')
figure
plot(t,y(:,5))
grid on
xlabel('t (sec)')
ylabel('Omega2 (rad/sec)')
title('Angular velocity of pitch')
figure
plot(t,y(:,6))
grid on
xlabel('t (sec)')
ylabel('Omega3 (rad/sec)')
title('Angular velocity of yaw')
delta = acos(cos(y(:,1)).*cos(y(:,2)));
figure
plot(t,delta)
grid on
xlabel('t (sec)')
ylabel('Delta')
title('Angle between vertical axes of ellipsoid and surface')
Appendix D: nemic.m
nemic.m
function ydot = nemic(t,y)
a = .2;
b = .03;
c = .02;
h = .01;
A = 2.0*(10^-4);
B = 1.6*(10^-3);
C = 1.7*(10^-3);
D = -2.0*(10^-5);
M = 1;
g = 9.81;
sigma = 0;
u1 = -cos(y(1))*sin(y(2));
u2 = sin(y(1));
u3 = cos(y(1))*cos(y(2));
eps = sqrt((a*u1)^2 + (b*u2)^2 + (c*u3)^2);
x1 = (a^2)*u1/eps;
x2 = (b^2)*u2/eps;
x3 = (c^2)*u3/eps;
u1t = y(6)*u2 - y(5)*u3;
u2t = y(4)*u3 - y(6)*u1;
u3t = y(5)*u1 - y(4)*u2;
epst = ((a^2)*u1*u1t + (b^2)*u2*u2t + (c^2)*u3*u3t)/eps;
x1t = ((a^2)*(eps*u1t -epst*u1))/eps^2;
x2t = ((b^2)*(eps*u2t -epst*u2))/eps^2;
x3t = ((c^2)*(eps*u3t -epst*u3))/eps^2;
v1 = y(5)*(h-x3)+ y(6)*x2;
v2 = -y(6)*x1 - y(4)*(h-x3);
v3 = -y(4)*x2 + y(5)*x1;
zeta1 = y(5)*(v3-x3t) - y(6)*(v2-x2t);
zeta2 = y(6)*(v1-x1t) - y(4)*(v3-x3t);
zeta3 = y(4)*(v2-x2t) - y(5)*(v1-x1t);
F1 = -sigma*y(4) + M*g*((x3-h)*u2 - x2*u3);
F2 = -sigma * y(5) + M * g * ((h - x3) * u1 + x1*u3);
F3 = -sigma * y(6) + M * g * (x2 * u1 - x1 * u2);
R1 = (D * y(4) + (B - C) * y(5)) * y(6);
R2 = ((C - A) * y(4) - D * y(5)) * y(6);
R3 = D * (y(5)^2 - y(4)^2) + (A - B) * y(4) * y(5);
I11 = A + M * (x2^2 + (h - x3)^2);
I22 = B + M * (x1^2 + (h - x3)^2);
I33 = C + M * (x1^2 + x2^2);
I12 = D - M * x1 * x2;
I23 = M * (h - x3) * x2;
I13 = M * (h - x3) * x1;
I31 = I13;
I32 = I23;
I21 = I12;
S1 = M * ((h - x3) * zeta2 + x2 * zeta3);
S2 = M * ((x3 - h) * zeta1 - x1 * zeta3);
S3 = M * (x1 * zeta2 - x2 * zeta1);
Q1 = F1 + R1 + S1;
Q2 = F2 + R2 + S2;
Q3 = F3 + R3 + S3;
g1 = [Q1 I12 I13; Q2 I22 I23; Q3 I32 I33];
g2 = [I11 Q1 I13; I21 Q2 I23; I31 Q3 I33];
g3 = [I11 I12 Q1; I21 I22 Q2; I31 I32 Q3];
g0 = [I11 I12 I13; I21 I22 I23; I31 I32 I33];
detg0 = det(g0);
y4 = det(g1)/detg0;
y5 = det(g2)/detg0;
y6 = det(g3)/detg0;
y1 = y(6) * sin(y(2)) + y(4) * cos(y(2));
y2 = (-y(6)* cos(y(2)) + y(4) * sin(y(2))) * tan(y(1)) + y(5);
y3 = (y(6) * cos(y(2)) - y(4) * sin(y(2))) * sec(y(1));
ydot = [y1;y2;y3;y4;y5;y6];