SHARE
TWEET

Untitled

a guest Jun 30th, 2011 113 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. Off[General::"spell"]
  2. Off[General::"spell1"]
  3.  
  4. MotiPoin[A_, B_, C0_, r0_, theta0_, b_, alpha_] :=
  5. Module[{q0, trif, K2, T, h, eq},
  6. q0 = (C0 r0 Tan[theta0])/B;
  7. trif = (2 \[Pi] B)/(r0 (A Cos[theta0] + C0 Sin[theta0] Tan[theta0]));
  8. K2 = (B q0)^2 + (C0 r0)^2;
  9. T = 1/2 (B q0^2 + C0 r0^2);
  10. h = Sqrt[(2 T)/K2];
  11. eq1 = Derivative[1][p][t] == ((B - C0) q[t] r[t])/A;
  12. eq2 = Derivative[1][q][t] == ((C0 - A) p[t] r[t])/B;
  13. eq3 = Derivative[1][r][t] == ((A - B) p[t] q[t])/C0;
  14. eq4 = Derivative[1][psi][t] == (Cos[phi[t]] q[t] + p[t] Sin[phi[t]])/Sin[theta[t]];
  15. eq5 = Derivative[1][phi][t] == r[t] - Cot[theta[t]] (Cos[phi[t]] q[t] + p[t] Sin[phi[t]]);
  16. eq6 = Derivative[1][theta][t] == p[t] Cos[phi[t]] - q[t] Sin[phi[t]];
  17. w1 = (Cos[phi[t]] Cos[psi[t]] - Sin[phi[t]] Sin[psi[t]] Cos[theta[t]]) p[t] - (Cos[psi[t]] Sin[phi[t]] + Cos[phi[t]] Cos[theta[t]] Sin[psi[t]]) q[t] + r[t] Sin[psi[t]] Sin[theta[t]];
  18. w2 = (Cos[psi[t]] Cos[theta[t]] Sin[phi[t]] + Cos[phi[t]] Sin[psi[t]]) p[t] + (Cos[phi[t]] Cos[psi[t]] Cos[theta[t]] - Sin[phi[t]] Sin[psi[t]]) q[t] - Cos[psi[t]] r[t] Sin[theta[t]];
  19. w3 = Cos[theta[t]] r[t] + Cos[phi[t]] q[t] Sin[theta[t]] + p[t] Sin[phi[t]] Sin[theta[t]];
  20. sol = NDSolve[{eq1, eq2, eq3, eq4, eq5, eq6, p[0] == 0, q[0] == q0, r[0] == r0, psi[0] == 0, phi[0] == 0, theta[0] == theta0}, {p, q, r, psi, phi, theta}, {t, 0, b trif}];
  21. {x, y} = Flatten[{-((w1 h)/w3), -((w2 h)/w3)} /. sol];
  22. z = x^2 + y^2;
  23. If[A < C0 < B || B < C0 < A, Goto[2], Goto[1]];
  24.  
  25. Label[1];
  26. m = FindMinimum[z, {t, 0, 0, b trif}];
  27. M = FindMinimum[-z, {t, 0, 0, b trif}];
  28. ra1 = Sqrt[m[[1]]];
  29. ra2 = Sqrt[-M[[1]]];
  30. Print["L'erpoloide è contenuta in una corona circolare"];
  31. Print["avente raggio interno ra1 e raggio esterno ra2"];
  32. Print["ra1=", ra1]; Print["ra2=", ra2];
  33. c1 = ParametricPlot[{ra1 Sin[u], ra1 Cos[u]}, {u, 0, 2 \[Pi]}, AspectRatio -> 1, DisplayFunction -> Identity, PlotStyle -> RGBColor[0.8669, 0.258, 0.227]];
  34. c2 = ParametricPlot[{ra2 Sin[u], ra2 Cos[u]}, {u, 0, 2 \[Pi]}, AspectRatio -> 1, DisplayFunction -> Identity, PlotStyle -> RGBColor[0.925, 0.140, 0.129]];
  35. Plot[Sqrt[z], {t, 0, b trif}, AxesLabel -> {"t", "ra"}];
  36. erp = ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> 1, PlotRange -> All, DisplayFunction -> Identity];
  37. Show[erp, c1, c2, DisplayFunction -> $DisplayFunction];
  38. Goto[3];
  39.  
  40. Label[2];
  41. Plot[Sqrt[z], {t, 0, b trif}, AxesLabel -> {"t", "ra"}];
  42. erp = ParametricPlot[{x, y}, {t, 0, b trif}, AspectRatio -> 1, PlotRange -> All];
  43.  
  44. Label[3];
  45. xp = p[t]/Sqrt[2 T] /. sol;
  46. yp = q[t]/Sqrt[2 T] /. sol;
  47. zp = r[t]/Sqrt[2 T] /. sol;
  48. X = (Cos[u] Sin[v])/Sqrt[A];
  49. Y = (Sin[u] Sin[v])/Sqrt[B];
  50. Z = Cos[v]/Sqrt[C0];
  51. el = ParametricPlot3D[{X, Y, Z}, {u, 0, 2 \[Pi]}, {v, 0, alpha}, LightSources -> {{{-1, -1, 3}, GrayLevel[0.999]}}, Boxed -> False, DisplayFunction -> Identity];
  52. pol = ParametricPlot3D[Evaluate[Flatten[{xp, yp, zp}] /. sol], {t, 0, b trif}, PlotPoints -> 200, DisplayFunction -> Identity];
  53. Show[el, pol, DisplayFunction -> $DisplayFunction];
  54. ]
  55.  
  56. MotiPoin[1,1.5,0.5,3,Pi/4,1.5,Pi/4]
  57. MotiPoin[1, 1.5, 0.5, 3, 0.01, 1.5, Pi/100]
  58. MotiPoin[0.5, 1.5, 1, -3, 0.01, 3.5, Pi]
  59. MotiPoin[1, 1, 1.5, 3, Pi/4, 2.5, Pi/2]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top