diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr index 5ceec336f9..2a606bb777 100644 --- a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr +++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr @@ -19,7 +19,7 @@ * It can be found here: https://inis.iaea.org/collection/NCLCollectionStore/_Public/19/019/19019134.pdf * Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated monitor "Wavelength(Ang) as a function of x_screw(mm)" will be a linear function. The monochromator functions properly. * -* %Example: x_screw=60 Detector: w_monitor_I=1.10881e-06 +* %Example: x_screw=60 Detector: w_monitor_I=0.0721 * %Scan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={ * 2007.08,2406.85,2820.95,3242.17,3659.02,4075.98,4485.76,4901.23, * 5324.84,5733.95,6156.24,6569.68,6979.66,7395.58,7805.61,8222.2, diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md index fce75b8183..b8458b2f9c 100644 --- a/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md +++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/README.md @@ -22,7 +22,7 @@ Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated ## Examples -- **Test: x_screw=60 Detector: w_monitor_I=1.10881e-06** +- **Test: x_screw=60 Detector: w_monitor_I=0.0721** - **TestScan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={21 values}** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md index 5eed2352af..7943269940 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/README.md @@ -18,9 +18,9 @@ You may DISCO_scan the grating monochromator with e.g.: ## Examples -- **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={21 values}** +- **TestScan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 -n1e6 Detector: Wavelength_I={21 values}** - -- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=2.08292e+11** +- **Test: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.445e+16** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr index 2174c71802..43a787d493 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr @@ -13,12 +13,12 @@ * %Description * DISCO beamline description here: https://hal.archives-ouvertes.fr/hal-01479318/document * You may DISCO_scan the grating monochromator with e.g.: -* %Scan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 Detector: Wavelength_I={ +* %Scan: mxrun SOLEIL_DISCO.instr -N21 DISCO_x_screw=20,103 DISCO_scan=1 -n1e6 Detector: Wavelength_I={ * 2019.39,2365.97,2837.41,3242.59,3655.92,4070.23,4487.58,4921.15, * 5306.77,5738.39,6151.91,6565.99,6979.47,7395.13,7786.17,8199.52, * 8631.53,9051.46,9483.21,9913.98,10299.5 } * -* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=2.08292e+11 +* %Example: DISCO_grazing_angle_one=22.5 -n 1e5 Detector: DISCO_second_HFP_I=1.445e+16 * * %Parameters * DISCO_grazing_angle_one: [deg] Mirrors M1/M2 rotation angle diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md index a7eb00befa..aa2adae939 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/README.md @@ -28,7 +28,7 @@ This model implements the XRD station for powders. ## Examples -- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=2.56e+09** +- **Test: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.341e+08** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr index c422ef5893..8984370585 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_MARS/SOLEIL_MARS.instr @@ -24,7 +24,7 @@ * * This model implements the XRD station for powders. * -* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=2.56e+09 +* %Example: MARS_E0=16.99 Detector: MARS_detector_diffraction_I=1.341e+08 * * %Parameters * MARS_E0: [keV] Central energy to be emitted by the source diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md index 276519732e..04c4c8f2ea 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/README.md @@ -48,7 +48,7 @@ Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000 ## Examples -- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=2e+09** +- **Test: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.61e+11** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr index 7faf684e6c..7cd24dd0c8 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_ROCK/SOLEIL_ROCK.instr @@ -44,7 +44,7 @@ * Manganese and chrome ROCK_scan: mxrun SOLEIL_ROCK.instr ROCK_E0=5.700,6.800 ROCK_scan=1 ROCK_sample_file=MnCr ROCK_cc=2 -N100 * Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr ROCK_E0=17.000 ROCK_cc=2 -n1e8 * -* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=2e+09 +* %Example: ROCK_E0=15.918 Detector: ROCK_fluo_monitor_I=1.61e+11 * * %Parameters * ROCK_E0: [keV] Energy to hit, i.e. selected by the channel-cut monochromator diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md index 0e88d31f6f..cb3d87ca31 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/README.md @@ -34,7 +34,7 @@ Position | Element ## Examples -- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=1.23968e+10** +- **Test: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.69e+10** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr index a179c44abb..f3e96f61a1 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SAMBA/SOLEIL_SAMBA.instr @@ -34,7 +34,7 @@ * 30 | Sample Stage (XAS configuration) * 35 | Sample Stage (SEXAFS configuration) * -* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=1.23968e+10 +* %Example: SAMBA_E0=15 Detector: SAMBA_det_fluo_I=5.69e+10 * * %Parameters * SAMBA_E0: [keV] Central energy diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md index 008b20f8b3..4ebd4abbea 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/README.md @@ -30,7 +30,7 @@ Position | Element ## Examples -- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=2.14522e+10** +- **Test: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.69e+09** ## Input parameters diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr index fcbb84bbd0..2ad3ecd670 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_SIRIUS/SOLEIL_SIRIUS.instr @@ -28,7 +28,7 @@ * 30 | Sample Stage (XAS configuration) * 35 | Sample Stage (SEXAFS configuration) * -* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=2.14522e+10 +* %Example: SIRIUS_E0=8 Detector: SIRIUS_def_loc_I=6.69e+09 * * %Parameters * SIRIUS_E0: [keV] Central energy (1.4 - 13 keV) diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md index 29095ae97c..143cc28999 100644 --- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md +++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/README.md @@ -17,7 +17,7 @@ Simply a bending magnet illuminating a reflection grating. ## Examples -- **Test: E0=1 Detector: psd_monitor_I=3.79378e+06** +- **Test: E0=1 -n 1e7 Detector: psd_monitor_I=3.53e+09** ## Input parameters diff --git a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr index 685f61520e..44005c09e7 100644 --- a/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr +++ b/mcxtrace-comps/examples/Tests_optics/Test_grating_reflect/Test_grating_reflect.instr @@ -15,7 +15,7 @@ * * Simply a bending magnet illuminating a reflection grating. * -* %Example: E0=1 Detector: psd_monitor_I=3.79378e+06 +* %Example: E0=1 -n 1e7 Detector: psd_monitor_I=3.53e+09 * * %Parameters * E0: [keV] Source's center of emitted energy spectrum diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md index 8a662d9e62..28987cad2b 100644 --- a/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md +++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/README.md @@ -17,7 +17,7 @@ This is a simple test-instrument for the bending magnet component. ## Examples -- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=2.93138e+13** +- **Test: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.212e+14** - **Test: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15** ## Input parameters diff --git a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr index 8133248629..29bd12c2c4 100644 --- a/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr +++ b/mcxtrace-comps/examples/Tests_sources/Test_BM/Test_BM.instr @@ -16,7 +16,7 @@ * %Description * This is a simple test-instrument for the bending magnet component. * -* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=2.93138e+13 +* %Example: Test_BM.instr SOURCE=0 e0=12.5 de=0.5 Detector: emon_I=9.212e+14 * %Example: Test_BM.instr SOURCE=1 e0=12.5 de=0.5 Detector: emon_I=1.28549e+15 * * %Parameters diff --git a/mcxtrace-comps/sources/Bending_magnet.comp b/mcxtrace-comps/sources/Bending_magnet.comp index 8be92132bc..2b5e1d58c6 100644 --- a/mcxtrace-comps/sources/Bending_magnet.comp +++ b/mcxtrace-comps/sources/Bending_magnet.comp @@ -13,11 +13,20 @@ * Date: May, 2013. * Version: 1.0 * Origin: DTU Physics +* Modified by: Antoine Padovani, 25th Oct 2026 : fix flux and sampling * * Model of a bending magnet source * * %Description -* A source model based on the derivation from B.D. Patterson, Am. J. Phys. 79, 1046 (2011); doi: 10.1119/1.3614033 +* A source model using the full angular distribution of synchrotron radiation from a bending magnet, +* eq. (1) of K.-J. Kim, X-Ray Data Booklet, Section 2.1 (https://xdb.lbl.gov/Section2/Sec_2-1.html). +* In the horizontal plane (psi=0) this reduces to eq. (28) of B.D. Patterson, Am. J. Phys. 79, 1046 (2011); doi: 10.1119/1.3614033 +* +* The ray weight p is the flux in [photons/s] emitted into the sampled solid angle and the energy +* band [E0-dE, E0+dE] (i.e. NOT per mrad^2 or per 0.1% bandwidth). The energy is sampled uniformly, so dE>0 +* is required. In lambda0/dlambda mode the corresponding energy interval is sampled uniformly in energy. +* With focus_xw/focus_yh the target window is sampled; without them, the horizontal/vertical emission angle +* is sampled uniformly within +-2.5 mrad, so the total flux is then the flux into a 5 mrad horizontal fan. * * Example: Bending_magnet( * E0 = 14, dE = 7, Ee = 2.75, @@ -46,44 +55,44 @@ DEFINE COMPONENT Bending_magnet -SETTING PARAMETERS (E0=0, dE=0, lambda0=0,dlambda=0, phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0) +SETTING PARAMETERS (E0=0, dE=0, lambda0=0,dlambda=0, + phase=0, randomphase=1, Ee=2.4, Ie=0.4, B=1.6, sigey=0, sigex=0, focus_xw=0, focus_yh=0, dist=1, gauss_t=0) -/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ +/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ SHARE %{ - // %include "read_table-lib" - - #ifndef MCCODE_BESSELKNU - #define MCCODE_BESSELKNU 1 - - #pragma acc routine seq - double - besselKnu (double nu, double x) { - const double h = 0.5; - double KK = 0, dK; - int r = 0; - const int maxiter = 1000; - KK = exp (-x) / 2.0; - dK = 1; - while (dK > DBL_EPSILON && r < maxiter) { + +#ifndef MCCODE_BESSELKNU +#define MCCODE_BESSELKNU 1 + +#pragma acc routine seq +double besselKnu(double nu, double x){ + const double h=0.5; + double KK=0,dK; + int r=0; + const int maxiter=1000; + KK=exp(-x)/2.0; + dK=1; + while (dK>DBL_EPSILON && r= maxiter) { - fprintf (stderr, "Warning: Maximum number of iterations exceeded in besselKnu(%g,%g).\n", nu, x); +#ifndef OPENACC + if (r>=maxiter) { + fprintf(stderr,"Warning: Maximum number of iterations exceeded in besselKnu(%g,%g).\n",nu,x); } - #endif - KK *= h; +#endif + KK*=h; return KK; } - #endif /*MCCODE_BESSELKNU*/ +#endif /*MCCODE_BESSELKNU*/ + +#ifndef M_SQRT1_2 +#define M_SQRT1_2 0.70710678118654752440 +#endif - #ifndef M_SQRT1_2 - #define M_SQRT1_2 0.70710678118654752440 - #endif %} DECLARE @@ -93,175 +102,208 @@ DECLARE double igamma; double kc; /*characteristic wavenumber of radiation from bending magnet*/ double s1x; - double s1y; /*beam's size at dist (convolution of sigex/sigey and igamma)*/ + double s1y; /*sigma of the gaussian sampling of the target window (gauss_t)*/ double p0; + double ctemult; + double muc; + double qe; %} INITIALIZE %{ - // fprintf(stderr,"Warning (%s): Bending_magnet is an experimental component - testing is ongoing\n",NAME_CURRENT_COMP); - - if (B <= 0 || Ee <= 0 || Ie <= 0) { - fprintf (stderr, "Error (%s): B, Ee, and Ie must all be >= 0. Found (%g %g %g). Aborting.\n", NAME_CURRENT_COMP, B, Ee, Ie); - exit (1); + if(B<=0 || Ee<=0 || Ie<=0 ){ + fprintf(stderr, "Error (%s): B, Ee, and Ie must all be > 0. Found (%g %g %g). Aborting.\n",NAME_CURRENT_COMP,B,Ee,Ie); + exit(1); } - - if (sigex < 0 || sigey < 0) { - fprintf (stderr, "Error (%s): sigex and sigey must be > 0. Negative beam size isn't meaningful. Aborting.\n", NAME_CURRENT_COMP); - exit (1); + if (!E0 && lambda0 > 0 && dlambda > 0 && lambda0-dlambda > 0) { + // compute equivalent E0 and dE + double E1 = (2 * M_PI / (lambda0+dlambda)) * K2E; + double E2 = (2 * M_PI / (lambda0-dlambda)) * K2E; + E0 = (E1 + E2)/2; + dE = fabs(E1 - E2)/2; } - if (dist <= 0) { - fprintf (stderr, "Error (%s): Target undefined.\n", NAME_CURRENT_COMP); + if (E0 <= 0 || dE <= 0 || E0-dE <= 0) { + fprintf (stderr, "Error (%s): E0 and dE>0 (or lambda0 and dlambda>0) need to be set, with E0-dE>0. Aborting.\n", NAME_CURRENT_COMP); exit (1); } + if (sigex <0 || sigey<0){ + fprintf(stderr, "Error (%s): sigex and sigey must be > 0. Negative beam size isn't meaningful. Aborting.\n",NAME_CURRENT_COMP); + exit(1); + } + if (dist<=0){ + fprintf(stderr,"Error (%s): Target undefined.\n",NAME_CURRENT_COMP); + exit(1); + } /*compute gamma*/ - gamma = (Ee * 1e9) / (MELECTRON / CELE * M_C * M_C); /*the extra CELE is to convert to eV*/ - gamma2 = gamma * gamma; - igamma = 1.0 / gamma; + gamma=(Ee*1e9)/(MELECTRON/CELE*M_C*M_C);/*the extra CELE is to convert to eV*/ + gamma2=gamma*gamma; + igamma=1.0/gamma; + + ctemult=(3.0*ALPHA/(4.0*M_PI*M_PI))*(gamma2)*(Ie/CELE); - // printf("Bending_magnet (%s): gamma=%g, divergence is 1/gamma=%g rad.\n",NAME_CURRENT_COMP,gamma,igamma); + //printf("Bending_magnet (%s): gamma=%g, divergence is 1/gamma=%g rad.\n",NAME_CURRENT_COMP,gamma,igamma); /*compute characteristic energy in keV*/ - double Ec = 0.665 * Ee * Ee * B; - // double Ec=1.5*gamma2*HBAR*CELE*B/MELECTRON *1e-3; /*check units on this one. The 1e-3 factor is because energy is assumed to be in keV*/ + double Ec=0.665*Ee*Ee*B; + //double Ec=1.5*gamma2*HBAR*CELE*B/MELECTRON *1e-3; /*check units on this one. The 1e-3 factor is because energy is assumed to be in keV*/ /*We normally do computations in k so use that for transfer*/ - kc = E2K * Ec; + kc=E2K*Ec; + + muc = kc*M_C/(2.0*M_PI); - s1x = sqrt (sigex * sigex + igamma * igamma * dist * dist); - s1y = sqrt (sigey * sigey + igamma * igamma * dist * dist); - p0 = 1.0 / mcget_ncount (); + if (focus_xw){ + s1x=focus_xw/2.0; + } + if (focus_yh){ + s1y=focus_yh/2.0; + } + p0=1.0/mcget_ncount(); + + qe = 1.0/(2.0*dE); // energy sampling density, dE > 0 checked above %} TRACE %{ - double xx, yy, x1, y1, z1; - double k, e, l; - double F1 = 1.0; - double dx, dy, dz; + double xx,yy,x1,y1,z1; + double k,e; + double dx,dy,dz; + double dl; + + const double angle=2.5e-3; /*without focus, psi and theta are sampled uniformly in (-angle,angle)*/ + double psi, theta, xcap, oneplusxsquared, chi, K23, K13; // initial source area - xx = randnorm (); - yy = randnorm (); - x = xx * sigex; - y = yy * sigey; - z = 0; + xx=randnorm(); + yy=randnorm(); + x=xx*sigex; + y=yy*sigey; + z=0; // Gaussian distribution at origin - p = p0; /*initial weight is p0*/ - if (E0) { - if (!dE) { - e = E0; - } else { - e = randpm1 () * dE + E0; - } - k = E2K * e; - } else if (lambda0) { - if (!dlambda) { - l = lambda0; - } else { - l = randpm1 () * dlambda + lambda0; - } - k = (2 * M_PI / l); - } + p=p0;/*initial weight is p0*/ + // E0 and dE are defined + e=randpm1()*dE + E0; + k=E2K*e; + + double mu = k*M_C/(2*M_PI); + double mumuc = mu/muc; + // targeted area calculation - if (focus_xw) { - if (!gauss_t) { + if (focus_xw){ + if (!gauss_t){ /*sample uniformly but adjust weight*/ - x1 = randpm1 () * focus_xw / 2.0; - p *= exp (-(x1 * x1) / (2.0 * s1x * s1x)); - } else { + x1=randpm1()*focus_xw/2.0; + theta = atan((x1-x)/dist); + p*=focus_xw; //sampling probability density is 1.0/focus_xw + p*= dist/(dist*dist + (x1-x)*(x1-x)); //take into account the fact that we are sampling a length, not an angle (and our flux is angular). So need to apply that correction also. + }else { do { - x1 = randnorm () * s1x; - } while (focus_xw != 0 && fabs (x1) > focus_xw / 2.0); - /*adjust for restricted sampling window*/ - p *= erf (focus_xw * 0.5 * M_SQRT1_2 / s1x); + x1=randnorm()*s1x; + }while (fabs(x1)>focus_xw/2.0); + theta = atan((x1-x)/dist); + /*adjust for restricted sampling window*/ + p*=erf(focus_xw*0.5*M_SQRT1_2/s1x)*s1x*sqrt(2.0*M_PI)*exp((x1*x1)/(2*s1x*s1x)); // 1/sampling probability density + p*= dist/(dist*dist + (x1-x)*(x1-x)); } - } else { - x1 = randnorm () * igamma; + }else{ + theta = randpm1()*angle; + p*=2*angle; // 1/sampling probability density + x1=x+dist*tan(theta); } - if (focus_yh) { - if (!gauss_t) { + if (focus_yh){ + if (!gauss_t){ /*sample uniformly but adjust weight*/ - y1 = randpm1 () * focus_yh / 2.0; - p *= exp (-(y1 * y1) / (2.0 * s1y * s1y)); - } else { + y1=randpm1()*focus_yh/2.0; + psi = atan((y1-y)/dist); + p*=focus_yh; //sampling probability density is 1.0/focus_yh + p*= dist/(dist*dist + (y1-y)*(y1-y)); + }else { do { - y1 = randnorm () * s1y; - } while (fabs (y1) > focus_yh / 2.0); - /*adjust for restricted sampling window*/ - p *= erf (focus_yh * 0.5 * M_SQRT1_2 / s1y); + y1=randnorm()*s1y; + }while (fabs(y1)>focus_yh/2.0); + psi = atan((y1-y)/dist); + /*adjust for restricted sampling window*/ + p*=erf(focus_yh*0.5*M_SQRT1_2/s1y)*s1y*sqrt(2.0*M_PI)*exp((y1*y1)/(2*s1y*s1y)); // 1/sampling probability density + p*= dist/(dist*dist + (y1-y)*(y1-y)); } - } else { - y1 = randnorm () * igamma; + }else{ + psi = randpm1()*angle; + p*=2*angle; // 1/sampling probability density + y1=y+dist*tan(psi); } - z1 = dist; - dx = x1 - x; - dy = y1 - y; - dz = sqrt (dx * dx + dy * dy + dist * dist); + + xcap = gamma*psi; + oneplusxsquared = 1+xcap*xcap; - kx = (k * dx) / dz; - ky = (k * dy) / dz; - kz = (k * dist) / dz; + chi = mumuc*pow(oneplusxsquared,1.5)*(1.0/2.0); + + K23 = besselKnu(0.666666666666666666666666667,chi); + K13 = besselKnu(0.333333333333333333333333334,chi); + + z1=z+dist; + dx=x1-x; + dy=y1-y; + dz=z1-z; + + dl = sqrt(dx*dx+dy*dy+dz*dz); + + kx=(k*dx)/dl; + ky=(k*dy)/dl; + kz=(k*dz)/dl; + + p*=ctemult*mumuc*mumuc*(oneplusxsquared*oneplusxsquared)*(K23*K23+(xcap*xcap/(oneplusxsquared))*K13*K13); - /*spectral strength of radiation is given by Patterson*/ - double k_kc = k / kc; - double K2_3 = besselKnu (0.666666666666666666666666667, k_kc * 0.5); - p *= 1.33e13 * Ee * Ee * Ie * k_kc * k_kc * K2_3 * K2_3; - // p*=ALPHA/(M_PI*M_PI)*gamma2*Ie/CELE* 1e-4 * 0.75 *k_kc*k_kc*K2_3*K2_3; + p*= 1.0/(qe*e); // energy sampling density 1.0/qe and 1/e coming from the physical flux /*randomly pick phase*/ - if (randomphase) { - phi = rand01 () * 2 * M_PI; - } else { - phi = phase; + if (randomphase){ + phi=rand01()*2*M_PI; + }else{ + phi=phase; } /*set polarization vector*/ - Ex = 0; - Ey = 0; - Ez = 0; + Ex=Ey=Ez=0; %} MCDISPLAY %{ - - double radius, D; - radius = 3.3 * Ee / B; - ; - D = 0.1; - double x0, x1, z0, z1; - const double phimin = -2.0 * DEG2RAD, phimax = 2.0 * DEG2RAD; - double phi = phimin, dphi; - dphi = (phimax - phimin) / 32; - while (phi < phimax) { - x0 = radius * (1 - cos (phi)); - x1 = radius * (1 - cos (phi + dphi)); - z0 = radius * sin (phi); - z1 = radius * sin (phi + dphi); - line (x0, 0.0, z0, x1, 0.0, z1); - phi += dphi; + double radius,D; + radius=3.3*Ee/B; + D=0.1; + double x0,x1,z0,z1; + const double phimin=-2.0*DEG2RAD, phimax=2.0*DEG2RAD; + double phi=phimin,dphi; + dphi=(phimax-phimin)/32; + while (phi