Thread: 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

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

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

I think quzah was being sarcastic

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