scattering project for Petnica Pi (Step by Step Photon 3d Simulation)

config.h

////// PROCESS //////
// distance to move per CR | length unit
#define STEP -5000
// How many photons to calc
#define SROUND 500000

////// PHYSICAL //////
// starting angel 1 (°)
#define ANGLE1 315
// starting angel 2(°)
#define ANGLE2 315
// number density of atmosphere (particles per L²)
#define NDENS 0.025e27 
// cross section (L²)
#define CSECT 10000
// wave length (nm)
#define WAVEL 450

////// GEOMETRY //////
#define HEIGHT 1000000
#define WIDTH 1000000

////// STOCHASTIC //////
// delta wave length
#define WAVED 200
// step for wave length
#define WAVES 50

// angle to hit at OFFSET
#define HANGLE 0
// delta for angel
#define ANGELD 180
// step of angel
#define ANGELS 20
// angle to hit at OFFSET
#define ANGELOFF 20

// how wide is hitting area
#define WATCHX 1000000
// offset of HA from left
#define OFFSETX 0
// how wide is hitting area
#define WATCHZ 1000000
// offset of HA from left
#define OFFSETZ 0

////// LOGGING //////
// log all
#define LOGa	0
// log all steps
#define LOGas	0
// log scattering on floor
#define LOGs	0
// log x, y per CR (csv)
#define LOGsz	0
// log x, z per CR (csv)
#define LOGsy	0
// log y, z per CR (csv)
#define LOGsx	0
// log x, y, z per CR (csv)
#define LOGsa	0
// log x, z for hits (csv)
#define LOGh	0
// log stats
#define LOGr	1
// log stats (csv)
#define LOGrc	1
// log chance data
#define LOGc	0
// log more info with stats
#define LOGar	0

Makefile

all: $(ALL)
	 cc scatter.c -O9 -lm -o scatter

scatter.c

/** 3D Modular Phonon Scattering Simulation - scatter.c
  * 
  * "config.h" contains all physical, geometrical, stochastical process
  * and logging information and values.
  **/

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

#include "config.h"

// simulate photon path with specified wave length and angle
void run(int lambda, double* X, double* Y, double* Z, double *A1, double *A2) {
	// One loop = One movement calculation
	while (1) {

		#if LOGas != 0
			printf("(%f,\t%f,\t%f)\t(%f,\t%f)\n",
				*X, *Z, *Y, *A1, *A2);
		#endif
		
		// calculate new photon position trigonometrically
		*X += STEP * sin(*A1) * cos(*A2);
		*Y += STEP * sin(*A1) * sin(*A2);
		*Z += STEP * cos(*A1);
		
		// Calculate news positions when overstepping (X || Z)-boundaries
		while (*X < 0)
			*X = WIDTH + *X;
		while (*X > WIDTH)
			*X = *X - WIDTH;
			
		while (*Y < 0)
			*Y = WIDTH + *Y;
		while (*Y > WIDTH)
			*Y = *Y - WIDTH;

		// Quit for overstepping Y-boundaries
		if (*Z < 0 || *Z > HEIGHT)
			return;
		
		#if LOGsy != 0
			printf("%f, %f\n", *X, *Z);
		#endif
		#if LOGsz != 0
			printf("%f, %f\n", *X, *Y);
		#endif
		#if LOGsx != 0
			printf("%f, %f\n", *Y, *Z);
		#endif
		#if LOGsa != 0
			printf("%f, %f, %f\n", *X, *Y, *Z);
		#endif
		
		// chan: Random chance from 0-1
		// comp: Chance of scattering depending on physical constants
		double chan = fmod(rand(), 1e2)*1e-2,
			comp = 1 - pow(M_E, (-1 *(NDENS * pow(M_E, -(*Z)/HEIGHT)) * 
					abs(STEP) * 33.4E-32 * pow((589.3/lambda), 4)));
					
		#if LOGc != 0
			printf("%d\t%f - %f\n",chan < comp, chan, comp);
		#endif
		
		// Calculate new angle
		if (chan < comp) {
			*A1 = fmod(rand()*1e-2, 360);
			*A2 = fmod(rand()*1e-2, 360);
		}
	}
}

int main() {
	srand(abs(time(NULL)*getpid())*10);
	
	// for counter and hit counter for X positions within
	int A1, A2, W, i;
	double x, y, z, a1, a2;
	
	#if LOGr != 0
		int hits = 0;
	#endif
	
	for (A1 = HANGLE; A1 <= HANGLE+ANGELD; A1 += ANGELS)
	for (A2 = HANGLE; A2 <= HANGLE+ANGELD; A2 += ANGELS) {
		for (W = WAVEL; W <= WAVEL+WAVED; W += WAVES) {
			#if LOGr != 0
				hits = 0;
			#endif
			for (i = 0; i < SROUND; i++) {
			
				x = (double) fmod(rand(), WIDTH);
				y = (double) fmod(rand(), WIDTH);
				z = HEIGHT;
				a1 = ANGLE1;
				a2 = ANGLE2;
				
				run(W, &x, &y, &z, &a1, &a2);
				
				#if LOGa != 0
					printf("(%f, %f, %f)\t(%f, %f)\n",
						x, y, z, a1, a2);
				#endif
				#if LOGs != 0
					printf("%f, %f\n", x, z);
				#endif
				#if LOGr != 0
					if (x > OFFSETX && x < OFFSETX + WATCHX
						&& y > OFFSETZ && y < OFFSETZ + WATCHZ
					#if LOGh != 0
						&& a1 > A1 && a1 < A1 + ANGELOFF
						&& a2 > A2 && a2 < A2 + ANGELOFF
					#endif
						) hits++;
					#if LOGar != 0
						printf("x: %f, z: %f, y: %f, a1: %f, a2: %f\n",
							x, y, z, a1, a2);
					#endif
				#endif
			}
	
			#if LOGr != 0
				#if LOGrc != 0
					printf("%d, %d, %d, %d\n",
						hits, W, A1, A2);
				#else
					printf("P(h): %f%%, h: %d, λ: %dnm, °₁: %d, °₂: %d\n",
						((double) hits/SROUND)*1e2, hits, W, A1, A2);
				#endif
			#endif
		}
		#if LOGr != 0 && LOGrc == 0
			printf("=====================================\n");
		#endif
	}
	return 0;
}