API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 18th, 2019
101
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.96 KB | None | 0 0
raw download clone embed print report
  1. mags = {{1.2, 0, 0, 1, 0, 0, -1}, {1.2, Pi/6, 0, 1, -Sin[Pi/6],
  2. 0, -Cos[Pi/6]}, {1.2, Pi/6, 2*Pi/3,
  3. 1, -Sin[Pi/6]*Cos[2*Pi/3], -Sin[Pi/6]*
  4. Sin[2*Pi/3], -Cos[Pi/6]}, {1.2, Pi/6, 4*Pi/3,
  5. 1, -Sin[Pi/6]*Cos[4*Pi/3], -Sin[Pi/6]*Sin[4*Pi/3], -Cos[Pi/6]}};
  6.  
  7. nf = Nearest[{{-0.082788, 0}, {0.041394,
  8. 0.071696}, {0.041394, -0.071696}} -> Automatic] /* First;
  9. eq = {tht''[
  10. t] - (Sin[tht[t]]*Cos[tht[t]]*
  11. phi'[t]^2 - (Sum[
  12. c*mags[[i,
  13. 4]]*((mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
  14. mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
  15. mags[[i, 7]]*
  16. Sin[tht[t]]) + (3*
  17. mags[[i,
  18. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  19. Cos[tht[t]]*Cos[phi[t]] +
  20. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
  21. Sin[phi[t]] -
  22. Cos[mags[[i, 2]]]*
  23. Sin[tht[
  24. t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
  25. mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
  26. mags[[i, 7]]*Cos[tht[t]]) -
  27. mags[[i,
  28. 1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
  29. Cos[mags[[i, 3]]] +
  30. mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
  31. mags[[i, 7]]*Cos[mags[[i, 2]]])) -
  32. 3*(mags[[i, 5]]*Cos[tht[t]]*Cos[phi[t]] +
  33. mags[[i, 6]]*Cos[tht[t]]*Sin[phi[t]] -
  34. mags[[i, 7]]*Sin[tht[t]])*(1 -
  35. mags[[i,
  36. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  37. Sin[tht[t]]*Cos[phi[t]] +
  38. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  39. Sin[phi[t]] +
  40. Cos[mags[[i, 2]]]*Cos[tht[t]])))/(((1 +
  41. mags[[i, 1]]^2 -
  42. 2*mags[[i,
  43. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  44. Sin[tht[t]]*Cos[phi[t]] +
  45. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  46. Sin[phi[t]] +
  47. Cos[mags[[i, 2]]]*Cos[tht[t]])))^1) -
  48. 15*mags[[i,
  49. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Cos[tht[t]]*
  50. Cos[phi[t]] +
  51. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Cos[tht[t]]*
  52. Sin[phi[t]] -
  53. Cos[mags[[i, 2]]]*
  54. Sin[tht[t]])*((mags[[i, 5]]*Sin[tht[t]]*
  55. Cos[phi[t]] +
  56. mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
  57. mags[[i, 7]]*Cos[tht[t]]) -
  58. mags[[i,
  59. 1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
  60. Cos[mags[[i, 3]]] +
  61. mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
  62. mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
  63. mags[[i,
  64. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  65. Sin[tht[t]]*Cos[phi[t]] +
  66. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  67. Sin[phi[t]] +
  68. Cos[mags[[i, 2]]]*Cos[tht[t]]))/((1 +
  69. mags[[i, 1]]^2 -
  70. 2*mags[[i,
  71. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  72. Sin[tht[t]]*Cos[phi[t]] +
  73. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  74. Sin[phi[t]] +
  75. Cos[mags[[i, 2]]]*Cos[tht[t]])))^2)/((1 +
  76. mags[[i, 1]]^2 -
  77. 2*mags[[i,
  78. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  79. Sin[tht[t]]*Cos[phi[t]] +
  80. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  81. Sin[phi[t]] +
  82. Cos[mags[[i, 2]]]*Cos[tht[t]])))^1.5, {i,
  83. 4}] + (m + M/2)*g*L*Sin[tht[t]] +
  84. k*(L^2)*tht'[t]/3)/((M/3 + m)*L^2)) == 0,
  85. phi''[t] - (-Sum[
  86. c*mags[[i,
  87. 4]]*((-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
  88. mags[[i, 6]]*Sin[tht[t]]*
  89. Cos[phi[t]]) + (3*
  90. mags[[i,
  91. 1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  92. Sin[tht[t]]*Sin[phi[t]] +
  93. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  94.  
  95. Cos[phi[
  96. t]])*(2*(mags[[i, 5]]*Sin[tht[t]]*Cos[phi[t]] +
  97. mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
  98. mags[[i, 7]]*Cos[tht[t]]) -
  99. mags[[i,
  100. 1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
  101. Cos[mags[[i, 3]]] +
  102. mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
  103. mags[[i, 7]]*Cos[mags[[i, 2]]])) -
  104. 3*(-mags[[i, 5]]*Sin[tht[t]]*Sin[phi[t]] +
  105. mags[[i, 6]]*Sin[tht[t]]*Cos[phi[t]])*(1 -
  106. mags[[i,
  107. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  108. Sin[tht[t]]*Cos[phi[t]] +
  109. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  110. Sin[phi[t]] +
  111. Cos[mags[[i, 2]]]*Cos[tht[t]])))/((1 +
  112. mags[[i, 1]]^2 -
  113. 2*mags[[i,
  114. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  115. Sin[tht[t]]*Cos[phi[t]] +
  116. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  117. Sin[phi[t]] +
  118. Cos[mags[[i, 2]]]*Cos[tht[t]]))^1) -
  119. 15*mags[[i,
  120. 1]]*(-Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
  121. Sin[phi[t]] +
  122. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  123. Cos[phi[t]])*((mags[[i, 5]]*Sin[tht[t]]*
  124. Cos[phi[t]] +
  125. mags[[i, 6]]*Sin[tht[t]]*Sin[phi[t]] +
  126. mags[[i, 7]]*Cos[tht[t]]) -
  127. mags[[i,
  128. 1]]*(mags[[i, 5]]*Sin[mags[[i, 2]]]*
  129. Cos[mags[[i, 3]]] +
  130. mags[[i, 6]]*Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]] +
  131. mags[[i, 7]]*Cos[mags[[i, 2]]]))*(1 -
  132. mags[[i,
  133. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  134. Sin[tht[t]]*Cos[phi[t]] +
  135. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  136. Sin[phi[t]] + Cos[mags[[i, 2]]]*Cos[tht[t]]))/(1 +
  137. mags[[i, 1]]^2 -
  138. 2*mags[[i,
  139. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*
  140. Sin[tht[t]]*Cos[phi[t]] +
  141. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  142. Sin[phi[t]] +
  143. Cos[mags[[i, 2]]]*Cos[tht[t]]))^2)/(1 +
  144. mags[[i, 1]]^2 -
  145. 2*mags[[i,
  146. 1]]*(Sin[mags[[i, 2]]]*Cos[mags[[i, 3]]]*Sin[tht[t]]*
  147. Cos[phi[t]] +
  148. Sin[mags[[i, 2]]]*Sin[mags[[i, 3]]]*Sin[tht[t]]*
  149. Sin[phi[t]] +
  150. Cos[mags[[i, 2]]]*Cos[tht[t]]))^1.5, {i,
  151. 4}]/(((M/3 + m)*L^2)*Sin[tht[t]]^2) -
  152. 2*tht'[t]*phi'[t]/Tan[tht[t]] -
  153. k*(L^2)*phi'[t]/(3*((M/3 + m)*L^2))) == 0, tht'[0] == 0,
  154. phi'[0] == 0.5, tht[0] == tht0, phi[0] == phi0};
  155.  
  156. getStateData[op_] :=
  157. First@NDSolve`ProcessEquations[eq /. op, {tht, phi}, t,
  158. Method -> "Adams"];
  159. op = {k -> .1, c -> .0004, M -> .05, m -> .01, L -> .3, g -> 9.8,
  160. tht0 -> 1, phi0 -> 3};
  161. sd = getStateData[op]
  162.  
  163. getBaSin[x_, y_] :=
  164. If[x*x + y*y <= 0.3*0.3,
  165. Module[{state = sd, sol},
  166. state = First[
  167. NDSolve`Reinitialize[
  168. sd, {tht[0] == ArcSin[((x^2 + y^2)^0.5)/0.3],
  169. phi[0] ==
  170. If[x != 0, If[x > 0, ArcTan[y/x], ArcTan[y/x] + Pi],
  171. If[y > 0, Pi/2, -Pi/2]]}]];
  172. NDSolve`Iterate[state, 21];
  173. sol = {tht[20], phi[20]} /. NDSolve`ProcessSolutions[state];
  174. nf[{0.3*Sin[sol[[1]]]*Cos[sol[[2]]],
  175. 0.3*Sin[sol[[1]]]*Sin[sol[[2]]]}]], 0]
  176.  
  177. plot1 = ArrayPlot[
  178. ParallelTable[
  179. getBaSin[x, y], {y, 0.3, -0.3, -0.1}, {x, -0.3, 0.3, 0.1}],
  180. ColorRules -> {0 -> White, 1 -> Red, 2 -> Green, 3 -> Blue}] //
  181. AbsoluteTiming
  182. plot2 = ListPlot[{{-0.082788, 0}, {0.041394,
  183. 0.071696}, {0.041394, -0.071696}}]
  184. Show[{plot1, plot2}]
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!