C-Programming: My Code & Problem Sheet

[Hoofbeat]

Thanks to everyone who offered to take a look at my problem sheet and my relevant code. I have scanned in the 3 relevant pages of my problem sheet and they can be seen at:

www.hoofbeat.tv/tempcomp.htm[/URL] (the files are too large to upload onto here)

Below is a copy of my original code (we've tried soooo many versions for the acceleration formula :rolleyes: ). NB. Incase you're not familiar with programming in Windows, the code "system (PAUSE)" is required when using a windows programmer to stop the new window closing once the program has been run.

=========

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*A program designed to plot the trajectory of a rocket on a mission to land on the moon*/

int main(int argc, char *argv[])
{
const float pi=3.1415; /*This defines the value of pi*/
const float u=0.0066; /*This defines the intial velocity of rocket and is equivalent to v0*/
const float angle=((89.9*pi)/180.0); /*This defines the angle of projection*/
const float xe=0, ye=0; /*These are the Earth's coords: (0, 0) - ie. where it's centre lies*/
const float xm=0, ym=225.7; /*These are the Moon's coords: (0, 225.7)*/
const float Mm=1.0, Me=83.3; /*These are the masses of the Moon & Earth Respectively*/
const float Rm=1.0, Re=3.7; /*These are the radii of the Moon & Earth Respectively*/
const float G=(9.63*pow(10, -7)); /*This is the Gravitational constant)*/

int dt=10; /*This is the step size delta.t*/
float xr=0, yr=3.7; /*These are the coords of the Rocket)*/
float vx, vy; /*These are the current velocities in both x & y directions*/
float ax=0, ay=0; /*These are the accelerations in both x & y directions*/
float rrx, rry; /*These are the components of the vector from Earth to Rocket*/
float rmx, rmy; /*These are the components of the vector from Earth to Moon*/
/*Do NOT confuse with Re & Rm which are radii*/

vx = u*cos(angle);
vy = u*sin(angle);

int t=0;
/*Using a while loop*/
while(pow((pow(xm-xr,2)) + (pow(ym-yr,2)),1/2)>=Rm)
/*When the distance to the centre of the moon < radius of the moon,*/
/*then the rocket has landed and the calculation can be stopped*/
{
t=t+dt; /*Time increases by 10 each time*/
rmx = xm - xe;
rmy = ym - ye;
rrx = xr - xe;
rry = yr - ye;
ax = (((-G*Me)*xr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(xr-xm))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in x-direction*/
ay = (((-G*Me)*yr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(yr-ym))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in y-direction*/
vx=vx+(dt*ax); /*Step 4 of Euler's method*/
vy=vy+(dt*ay); /*Step 4 of Euler's method*/
xr=xr+(dt*vx); /*Step 5 of Euler's method*/
yr=yr+(dt*vy); /*Step 5 of Euler's method*/
printf("Acceleration x: %f\n", ax);
printf("Acceleration y: %f\n", ay);
printf("Time: %d\n", t);
printf("X velocity: %f\n", vx);
printf("Y velocity: %f\n", vy);
printf("X coord: %f\n", xr);
printf("Y coord: %f\n", yr);
system("PAUSE");
}

system("PAUSE"); /*This code is specific to a Windows system*/
return 0;
}

======================

Thanks once again

 

[cronxeh]

ax = (((-G*Me)*xr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(xr-xm))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in x-direction*/

why does it say..
(pow(((rrx*rrx)+(rry*rry)),(3/2)))) ?

 

[NateTG]

A suggestion on presentation: use the [ code ] tags since they make code much more legible because it preserves indentation.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*A program designed to plot the trajectory of a rocket on a mission to land on the moon*/

int main(int argc, char *argv[])
{
  const float pi=3.1415; /*This defines the value of pi*/
  const float u=0.0066; /*This defines the intial velocity of rocket and is equivalent to v0*/
  const float angle=((89.9*pi)/180.0); /*This defines the angle of projection*/
  const float xe=0, ye=0; /*These are the Earth's coords: (0, 0) - ie. where it's centre lies*/
  const float xm=0, ym=225.7; /*These are the Moon's coords: (0, 225.7)*/
  const float Mm=1.0, Me=83.3; /*These are the masses of the Moon & Earth Respectively*/
  const float Rm=1.0, Re=3.7; /*These are the radii of the Moon & Earth Respectively*/
  const float G=(9.63*pow(10, -7)); /*This is the Gravitational constant)*/
  
  int dt=10; /*This is the step size delta.t*/
  float xr=0, yr=3.7; /*These are the coords of the Rocket)*/
  float vx, vy; /*These are the current velocities in both x & y directions*/
  float ax=0, ay=0; /*These are the accelerations in both x & y directions*/
  float rrx, rry; /*These are the components of the vector from Earth to Rocket*/
  float rmx, rmy; /*These are the components of the vector from Earth to Moon*/
                         /*Do NOT confuse with Re & Rm which are radii*/ 
  
  vx = u*cos(angle);
  vy = u*sin(angle);
  
  int t=0;
  /*Using a while loop*/
   while(pow((pow(xm-xr,2)) + (pow(ym-yr,2)),1/2)>=Rm)
  /*When the distance to the centre of the moon < radius of the moon,*/
   /*then the rocket has landed and the calculation can be stopped*/
    {
        t=t+dt; /*Time increases by 10 each time*/
        rmx = xm - xe;                
        rmy = ym - ye;
        rrx = xr - xe;
        rry = yr - ye;
        ax = (((-G*Me)*xr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(xr-xm))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in x-direction*/
        ay = (((-G*Me)*yr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(yr-ym))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in y-direction*/
        vx=vx+(dt*ax); /*Step 4 of Euler's method*/
        vy=vy+(dt*ay); /*Step 4 of Euler's method*/
        xr=xr+(dt*vx); /*Step 5 of Euler's method*/
        yr=yr+(dt*vy); /*Step 5 of Euler's method*/
        printf("Acceleration x: %f\n", ax);
        printf("Acceleration y: %f\n", ay);        
        printf("Time: %d\n", t);
        printf("X velocity: %f\n", vx);
        printf("Y velocity: %f\n", vy);
        printf("X coord: %f\n", xr);
        printf("Y coord: %f\n", yr);
        system("PAUSE");
    }
  
  system("PAUSE");	/*This code is specific to a Windows system*/
  return 0;
}

I don't know if your code compiles, but some things you might consider:

ax = (((-G*Me)*xr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(xr-xm))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in x-direction*/        
ay = (((-G*Me)*yr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(yr-ym))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in y-direction*/

Is pretty hard to read. and not actually all that efficient considering that you're doing the exponentiation (your most expensive operation) multiple times. Consider using de and dm - variables for the distance to the Earth and moon respectively instead:

de = sqrt(rrx*rrx+rry*rry);
dm = sqrt((rrx-rmx)*(rrx-rmx) + (rry-rmy)*(rry-rmy));

A good compiler is likely to recognize squares and whatnot and perform the optimization inteself, but I find that more readable than, say,

de = pow(pow(rrx,2)+pow(rry,2),1/2);

Then you can write the acceleration equations (assuming Neutonian gravity):

de = sqrt(rrx*rrx+rry*rry);
dm = sqrt((rrx-rmx)*(rrx-rmx) + (rry-rmy)*(rry-rmy));
ax=
 G*Me/(de*de)*(-rrx)/de +
 G*Mm/(dm*dm)*(rmx-rrx)/dm
ay=
 G*Me/(de*de)*(-rry)/de +
 G*Mm/(dm*dm)*(rmy-rry)/dm

Which is much easier to read, and certainly no less efficient.

Of course, the equations, as you have them written, appear to correspond to:

ax=
 G*Me/sqrt(de*de*de) +
 G*Mm/sqrt(dm*dm*dm)
ay=
 G*Me/sqrt(de*de*de) +
 G*Mm/sqrt(dm*dm*dm)

Which is incorrect unless you're using an unusual force of gravity.

Using the radius of the moon also provides a much nicer interface to the loop clause.

Which brings me to the next issue:
Although you may expect the rocket to make it, it's entirely possible that the rocket never reaches the moon. You should really include a maximum time in your while clause.

In fact, I would recommend something more along the lines of:

  for(t=0;t<time_max;t+=time_increment) {
.
.
.
  if(dm < rm) {
     printf("The rocket has hit the moon\n");
     break;
  }
  if(de < re) {
     printf("The rocket has hit the earth\n");
     break;
   }
.
.
.
}
printf ("The position of the rocket is %f,%f at time %d \n", rrx,rry,t);

 

[Hoofbeat]

Wow! Thanks for the tips NateTG. I'll try re-doing my code tomorrow after my lectures and hope that it works - you've certainly given me a lot to think about (I'd already decided to change the while loop to a for loop, but your for loop is even better than mine!).

Cronxeh - I was using the power of 3/2 as a way of simplifying the number of equations I had to do, since I was considering the moduli of vectors etc.

If anyone else has any suggestions as to what may not be allowing my program to work (it returns a value of zero for both accelerations at every iteration) then they'll be gladly appreciated. As I said, I don't have time to do it now (I need to finish some problem sheets) but I'll try it tomorrow and no doubt be back asking more questions.

Thanks once again, you've given me hope that I may be able to do this program!

 

[chronon]

When I tried your program I found that the accelerations weren't actually zero, they were just small, especially in the x-direction, but that is to be expected as the rocket starts off almost vertical, so there won't be much sideways force.

Also the rocket seemed underpowered as it began to fall back towards the centre of the Earth.

I would advise you to do some 'back of the envelope' calculations to see what sort of numbers you would expect to get, and then you can see whether the results of your program agree with these.

 

[NateTG]

This is in response to a PM.

1. Thanks for the tips on the layout of the acceleration formulae. Firstly, just want to check that you think my interpretation of the acceleration is correct, in that even when we are just considering the acc in x-direction that we have to consider the y-components of the position vector in working out the contribution of the gravitational force. You seem to suggest this with your formulae, which is fine, it's just that the demonstrators in the computing lab were not 100% convinced that this was correct (possibly as my code wasn't working and they felt this was where the problem lay).

The distance between two objects is affected by both the x and y coordinates so both components must affect gravity.

2. When I use floats, the accelerations appear almost zero, so I have changed it so that I use 'long floats'. However, no more decimal places occur in the values unless I specifically specify that I want 15d.p, eg. "%.15lf". Is there away around this? I want the values to the most decimal places as is feasible.

"%.15lf" is really only a formatting thing and has very little to do with the way that the number is stored internally. A long float should have more that adequate precision for your needs.

3. I'm not entirely sure how to use your for-loop. I understand that we're considering the calculation to go on for a certain length of time (my demonstrators recommended 10000) as the rocket may return to the Earth before it even reaches the moon. Am I correct in assuming that the if function goes inside the loop and merely terminates the loop if the rocket reaches the moon or returns to the Earth during the iteration, to stop unnecessary calculations?

That's roughly it. Actually I was concerned about the possibility of the rocket getting shot into space and never hitting anything. If I had formatted my code a little better it might be clearer that the if statements should be inside the loop.

Also, do I insert the calculations before the "if" function?

Yes. (Otherwise de and dm would not be set before you used them in the comparison the first time.) There are other options, like putting the conditions inside the for statement but they don't work quite as nicely here.

P.S.
Consider doing something like:

if(0==time%1000) {
  printf("Acceleration x: %f\n", ax);
  printf("Acceleration y: %f\n", ay);
  printf("Time: %d\n", t);
  printf("X velocity: %f\n", vx);
  printf("Y velocity: %f\n", vy);
  printf("X coord: %f\n", xr);
  printf("Y coord: %f\n", yr);
  system("PAUSE");
}

So that you only get a pause after a bunch of iterations.

 

[Hoofbeat]

"%.15lf" is really only a formatting thing and has very little to do with the way that the number is stored internally. A long float should have more that adequate precision for your needs.

That's what I always thought, but if I just use "%lf" then I always get a value of zero for the acceleration as it only gives me a value 0.000000 (ie. 6dp). Any ideas what could be causing this?

 

[NateTG]

That's what I always thought, but if I just use "%lf" then I always get a value of zero for the acceleration as it only gives me a value 0.000000 (ie. 6dp). Any ideas what could be causing this?

It has to do with the formatting string. You could use "%g" instead of "%f" to get scienfic notation.

 

[Hoofbeat]

Thanks once again to everyone who offered some advice. I popped down to computing labs earlier and ran my program and plotted a graph of the trajectory and got the correct curve and discovered that my rocket did actually land on the moon (I increased the maximum time).

So now I'm on to the second part of my project. I now have to change the coords of the rocket so that it's launched from (3.7,0) and calculate the angle required to get the rocket to land on the moon. I spoke to someone who's already done this project and apparantely he created a loop within a loop that incremented the angle until the required angle was found. Is this the optimum way of doing the program?

Thanks once again

 

[NateTG]

.
So now I'm on to the second part of my project. I now have to change the coords of the rocket so that it's launched from (3.7,0) and calculate the angle required to get the rocket to land on the moon. I spoke to someone who's already done this project and apparantely he created a loop within a loop that incremented the angle until the required angle was found. Is this the optimum way of doing the program?

It's a brute force solution, but it apparently works. Depending on what you mean by optimum, it may, or may not be an optimal solution. It certainly does not provide optimal compute time in the general case.