Guest User

Untitled

a guest
Apr 24th, 2018
70
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.57 KB | None | 0 0
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. const float G=-6.673E-11;//stores the gravitational constant
  6.  
  7. double funcx(double M, double x, double y);//the functions for the velocities
  8. double funcy(double M, double x, double y);
  9.  
  10. int main()
  11. {
  12. /*defining variables*/
  13. FILE *path;
  14. double x,y,t,T,dt,vx,vy,M,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y,k1vx,k2vx,k3vx,k4vx,k1vy,k2vy,k3vy,k4vy,KE;
  15.  
  16. printf("Calculating a basic orbit of a moving body\n");
  17.  
  18. /*setting values for independant variables*/
  19. printf("\nEnter the mass of the large object: ");
  20. scanf("%lf",&M);
  21. printf("\nEnter the initial coordinates of the moving object\n");
  22. printf("x = ");
  23. scanf("%lf",&x);
  24. printf("y = ");
  25. scanf("%lf",&y);
  26. printf("\nEnter the initial velocities of the moving object\n");
  27. printf("x velocity = ");
  28. scanf("%lf",&vx);
  29. printf("y velocity = ");
  30. scanf("%lf",&vy);
  31. printf("\nEnter the time over which to map the orbit: ");
  32. scanf("%lf",&T);
  33. printf("Enter the size of the time interval to calculate over: ");
  34. scanf("%lf",&dt);
  35.  
  36. path=fopen("path.txt","w");
  37.  
  38. for(t=0;t<T;t+=dt)//main loop for finding values
  39. {
  40. /*defining runge-kutta k values for x direction*/
  41. k1vx=funcx(M,x,y);//vx and x are interdependent so have to be calculated together
  42. k1x=vx;
  43. k2vx=funcx(M,(x+(k1x*dt/2)),y);
  44. k2x=vx*k1vx*(dt/2);
  45. k3vx=funcx(M,(x+(k2x*dt/2)),y);
  46. k3x=vx*k2vx*(dt/2);
  47. k4vx=funcx(M,(x+(k3x*dt)),y);
  48. k4x=vx*k3vx*dt;
  49.  
  50. /*defining runge-kutta k values for x direction*/
  51. k1vy=funcy(M,x,y);//vy and y are interdependent so have to be calculated together
  52. k1y=vy;
  53. k2vy=funcy(M,x,(y+(k1y*dt/2)));
  54. k2y=vy*k1vy*(dt/2);
  55. k3vy=funcy(M,x,(y+(k2y*dt/2)));
  56. k3y=vy*k2vy*(dt/2);
  57. k4vy=funcy(M,x,(y+(k3y*dt)));
  58. k4y=vy*k3vy*dt;
  59.  
  60. /*calculating from runge-kutta for v*/
  61. vx+=(dt/6)*(k1vx+(2*k2vx)+(2*k3vx)+k4vx);
  62. vy+=(dt/6)*(k1vy+(2*k2vy)+(2*k3vy)+k4vy);
  63.  
  64. /*distance dependence was already incorporated in finding the v values*/
  65. x+=dt*vx;
  66. y+=dt*vy;
  67.  
  68. KE=0.5*((vx*vx)+(vy*vy));//not true KE, indicates conservation over time
  69. fprintf(path,"%f%c%f%c%f%c%f%c%f%c%f\n",t,9,x,9,y,9,vx,9,vy,9,KE);
  70. }
  71.  
  72. return 0;
  73. }
  74. double funcx(double M, double x, double y)
  75. {
  76. return (M*G*x)/pow((x*x)+(y*y),3/(double)2);
  77. }
  78. double funcy(double M, double x, double y)
  79. {
  80. return (M*G*y)/pow((x*x)+(y*y),3/(double)2);
  81. }
Add Comment
Please, Sign In to add comment