A few issues with my first (physics based) program

[jryckman3]

Hello all,
I am doing this project for a class and it is due in a few days.. This is the first "easy" part and there is a much harder second part that I am really unsure about.... but Being a complete newbie to programming in general, I searched all over the web and found this site to be very helpful with the tutorial as well as other sites with example problems that I followed in some parts to try to get me along my way. I believe I am setting up everything correctly, however it keeps spitting out zero when I'm supposed to get a decimal number. I'm trying to be as optimistic as possible, and this certainly has not been a walk in the park for me.

The program involves setting up an isotropic scattering source at a position of 80 cm above the axis shooting into a cylindrical phantom with a height and diameter of 30. The top of the cylinder is at the origin, and it is surrounding the z axis. What I'm trying to do is decide what fraction of the initial photons (which is 10,000, given by the variable "ratio") crosses into the phantom. For some reason, I keep getting 0 as my output! I think I am close to finding an answer, however, I could be very far off in terms of the code. It's hard when it is a very small class and none of my friends have had any experience with C++.

Please find my attached (copied and pasted, rather :P)".cpp" file.

Oh, and the subroutine portion is obviously not my work! I included it so that the whole file is posted. In fact, I found some bits of the code online from different sources (such as for the isotropic source), as looking at examples has helped out a bunch.

#include <math.h>
#include <stdio.h>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <string>
using namespace std;

#define InitRandomGen    (double) RandomGen(0, 1, NULL)
/* Initializes the seed for the random number generator. */     
#define RandomNum        (double) RandomGen(1, 0, NULL)
/* Calls for a random number from the randum number generator. */


/* DECLARE FUNCTION */
double RandomGen(char Type, long Seed, long *Status);  
/* Random number generator */


main() {
	double	x, y, z;    /* photon position */
	double	u, v, w; /* photon trajectory as cosines */
	double	s;          /* step sizes [cm] */
	double	psi;        /* azimuthal angle */
	double	i_photon;   /* current photon */
	short   photon_status;  /* flag = ALIVE=1 or DEAD=0 */
	double ratio;
	double count=0;
	double PI=3.1415926;
	double phantomradius=15;
	double ALIVE=1;
	double DEAD=0;
	double costheta;
	double sintheta;
	
	
	
	/* other variables */
	double	Nphotons;   /* number of photons in simulation */
	double	r;          /* radial position */
	
	Nphotons    = 10000; /* set nuber of photons in simulation */
	ratio = count/Nphotons;
	i_photon = 0;
	InitRandomGen;
	/*RUN*/
	do {
		i_photon += 1;	/* increment photon count */
		photon_status = ALIVE;      /* Launch an ALIVE photon */
		
		x = 0;                      /* Set photon position to z=80 */
		y = 0;
		z = 80;
		
		/* Randomly set photon trajectory to yield an isotropic source. */
		costheta = 2.0*RandomNum - 1.0;   
		sintheta = sqrt(1.0 - costheta*costheta);
		psi = 2*PI*RandomNum; 
		w = costheta;
		double wsqr = w*w;
		double findp = 1-wsqr;
		double p = sqrt(findp);
		u = p*cos(psi);
		v = p*sin(psi);
		
		do {
			s = .01;          /* Step size */
			x += s * u;                        /* Update position */
			y += s * v;
			z += s * w;
			
			/* CHECK to see if photon is alive or dead */
			double xyhit;
			if( x*x + y*y <= phantomradius*phantomradius);
			{xyhit = 1;}
			double zhit; 
			if (z<=0);
			{zhit = 1;}
			if (w>0);
			{
				photon_status=DEAD; /* If photon is moving in positive Z direction, then its importance is 0 */
			break;
			}
			if (r>phantomradius);
			{
				photon_status=DEAD; /* If radial location of photon is greater than 15, then its importance is 0 */
			break;
			}
			double count;
			if (xyhit == 1 && zhit == 1) {count++, photon_status=DEAD;}
		}
		while (photon_status == ALIVE);
	} /* end RUN */
	while (i_photon < Nphotons);
	string mystring="Hits in phantom over initial number of photons: ";
	cout << mystring;
	cout << ratio;
	cout << '\n';
	return 0;
}
/* end of main */

/* SUBROUTINES */

/**************************************************************************
 *	RandomGen
 *      A random number generator that generates uniformly
 *      distributed random numbers between 0 and 1 inclusive.
 *      The algorithm is based on:
 *      W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P.
 *      Flannery, "Numerical Recipes in C," Cambridge University
 *      Press, 2nd edition, (1992).
 *      and
 *      D.E. Knuth, "Seminumerical Algorithms," 2nd edition, vol. 2
 *      of "The Art of Computer Programming", Addison-Wesley, (1981).
 *
 *      When Type is 0, sets Seed as the seed. Make sure 0<Seed<32000.
 *      When Type is 1, returns a random number.
 *      When Type is 2, gets the status of the generator.
 *      When Type is 3, restores the status of the generator.
 *
 *      The status of the generator is represented by Status[0..56].
 *
 *      Make sure you initialize the seed before you get random
 *      numbers.
 ****/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC 1.0E-9

double RandomGen(char Type, long Seed, long *Status){
	static long i1, i2, ma[56];   /* ma[0] is not used. */
	long        mj, mk;
	short       i, ii;
	
	if (Type == 0) {              /* set seed. */
		mj = MSEED - (Seed < 0 ? -Seed : Seed);
		mj %= MBIG;
		ma[55] = mj;
		mk = 1;
		for (i = 1; i <= 54; i++) {
			ii = (21 * i) % 55;
			ma[ii] = mk;
			mk = mj - mk;
			if (mk < MZ)
				mk += MBIG;
			mj = ma[ii];
		}
		for (ii = 1; ii <= 4; ii++)
			for (i = 1; i <= 55; i++) {
				ma[i] -= ma[1 + (i + 30) % 55];
				if (ma[i] < MZ)
					ma[i] += MBIG;
			}
		i1 = 0;
		i2 = 31;
	} else if (Type == 1) {       /* get a number. */
		if (++i1 == 56)
			i1 = 1;
		if (++i2 == 56)
			i2 = 1;
		mj = ma[i1] - ma[i2];
		if (mj < MZ)
			mj += MBIG;
		ma[i1] = mj;
		return (mj * FAC);
	} else if (Type == 2) {       /* get status. */
		for (i = 0; i < 55; i++)
			Status[i] = ma[i + 1];
		Status[55] = i1;
		Status[56] = i2;
	} else if (Type == 3) {       /* restore status. */
		for (i = 0; i < 55; i++)
			ma[i + 1] = Status[i];
		i1 = Status[55];
		i2 = Status[56];
	} else
		puts("Wrong parameter to RandomGen().");
	return (0);
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC

[mordekai]

I'm not sure, if the problem is really that simple, but:

You never update the variable ratio. In line 43 it is set to 0 (because count is zero) and in line 97 it is outputted and in between it has never been assigned another value. Probably you need the command ratio=count/NPhotons; after your loops...

[jryckman3]

Thanks for the response, mordekai. I've moved the ratio=count/Nphotons after the loops and am still getting zero... Is there anything else that might be a problem? I'm pretty sure either A) the loop does not work at all, ie every particle is named dead through the first loop and not counted B) there is something wrong with the isotropic point source or C) The random number generator that I have used is not working and needs to be redefined in my own words.. I will try to see if I can put together a new random number generator and see if that gets me anywhere.

[Disch]

I just skimmed for things which might be a problem. This stuck out, but I don't know if this is what's causing your troubles:

double ALIVE=1;
double DEAD=0;
//...
short   photon_status;  /* flag = ALIVE=1 or DEAD=0 */
//...
while (photon_status == ALIVE);

This is a no-no. Floating point variables (such as doubles) are lossy. photon_status == ALIVE will only return true if photon_status matches ALIVE exactly which might not be the case because of rounding or truncation issues involved in floating point computations. For example... 1.0 might not be exactly 1.0, but might be something like 1.000000000001. Because of this you should never use floating points for this kind of check, and should instead stick to lossless integers (or bools, since this is a boolean check).

Casting between short/double probably isn't helping either.