Commit: 3a89a23738d3e253f9d93c2943f8748a479861e9
Parent: 4e8459bb31f45737d053b973d2a786d15033faef
Author: Randy Palamar
Date: Tue, 26 Mar 2024 22:15:16 -0600
rewrite in c
Diffstat:
| D | Photon.cpp | | | 223 | ------------------------------------------------------------------------------- |
| D | Photon.h | | | 34 | ---------------------------------- |
| M | build.sh | | | 5 | +++-- |
| D | main.cpp | | | 133 | ------------------------------------------------------------------------------- |
| A | mcml.c | | | 470 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| D | utils.cpp | | | 54 | ------------------------------------------------------ |
| D | utils.h | | | 60 | ------------------------------------------------------------ |
7 files changed, 473 insertions(+), 506 deletions(-)
diff --git a/Photon.cpp b/Photon.cpp
@@ -1,223 +0,0 @@
-// Photon.cpp: implementation of the Photon class.
-//
-//////////////////////////////////////////////////////////////////////
-
-#include "Photon.h"
-#include "utils.h"
-
-#include <math.h>
-/////////////////////////////////////////
-
-
-double Photon::Rsp=0.0;
-/////////////////////////////////////////////////////////////////////////
-
-
-//////////////////////////////////////////////////////////////////////
-// Construction/Destruction
-//////////////////////////////////////////////////////////////////////
-
-Photon::Photon()
-{
- w=1.0;
- nScatter=0;
- dead=false;
- posx=0.0;
- posy=0.0;
- posz=0.0;
-////////incident direction
- ux=sin(th_in)*n0/n;
- uy=0.0;
- uz=sqrt(1-ux*ux);
-//////////////////////////////////////
- s=0.0;
- db=0.0;
-}
-
-Photon::~Photon()
-{
-
-}
-
-void Photon::Launch()
-{
- double rsp, sinthi, sintht, costhi, costht, raux1, raux2;
- sinthi=sin(th_in);
- costhi=cos(th_in);
- sintht=ux;
- costht=uz;
- raux1=sinthi*costht-costhi*sintht;
- raux1=raux1*raux1;
- raux2=sinthi*costht+costhi*sintht;
- raux2=raux2*raux2;
- rsp=0.5*(raux1/raux2+raux1*(1-raux2)/raux2/(1-raux1));
- Rsp+=w*rsp;
- w=w*(1-rsp);
-}
-
-void Photon::CalStep()
-{
- if (s < ZERO) {
- double ran = random_uniform();
- s = -log(ran + EPS);
- }
-}
-
-void Photon::Caldb()
-{
- if (uz<-ZERO) //travel up to reach upper layer
- {
- db=-posz/uz;
- }
- else if (uz>ZERO) //travel down to reach bottom layer
- {
- db=(d-posz)/uz;
- }
- else
- {
- db=BIGVAL; // A big value
- }
-}
-
-bool Photon::HitBoundary()
-{
- if(db*mut<=s)
- return true;
- else
- return false;
-}
-
-void Photon::Move(double ds) //ds has dimension in cm
-{
- posx+=ds*ux;
- posy+=ds*uy;
- posz+=ds*uz;
- s=s-ds*mut;
-}
-
-void Photon::Absorb() // calculate absorption
-{
- double dw;
- dw=w*mua/mut;//3.24
- w-=dw;
- if (w<0)
- {
- dead=true;
- w=0.0;
- dw=w+dw; //dw equals the small weight dw=(w-dw)+dw
- }
-}
-
-void Photon::Scatter() // calculate next propagation direction
-{
- if (dead)
- return;
-
- double costheta, fei, ran;
- ran = random_uniform();
- if (g != 0) {
- double aa = (1 - g * g) / (1 - g + 2 * g * ran);
- costheta = (1 + g * g - aa * aa) / (2 * g);
- if (costheta < -1)
- costheta = -1;
- else if (costheta > 1)
- costheta = 1;
- } else {
- costheta = 2 * ran - 1;
- } // (3.28)
-
- double sintheta, sinfei, cosfei;
- sintheta = sqrt(1 - costheta * costheta);
- ran = random_uniform();
- fei = 2 * PI * ran; // (3.29)
- cosfei = cos(fei);
- if (fei < PI) {
- sinfei = sqrt(1 - cosfei * cosfei);
- } else {
- sinfei = -sqrt(1 - cosfei * cosfei);
- }
-
- double uux, uuy, uuz;
- if ((SIGN(uz) * uz) <= (1 - ZERO)) {
- uux = sintheta * (ux * uz * cosfei - uy * sinfei)
- / sqrt(1 - uz * uz) + ux * costheta;
- uuy = sintheta * (uy * uz * cosfei + ux * sinfei)
- / sqrt(1 - uz * uz) + uy * costheta;
- uuz = -sintheta * cosfei * sqrt(1 - uz * uz) + uz * costheta;
- // (3.24)
- } else {
- // close to normal propagation
- uux = sintheta * cosfei;
- uuy = sintheta * sinfei;
- uuz = SIGN(uz) * costheta;
- } // (3.30)
- ux = uux;
- uy = uuy;
- uz = uuz;
- nScatter++;
-}
-
-void Photon::ReflectTransmit() // deal with photon across boundary
-{
- double sinalphai,sinalphat,ralphai;
- sinalphai=sqrt(1-(uz)*(uz));
- sinalphat=sinalphai*n/n0; //3.35
- if (sinalphat<ZERO) // close to normal
- {
- ralphai=(n-n0)/(n+n0);
- ralphai=ralphai*ralphai;
- }
- else if (sinalphat>1)// Total reflection
- {
- ralphai=1;
- }
- else
- {
- double A,B;
- A=sinalphai*sqrt(1-sinalphat*sinalphat)-sinalphat*sqrt(1-sinalphai*sinalphai);
- A=A*A;
- B=sinalphai*sqrt(1-sinalphat*sinalphat)+sinalphat*sqrt(1-sinalphai*sinalphai);
- B=B*B;
- ralphai=0.5*(2*A-A*A-A*B)/(B*(1-A));//(3.36)
- }
- double ran = random_uniform();
- if (ran<=ralphai) //rebound to current layer
- {
- uz=-uz;
- }
- else //cross the boundary
- {
- if (uz<-ZERO) //hit upper doundary
- {
- if(posx>-x_offset&&posx<x_offset&&posy>-y_offset&&posy<y_offset)
- {
- u32 xoff = (posx + x_offset) / dx * Rd_xy.Ny;
- u32 yoff = (posy + y_offset) / dy;
- Rd_xy.b[xoff + yoff] += w;
- }
- dead=true;
- }
- else //hit bottom boundary
- {
- dead=true;
- }
- }
-}
-
-void Photon::Termination() //Russian roulette survial test
-{
- if (dead)
- return;
-
- if (w > 1e-4)
- return;
-
- double m = 10;
- double e = random_uniform();
- if (m * e > 1) {
- w = 0;
- dead = true;
- } else {
- w *= m;
- }
-}
diff --git a/Photon.h b/Photon.h
@@ -1,34 +0,0 @@
-// Photon.h: interface for the Photon class.
-// describe the behaviors of photons in Tissue
-//////////////////////////////////////////////////////////////////////
-
-#ifndef _PHOTON_H
-#define _PHOTON_H
-
-class Photon
-{
-public:
- double posx, posy, posz; // current position of photon
- double ux, uy, uz; // current direction of photon
- int nScatter; // # of scattering events experienced
- double w; // weight
- double s; // s is the remaining dimensionless step-size
- double db; // the distance from boundary
- bool dead; // True if photon is dead, otherwise False
-
- static double Rsp; // specular reflectance
-
- Photon();
- virtual ~Photon();
- void Launch(); // calculate specular reflectance at the surface
- void CalStep(); // calculate dimensionless step length
- void Caldb(); // calculate the travel length to reach the boundary
- bool HitBoundary(); // True if hitboudary
- void Move(double ds); // move photon towards boundary
- void Absorb(); // calculate absorption
- void Scatter(); // calculate next propagation direction
- void ReflectTransmit(); // deal with photon across boundary
- void Termination(); // Russian roulette survial test
-};
-
-#endif // _PHOTON_H
diff --git a/build.sh b/build.sh
@@ -1,5 +1,6 @@
#!/bin/sh
+set -x
-srcs="main.cpp Photon.cpp utils.cpp"
+srcs="mcml.c"
-c++ -O3 -Wall -march=native $srcs -o mcml
+clang -O3 -Wall -march=native $srcs -o mcml
diff --git a/main.cpp b/main.cpp
@@ -1,133 +0,0 @@
- /* Hybrid-MC by Li Li, Last modified 03/07
-refrence: Wang, "Biomedical optics", chpater3&6 */
-
-
-/* Programming notes
-1. Use Cartesian coordinates.
-2. Incident beam in x-z plane.
-*/
-
-#include <iostream>
-#include <fstream>
-#include <string>
-
-#define GLOBALORIGIN
-#include "utils.h"
-#undef GLOABLORIGIN
-
-#include "Photon.h"
-
-
-using namespace std;
-////////////////////////////////////////////////
-int main(void)
-{
-////////Input Interface////////////////////
- cout<<"*****************************************"<<endl;
- cout<<"Hybrid Monte Carlo for Homogeneous Slab with Oblique incidence"<<endl;
- cout<<"By Li Li, 04/2007"<<endl;
- cout<<"Refrence: Wang, Biomedical optics, chpater3&6"<<endl;
- cout<<"*****************************************"<<endl;
- cout<<endl;
-
- random_init();
-
- cout<<"File name(Input:*.mci; Output:*.mco): ";
- string sFilename, sInfilename, sOutfilename;
- cin>>sFilename;
- sInfilename=sFilename+".mci";
- sOutfilename=sFilename+".mco";
- ifstream Infile;
- Infile.open(sInfilename.c_str());
- if(Infile.fail())
- {
- cerr<<"Error! Can't open input file. Err Code 1"<<endl;
- exit(1);
- }
- // double P0, R; //Incident energy/power P0, Beam diameter R
- // Infile>>P0>>R;
- int nNum_photons; //# of photons used
- Infile>>nNum_photons;
- Infile>>th_in;
- th_in=th_in/180*PI;
- Infile>>dx>>dy>>Nx>>Ny;
- x_offset=dx*Nx/2;
- y_offset=dy*Ny/2;
- Infile>>n>>mua>>mus>>g>>d;
- mut=mua+mus;
- Infile.close();
-/////////////////////////////////////////////
-
-////////for timing///////////////////////////
- time_t start,end; //simulation beginning and finishing time
- time(&start);
-/////////////////////////////////////////////
- int tmp1;
-
- alloc_mat(&Rd_xy, Nx, Ny);
-
-//////////Propagate photons////////////////////////////////
- for (tmp1=1;tmp1<=nNum_photons;tmp1++)
- {
- Photon photon;
- photon.Launch();
- do
- {
- photon.CalStep();
- photon.Caldb();
- if(photon.HitBoundary())
- {
- photon.Move(photon.db);
- photon.ReflectTransmit();
- }
- else
- {
- photon.Move(photon.s/mut);
- photon.Absorb();
- photon.Scatter();
- }
- photon.Termination();
- } while(!photon.dead);
- if (tmp1%(nNum_photons/10)==0)
- {
- cout<<tmp1<<" photons are done!"<<endl;
- }
- }
-
- time(&end);
- int CompuTime;
- CompuTime=end-start;
- cout<<"This simulation take "<<CompuTime<<" seconds!"<<endl;
-/////////////////////////////////////////////
-
-
-////////Write data before convolution////////
- ofstream f_out(sOutfilename.c_str());
- f_out<<"/////////////Grid coordinates///////////"<<endl;
- f_out<<"x (cm): ";
- for (tmp1=0;tmp1<Nx;tmp1++)
- {
- f_out<<(tmp1+0.5)*dx-x_offset<<' ';
- }
- f_out<<endl;
- f_out<<"y (cm): ";
- for (tmp1=0;tmp1<Ny;tmp1++)
- {
- f_out<<(tmp1+0.5)*dy-y_offset<<' ';
- }
- f_out<<endl;
- f_out<<"/////////////////////////////////////////"<<endl;
-
-
- f_out<<"///////////////Rd_xy/////////////////////"<<endl;
- double scale=nNum_photons*dx*dy;
- double *b = Rd_xy.b;
- for (u32 i = 0; i < Rd_xy.Nx; i++) {
- for (u32 j = 0; j < Rd_xy.Ny; j++)
- f_out << b[j] / scale << ' ';
- f_out << endl;
- b += Rd_xy.Ny;
- }
- f_out<<endl;
- f_out<<"/////////////////////////////////////////"<<endl;
-}
diff --git a/mcml.c b/mcml.c
@@ -0,0 +1,470 @@
+/* Oblique Monte Carlo
+ * Written by Randy Palamar - March 2024
+ * Based on: Hybrid-MC by Li Li
+ * Refrence: Wang, "Biomedical optics", chpater 3 & 6
+ */
+
+/* Details:
+ * - Uses Cartesian coordinates.
+ * - Beam is incident in x-z plane.
+ */
+
+#include <math.h>
+#include <errno.h>
+#include <fcntl.h>
+#include <stdarg.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <unistd.h>
+
+#define SGN(x) ((x) >= 0 ? 1 : -1)
+#define ABS(x) ((x) >= 0 ? x : -x)
+
+#define IS_SPACE(x) ((x) == '\n' || (x) == '\r' || (x) == '\t' || (x) == ' ')
+
+#define ZERO 1.0e-6
+#define BIGVAL 1.0e6
+#define EPS 1.0e-9
+
+typedef uint64_t u64;
+typedef uint32_t u32;
+typedef uint8_t u8;
+typedef double f64;
+typedef ptrdiff_t size;
+
+typedef struct { u8 *data; size len; } s8;
+
+typedef struct {
+ f64 x, y ,z;
+} Vec3;
+
+typedef struct {
+ u32 Nx, Ny;
+ f64 *b;
+} Mat2;
+
+typedef struct {
+ Vec3 pos;
+ Vec3 dir;
+ f64 w;
+ u32 n_scatters;
+ u32 dead;
+} Photon;
+
+/* global data that will not be modified after startup */
+static struct {
+ f64 dx, dy;
+ f64 xoff, yoff;
+ u32 Nx, Ny;
+ u32 N_photons;
+
+ f64 mu_a, mu_s, mu_t;
+ f64 g, d;
+
+ f64 n, n0;
+ f64 theta_i;
+} gctx;
+
+/* these will be modified; FIXME: multithreading */
+static Mat2 Rd_xy;
+static u64 rand_state[2];
+
+static void
+die(const char *fmt, ...)
+{
+ va_list ap;
+
+ va_start(ap, fmt);
+ vfprintf(stderr, fmt, ap);
+ va_end(ap);
+
+ exit(1);
+}
+
+static s8
+read_file(const char *fname)
+{
+ int fd;
+ s8 ret;
+
+ if ((fd = open(fname, O_RDONLY)) < 0)
+ die("open: %s\n", fname);
+
+ ret.len = lseek(fd, 0L, SEEK_END);
+ ret.data = malloc(ret.len);
+ if (ret.data == NULL)
+ die("malloc\n");
+ lseek(fd, 0L, SEEK_SET);
+
+ size rlen;
+ do {
+ rlen = read(fd, ret.data, ret.len);
+ } while (rlen == -1 && errno == EINTR);
+ close(fd);
+
+ if (rlen != ret.len)
+ die("read: %s\n", fname);
+
+ return ret;
+}
+
+static s8
+split_at_space(s8 s)
+{
+ s8 ret;
+
+ /* skip leading white space */
+ while (IS_SPACE(*s.data) && s.len) {
+ s.data++;
+ s.len--;
+ }
+ ret.data = s.data;
+
+ while (!IS_SPACE(*s.data) && s.len) {
+ s.data++;
+ s.len--;
+ }
+ s.data[0] = 0;
+ ret.len = s.data - ret.data;
+ return ret;
+}
+
+/* fills in global ctx based on input */
+static void
+load_input(const char *name)
+{
+ s8 f = read_file(name);
+ s8 s = split_at_space(f);
+
+ char *e;
+ #define ADVANCE_S(_s) do {\
+ _s.data = (u8 *)e + 1;\
+ _s.len = f.len - (f.data - s.data);\
+ _s = split_at_space(_s);\
+ } while (0);
+
+ errno = 0;
+ gctx.N_photons = (u64)strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.theta_i = strtod((char *)s.data, &e) * M_PI / 180;
+ ADVANCE_S(s);
+ gctx.dx = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.dy = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.Nx = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.Ny = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.n = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.mu_a = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.mu_s = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.g = strtod((char *)s.data, &e);
+ ADVANCE_S(s);
+ gctx.d = strtod((char *)s.data, &e);
+
+ gctx.xoff = gctx.dx * gctx.Nx / 2;
+ gctx.yoff = gctx.dy * gctx.Ny / 2;
+ gctx.mu_t = gctx.mu_a + gctx.mu_s;
+
+ gctx.n0 = 1.0;
+
+ #undef ADVANCE_S
+
+ if (errno != 0)
+ die("strtod: errno = %d\n", errno);
+
+ free(f.data);
+}
+
+static void
+random_init(void)
+{
+ struct timespec ts;
+ clock_gettime(CLOCK_REALTIME, &ts);
+ rand_state[0] = (intptr_t)&printf ^ ts.tv_sec
+ ^ ((u64)ts.tv_nsec * 0xAC5533CD);
+ rand_state[1] = (intptr_t)&malloc ^ ts.tv_sec
+ ^ ((u64)ts.tv_nsec * 0xAC5533CD);
+}
+
+static u64
+xoroshiro128plus(u64 s[2])
+{
+ u64 s0 = s[0];
+ u64 s1 = s[1];
+ u64 result = s0 + s1;
+ s1 ^= s0;
+ s[0] = ((s0 << 24) | (s0 >> 40)) ^ s1 ^ (s1 << 16);
+ s[1] = (s1 << 37) | (s1 >> 27);
+ return result;
+}
+
+static f64
+random_uniform(void)
+{
+ return xoroshiro128plus(rand_state) / (f64)UINT64_MAX;
+}
+
+static void
+alloc_mat2(Mat2 *m, u32 x, u32 y)
+{
+ m->b = calloc(x * y, sizeof(f64));
+ if (m->b == NULL)
+ die("calloc\n");
+ m->Nx = x;
+ m->Ny = y;
+}
+
+static void
+launch_photon(Photon *p)
+{
+ p->pos = (Vec3){0, 0, 0};
+ p->dir.x = sin(gctx.theta_i) * gctx.n0 / gctx.n;
+ p->dir.y = 0;
+ p->dir.z = sqrt(1 - p->dir.x * p->dir.x);
+ p->n_scatters = 0;
+ p->dead = 0;
+
+ f64 sin_ti = sin(gctx.theta_i);
+ f64 cos_ti = cos(gctx.theta_i);
+ f64 sin_tt = p->dir.x;
+ f64 cos_tt = p->dir.z;
+ f64 raux1 = sin_ti * cos_tt - cos_ti * sin_tt;
+ raux1 *= raux1;
+ f64 raux2 = sin_ti * cos_tt + cos_ti * sin_tt;
+ raux2 *= raux2;
+ f64 rsp = raux1 / raux2 + raux1 * (1 - raux2) / raux2 / (1 - raux1);
+ rsp *= 0.5;
+ p->w = (1.0 - rsp);
+}
+
+static void
+move_photon(Photon *p, f64 s)
+{
+ p->pos.x += s * p->dir.x;
+ p->pos.y += s * p->dir.y;
+ p->pos.z += s * p->dir.z;
+}
+
+static void
+absorb_photon(Photon *p)
+{
+ f64 dw = p->w * gctx.mu_a / gctx.mu_t; /* eq 3.24 */
+ p->w -= dw;
+ if (p->w < 0)
+ p->dead = 1;
+}
+
+static void
+reflect_or_transmit_photon(Photon *p)
+{
+ f64 sin_ai = sqrt(1 - p->dir.z * p->dir.z);
+ f64 sin_at = sin_ai * gctx.n / gctx.n0; /* eq. 3.35 */
+
+ f64 r_ai;
+ if (sin_at < ZERO) { /* roughly normal */
+ r_ai = (gctx.n - gctx.n0) / (gctx.n + gctx.n0);
+ r_ai *= r_ai;
+ } else if (sin_at > 1.0) { /* TIR */
+ r_ai = 1;
+ } else {
+ f64 A = sin_ai * sqrt(1 - sin_at * sin_at)
+ - sin_at * sqrt(1 - sin_ai * sin_ai);
+ A *= A;
+ f64 B = sin_ai * sqrt(1 - sin_at * sin_at)
+ + sin_at * sqrt(1 - sin_ai * sin_ai);
+ B *= B;
+ r_ai = 0.5 * (2 * A - A * A - A * B) / (B * (1 - A)); /* eq 3.36 */
+ }
+
+ f64 r = random_uniform();
+ if (r <= r_ai) {
+ /* rebound to current layer */
+ p->dir.z = -p->dir.z;
+ } else {
+ /* cross the boundary */
+ if (p->dir.z < -ZERO) {
+ /* hit upper boundary. record reflactance if
+ * we are in bounds of output region */
+ if (p->pos.x > -gctx.xoff && p->pos.x < gctx.xoff &&
+ p->pos.y > -gctx.yoff && p->pos.y < gctx.yoff) {
+ u32 ri = (p->pos.x + gctx.xoff) / gctx.dx;
+ u32 ci = (p->pos.y + gctx.yoff) / gctx.dy;
+ Rd_xy.b[ri * Rd_xy.Ny + ci] += p->w;
+ }
+ }
+ p->dead = 1;
+ }
+}
+
+static void
+scatter_photon(Photon *p)
+{
+ if (p->dead)
+ return;
+
+ f64 cos_t, fei;
+ f64 g = gctx.g;
+ if (g != 0) {
+ f64 r = random_uniform();
+ f64 aa = (1 - g * g) / (1 - g + 2 * g * r);
+ cos_t = (1 + g * g - aa * aa) / (2 * g);
+ if (cos_t < -1)
+ cos_t = -1;
+ else if (cos_t > 1)
+ cos_t = 1;
+ } else {
+ cos_t = 2 * random_uniform() - 1;
+ } /* eq. (3.28) */
+
+ f64 sin_t, sin_fei, cos_fei;
+ sin_t = sqrt(1 - cos_t * cos_t);
+ fei = 2 * M_PI * random_uniform(); /* eq. (3.29) */
+ cos_fei = cos(fei);
+ sin_fei = sin(fei);
+
+ f64 uux, uuy, uuz;
+ if (ABS(p->dir.z) <= (1 - ZERO)) {
+ /* eq. (3.24) */
+ Vec3 u = p->dir;
+ uux = sin_t * (u.x * u.z * cos_fei - u.y * sin_fei)
+ / sqrt(1 - u.z * u.z) + u.x * cos_t;
+ uuy = sin_t * (u.y * u.z * cos_fei + u.x * sin_fei)
+ / sqrt(1 - u.z * u.z) + u.y * cos_t;
+ uuz = -sin_t * cos_fei * sqrt(1 - u.z * u.z) + u.z * cos_t;
+ } else {
+ /* close to normal propagation */
+ /* eq. (3.30) */
+ uux = sin_t * cos_fei;
+ uuy = sin_t * sin_fei;
+ uuz = SGN(p->dir.z) * cos_t;
+ }
+ p->dir.x = uux;
+ p->dir.y = uuy;
+ p->dir.z = uuz;
+ p->n_scatters++;
+
+}
+
+static void
+check_photon_life(Photon *p)
+{
+ if (p->dead)
+ return;
+
+ if (p->w > 1e-4)
+ return;
+
+ f64 m = 10;
+ f64 e = random_uniform();
+ if (m * e > 1) {
+ p->dead = 1;
+ } else {
+ p->w *= m;
+ }
+}
+
+static f64
+next_step(f64 s)
+{
+ if (s < ZERO) {
+ f64 r = random_uniform();
+ s = -log(r + EPS);
+ }
+ return s;
+}
+
+static f64
+step_towards_boundary(Photon *p)
+{
+ f64 db;
+ if (p->dir.z < -ZERO) /* travel up */
+ db = -p->pos.z / p->dir.z;
+ else if (p->dir.z > ZERO) /* travel down */
+ db = (gctx.d - p->pos.z) / p->dir.z;
+ else
+ db = BIGVAL;
+ return db;
+}
+
+static void
+simulate_photon(Photon *p)
+{
+ f64 step = 0, boundary_dist = 0;
+
+ launch_photon(p);
+ do {
+ step = next_step(step);
+ boundary_dist = step_towards_boundary(p);
+ if (boundary_dist * gctx.mu_t <= step) {
+ move_photon(p, boundary_dist);
+ step -= boundary_dist * gctx.mu_t;
+ reflect_or_transmit_photon(p);
+ } else {
+ move_photon(p, step / gctx.mu_t);
+ absorb_photon(p);
+ scatter_photon(p);
+ }
+ check_photon_life(p);
+ } while (!p->dead);
+}
+
+int
+main(int argc, char *argv[])
+{
+ if (argc != 3)
+ die("usage: %s input_file output_file\n", argv[0]);
+
+ load_input(argv[1]);
+ char *outfile = argv[2];
+
+ random_init();
+
+ alloc_mat2(&Rd_xy, gctx.Nx, gctx.Ny);
+
+ /* Propagate Photons */
+ time_t tstart, tend;
+ time(&tstart);
+ Photon p;
+ for (u32 i = 1; i <= gctx.N_photons; i++) {
+ simulate_photon(&p);
+ if (i % (gctx.N_photons / 10) == 0)
+ printf("[%u/%u] photons done!\n", i, gctx.N_photons);
+ }
+ time(&tend);
+ printf("Simulation took: %ld [s]\n", tend - tstart);
+
+
+ FILE *fd = fopen(outfile, "w");
+ if (fd == NULL)
+ die("rip can't open output file: %s\n", outfile);
+
+ fputs("/////////////Grid coordinates///////////\nx (cm): ", fd);
+ for (u32 i = 0; i < gctx.Nx; i++)
+ fprintf(fd, "%e ", (i + 0.5) * gctx.dx - gctx.xoff);
+ fputs("\ny (cm): ", fd);
+ for (u32 i = 0; i < gctx.Nx; i++)
+ fprintf(fd, "%e ", (i + 0.5) * gctx.dy - gctx.yoff);
+ fputs("\n/////////////////////////////////////////\n", fd);
+ fputs("///////////////Rd_xy/////////////////////\n", fd);
+
+ f64 scale = gctx.N_photons * gctx.dx * gctx.dy;
+ f64 *b = Rd_xy.b;
+ for (u32 i = 0; i < Rd_xy.Nx; i++) {
+ for (u32 j = 0; j < Rd_xy.Ny; j++)
+ fprintf(fd, "%e ", b[j] / scale);
+ fputc('\n', fd);
+ b += Rd_xy.Ny;
+ }
+ fputs("\n/////////////////////////////////////////\n", fd);
+ fclose(fd);
+
+ return 0;
+}
diff --git a/utils.cpp b/utils.cpp
@@ -1,54 +0,0 @@
-#include <stdint.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <time.h>
-
-#include "utils.h"
-
-static u64 rand_state[2];
-
-void
-random_init(void)
-{
- struct timespec ts;
- clock_gettime(CLOCK_REALTIME, &ts);
- rand_state[0] = (intptr_t)&printf ^ ts.tv_sec
- ^ ((u64)ts.tv_nsec * 0xAC5533CD);
- rand_state[1] = (intptr_t)&malloc ^ ts.tv_sec
- ^ ((u64)ts.tv_nsec * 0xAC5533CD);
-}
-
-static u64
-xoroshiro128plus(u64 s[2])
-{
- u64 s0 = s[0];
- u64 s1 = s[1];
- u64 result = s0 + s1;
- s1 ^= s0;
- s[0] = ((s0 << 24) | (s0 >> 40)) ^ s1 ^ (s1 << 16);
- s[1] = (s1 << 37) | (s1 >> 27);
- return result;
-}
-
-double
-random_uniform(void)
-{
- return xoroshiro128plus(rand_state) / (double)UINT64_MAX;
-}
-
-void
-alloc_mat(struct Mat *m, int x, int y)
-{
- if (x < 0 || y < 0) {
- fputs("alloc_mat: dimensions must be positive\n", stderr);
- exit(1);
- }
-
- m->b = (double *)calloc(x * y, sizeof(double));
- if (m->b == NULL) {
- fputs("malloc\n", stderr);
- exit(1);
- }
- m->Nx = x;
- m->Ny = y;
-}
diff --git a/utils.h b/utils.h
@@ -1,60 +0,0 @@
-/* Providing global variables and functions*/
-#include <math.h>
-#include <stdint.h>
-
-typedef uint64_t u64;
-typedef uint32_t u32;
-
-//global constants
-#define PI M_PI
-
-#ifndef ZERO
- #define ZERO 1.0e-6 //for small deflection angle
-#endif
-
-#ifndef n0
- #define n0 1.0 //refractive index of ambient medium
-#endif
-
-#ifndef BIGVAL
- #define BIGVAL 1e6//big value
-#endif
-
-#define EPS 1.0e-9
-
-//global variables
-#ifdef GLOBALORIGIN
- #define GLOBAL
-#else
- #define GLOBAL extern
-#endif
-GLOBAL double dx;
-GLOBAL double dy;
-GLOBAL int Nx;
-GLOBAL int Ny;
-GLOBAL double zc;
-GLOBAL double n;
-GLOBAL double mua;
-GLOBAL double mus;
-GLOBAL double mut;
-GLOBAL double g;
-GLOBAL double d;
-
-GLOBAL double th_in; // incident angle in degree
-
-GLOBAL double x_offset;
-GLOBAL double y_offset;
-
-struct Mat {
- u32 Nx, Ny;
- double *b;
-};
-
-GLOBAL struct Mat Rd_xy;
-
-#define SIGN(x) ((x) >= 0 ? 1 : -1)
-
-void random_init(void);
-double random_uniform(void);
-
-void alloc_mat(struct Mat *, int, int);