Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- const float G=-6.673E-11;//stores the gravitational constant
- double funcx(double M, double x, double y);//the functions for the velocities
- double funcy(double M, double x, double y);
- int main()
- {
- /*defining variables*/
- FILE *path;
- 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;
- printf("Calculating a basic orbit of a moving body\n");
- /*setting values for independant variables*/
- printf("\nEnter the mass of the large object: ");
- scanf("%lf",&M);
- printf("\nEnter the initial coordinates of the moving object\n");
- printf("x = ");
- scanf("%lf",&x);
- printf("y = ");
- scanf("%lf",&y);
- printf("\nEnter the initial velocities of the moving object\n");
- printf("x velocity = ");
- scanf("%lf",&vx);
- printf("y velocity = ");
- scanf("%lf",&vy);
- printf("\nEnter the time over which to map the orbit: ");
- scanf("%lf",&T);
- printf("Enter the size of the time interval to calculate over: ");
- scanf("%lf",&dt);
- path=fopen("path.txt","w");
- for(t=0;t<T;t+=dt)//main loop for finding values
- {
- /*defining runge-kutta k values for x direction*/
- k1vx=funcx(M,x,y);//vx and x are interdependent so have to be calculated together
- k1x=vx;
- k2vx=funcx(M,(x+(k1x*dt/2)),y);
- k2x=vx*k1vx*(dt/2);
- k3vx=funcx(M,(x+(k2x*dt/2)),y);
- k3x=vx*k2vx*(dt/2);
- k4vx=funcx(M,(x+(k3x*dt)),y);
- k4x=vx*k3vx*dt;
- /*defining runge-kutta k values for x direction*/
- k1vy=funcy(M,x,y);//vy and y are interdependent so have to be calculated together
- k1y=vy;
- k2vy=funcy(M,x,(y+(k1y*dt/2)));
- k2y=vy*k1vy*(dt/2);
- k3vy=funcy(M,x,(y+(k2y*dt/2)));
- k3y=vy*k2vy*(dt/2);
- k4vy=funcy(M,x,(y+(k3y*dt)));
- k4y=vy*k3vy*dt;
- /*calculating from runge-kutta for v*/
- vx+=(dt/6)*(k1vx+(2*k2vx)+(2*k3vx)+k4vx);
- vy+=(dt/6)*(k1vy+(2*k2vy)+(2*k3vy)+k4vy);
- /*distance dependence was already incorporated in finding the v values*/
- x+=dt*vx;
- y+=dt*vy;
- KE=0.5*((vx*vx)+(vy*vy));//not true KE, indicates conservation over time
- 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);
- }
- return 0;
- }
- double funcx(double M, double x, double y)
- {
- return (M*G*x)/pow((x*x)+(y*y),3/(double)2);
- }
- double funcy(double M, double x, double y)
- {
- return (M*G*y)/pow((x*x)+(y*y),3/(double)2);
- }
Add Comment
Please, Sign In to add comment