Math is simply 0

OK, I admit right now that I am wrong.

That said, I can find absolutly no reason when this should be. The function excerpted below should compute the Force F force each set of 2 particels in a system. Yet, no matter how many ways I've tried it, the line

Code:

---

F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));

---

always returns 0! I tried it with "G" instead of the 00....0667, with G as both a #declare and a local variable. I tried it with and without the type casts. I simple do not understand.

Here is the function:

Code:

---

void ComputeForces( void )
{
  int i,j,fxsign,fysign; //signs: 0 +, 1 -
  double r,F,theta,dx,dy;
        float G=0.0000000000667;

  ForceSum.i=0;
  ForceSum.j=0;

  for (i=0;i<NUMPARTS;i++)
  {
    for (j=0;j<NUMPARTS;j++)
    {
      if (j!=i)
      {
                                fxsign=0;
                                fysign=0;
                                if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
                                if ((Particles[j].y - Particles[i].y) < 0) fysign++;
                                dx = abs(Particles[i].x - Particles[j].x);
        dy = abs(Particles[i].y - Particles[j].y);
        theta = atan(dy/dx); //theta in radians
        r = sqrt((dx * dx) + (dy * dy)); //distance
                                F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
                                printf("F= %f, G= %f m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
        //dx = cos(rad(theta)) * r;
        //dy = sin(rad(theta)) * r;
                                dx = cos(rad(theta)) * F;
        dy = sin(rad(theta)) * F;
                                //if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,dx,dy,theta); //###DEBUG
                                if (fxsign==0)
                                        ForceSum.i += dx;
                                else
                                        ForceSum.i -= dx;

                                if (fysign==0)
                ForceSum.j += dy;
                                else
                                        ForceSum.j -= dy;
                                //if (i==2) printf("forcesum %f, %f\n",ForceSum.i,ForceSum.j); //###DEBUG
      }
    }
    MoveForce[i].i = ForceSum.i;
    MoveForce[i].j = ForceSum.j;
  }
}

---

And here is the whole source:

Code:

---

/* nbody.c by Jack Carrozzo */

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

#define NUMPARTS 3 //how many particles
#define GRID_X  600
#define GRID_Y  400
#define MASS    1000
#define TIMESTEP 0.04 //seconds per step
#define rad(d)  (4.1415926 * d) / 180
#define DEBUG    0

struct PartType
{
  int    x; //x coord for body
  int    y; //y coord
  int    m; //mass
  double vi;//i-hat of velocity vector
  double vj;//j-hat "                "
};

struct ForceType
{
  double i;
  double j;
};

struct PartType Particles[NUMPARTS];
struct ForceType ForceSum;
struct ForceType MoveForce[NUMPARTS]; //where forces are stored until move()
int ImageGrid[GRID_X][GRID_Y];
int Steps,ERR;

void InitParticles( void )
{
  int i,j;

        srand(time()*getpid());
  for (i=0;i<NUMPARTS;i++)
  {
    Particles[i].x = ((double)rand()/((double)pow(2,31))) * GRID_X;
    Particles[i].y = ((double)rand()/((double)pow(2,31))) * GRID_Y;
    Particles[i].m = MASS;
    Particles[i].vi = 0;
    Particles[i].vj = 0;
  }

        for (i=0;i<GRID_X;i++)
        {
                for (j=0;j<GRID_Y;j++)
                {
                        ImageGrid[i][j]=0;
                }
        }
}

void ComputeForces( void )
{
  int i,j,fxsign,fysign; //signs: 0 +, 1 -
  double r,F,theta,dx,dy;
        float G=0.0000000000667;

  ForceSum.i=0;
  ForceSum.j=0;

  for (i=0;i<NUMPARTS;i++)
  {
    for (j=0;j<NUMPARTS;j++)
    {
      if (j!=i)
      {
                                fxsign=0;
                                fysign=0;
                                if ((Particles[j].x - Particles[i].x) < 0) fxsign++;
                                if ((Particles[j].y - Particles[i].y) < 0) fysign++;
                                dx = abs(Particles[i].x - Particles[j].x);
        dy = abs(Particles[i].y - Particles[j].y);
        theta = atan(dy/dx); //theta in radians
        r = sqrt((dx * dx) + (dy * dy)); //distance
                                F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
                                printf("F= %f, G= %f m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);
        //dx = cos(rad(theta)) * r;
        //dy = sin(rad(theta)) * r;
                                dx = cos(rad(theta)) * F;
        dy = sin(rad(theta)) * F;
                                //if (i==2) printf("compareing to %d, r - %f, force x %f, force y %f, theta %f\n",j,r,dx,dy,theta); //###DEBUG
                                if (fxsign==0)
                                        ForceSum.i += dx;
                                else
                                        ForceSum.i -= dx;

                                if (fysign==0)
                ForceSum.j += dy;
                                else
                                        ForceSum.j -= dy;
                                //if (i==2) printf("forcesum %f, %f\n",ForceSum.i,ForceSum.j); //###DEBUG
      }
    }
    MoveForce[i].i = ForceSum.i;
    MoveForce[i].j = ForceSum.j;
  }
}

void MoveParticles( void )
{
  int i;
  double A; //acceleration

  for (i=0;i<NUMPARTS;i++)
  {
    //F=MA, A=F/M
    A = (MoveForce[i].i / Particles[i].m);
    Particles[i].vi += (A * TIMESTEP);
    A = (MoveForce[i].j / Particles[i].m);
                //###DEBUG!!!
                A *= 25;
    Particles[i].vj += (A * TIMESTEP);

                Particles[i].x += (int)(Particles[i].vi);
                Particles[i].y += (int)(Particles[i].vj);

                //if (i==2) printf("part 2, y %f, f %f, %f.\n",Particles[i].vj,MoveForce[i].i,MoveForce[i].j);
  }
  Steps++;
}

void DataOut( void )
{
  int i,j;
        char filename[100];

        if (DEBUG)
        {
                printf("\n  [+] Time: %d [+]\n",Steps);
                for (i=0;i<NUMPARTS;i++)
                {
                        printf("----------------------\n");
                        printf("Particle Num    : %d\n",i);
                        printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
                        printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
                        printf("----------------------\n");
                }
        }

        CreateBitmap(GRID_X,GRID_Y);

        for (i=0;i<NUMPARTS;i++)
  {
    for (j=0;j<NUMPARTS;j++)
    {
                        SetPixel(i, j, 0, 0, 0); //blackout image
                }
        }

        for (i=0;i<NUMPARTS;i++)
        {
                if (Particles[i].x >= 0 && Particles[i].y >= 0 && Particles[i].x <= GRID_X && Particles[i].y <= GRID_Y)
                {
                        SetPixel((Particles[i].x), (Particles[i].y), 255, 255, 255);
                        SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 255, 255);
                        SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 255, 255, 255);
                        SetPixel((Particles[i].x), (Particles[i].y -1), 255, 255, 255);

                        /*if (i==0)
                        {
                                SetPixel((Particles[i].x), (Particles[i].y), 255, 100, 100);
                                SetPixel((Particles[i].x + 1), (Particles[i].y), 255, 100, 100);
                                SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 255, 100, 100);
                                SetPixel((Particles[i].x), (Particles[i].y -1), 255, 100, 100);
                        }
                        if (i==1)
                        {
                                SetPixel((Particles[i].x), (Particles[i].y), 100, 255, 100);
                                SetPixel((Particles[i].x + 1), (Particles[i].y), 100, 255, 100);
                                SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 100, 255, 100);
                                SetPixel((Particles[i].x), (Particles[i].y -1), 100, 255, 100);
                        }
                        if (i==2)
                        {
                                SetPixel((Particles[i].x), (Particles[i].y), 100, 100, 255);
                                SetPixel((Particles[i].x + 1), (Particles[i].y), 100, 100, 255);
                                SetPixel((Particles[i].x + 1), (Particles[i].y - 1), 100, 100, 255);
                                SetPixel((Particles[i].x), (Particles[i].y -1), 100, 100, 255);
                        }*/
                }
        }

        sprintf(filename,"img/step%d.bmp",Steps);

        WriteBitmapToFile(&filename);

        DeleteBitmap();
}

//JACK: then take ImageGrid to setpoint(), REMEMBER TO zero IMAGE GRID!!



//go to bottom to see main

// bitmap.c
// --------
// author: Levi Lansing
// email: [email protected]
// last updated: 10/15/03
//
// DESCRIPTION:
// - a simple set of bitmap opperations
// - outputs a standard windows bitmap file

//#include <stdio.h>

typedef int BOOL;
#define TRUE        1
#define FALSE        0

typedef unsigned short WORD;
typedef unsigned long DWORD;
typedef unsigned char BYTE;


//bitmap headers
typedef struct tag_BITMAPFILEHEADER {
        WORD    bfType;
        DWORD      bfSize;
        WORD    bfReserved1;
        WORD    bfReserved2;
        DWORD      bfOffBits;
} BITMAPFILEHEADER;

typedef struct tag_BITMAPINFOHEADER{
        DWORD      biSize;
        DWORD      biWidth;
        DWORD      biHeight;
        WORD    biPlanes;
        WORD    biBitCount;
        DWORD      biCompression;
        DWORD      biSizeImage;
        DWORD      biXPelsPerMeter;
        DWORD      biYPelsPerMeter;
        DWORD      biClrUsed;
        DWORD      biClrImportant;
        DWORD      bmiColors;
} BITMAPINFOHEADER;



//some global variables to keep track of the bitmap info
//don't mess with this stuff
DWORD b_height;
DWORD b_width;
DWORD b_widthbytes;
BYTE* b_pBits = NULL;

// CALL THIS FIRST!
// creates a bitmap and allocates the memory for it
void CreateBitmap(int width,int height)
{
        if(b_pBits)
                free(b_pBits);

        b_width = width;
        b_height = height;
        b_widthbytes = b_width*3;
        if(b_width%4)
                b_widthbytes += 4-(b_width%4);
        b_pBits = (BYTE*)malloc(b_widthbytes*height);
}

//frees the allocated memory for the bitmap
void DeleteBitmap()
{
        if(b_pBits)
                free(b_pBits);
        b_pBits = NULL;
}

//sets the color of a single pixel in the bitmap
// colors red,green,and blue accept values between 0 (dark) and 255 (light) only!
void SetPixel(int x, int y, BYTE red, BYTE green, BYTE blue)
{
        if(!b_pBits)
        {
                printf("Error: you forgot to create the bitmap first!\n");
                return;
        }
        if(x>=b_width || y>=b_height || x<0 || y<0)
        {
                printf("pixel (%d,%d) is out of bounds- bitmap is only %dx%d!\n",x,y,b_width,b_height);
                return;
        }

        BYTE* pBits = b_pBits + y*b_widthbytes + x*3;
        pBits[0] = (BYTE) blue;
        pBits[1] = (BYTE) green;
        pBits[2] = (BYTE) red;
}

BOOL WriteBitmapToFile(char* filename)
{
        if(!b_pBits)
        {
                printf("Error: you forgot to create the bitmap first!\n");
                return;
        }

        BITMAPFILEHEADER bmfh;
        BITMAPINFOHEADER bmi;
        int i;

        FILE* fp = fopen(filename,"wb");
        if(!fp)
                return FALSE;

        bmfh.bfType = 0x4D42;
        bmfh.bfSize = 40+14+b_widthbytes*b_height;
        bmfh.bfReserved1 = 0;
        bmfh.bfReserved2 = 0;
        bmfh.bfOffBits = 14+40;

        //write the file header
        fwrite(&bmfh.bfType,2,1,fp);
        fwrite(&bmfh.bfSize,4,1,fp);
        fwrite(&bmfh.bfReserved1,2,1,fp);
        fwrite(&bmfh.bfReserved2,2,1,fp);
        fwrite(&bmfh.bfOffBits,4,1,fp);

        bmi.biSize = 40;
        bmi.biBitCount = 24;
        bmi.biClrImportant = 0;
        bmi.biClrUsed = 0;
        bmi.biCompression = 0;
        bmi.biHeight = b_height;
        bmi.biPlanes = 1;
        bmi.biSizeImage = b_width*b_height*3;
        bmi.biWidth = b_width;
        bmi.biXPelsPerMeter = 3780;
        bmi.biYPelsPerMeter = 3780;
        bmi.bmiColors = 0;

        //write the bitmap header
        fwrite(&bmi.biSize,40,1,fp);

        //bitmap files are upside down! flip it over!
        BYTE* pBits = b_pBits + (b_height-1)*b_widthbytes;
        for(i=0; i<b_height; i++)
        {
                fwrite((void*)pBits,b_widthbytes,1,fp);
                pBits = pBits - b_widthbytes;
        }
        fclose(fp);

        return TRUE;
}

int main( int argc, char *argv[] )
{
        int i, iters;

        printf("Simple NBODY solver by Jack Carrozzo.\n\n");

        if (argc!=2)
        {
                printf("Usage: ./nbody <iterations\n");
                exit(1);
        }

        iters = atoi(argv[1]);

        InitParticles();
        printf("\n  [+] Time: %d [+]\n",Steps);
                for (i=0;i<NUMPARTS;i++)
                {
                        printf("----------------------\n");
                        printf("Particle Num    : %d\n",i);
                        printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
                        printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
                        printf("----------------------\n");
                }

        DataOut();

        for (i=0;i<iters;i++)
        {
                ComputeForces();
                MoveParticles();
                DataOut();
        }

        //Final
        printf("\n  [+] Time: %d [+]\n",Steps);
                for (i=0;i<NUMPARTS;i++)
                {
                        printf("----------------------\n");
                        printf("Particle Num    : %d\n",i);
                        printf("Particle Position: (%d,%d)\n",Particles[i].x,Particles[i].y);
                        printf("Particle Velocity: (%f,%f)\n",Particles[i].vi,Particles[i].vj);
                        printf("----------------------\n");
                }

        printf("Done.\n");

  return 0;
}

---

Thank you very much, sorry I'm stupid,

-Jack Carrozzo
jack [{@}] crepinc.com

Quote:

---

Originally Posted by crepincdotcom

OK, I admit right now that I am wrong.

That said, I can find absolutly no reason when this should be. The function excerpted below should compute the Force F force each set of 2 particels in a system. Yet, no matter how many ways I've tried it, the line

Code:

---

F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));

---

always returns 0!

---

Maybe it's not zero. Try using %e or %g for the print statement:

Code:

---


  printf("F= %e, G= %e m1= %d, m2= %d, r= %f\n",F,G,Particles[i].m,Particles[j].m,r);

---

Regards,

Dave

As a test try setting G to a larger number, say 1.

Off the cuff without really looking I wouldn't be surpised if the number was too small which grav forces are between particles (unless the particals are say the earth and moon)

float doesn't "see" precision beyond 6 digits. limits.h will show you FLT_DIG which tells you how many digits of precision float variables have on your system.

use double.

He is using a double Jim

Quote:

---

Originally Posted by Thantos

He is using a double Jim

---

But the he's using %f to print numbers whose magnitudes are are less than 1.0e-6. So it shows 0.000000.


with %e everywhere:

Code:

---

F= 3.268262e-10, G= 6.670000e-11 m1= 1000, m2= 1000, r= 451.756572

  [+] Time: 2 [+]
----------------------
Particle Num    : 0
Particle Position: (328,185)
Particle Velocity: (4.218655e-13,-8.258777e-14)
----------------------
----------------------
Particle Num    : 1
Particle Position: (38,375)
Particle Velocity: (4.923984e-13,-1.053544e-13)
----------------------
----------------------
Particle Num    : 2
Particle Position: (428,147)
Particle Velocity: (-2.067952e-29,4.846761e-30)

---

Regards,

Dave

Well that has other problems also since the exponent part of a double has more bits then a float. So basically printf will be decoding it improperly

Quote:

---

Originally Posted by Thantos

Well that has other problems also since the exponent part of a double has more bits then a float. So basically printf will be decoding it improperly

---

printf() works with %e, %f, %g whether it's given a double or float (since float arguments are always promoted to double whenever printf() is called).

And I don't propose to determine whether the answers are correct; I just wanted to point out how he can get some output from his program so that he can check the results.

Regards,

Dave

Thanks Guys, so that fixes the displayed problem. But somewhere in the I still have type issues, because even after several iterations, all the ForceSums are still zero or my personal favorite, negative zero. Any ideas?

By the way Thantos, that value for G is the gravitational constant, and when multiplyed by the masses of the 2 bodies it becomes sufficient enough to putt the bodies together.

Thanks again,

-Jack Carrozzo
jack [{@}] crepinc.com

You mean to say that when you keep multiplying a number by a decimal less than one, it eventually goes to zero?! No! Tell me it's not so!

Quzah.

Quote:

---

Originally Posted by quzah

You mean to say that when you keep multiplying a number by a decimal less than one, it eventually goes to zero?! No! Tell me it's not so!

Quzah.

---

That wasn't the case for the original problem (see my previous post), but it can happen. It's called floating point underflow:

Code:

---

#include <stdio.h>

int main()
{
  double x = 1.0;

  int i;

  for (i = 0; i < 10000; i++) {
    x *= 0.1;
    if (x == 0.0) {
      printf("for i = %d, x = %e\n", i, x);
      break;
    }
  }
  printf("After the loop, i = %d, x = %e\n", i, x);
  return 0;
}

---

Regards,

Dave

I think quzah was being sarcastic

Quote:

---

Originally Posted by Thantos

I think quzah was being sarcastic

---

Ok by me, but I can see how the original poster might have missed it. (By the way, if you use x = -1.0 in my example, you do get to see negative zero.)

Regards,

Dave

So to recap, if I simply use a double rather than float it should work?

Sorry to not have figured this out myself.

Thanks guys, interesting discussion too ;-)

> 0.0000000000667
As you're probably aware now, floating point numbers have a number of decimal digits of precision (in the case of doubles, it's about 15).

In this particular number, you're using just 3 decimal digits (all good so far).

> (double)(Particles[i].m) * (double)(Particles[j].m)))
Imagine the partial result of this is huge, I mean something like
123450000000000.000

Now here comes the grief
> 0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m)))
In attempting to normalise your small number to be on the same level as your large number (so the calculation can be done), it gets moved more than 15 decimal digits away from its original value, and as a result, winds up being zero!.

> F = (((0.0000000000667 * (double)(Particles[i].m) * (double)(Particles[j].m))) / (r * r));
Since it's just a series of multiplies and divides, then rearranging the expression might keep more significant digits overlapping, thus giving a better answer.
Say something like
F = (((0.0000000000667 * (double)(Particles[i].m) / r * (double)(Particles[j].m) / r));


http://cch.loria.fr/documentation/IE...M/goldberg.pdf