Dear Qianqian,
In the MCX code is there an option available for a spot focused Gaussian beam similar to;or?
hi Gijs,
mcx supports collimated Gaussian beam (srctype='gaussian'), convergent/divergent gaussian beam (srcdir(4)>0 or <0), or angular Gaussian beam (srctype='zgaussian').
For example, collimated Gaussian beam
a point-focusing Gaussian beam
the focused gaussian beam output screenshot is attached.
Qianqian
Best regards,Gijs Buist--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/35508d2d-2155-4f17-8f7b-d3fcf7cd55e5n%40googlegroups.com.
Hi Qianqian,
Does the focussed Gaussian of MCX work via method b) in the figure below?
I wish to simulated a setup with geometry b) where the focus inside the medium.

Best regards,
Gijs Buist
--
You received this message because you are subscribed to a topic in the Google Groups "mcx-users" group.
To unsubscribe from this topic, visit
https://groups.google.com/d/topic/mcx-users/wauKd1IbEJE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to
mcx-users+...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/mcx-users/f747312e-eb93-9f2a-fb02-cc27fb738bdf%40neu.edu.
Hi Qianqian,
Does the focussed Gaussian of MCX work via method b) in the figure below?
I wish to simulated a setup with geometry b) where the focus inside the medium.
hi Gijs, the one that I showed earlier is method (a). true gaussian beam is not yet supported in mcx.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/d991bbc71f814397a91026933e194ba3%40vu.nl.
Hi Qianqian,
I have written code to determine the start positions and directions of photons to be launched in order to create a true Gaussian beam.
It works by realizing that the Gaussian beam is a hyperboloid of one sheet and that such shapes can be fully described by straight lines.
And needs scrparam1=(w,zf,zr) (radius beam waist at focus, z-position focus, rayleigh range) as input.
I now wish to add this code to the MCX code but am not sure on how exactly.
I believe it should look somewhat like this:
But I am not sure how integrate it with the rest of the MCX code.
Could you tell me how I get xi,yi and zi to be the start position of the photon in the MCX code, as well as how to make vx,vy,vz the directional components?
Many thanks in advance.
Gijs Buist
Hi Qianqian,
I have written code to determine the start positions and directions of photons to be launched in order to create a true Gaussian beam.
It works by realizing that the Gaussian beam is a hyperboloid of one sheet and that such shapes can be fully described by straight lines.
And needs scrparam1=(w,zf,zr) (radius beam waist at focus, z-position focus, rayleigh range) as input.
hi Gijs, thanks for contributing.
currently, all mcx sources can be pointed to any direction, but reading your below code, it seems it is only oriented along the +z direction.
so, in order for me to integrate this code, I would like to understand this better and generalize it for all directions.
can you send me offline a drawing/slides showing how is this source implemented? I just want to get some idea how is this
different from the existing point-focusing Gaussian and how to add it with minimal lines of codes. please email that drawing/schematic
in a reply or to my email directly.
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/007f1be2ee53408b98225afeaa234d70%40vu.nl.
Hi Qianqian,
I indeed started with the simple case of propagating along the z-direction.
Attached are slides explaining the math.
Regards,
Gijs Buist
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/1ea3c320-fb53-4634-cb68-5a7c46d28eea%40neu.edu.
Hi Qianqian,
I indeed started with the simple case of propagating along the z-direction.
Attached are slides explaining the math.
hi Gijs, thank you for your detailed slides! well done and thanks for sharing with the list.
I have a busy schedule today but will take a look during the weekend. I may likely follow up and schedule a meeting with you next week and learn more about your implementation.
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/ca747c60a57d4a59a8a2ad6394d292ea%40vu.nl.
Hi Qianqian,
If you have any questions I would be happy to schedule a meeting.
Additionally I have been thinking on how to create the gaussians for any direction and I believe the best way is by rotation and translation of the ‘standard’ z-propagation centred at the origin result.
Regards,
Gijs
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/417b54c8-8708-eb0e-8317-c51527b3b7ad%40neu.edu.
Hi Qianqian,
If you have any questions I would be happy to schedule a meeting.
Additionally I have been thinking on how to create the gaussians for any direction and I believe the best way is by rotation and translation of the ‘standard’ z-propagation centred at the origin result.
that's what I thought too. we've already have a 3-D vector rotation function and have been using it for launching photons
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L902
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L1201-L1203
can you update your code and add rotation to allow it to point to srcdir?
I am CCing my student Shijie Yan who can answer questions about our part of the code.
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/ebb031604fab491dbcb7f327ee70fb33%40vu.nl.
Dear,
I wrote and tested a first draft of code in Python that should give a true gaussian for arbitrary source position and direction (see below).
First the position and direction for source at (0,0,0) with dir (0,0,1) is calculated and then based on the desired direction a rotation around the y-axis and z-axis (in that order) is preformed on the direction and position.
Then the position is shifted by the desired source position.
I did not fully understand the rotatevector of MCX so the code still needs to be translated to .cu
Best,
Gijs
PYTHON CODE
###########################################################
Import numpy as np
# INPUT PARAMETERS NECESSARY FOR PHOTON LAUNCH OF TRUE GAUSSIAN FROM ARBITRARY POSITION AND ARBITRARY BEAM DIRECTION
#beam parameters
w0 = 1 #radius beam waist 1/e reduction E-field (1/e2 for intensity) at focus
df = -1 # distance between launch plane and focus
zr = 2 #rayleigh range
# direction and position beam for rotation and shift respectively
beam_dir = np.array([1,0,0]) # (x,y,z) beam direction
source_pos = np.array([2,3,0]) # (x,y,z) position center beam source
# GET SOLUTION FOR SOURCE AT ORIGIN PROPAGATING IN +Z DIRECTION
rng = np.random.default_rng()
num_phot=1 # In MCX core it should be a single photon launch function, here arrays of num_phot length
ran1 = rng.uniform(low=0,high=1.0, size=num_phot) # for distance from beam center
ran2 = rng.uniform(low=0,high=1.0, size=num_phot) # for azimuthal angle
#r_samps = r = w0*np.sqrt(-0.5*np.log(ran1)) # no sampling to test rotations reliably
r_samps = np.zeros(num_phot)+1.0 # all at r=1.0
thet_samps = 2*np.pi*ran2 # azimuthal angle of photon before rotation
x0 = r_samps*np.cos(thet_samps) # x-position at focus
y0 = r_samps*np.sin(thet_samps) # y-position at focus
xi,yi = entry_cor(r_samps,thet_samps,zr,df) #positions at launch plane
zi = np.zeros(len(xi)) # launch at plane z=0
vx,vy,vz = dir_zf0(r_samps,thet_samps,zr) # components directional vector photon before rotation
# GET (COSINES AND SINES OF) ROTATION ANGLES TO GO FROM (0,0,1) TO BEAM_DIR ([0],[1],[2])
# Rotation of phi around y-axis and rotation of theta around z-axis in that order
# as rotations around z-axis do not do anything for (0,0,1) vector
# theta is counter-clockwise angle with x-axis (azimuthal angle) (from +x axis to +y axis)
# phi is angle with z-axis towards positive x-axis (zenith) (+z axis to +x axis)
# to prevent division by zero
if beam_dir[2]==1:
0 #no rotations necessary
else:
tmp0 = 1/np.sqrt(beam_dir[0]**2 + beam_dir[1]**2) #1/length in xy plane
tmp1 = 1/np.sqrt(beam_dir[0]**2 + beam_dir[1]**2+beam_dir[2]**2) #1/total_length
ctheta = beam_dir[0]*tmp0
stheta = beam_dir[1]*tmp0
cphi = beam_dir[2]*tmp1
sphi = np.sqrt(beam_dir[0]**2+beam_dir[1]**2)*tmp1
# ROTATE LAUNCH POSITION AND DIRECTION
# first rotation around y-axis then around z-axis.
# Rotating first around z-axis does not do anything due to symmetry
# Replace with the rotatevector of mcx if possible
# rotated launch coordinates, zi part left out as zi=0
xrot = (xi*ctheta*cphi - yi*stheta) #+zi*ctheta*sphi
yrot = (xi*stheta*cphi + yi*ctheta) #+ zi*stheta*sphi
zrot = (-xi*sphi) #+ yi*0 + zi*cphi
# rotated directional vector
vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi
vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi
vzrot = -vx*sphi + vy*0 + vz*cphi
# SHIFT LAUNCH POSITION OF PHOTON WITH SOURCE_POS ([0],[1],[2])
xrot_shift = xrot + source_pos[0]
yrot_shift = yrot + source_pos[1]
zrot_shift = zrot + source_pos[2]
# RETURN FINAL LAUNCH POSITION AND DIRECTION OF PHOTON
phot_launch_pos = np.array([xrot_shift,yrot_shift,zrot_shift])
phot_launch_dir = np.array([vxrot,vyrot,vzrot])
#########################################################################################
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/5e56891c-f30b-4d1b-802d-28ede58f6f58%40neu.edu.
Hi,
I believe that the code for a focused Gaussian with arbitrary direction and position should look something like the code below.
However, I am not sure if the naming is consistent with the rest of the MCX code nor what to do with the resulting coordinates and directions.
Regards,
Gijs
case(MCX_SRC_TRUE_GAUSSIAN): { // Gaussian-beam with spot focus
//hyperboloid of one sheet, doubly ruled surface
// scrparam1=(w,df,zr) beam waist, distance from source to focus, rayleigh range
// srcdir = (dx,dy,dz) direction of gausssian beam
// srcpos = (x,y,z) location of the center of the gaussian source
// Sampling for srcdir = (0,0,1) and srcpos = (0,0,0)
float sphi, cphi;
float phi=TWO_PI*rand_uniform01(t);
sincosf(phi,&sphi,&cphi);
float r; // sampled distance from z-axis at focus
float x0; // sampled x coordinate photon at focus
float y0; // sampled y coordinate photon at focus
r=sqrtf(-0.5f*logf(rand_uniform01(t)))*gcfg->srcparam1.x;
x0=r*cphi;
y0=r*sphi;
float xi; // x position photon at source
float yi; // y position photon at source
float zi=0; // z position photon at source
float t = -gcfc->srcparam1.y/gcfc->scrparam1.z; //parameter to generate photon path from coordinates at focus (depends focal distance and rayleigh range)
xi= r*(cphi-t*sphi);
yi= r*(sphi+tcphi);
float l; //use rsqrtf instead to calculate 1/l ? used to normalize directional vector
float vx; // normalized x part of directional vector photon
float vy; // normalized y part of directional vector photon
float vz; // normalized z part of directional vector photon
l=sqrtf(r*r+gcfc->srcparam1.z*gcfc->srcparam1.z); //use rsqrtf instead to calculate 1/l?
vx=-r*sphi/l;
vy=r*cphi/l;
vz=gcfc->srcparam1.z/l;
// rotate photon direction and position if srcdir!=(0,0,1)
// rotatevector function updates photon direction so cannot be used for position
if( v->z>-1.f+EPS && v->z<1.f-EPS ) {
float tmp0=1.f-srcdir->z*srcdir->z; // if srcdir unit length then this is transverse direction squared
float ctheta = srcdir->x*rsqrt(tmp0); // cos(phi), phi counter clockwise angle with x-axis (rotation of phi around z-axis)
float stheta = srcdir->y*rsqrt(tmp0);
float cphi = srcdir->z; // only in case of unit vector srcdir, theta angle with z-axis (rotation towards +x)
float sphi = sqrt(tmp0); // sqrt(dx*dx+dy*dy)/1
// position photon after rotating (xi,yi,zi=0)=>(xrot,yrot,zrot)
float xrot = (xi*ctheta*cphi - yi*stheta)
float yrot = (xi*stheta*cphi + yi*ctheta)
float zrot = (-xi*sphi)
// add the shifted source position to the photon
xrot+=srcpos->x
yrot+=scrpos->y
zrot+=srcpos->z
// rotate direction of photon
// could be done by rotatevector function?
float vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi
float vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi
float vzrot = -vx*sphi + vy*0 + vz*cphi
// photon position (xrot,yrot,zrot) now need to be changed in the correct object
// photon direction (vxrot,vyrott,vzrot) now need to be changed in the correct object
}
// no rotation happened as srcdir==(0,0,1)
// photon position (xi,yi,zi) now need to be changed in the correct object
// photon direction (vx,vy,vz) now need to be changed in the correct object
hi Gijs,
thanks a lot for the patch. unfortunately my coding time becomes
extremely limited as my dept service load eats most of my time.
can only answer short questions, but for testing, I am CCing my
student Shijie -
Shijie, can you please test Gijs's patch? I think this will be a great addition to mcx! please read the slides that Gijs sent previously to understand the algorithm. if everything works well, please create a PR on github. use new sourcetype string 'hyperboloid' and MCX_SRC_HYPERBOLOID_GAUSSIAN macro.
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/e2a02fb392414ee19029c5932fd35689%40vu.nl.
Hi Shijie,
Attached you will find a complete python script for the hyperboloid source distribution.
As for your questions;
Best,
Gijs
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB7141222327175A623DBE7D2F9CAF9%40PH0PR06MB7141.namprd06.prod.outlook.com.
Hi Shijie,
There should be multiple figures after the one you mentioned. I usually use Spyder to execute the code which gives all of the figures, if the Python command is used the next figure will appear after the previous figure is closed. Therefore a different figure should appear after closing “unrotated launch positions”.
Additionally about the rotatevector in MCX, if I understand correctly it updates the direction vector based the scattering angles (which are local properties).
In that case I do not think it would be useful for rotating the Gaussian beam as you do not know the local rotation (scattering angles) but the global rotation of the whole beam. I therefore determine the (global) sines and cosines of the rotation from (0,0,1) to an arbitrary vector of unit length. After determining those once they can be used in the rotation matrix RzRy (see Wikipedia) to rotate all the direction vectors. Finally as the hyperboloid is an extended object all the starting coordinates have to be rotated too.
Would something work in the MCX code?
case(MCX_SRC_TRUE_GAUSSIAN): { // Gaussian-beam with spot focus
//hyperboloid of one sheet, doubly ruled surface
// scrparam1=(w,df,zr) beam waist, distance from source to focus, rayleigh range
// srcdir = (dx,dy,dz) direction of gausssian beam
// srcpos = (x,y,z) location of the center of the gaussian source
// Sampling for srcdir = (0,0,1) and srcpos = (0,0,0)
float sphi, cphi;
float phi=TWO_PI*rand_uniform01(t);
sincosf(phi,&sphi,&cphi);
float r; // sampled distance from z-axis at focus
r=sqrtf(-0.5f*logf(rand_uniform01(t)))*gcfg->srcparam1.x;
float xi; // x position photon at source
float yi; // y position photon at source
float zi=0; // z position photon at source
float t = -gcfc->srcparam1.y/gcfc->scrparam1.z; //parameter to generate photon path from coordinates at focus (depends focal distance and rayleigh range)
xi= r*(cphi-t*sphi);
yi= r*(sphi+t*cphi);
float l; //use rsqrtf instead to calculate 1/l ? used to normalize directional vector
float vx; // normalized x part of directional vector photon
float vy; // normalized y part of directional vector photon
float vz; // normalized z part of directional vector photon
l=rsqrtf(r*r+gcfc->srcparam1.z*gcfc->srcparam1.z); //use rsqrtf instead to calculate 1/l?
vx=-r*sphi*l;
vy=r*cphi*l;
vz=gcfc->srcparam1.z*l;
// calculate global rotation matrix and rotate photon direction and position if srcdir!=(0,0,1)
if( v->z>-1.f+EPS && v->z<1.f-EPS ) {
// here the sines and cosines for the rotation matrix are determined
// these are global properties and thus the same for all photons
// the way it is now it is calculated every time but maybe there is a good way to change that
float tmp0=1.f-srcdir->z*srcdir->z; // if srcdir unit length then this is transverse direction squared
float ctheta = srcdir->x*rsqrt(tmp0); // cos(phi), phi counter clockwise angle with x-axis (rotation of phi around z-axis)
float stheta = srcdir->y*rsqrt(tmp0);
float cphi = srcdir->z; // only in case of unit vector srcdir, theta angle with z-axis (rotation towards +x)
float sphi = sqrt(tmp0); // sqrt(dx*dx+dy*dy)/1
// here the rotations based on the global rotation matrix are done
// for the photon launch coordinates and directional vector
// position photon after rotating (xi,yi,zi=0)=>(xrot,yrot,zrot)
float xrot = (xi*ctheta*cphi - yi*stheta)
float yrot = (xi*stheta*cphi + yi*ctheta)
float zrot = (-xi*sphi)
// add the shifted source position to the photon
xrot+=srcpos->x
yrot+=scrpos->y
zrot+=srcpos->z
// rotate direction of photon
float vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi
float vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi
float vzrot = -vx*sphi + vy*0 + vz*cphi
// photon position (xrot,yrot,zrot) now need to be changed in the correct object
// photon direction (vxrot,vyrott,vzrot) now need to be changed in the correct object
}
// no rotation happened as srcdir==(0,0,1)
// photon position (xi,yi,zi) now need to be changed in the correct object
// photon direction (vx,vy,vz) now need to be changed in the correct object
}
Best,
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71411B0B98D698A327D786D09CB09%40PH0PR06MB7141.namprd06.prod.outlook.com.
Hi Shijie,
Thanks for your response and if you have any questions about the code I am happy to answer them. As for the cross-section of the hyperboloid, this should indeed remain circular but I have not extensively tested this. Just looking at the 3d plots seems to be of no help as the perspective makes gauging the precise 3d shape difficult. To check this robustly the distance to the optical axis (at the same t-parameter) needs to be calculated for every photon. At the bottom of the attached script this is done for a beam_dir = (0,0,1), (0,1,0) and (1,0,0) (different lines need to be (un)commented for different beams). For arbitrary beam_dir the coordinates of the optical axis for all the t-parameter values need to be determined for which I have currently no function.
I hope this clears things up a bit.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71415442A1B2769E7F3EAA679CB59%40PH0PR06MB7141.namprd06.prod.outlook.com.
Hi Shijie,
Thank you for your work in implementing the hypberboloid gaussian. At the moment I cannot compile from source with my Windows pc’s so testing the github code will have to wait.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71411193D328A7A40CBB70059CB69%40PH0PR06MB7141.namprd06.prod.outlook.com.
Hi Gijs,
We are merging the code with the main branch (https://github.com/fangq/mcx/pull/127/commits). Once it is done, the new feature should be available in the MCX nightly build.
hi Gijs, I have accepted Shijie's pull request, and manually built the windows/linux/mac binaries. feel free to go ahead and download it from our nightly built site and see if it works as expected.
thanks again for making this useful feature available on mcx.
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/PH0PR06MB71416E4D0CC538DC3B92917A9CB89%40PH0PR06MB7141.namprd06.prod.outlook.com.
Hi Qianqian, Shijie
I have installed the nightlybuild and the (true) hyberboloid gaussian seems to work as expected.
Thanks for implementing it in MCX.
Best,
Gijs
From: mcx-...@googlegroups.com <mcx-...@googlegroups.com>
On Behalf Of Qianqian Fang
Sent: vrijdag 15 oktober 2021 05:57
To: mcx-...@googlegroups.com; Shijie Yan <yan....@northeastern.edu>
Subject: Re: [mcx-users] spot focused Gaussian source distribution
On 10/14/21 5:26 PM, Shijie Yan wrote:
Hi Gijs,
We are merging the code with the main branch (https://github.com/fangq/mcx/pull/127/commits). Once it is done, the new feature should be available in the MCX nightly build.
hi Gijs, I have accepted Shijie's pull request, and manually built the windows/linux/mac binaries. feel free to go ahead and download it from our nightly built site and see if it works as expected.
thanks again for making this useful feature available on mcx.
Qianqian
Thanks again for your contribution,
Shijie
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/1fe91fcf-7836-2fc4-454e-adea554546a9%40neu.edu.
Hi Qianqian, Shijie
I have installed the nightlybuild and the (true) hyberboloid gaussian seems to work as expected.
Thanks for implementing it in MCX.
glad to hear! thank you for the initial patch and Shijie for
implementing/debugging!
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/3af89fb87adf48e3a8be4aad4b861dd3%40vu.nl.
hi Gijs,
in an attempt of reducing register counts in mcx kernels, I
optimized the hyperboloid gaussian source implementation by
reusing some existing registers. see
https://github.com/fangq/mcx/commit/7b6e0e065ebda62952c1a979f97080d2ad7a72c6
can you please download the latest nightly build version (https://mcx.space/nightly/, or recompile on your side) and see if the source still functions as expected?
thanks
Qianqian
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/3af89fb87adf48e3a8be4aad4b861dd3%40vu.nl.