oblique_mc

Monte Carlo in Single Layer Biological Tissue
git clone anongit@rnpnr.xyz:oblique_mc.git
Log | Files | Refs | Feed | README | LICENSE

Commit: 4e8459bb31f45737d053b973d2a786d15033faef
Parent: 2a683bd8546d85140a2feae4df4d021e772f43bf
Author: Randy Palamar
Date:   Tue, 26 Mar 2024 13:36:41 -0600

delete template garbage

Diffstat:
MPhoton.cpp | 5+++--
MPhoton.h | 1-
Mmain.cpp | 30+++++++++---------------------
Mutils.cpp | 21+++++++++++++++++++--
Mutils.h | 112++++++++++++-------------------------------------------------------------------
5 files changed, 47 insertions(+), 122 deletions(-)

diff --git a/Photon.cpp b/Photon.cpp @@ -9,7 +9,6 @@ ///////////////////////////////////////// -double** Photon::Rd_xy=NULL; double Photon::Rsp=0.0; ///////////////////////////////////////////////////////////////////////// @@ -192,7 +191,9 @@ void Photon::ReflectTransmit() // deal with photon across boundary { if(posx>-x_offset&&posx<x_offset&&posy>-y_offset&&posy<y_offset) { - Rd_xy[int((posx+x_offset)/dx)][int((posy+y_offset)/dy)]+=w; + u32 xoff = (posx + x_offset) / dx * Rd_xy.Ny; + u32 yoff = (posy + y_offset) / dy; + Rd_xy.b[xoff + yoff] += w; } dead=true; } diff --git a/Photon.h b/Photon.h @@ -16,7 +16,6 @@ public: double db; // the distance from boundary bool dead; // True if photon is dead, otherwise False - static double **Rd_xy; static double Rsp; // specular reflectance Photon(); diff --git a/main.cpp b/main.cpp @@ -7,7 +7,6 @@ refrence: Wang, "Biomedical optics", chpater3&6 */ 2. Incident beam in x-z plane. */ - #include <iostream> #include <fstream> #include <string> @@ -19,6 +18,7 @@ refrence: Wang, "Biomedical optics", chpater3&6 */ #include "Photon.h" +using namespace std; //////////////////////////////////////////////// int main(void) { @@ -62,13 +62,9 @@ int main(void) time_t start,end; //simulation beginning and finishing time time(&start); ///////////////////////////////////////////// - int tmp1, tmp2; - -////////initiate recording matrix//////////// - - double** Rd_xy=alloc2D(Nx,Ny,0.0); + int tmp1; - Photon::Rd_xy=Rd_xy; + alloc_mat(&Rd_xy, Nx, Ny); //////////Propagate photons//////////////////////////////// for (tmp1=1;tmp1<=nNum_photons;tmp1++) @@ -125,21 +121,13 @@ int main(void) f_out<<"///////////////Rd_xy/////////////////////"<<endl; double scale=nNum_photons*dx*dy; - for (tmp1=0; tmp1<Nx; tmp1++) - { - for (tmp2=0; tmp2<Ny; tmp2++) - { - f_out<<Rd_xy[tmp1][tmp2]/scale<<' '; - } - f_out<<endl; + 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; - - - - -////////deallocate recording matrix////////// - dealloc2D(Nx,Rd_xy); -///////////////////////////////////////////// } diff --git a/utils.cpp b/utils.cpp @@ -1,8 +1,8 @@ #include <stdint.h> +#include <stdlib.h> +#include <stdio.h> #include <time.h> -typedef uint64_t u64; - #include "utils.h" static u64 rand_state[2]; @@ -35,3 +35,20 @@ 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,13 +1,12 @@ /* Providing global variables and functions*/ #include <math.h> -#include <new> -#include <iostream> -#include <time.h> +#include <stdint.h> + +typedef uint64_t u64; +typedef uint32_t u32; //global constants -#ifndef PI - #define PI 3.14159265 -#endif +#define PI M_PI #ifndef ZERO #define ZERO 1.0e-6 //for small deflection angle @@ -31,8 +30,8 @@ #endif GLOBAL double dx; GLOBAL double dy; -GLOBAL int Nx; -GLOBAL int Ny; +GLOBAL int Nx; +GLOBAL int Ny; GLOBAL double zc; GLOBAL double n; GLOBAL double mua; @@ -41,100 +40,21 @@ 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; -//generate random data -void random_init(void); -double random_uniform(void); +struct Mat { + u32 Nx, Ny; + double *b; +}; -#define SIGN(x) ((x) >= 0 ? 1 : -1) +GLOBAL struct Mat Rd_xy; +#define SIGN(x) ((x) >= 0 ? 1 : -1) -////////allocte dynamic array///////////////////////// -using namespace std; -template <class T> -T* alloc1D(int dim1, T t) -{ - T* p; - try - { - p=new T[dim1]; - } - catch(bad_alloc) - { - cerr<<"Fail to allocate 1D array!"<<endl; - exit(1); - } - for (int tmp1=0;tmp1<dim1;tmp1++) - { - p[tmp1]=t; - } - return p; -} -template <class T> -void dealloc1D(T* p) -{ - delete [] p; -} -/////////// -template <class T> -T** alloc2D(int dim1,int dim2,T t) -{ - T** p; - try - { - p=new T*[dim1]; - for (int tmp1=0;tmp1<dim1;tmp1++) - { - p[tmp1]=alloc1D(dim2,t); - } - } - catch (bad_alloc) - { - cerr<<"Fail to allocate 2D array!"; - } - return p; -} -template <class T> -void dealloc2D(int dim1,T** p) -{ - for (int tmp1=0;tmp1<dim1;tmp1++) - { - dealloc1D(p[tmp1]); - } - delete [] p; -} -///////// -template <class T> -T*** alloc3D(int dim1,int dim2,int dim3,T t) -{ - T*** p; - try - { - p=new T**[dim1]; - for (int tmp1=0;tmp1<dim1;tmp1++) - { - p[tmp1]=alloc2D(dim2,dim3,t); - } - } - catch (bad_alloc) - { - cerr<<"Fail to allocate 3D array!"; - } - return p; -} -template <class T> -void dealloc3D(int dim1,int dim2,T*** p) -{ - for (int tmp1=0;tmp1<dim1;tmp1++) - { - dealloc2D(dim2,p[tmp1]); - } - delete [] p; -} -//////////////////////////////////////////////// +void random_init(void); +double random_uniform(void); +void alloc_mat(struct Mat *, int, int);