scattering project for Petnica Pi (Monte Carlo Radiative Transfer Simulation)

config.h

////// PROCESS //////
// How many photons to calc
#define SROUND 1e6

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

////// GEOMETRY //////
#define HEIGHT 80000

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

// angle to hit at OFFSET
#define HANGLE1 0
// delta for angel
#define ANGELD1 90 
// step of angel
#define ANGELS1 5
// angle to hit at OFFSET
#define ANGELOFF1 5

// angle to hit at OFFSET
#define HANGLE2 0
// delta for angel
#define ANGELD2 360
// step of angel
#define ANGELS2 5
// angle to hit at OFFSET
#define ANGELOFF2 5

////// LOGGING //////
// log only right angle
#define LOGh	1
// 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
// log z over time
#define LOGz	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"

double sigma(int l) {
	return 3.4E-32 * pow((589.3/l), 4);
}

double ndens(double h) {
	return NDENS * pow(M_E, -(h)/HEIGHT);
}

int c = 0;

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

		double tau = -log(1-fmod(rand(), 1e5)*1e-5);
				
		#if LOGc != 0
			printf("%d\t%f - %f\n",tau < taumax, chan, comp);
		#endif
		
		*Z = tau*HEIGHT/taumax;
				
		#if LOGz != 0
			printf("%d, %f\n", c++, *Z);
		#endif
		
		// Quit for overstepping Y-boundaries
		if (*Z < 0 || *Z > HEIGHT)
			return;
		
		// Calculate new angle
		if (tau < taumax) {
			*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 z, a1, a2;
	
	#if LOGr != 0
		int hits = 0;
	#endif
	
	for (A1 = HANGLE1; A1 <= HANGLE1+ANGELD1; A1 += ANGELS1)
	for (A2 = HANGLE2; A2 <= HANGLE2+ANGELD2; A2 += ANGELS2) {
		for (W = WAVEL; W <= WAVEL+WAVED; W += WAVES) {
			#if LOGr != 0
				hits = 0;
			#endif
			for (i = 0; i < SROUND; i++) {
			
				z = 0;
				a1 = ANGLE1;
				a2 = ANGLE2;
				
				run(NDENS * sigma(W) * HEIGHT, &z, &a1, &a2);

				#if LOGr != 0
					if (1
					#if LOGh != 0
						&& a1 > A1 && a1 < A1 + ANGELOFF1
						&& a2 > A2 && a2 < A2 + ANGELOFF2
					#endif
						) hits++;
					#if LOGar != 0
						printf("z: %f, a1: %f, a2: %f\n",
							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;
}