Monte Carlo simulations for x-ray shielding

main.c

typedef struct {
	Vec3 position;
	Vec3 direction;
	B32 alive;
} PhotonState;

int main(void) {
	RNG* rng = rng_create();

	// Source: https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z82.html
	F32 pb_mass_attn = 2.014; // cm^2/g from NIST tables at 0.150MeV (150keV)
	F32 pb_density = 11.34; // g/cm^3
	F32 pb_linear_attn = pb_mass_attn * pb_density;

	// Linear Attenuation Coefficients for "component breakdowns":
	//   - Photoelectric Absorption at 150keV
	//   - Compton Scattering at 150keV
	//
	// Source: https://physics.nist.gov/PhysRefData/Xcom/html/xcom1-t.html
	F32 pb_mass_attn_scatter = 1.049e-1f; // cm^2/g
	F32 pb_mass_attn_absorbs = 9.484e-2f; // cm^2/g

	F32 pb_linear_attn_scatter = pb_mass_attn_scatter * pb_density; // 1/cm
	F32 pb_linear_attn_absorbs = pb_mass_attn_absorbs * pb_density; // 1/cm

	// scene
	F32 mat_beg = 0.0000f;
	F32 mat_end = 0.3175f;
	
	F32 sum = 0;
	U64 trials = 1000;
	for (U64 trial = 0; trial < trials; trial += 1) {
		U64 escaped = 0;
		U64 num_photons = 140000;
		for (U64 i = 0; i < num_photons; i += 1) {
			PhotonState photon = {
				.position = vec3(0.0f, 0.0f, 0.0f),
				.direction = vec3(0.0f, +1.0f, 0.0f),
				.alive = true,
			};
			while (photon.alive) {
				// "free path length"
				F32 dd = -logf(rng_next_f32(rng)) / pb_linear_attn;
				Vec3 dpos = vec3_scale(vec3_normalize(photon.direction), dd);
				Vec3 pos2 = vec3_add(photon.position, dpos);

				// Check escaped (infinitely wide, infinitely tall)
				if ((photon.direction.y > 0 && pos2.y > mat_end) ||
					(photon.direction.y < 0 && pos2.y < mat_beg)) {
					escaped += 1;
					photon.alive = false;
				}
				else {
					photon.position = pos2;
					F32 p_absorbed = pb_linear_attn_absorbs / (pb_linear_attn_absorbs + pb_linear_attn_scatter);
					// Absorption
					if (rng_next_f32(rng) < p_absorbed) {
						photon.alive = false;
					}
					// Scattering (change direction)
					else {
						Vec3 ddir = vec3(rng_next_f32(rng), rng_next_f32(rng), rng_next_f32(rng));
						Vec3 dir2 = vec3_normalize(ddir);
						photon.direction = dir2;
					}
				}
			}
		}
		printf("escaped=%ld of %ld (%0.3f%%, trial=%ld of %ld)\n", escaped, num_photons, (F32)escaped / (F32)num_photons * 100.0f, trial, trials);
		sum += (F32)escaped / (F32)num_photons;
	}
	printf("escape rate (avg): %0.3f%%\n", sum / (F32)trials * 100.0f);
}