Guest User

Untitled

a guest
May 23rd, 2018
82
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.96 KB | None | 0 0
  1. import std.stdio; [59/187]
  2. // import std.math;
  3. import core.stdc.stdio;
  4. import core.stdc.math;
  5.  
  6. enum pi = 3.141592653589793;
  7. enum solar_mass = (4 * pi * pi);
  8. enum days_per_year = 365.24;
  9.  
  10. struct Planet {
  11. double x, y, z;
  12. double vx, vy, vz;
  13. double mass;
  14. }
  15.  
  16. void advance(Planet[] bodies, double dt)
  17. {
  18. foreach (i; 0 .. bodies.length) {
  19. auto b = &(bodies[i]);
  20. foreach (j; i+1 .. bodies.length) {
  21. auto b2 = &(bodies[j]);
  22. auto dx = b.x - b2.x;
  23. auto dy = b.y - b2.y;
  24. auto dz = b.z - b2.z;
  25. auto distance = sqrt(dx * dx + dy * dy + dz * dz);
  26. auto mag = dt / (distance * distance * distance);
  27. b.vx -= dx * b2.mass * mag;
  28. b.vy -= dy * b2.mass * mag;
  29. b.vz -= dz * b2.mass * mag;
  30. b2.vx += dx * b.mass * mag;
  31. b2.vy += dy * b.mass * mag;
  32. b2.vz += dz * b.mass * mag;
  33. }
  34. }
  35. foreach (i; 0 .. bodies.length) {
  36. auto b = &(bodies[i]);
  37. b.x += dt * b.vx;
  38. b.y += dt * b.vy;
  39. b.z += dt * b.vz;
  40. }
  41. }
  42.  
  43. double energy(Planet[] bodies)
  44. {
  45. double e = 0.0;
  46. foreach (i; 0 .. bodies.length) {
  47. auto b = &(bodies[i]);
  48. e += 0.5 * b.mass * (b.vx * b.vx + b.vy * b.vy + b.vz * b.vz);
  49. foreach (j; i+1 .. bodies.length) {
  50. auto b2 = &(bodies[j]);
  51. auto dx = b.x - b2.x;
  52. auto dy = b.y - b2.y;
  53. auto dz = b.z - b2.z;
  54. auto distance = sqrt(dx * dx + dy * dy + dz * dz);
  55. e -= (b.mass * b2.mass) / distance;
  56. }
  57. }
  58. return e;
  59. }
  60.  
  61. void offset_momentum(Planet[] bodies)
  62. {
  63. double px = 0.0, py = 0.0, pz = 0.0;
  64. foreach (i; 0 .. bodies.length) {
  65. px += bodies[i].vx * bodies[i].mass;
  66. py += bodies[i].vy * bodies[i].mass;
  67. pz += bodies[i].vz * bodies[i].mass;
  68. }
  69. bodies[0].vx = - px / solar_mass;
  70. bodies[0].vy = - py / solar_mass;
  71. bodies[0].vz = - pz / solar_mass;
  72. }
  73.  
  74. enum NBODIES = 5;
  75.  
  76. Planet[NBODIES] bodies = [
  77. { /* sun */
  78. x: 0.0, y: 0.0, z: 0.0, vx: 0.0, vy: 0.0, vz: 0.0, mass:solar_mass
  79. },
  80. { /* jupiter */
  81. 4.84143144246472090e+00,
  82. -1.16032004402742839e+00,
  83. -1.03622044471123109e-01,
  84. 1.66007664274403694e-03 * days_per_year,
  85. 7.69901118419740425e-03 * days_per_year,
  86. -6.90460016972063023e-05 * days_per_year,
  87. 9.54791938424326609e-04 * solar_mass
  88. },
  89. { /* saturn */
  90. 8.34336671824457987e+00,
  91. 4.12479856412430479e+00,
  92. -4.03523417114321381e-01,
  93. -2.76742510726862411e-03 * days_per_year,
  94. 4.99852801234917238e-03 * days_per_year,
  95. 2.30417297573763929e-05 * days_per_year,
  96. 2.85885980666130812e-04 * solar_mass
  97. },
  98. { /* uranus */
  99. 1.28943695621391310e+01,
  100. -1.51111514016986312e+01,
  101. -2.23307578892655734e-01,
  102. 2.96460137564761618e-03 * days_per_year,
  103. 2.37847173959480950e-03 * days_per_year,
  104. -2.96589568540237556e-05 * days_per_year,
  105. 4.36624404335156298e-05 * solar_mass
  106. },
  107. { /* neptune */
  108. 1.53796971148509165e+01,
  109. -2.59193146099879641e+01,
  110. 1.79258772950371181e-01,
  111. 2.68067772490389322e-03 * days_per_year,
  112. 1.62824170038242295e-03 * days_per_year,
  113. -9.51592254519715870e-05 * days_per_year,
  114. 5.15138902046611451e-05 * solar_mass
  115. }
  116. ];
  117.  
  118. void main(string[] args)
  119. {
  120. import std.conv;
  121. int n = to!(int)(args[1]);
  122.  
  123. offset_momentum(bodies);
  124. printf("%.9f\n", energy(bodies));
  125. foreach (i; 0 .. n)
  126. advance(bodies, 0.01);
  127. printf("%.9f\n", energy(bodies));
  128. }
Add Comment
Please, Sign In to add comment