From a19c880b6e188e3d6d70e55ece37e040f040aad7 Mon Sep 17 00:00:00 2001 From: Peter Willendrup Date: Fri, 20 Feb 2026 09:19:41 +0100 Subject: [PATCH] Format McXtrace contrib --- mcxtrace-comps/contrib/Attenuating_mask.comp | 84 +- mcxtrace-comps/contrib/Bragg_crystal_BC.comp | 279 +++--- .../contrib/Bragg_crystal_bent_BC.comp | 504 +++++----- .../contrib/Bragg_crystal_simple.comp | 30 +- mcxtrace-comps/contrib/Detector_pn.comp | 133 +-- mcxtrace-comps/contrib/Laue_crystal_BC.comp | 383 ++++---- .../contrib/Mirror_toroid_pothole.comp | 359 +++---- mcxtrace-comps/contrib/PSD_monitor_rad.comp | 72 +- mcxtrace-comps/contrib/SAXSCurve.comp | 76 +- mcxtrace-comps/contrib/SAXSCylinders.comp | 184 ++-- .../contrib/SAXSEllipticCylinders.comp | 208 ++-- mcxtrace-comps/contrib/SAXSLiposomes.comp | 300 +++--- mcxtrace-comps/contrib/SAXSNanodiscs.comp | 310 +++--- mcxtrace-comps/contrib/SAXSNanodiscsFast.comp | 375 ++++---- .../contrib/SAXSNanodiscsWithTags.comp | 353 +++---- .../contrib/SAXSNanodiscsWithTagsFast.comp | 421 +++++---- mcxtrace-comps/contrib/SAXSPDB.comp | 276 +++--- mcxtrace-comps/contrib/SAXSPDBFast.comp | 887 +++++++++--------- mcxtrace-comps/contrib/SAXSQMonitor.comp | 263 +++--- mcxtrace-comps/contrib/SAXSShells.comp | 184 ++-- mcxtrace-comps/contrib/SAXSSpheres.comp | 152 +-- 21 files changed, 2925 insertions(+), 2908 deletions(-) diff --git a/mcxtrace-comps/contrib/Attenuating_mask.comp b/mcxtrace-comps/contrib/Attenuating_mask.comp index 3543b0ed5d..10d685230d 100644 --- a/mcxtrace-comps/contrib/Attenuating_mask.comp +++ b/mcxtrace-comps/contrib/Attenuating_mask.comp @@ -52,58 +52,62 @@ SHARE DECLARE %{ - double rho; - t_Table att_table; + double rho; + t_Table att_table; %} INITIALIZE %{ - int status=0; - if((status=Table_Read(&(att_table),att_file,0))==-1){ // load the attenuation crossection - fprintf(stderr,"Error: Could not parse file \"%s\" in COMP %s\n",att_file,NAME_CURRENT_COMP); exit(-1); - } - /* checking the header for retrieving density */ - char **header_parsed; - header_parsed = Table_ParseHeader(att_table.header,"rho",NULL); - if(header_parsed[0]){rho = strtod(header_parsed[0],NULL);} - else{fprintf(stderr,"Warning(%s): %s not found in header of %s, set to 1\n",NAME_CURRENT_COMP,"rho",att_file); rho = 1;} - printf("Rho = %1.3f.\n\n",rho); + int status = 0; + if ((status = Table_Read (&(att_table), att_file, 0)) == -1) { // load the attenuation crossection + fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", att_file, NAME_CURRENT_COMP); + exit (-1); + } + /* checking the header for retrieving density */ + char** header_parsed; + header_parsed = Table_ParseHeader (att_table.header, "rho", NULL); + if (header_parsed[0]) { + rho = strtod (header_parsed[0], NULL); + } else { + fprintf (stderr, "Warning(%s): %s not found in header of %s, set to 1\n", NAME_CURRENT_COMP, "rho", att_file); + rho = 1; + } + printf ("Rho = %1.3f.\n\n", rho); %} TRACE %{ - int hit_sample = 0; - double d0, d1; - - hit_sample = box_intersect(&d0, &d1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); // check intersection with the component - - if(hit_sample){ - PROP_DL(d0); - double cosleftx = cos(2*PI*x/blocks_xdist); - double cosrightx = cos(PI*blocks_xwidth/blocks_xdist); - double coslefty = cos(2*PI*y/blocks_ydist); - double cosrighty = cos(PI*blocks_yheight/blocks_ydist); - double k = sqrt(kx*kx+ky*ky+kz*kz); - double E = K2E*k; - if(holed_mask==0){ // mask composed of blocking disks - if((cosleftx>cosrightx || cosleftx==cosrightx) && (coslefty>cosrighty || coslefty==cosrighty)){ //absorbed - double mul = 1e2*(d1-d0)*rho*Table_Value(att_table,E,5); - p*=exp(-mul); - } //else transmitted - } - else{ // mask composed of apertures on a grating - if(cosleftx cosrightx || cosleftx == cosrightx) && (coslefty > cosrighty || coslefty == cosrighty)) { // absorbed + double mul = 1e2 * (d1 - d0) * rho * Table_Value (att_table, E, 5); + p *= exp (-mul); + } // else transmitted + } else { // mask composed of apertures on a grating + if (cosleftx < cosrightx || coslefty < cosrighty) { // absorbed + double mul = 1e2 * (d1 - d0) * rho * Table_Value (att_table, E, 5); + p *= exp (-mul); + } // else transmitted + } + PROP_DL (d1 - d0); + } %} MCDISPLAY %{ - box(0,0,0,xwidth,yheight,zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/Bragg_crystal_BC.comp b/mcxtrace-comps/contrib/Bragg_crystal_BC.comp index 3e0b632a09..942e0f28d6 100644 --- a/mcxtrace-comps/contrib/Bragg_crystal_BC.comp +++ b/mcxtrace-comps/contrib/Bragg_crystal_BC.comp @@ -92,160 +92,163 @@ DECLARE double f_rel; double f_nt; t_Table m_t; - t_Table f0_t; + t_Table f0_t; %} INITIALIZE %{ - int status; - if (material){ - if ((status=Table_Read(&(m_t),material,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,material); - exit(-1); - } - char **header_parsed; - header_parsed=Table_ParseHeader(m_t.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL); - if(header_parsed[2]){rho=strtod(header_parsed[2],NULL);} - if(header_parsed[0]){Z=strtod(header_parsed[0],NULL);} - if(header_parsed[1]){At=strtod(header_parsed[1],NULL);} - }else{ - fprintf(stderr,"Error(%s): No material file specified\n",NAME_CURRENT_COMP); + int status; + if (material) { + if ((status = Table_Read (&(m_t), material, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, material); + exit (-1); } - if(form_factors){ - if ((status=Table_Read(&(f0_t),form_factors,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,form_factors); - exit(-1); - } + char** header_parsed; + header_parsed = Table_ParseHeader (m_t.header, "Z", "A[r]", "rho", "Z/A", "sigma[a]", NULL); + if (header_parsed[2]) { + rho = strtod (header_parsed[2], NULL); } - - if(verbose) - printf("INFO (%s): initialized\n",NAME_CURRENT_COMP); + if (header_parsed[0]) { + Z = strtod (header_parsed[0], NULL); + } + if (header_parsed[1]) { + At = strtod (header_parsed[1], NULL); + } + } else { + fprintf (stderr, "Error(%s): No material file specified\n", NAME_CURRENT_COMP); + } + if (form_factors) { + if ((status = Table_Read (&(f0_t), form_factors, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, form_factors); + exit (-1); + } + } + + if (verbose) + printf ("INFO (%s): initialized\n", NAME_CURRENT_COMP); %} TRACE %{ - double E; // (keV) x-ray energy - double K; // length of k-vector - double tin; // 'time' of intersection of ray with y=0 plane (which include the crystal surface) - double x_int,y_int,z_int; // intersection with the y=0 plane - double dist; // distance from position at t=0 to the y=0 plane - double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). - double R; // Reflectivity values calculated by Mx_DarwinReflectivityBC() function for each incoming photon - - /* get the photon's kvector and energy */ - K=sqrt(kx*kx+ky*ky+kz*kz); - E = K2E*K; /* use built-in constants for consistency */ - /* make unit vector in the direction of k :*/ - double k0hat[3] = {kx/K, ky/K, kz/K}; // unit vector in the direction of k-vector. - double kk0[3]={kx,ky,kz}; - - /*intersection calculation*/ - tin = -y/k0hat[1]; - if (tin > 1e-8){ - /* check whether our intersection lies within the boundaries of the crystal*/ - x_int=x+k0hat[0]*tin; - y_int=y+k0hat[1]*tin; - z_int=z+k0hat[2]*tin; - - if (fabs(x_int)<=width/2 && fabs(z_int)<=length/2){ - dist=sqrt(SQR(x-x_int)+SQR(y-y_int)+SQR(z-z_int)); - PROP_DL(dist); /* now the photon is on the crystal surface, ready to be reflected... */ - SCATTER; - double d=cbrt(V)/(sqrt(h*h+k*k+l*l));/*this is valid only for cubic structures*/ - f00 = Z; - f0h = Table_Value(f0_t,1/(2*d),Z); - fp = Table_Value(m_t,E,1)-Z; - fpp = Table_Value(m_t,E,2); - - double nhat[3]={0,1,0}; // crystal surface normal is yhat, and is set up by mcxtrace to be into the crystal - - double alpha[3]={alphax, alphay, alphaz}; - double Rsig, Rpi; - double kh[3], sig_axis[3], pi_axis[3]; - int fail; - - double complex chi0, chih; - double k0mag, hscale, thetaB; - - Mx_CubicCrystalChi(&chi0, &chih, &k0mag, &hscale, &thetaB, - f00, f0h, fp, fpp, V, h, k, l, - debye_waller_B, E, - crystal_type,structure_factor_scale_r,structure_factor_scale_i); - if(thetaB==0){ - if(verbose) - fprintf(stderr,"WARNING (%s): reflection [%d %d %d] is inaccessible for E= %g keV. Terminating photon,\n",NAME_CURRENT_COMP,h,k,l,E); - ABSORB; - } - fail=Mx_DarwinReflectivityBC(&Rsig, &Rpi, kh, // these are the return values - k0hat, nhat, alpha, - chi0, chih, chih, k0mag, hscale, thetaB - ); - - cross(sig_axis, k0hat, kh, 1); /* kin x kout is sigma direction */ - /* sig_axis is the sigma direction, as returned by Bragg_Geometry inside DarwinReflectivity_BC */ - /* pi is a vector perpendicular to k_in and sig i.e. the direction of pi polarization incoming */ - cross(pi_axis, k0hat, sig_axis, 1); - - /* update outgoing direction vector kx, ky, kz = kh, fully scaled outgoing k vector */ - kx=kh[0]; ky=kh[1]; kz=kh[2]; - - /* resolve incoming polarization into sig and pi bits, and scale by sqrt(reflectivity) which is amplitude scale */ - double Esig=(Ex*sig_axis[0]+Ey*sig_axis[1]+Ez*sig_axis[2]); - double Epi= (Ex* pi_axis[0]+Ey* pi_axis[1]+Ez* pi_axis[2]); - if(Esig==0 && Epi==0) { /* someone didn't set the polarization direction; set it now to a random value and it will propagate */ - double psi=rand01()*PI/2; - Esig=cos(psi); Epi=sin(psi); - } - Esig=Esig*sqrt(Rsig); - Epi=Epi*sqrt(Rpi); - R=Esig*Esig+Epi*Epi; /* projected reflectivity, squared back to intensity */ - - double pi1_axis[3]; - /* pi1 is now a vector perpendicular to k_out and sig i.e. the direction of pi polarization outgoing */ - cross(pi1_axis, kh, sig_axis,1); - - /* a linear combination of these is still perpendicular to k, but has the correct polarization weighting */ - Ex=Epi*pi1_axis[0]+Esig*sig_axis[0]; - Ey=Epi*pi1_axis[1]+Esig*sig_axis[1]; - Ez=Epi*pi1_axis[2]+Esig*sig_axis[2]; - NORM(Ex, Ey, Ez); - -#ifdef MCDEBUG - fprintf(stderr,"Crystal: %s: Rsig=%.4g Rpi=%.4g " - " k0=(%.8f, %.8f, %.8f), nhat=(%.8f, %.8f, %.8f) alpha=(%.8f, %.8f, %.8f), " - " kout=(%.8f, %.8f, %.8f) " - " axis=(%.8f, %.8f, %.8f) E=(%.4f, %.4f, %.4f) \n", - NAME_CURRENT_COMP, Rsig, Rpi, - kk0[0], kk0[1], kk0[2], nhat[0], nhat[1], nhat[2], - alpha[0], alpha[1], alpha[2], - kh[0], kh[1], kh[2], - sig_axis[0], sig_axis[1], sig_axis[2], Ex, Ey, Ez); -#endif - /* apply Darwin reflectivity if not is supplied from outside*/ - if (!R0){ - p*=R; - }else{ - p*=R0; - } - /*catch dead rays*/ - if (p==0) ABSORB; - } else { - RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); - } + double E; // (keV) x-ray energy + double K; // length of k-vector + double tin; // 'time' of intersection of ray with y=0 plane (which include the crystal surface) + double x_int, y_int, z_int; // intersection with the y=0 plane + double dist; // distance from position at t=0 to the y=0 plane + double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). + double R; // Reflectivity values calculated by Mx_DarwinReflectivityBC() function for each incoming photon + + /* get the photon's kvector and energy */ + K = sqrt (kx * kx + ky * ky + kz * kz); + E = K2E * K; /* use built-in constants for consistency */ + /* make unit vector in the direction of k :*/ + double k0hat[3] = { kx / K, ky / K, kz / K }; // unit vector in the direction of k-vector. + double kk0[3] = { kx, ky, kz }; + + /*intersection calculation*/ + tin = -y / k0hat[1]; + if (tin > 1e-8) { + /* check whether our intersection lies within the boundaries of the crystal*/ + x_int = x + k0hat[0] * tin; + y_int = y + k0hat[1] * tin; + z_int = z + k0hat[2] * tin; + + if (fabs (x_int) <= width / 2 && fabs (z_int) <= length / 2) { + dist = sqrt (SQR (x - x_int) + SQR (y - y_int) + SQR (z - z_int)); + PROP_DL (dist); /* now the photon is on the crystal surface, ready to be reflected... */ + SCATTER; + double d = cbrt (V) / (sqrt (h * h + k * k + l * l)); /*this is valid only for cubic structures*/ + f00 = Z; + f0h = Table_Value (f0_t, 1 / (2 * d), Z); + fp = Table_Value (m_t, E, 1) - Z; + fpp = Table_Value (m_t, E, 2); + + double nhat[3] = { 0, 1, 0 }; // crystal surface normal is yhat, and is set up by mcxtrace to be into the crystal + + double alpha[3] = { alphax, alphay, alphaz }; + double Rsig, Rpi; + double kh[3], sig_axis[3], pi_axis[3]; + int fail; + + double complex chi0, chih; + double k0mag, hscale, thetaB; + + Mx_CubicCrystalChi (&chi0, &chih, &k0mag, &hscale, &thetaB, f00, f0h, fp, fpp, V, h, k, l, debye_waller_B, E, crystal_type, structure_factor_scale_r, + structure_factor_scale_i); + if (thetaB == 0) { + if (verbose) + fprintf (stderr, "WARNING (%s): reflection [%d %d %d] is inaccessible for E= %g keV. Terminating photon,\n", NAME_CURRENT_COMP, h, k, l, E); + ABSORB; + } + fail = Mx_DarwinReflectivityBC (&Rsig, &Rpi, kh, // these are the return values + k0hat, nhat, alpha, chi0, chih, chih, k0mag, hscale, thetaB); + + cross (sig_axis, k0hat, kh, 1); /* kin x kout is sigma direction */ + /* sig_axis is the sigma direction, as returned by Bragg_Geometry inside DarwinReflectivity_BC */ + /* pi is a vector perpendicular to k_in and sig i.e. the direction of pi polarization incoming */ + cross (pi_axis, k0hat, sig_axis, 1); + + /* update outgoing direction vector kx, ky, kz = kh, fully scaled outgoing k vector */ + kx = kh[0]; + ky = kh[1]; + kz = kh[2]; + + /* resolve incoming polarization into sig and pi bits, and scale by sqrt(reflectivity) which is amplitude scale */ + double Esig = (Ex * sig_axis[0] + Ey * sig_axis[1] + Ez * sig_axis[2]); + double Epi = (Ex * pi_axis[0] + Ey * pi_axis[1] + Ez * pi_axis[2]); + if (Esig == 0 && Epi == 0) { /* someone didn't set the polarization direction; set it now to a random value and it will propagate */ + double psi = rand01 () * PI / 2; + Esig = cos (psi); + Epi = sin (psi); + } + Esig = Esig * sqrt (Rsig); + Epi = Epi * sqrt (Rpi); + R = Esig * Esig + Epi * Epi; /* projected reflectivity, squared back to intensity */ + + double pi1_axis[3]; + /* pi1 is now a vector perpendicular to k_out and sig i.e. the direction of pi polarization outgoing */ + cross (pi1_axis, kh, sig_axis, 1); + + /* a linear combination of these is still perpendicular to k, but has the correct polarization weighting */ + Ex = Epi * pi1_axis[0] + Esig * sig_axis[0]; + Ey = Epi * pi1_axis[1] + Esig * sig_axis[1]; + Ez = Epi * pi1_axis[2] + Esig * sig_axis[2]; + NORM (Ex, Ey, Ez); + + #ifdef MCDEBUG + fprintf (stderr, + "Crystal: %s: Rsig=%.4g Rpi=%.4g " + " k0=(%.8f, %.8f, %.8f), nhat=(%.8f, %.8f, %.8f) alpha=(%.8f, %.8f, %.8f), " + " kout=(%.8f, %.8f, %.8f) " + " axis=(%.8f, %.8f, %.8f) E=(%.4f, %.4f, %.4f) \n", + NAME_CURRENT_COMP, Rsig, Rpi, kk0[0], kk0[1], kk0[2], nhat[0], nhat[1], nhat[2], alpha[0], alpha[1], alpha[2], kh[0], kh[1], kh[2], sig_axis[0], + sig_axis[1], sig_axis[2], Ex, Ey, Ez); + #endif + /* apply Darwin reflectivity if not is supplied from outside*/ + if (!R0) { + p *= R; + } else { + p *= R0; + } + /*catch dead rays*/ + if (p == 0) + ABSORB; + } else { + RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); } - -#ifdef MCDEBUG - fflush(stderr); - fflush(stdout); // make sure all error messages are in order! -#endif + } + #ifdef MCDEBUG + fflush (stderr); + fflush (stdout); // make sure all error messages are in order! + #endif %} MCDISPLAY %{ - magnify(""); - rectangle("xz",0,0,0,width,length); // display crystal outline - line(0,0,0,alphax*(width+length)/2.0, alphay*(width+length)/2.0, alphaz*(width+length)/2.0); // display alpha vector + magnify (""); + rectangle ("xz", 0, 0, 0, width, length); // display crystal outline + line (0, 0, 0, alphax*(width + length) / 2.0, alphay*(width + length) / 2.0, alphaz*(width + length) / 2.0); // display alpha vector %} END diff --git a/mcxtrace-comps/contrib/Bragg_crystal_bent_BC.comp b/mcxtrace-comps/contrib/Bragg_crystal_bent_BC.comp index 2318fa4c97..e4a51ca352 100644 --- a/mcxtrace-comps/contrib/Bragg_crystal_bent_BC.comp +++ b/mcxtrace-comps/contrib/Bragg_crystal_bent_BC.comp @@ -116,268 +116,296 @@ DECLARE INITIALIZE %{ - int status; - if (material){ - if ((status=Table_Read(&(m_t),material,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,material); - exit(-1); - } - char **header_parsed; - header_parsed=Table_ParseHeader(m_t.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL); - if(header_parsed[2]){rho=strtod(header_parsed[2],NULL);} - if(header_parsed[0]){Z=strtod(header_parsed[0],NULL);} - if(header_parsed[1]){At=strtod(header_parsed[1],NULL);} - }else{ - fprintf(stderr,"Error(%s): No material file specified\n",NAME_CURRENT_COMP); + int status; + if (material) { + if ((status = Table_Read (&(m_t), material, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, material); + exit (-1); } - if(form_factors){ - if ((status=Table_Read(&(f0_t),form_factors,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,form_factors); - exit(-1); - } + char** header_parsed; + header_parsed = Table_ParseHeader (m_t.header, "Z", "A[r]", "rho", "Z/A", "sigma[a]", NULL); + if (header_parsed[2]) { + rho = strtod (header_parsed[2], NULL); + } + if (header_parsed[0]) { + Z = strtod (header_parsed[0], NULL); + } + if (header_parsed[1]) { + At = strtod (header_parsed[1], NULL); + } + } else { + fprintf (stderr, "Error(%s): No material file specified\n", NAME_CURRENT_COMP); + } + if (form_factors) { + if ((status = Table_Read (&(f0_t), form_factors, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, form_factors); + exit (-1); } + } - a2inv=(x_a)?1/(x_a*x_a):0; /* 0 really means infinity for x direction */ - b2inv=(y_b)?1/(y_b*y_b):0; - c2inv=(z_c)?1/(z_c*z_c):0; + a2inv = (x_a) ? 1 / (x_a * x_a) : 0; /* 0 really means infinity for x direction */ + b2inv = (y_b) ? 1 / (y_b * y_b) : 0; + c2inv = (z_c) ? 1 / (z_c * z_c) : 0; - l_a2inv=(lattice_x_a)?1/(lattice_x_a*lattice_x_a):0; /* 0 really means infinity for x direction */ - l_b2inv=(lattice_y_b)?1/(lattice_y_b*lattice_y_b):0; - l_c2inv=(lattice_z_c)?1/(lattice_z_c*lattice_z_c):0; + l_a2inv = (lattice_x_a) ? 1 / (lattice_x_a * lattice_x_a) : 0; /* 0 really means infinity for x direction */ + l_b2inv = (lattice_y_b) ? 1 / (lattice_y_b * lattice_y_b) : 0; + l_c2inv = (lattice_z_c) ? 1 / (lattice_z_c * lattice_z_c) : 0; - cos_alpha0=cos(alpha); - sin_alpha0=sin(alpha); + cos_alpha0 = cos (alpha); + sin_alpha0 = sin (alpha); - if (verbose){ - printf("INFO (%s): curve=(%.5f %.5f %.5f) lcurve=(%.3f %.5f %.5f)\n\n",NAME_CURRENT_COMP, - a2inv, b2inv, c2inv, l_a2inv, l_b2inv, l_c2inv); - printf("INFO (%s): initialized\n",NAME_CURRENT_COMP); - } + if (verbose) { + printf ("INFO (%s): curve=(%.5f %.5f %.5f) lcurve=(%.3f %.5f %.5f)\n\n", NAME_CURRENT_COMP, a2inv, b2inv, c2inv, l_a2inv, l_b2inv, l_c2inv); + printf ("INFO (%s): initialized\n", NAME_CURRENT_COMP); + } %} TRACE %{ - double E; // (keV) x-ray energy - double K; // length of k-vector - double kxu,kyu,kzu; // unit vector in the direction of k-vector. - double x_int,y_int,z_int; // intersection with the y=0 plane - double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). - //double Thetain; // (rad) angle between the crystal surface and the incident ray - //double Theta0; // (rad) angle between the Bragg planes and the incident ray - //double Thetah; // (rad) angle between the Bragg planes and the reflected ray - //double DeltaTheta0; // (rad) the center of the reflectivity curve is at asin(n*lambda/(2*d)) + DeltaTheta0 - //double Rpi, Rsig; - double R; // Reflectivity value calculated by Mx_DarwinReflectivity() function for each incoming photon - - /* get the photon's kvector and energy */ - K=sqrt(kx*kx+ky*ky+kz*kz); - E = K2E*K; /* use built-in constants for consistency */ - /* make unit vector in the direction of k :*/ - kxu = kx; kyu = ky; kzu = kz; - NORM(kxu,kyu,kzu); - double k0hat[3]={kxu,kyu,kzu}; - /* printf("incoming kx,ky,kz, Ex, Ey, Ez, k.E: %f %f %f %g %g %g %g\n", kx,ky,kz,Ex,Ey,Ez, kxu*Ex+kyu*Ey+kzu*Ez); */ - + double E; // (keV) x-ray energy + double K; // length of k-vector + double kxu, kyu, kzu; // unit vector in the direction of k-vector. + double x_int, y_int, z_int; // intersection with the y=0 plane + double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). + // double Thetain; // (rad) angle between the crystal surface and the incident ray + // double Theta0; // (rad) angle between the Bragg planes and the incident ray + // double Thetah; // (rad) angle between the Bragg planes and the reflected ray + // double DeltaTheta0; // (rad) the center of the reflectivity curve is at asin(n*lambda/(2*d)) + DeltaTheta0 + // double Rpi, Rsig; + double R; // Reflectivity value calculated by Mx_DarwinReflectivity() function for each incoming photon + + /* get the photon's kvector and energy */ + K = sqrt (kx * kx + ky * ky + kz * kz); + E = K2E * K; /* use built-in constants for consistency */ + /* make unit vector in the direction of k :*/ + kxu = kx; + kyu = ky; + kzu = kz; + NORM (kxu, kyu, kzu); + double k0hat[3] = { kxu, kyu, kzu }; + /* printf("incoming kx,ky,kz, Ex, Ey, Ez, k.E: %f %f %f %g %g %g %g\n", kx,ky,kz,Ex,Ey,Ez, kxu*Ex+kyu*Ey+kzu*Ez); */ + + #ifdef MCDEBUG + printf ("%s: starting trace, k0hat=(%.5f %.5f %.5f)\n", NAME_CURRENT_COMP, k0hat[0], k0hat[1], k0hat[2]); + #endif + + /* this intersection code copied from Mirror_elliptic.comp, with coordinates modified */ + double A, B, C, xt, yt, zt; + double t0, t1; + /*an offset to the mirror parameters perhaps*/ + + xt = x; + zt = z; + /*the reference point is on the ellipsoid surface such that the ellipsoid mass lies on the positive y-side of the zy-plane*/ + yt = y - y_b; + + C = xt * xt * a2inv + yt * yt * b2inv + zt * zt * c2inv - 1; + B = 2 * (kxu * xt * a2inv + kyu * yt * b2inv + kzu * zt * c2inv); + A = kxu * kxu * a2inv + kyu * kyu * b2inv + kzu * kzu * c2inv; + + if (solve_2nd_order (&t0, &t1, A, B, C)) { + double xx0, xx1, yy0, yy1, zz0, zz1; /* we will have to tentatively propagate twice to see which surface we hit */ + xx0 = x + kxu * t0; + yy0 = y + kyu * t0; + zz0 = z + kzu * t0; + xx1 = x + kxu * t1; + yy1 = y + kyu * t1; + zz1 = z + kzu * t1; + + /*Check if we hit the mirror and whether the hit it is in front of the ray. + * This does not account for mirror curvature */ + int hit0 = (fabs (xx0) < width / 2.0) && (fabs (zz0) < length / 2.0) && t0 > 0; + int hit1 = (fabs (xx1) < width / 2.0) && (fabs (zz1) < length / 2.0) && t1 > 0; + int doit = hit0 || hit1; + if (hit0 && !hit1) + PROP_DL (t0); /* only one intersection actually on mirror */ + else if (hit1 && !hit0) + PROP_DL (t1); /* other intersection */ + else if (hit0 && hit1) { /* both, take first strike (which may be back of mirror) */ + PROP_DL (t0 < t1 ? t0 : t1); + } else { + RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); + } #ifdef MCDEBUG - printf("%s: starting trace, k0hat=(%.5f %.5f %.5f)\n", NAME_CURRENT_COMP, - k0hat[0], k0hat[1], k0hat[2]); + printf ("%s: solving for crystal hit, l0=%f, l1=%f, x=(%f %f %f) doit=%d, hit0=%d, hit1=%d\n", NAME_CURRENT_COMP, t0, t1, x, y, z, doit, hit0, hit1); #endif - - /* this intersection code copied from Mirror_elliptic.comp, with coordinates modified */ - double A,B,C, xt, yt, zt; - double t0,t1; - /*an offset to the mirror parameters perhaps*/ - - xt=x; - zt=z; - /*the reference point is on the ellipsoid surface such that the ellipsoid mass lies on the positive y-side of the zy-plane*/ - yt=y-y_b; - - C=xt*xt*a2inv + yt*yt*b2inv + zt*zt*c2inv -1; - B=2*(kxu*xt*a2inv + kyu*yt*b2inv + kzu*zt*c2inv); - A=kxu*kxu*a2inv + kyu*kyu*b2inv + kzu*kzu*c2inv; - - if(solve_2nd_order(&t0,&t1,A,B,C)){ - double xx0, xx1, yy0, yy1, zz0, zz1; /* we will have to tentatively propagate twice to see which surface we hit */ - xx0=x+kxu*t0; yy0=y+kyu*t0; zz0=z+kzu*t0; - xx1=x+kxu*t1; yy1=y+kyu*t1; zz1=z+kzu*t1; - - /*Check if we hit the mirror and whether the hit it is in front of the ray. - * This does not account for mirror curvature */ - int hit0=(fabs(xx0)0; - int hit1=(fabs(xx1)0; - int doit=hit0 || hit1; - if(hit0 && !hit1) PROP_DL(t0); /* only one intersection actually on mirror */ - else if (hit1 && !hit0) PROP_DL(t1); /* other intersection */ - else if (hit0 && hit1) { /* both, take first strike (which may be back of mirror) */ - PROP_DL(t0(2*k*sin(theta_bragg)-dQ) && Q<(2*k*sin(theta_bragg)+dQ)) - {kz=-kz; + PROP_Z0; // MOVING TOWARDS THE CRYSTAL. + if (inside_rectangle (x, y, xwidth, yheight)) { + if (Q > (2 * k * sin (theta_bragg) - dQ) && Q < (2 * k * sin (theta_bragg) + dQ)) { + kz = -kz; SCATTER; - }else{ + } else { ABSORB; } } - %} MCDISPLAY %{ /* A bit ugly; hard-coded dimensions. */ - magnify(""); - line(-xwidth/2,-yheight/2,0,xwidth/2,-yheight/2,0); - line(-xwidth/2,yheight/2,0,xwidth/2,yheight/2,0); - line(-xwidth/2,-yheight/2,0,-xwidth/2,yheight/2,0); - line(xwidth/2,-yheight/2,0,xwidth/2,yheight/2,0); - + magnify (""); + line (-xwidth / 2, -yheight / 2, 0, xwidth / 2, -yheight / 2, 0); + line (-xwidth / 2, yheight / 2, 0, xwidth / 2, yheight / 2, 0); + line (-xwidth / 2, -yheight / 2, 0, -xwidth / 2, yheight / 2, 0); + line (xwidth / 2, -yheight / 2, 0, xwidth / 2, yheight / 2, 0); %} END diff --git a/mcxtrace-comps/contrib/Detector_pn.comp b/mcxtrace-comps/contrib/Detector_pn.comp index 13734db231..a6df5b582a 100644 --- a/mcxtrace-comps/contrib/Detector_pn.comp +++ b/mcxtrace-comps/contrib/Detector_pn.comp @@ -73,103 +73,104 @@ DECLARE INITIALIZE %{ - int status=0; - - if(!xwidth || !yheight){ - fprintf(stderr,"%s: Detector has zero effective area\n",NAME_CURRENT_COMP); - exit(0); + int status = 0; + + if (!xwidth || !yheight) { + fprintf (stderr, "%s: Detector has zero effective area\n", NAME_CURRENT_COMP); + exit (0); } - xmax=xwidth/2.0; - xmin=-xmax; - ymax=yheight/2.0; - ymin=-ymax; - - //t_Table T; - if ( (status=Table_Read(&T,material_datafile,0))==-1){ - fprintf(stderr,"Error: Could not parse file \"%s\" in COMP %s\n",material_datafile,NAME_CURRENT_COMP); - exit(-1); + xmax = xwidth / 2.0; + xmin = -xmax; + ymax = yheight / 2.0; + ymin = -ymax; + + // t_Table T; + if ((status = Table_Read (&T, material_datafile, 0)) == -1) { + fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", material_datafile, NAME_CURRENT_COMP); + exit (-1); + } + char** header_parsed; + header_parsed = Table_ParseHeader (T.header, "Z", "A[r]", "rho", "sigma[a]"); + // if (!At) At=strtod(header_parsed[1],NULL); + // if (!Z) Z=strtol(header_parsed[0],NULL,10); + if (!rho) + rho = strtod (header_parsed[2], NULL); + + if (xwidth > 0) { + xmax = xwidth / 2; + xmin = -xmax; + } + if (yheight > 0) { + ymax = yheight / 2; + ymin = -ymax; } - char **header_parsed; - header_parsed=Table_ParseHeader(T.header,"Z","A[r]","rho","sigma[a]"); - //if (!At) At=strtod(header_parsed[1],NULL); - //if (!Z) Z=strtol(header_parsed[0],NULL,10); - if (!rho) rho=strtod(header_parsed[2],NULL); - - if (xwidth > 0) { xmax = xwidth/2; xmin = -xmax; } - if (yheight > 0) { ymax = yheight/2; ymin = -ymax; } if ((xmin >= xmax) || (ymin >= ymax)) { - printf("Detector_pn: %s: Null detection area !\n" - "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", - NAME_CURRENT_COMP); - exit(0); + printf ("Detector_pn: %s: Null detection area !\n" + "ERROR (xwidth,yheight,xmin,xmax,ymin,ymax). Exiting", + NAME_CURRENT_COMP); + exit (0); } - PSD_N = create_darr2d(nx, ny); - PSD_p = create_darr2d(nx, ny); - PSD_p2 = create_darr2d(nx, ny); - + PSD_N = create_darr2d (nx, ny); + PSD_p = create_darr2d (nx, ny); + PSD_p2 = create_darr2d (nx, ny); + // Use instance name for monitor output if no input was given - if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP); + if (!strcmp (filename, "\0")) + sprintf (filename, "%s", NAME_CURRENT_COMP); %} TRACE %{ - double alpha,e,k,mu; - double l0,l1; - int i,j; - if (box_intersect(&l0,&l1,x,y,z,kx,ky,kz,xwidth,yheight,zdepth)){ - PROP_DL(l0); + double alpha, e, k, mu; + double l0, l1; + int i, j; + if (box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth)) { + PROP_DL (l0); /*table interpolation*/ - k=sqrt(kx*kx+ky*ky+kz*kz); - e=k*K2E; - mu=Table_Value(T,e,1)*rho*1e2; - - l1-=l0; - p*=(1-exp(-mu*l1)); - - //photon detected in surface pixel - if (x>xmin && xymin && y xmin && x < xmax && y > ymin && y < ymax) { + i = floor ((x - xmin) * nx / (xmax - xmin)); + j = floor ((y - ymin) * ny / (ymax - ymin)); + #pragma acc atomic PSD_N[i][j] = PSD_N[i][j] + 1; -#pragma acc atomic + #pragma acc atomic PSD_p[i][j] = PSD_p[i][j] + p; -#pragma acc atomic - PSD_p2[i][j] = PSD_p2[i][j] + p*p; + #pragma acc atomic + PSD_p2[i][j] = PSD_p2[i][j] + p * p; SCATTER; } if (restore_xray) { - RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); + RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); } } %} SAVE %{ - DETECTOR_OUT_2D( - "Detector_pn", - "X position [cm]", - "Y position [cm]", - xmin*100.0, xmax*100.0, ymin*100.0, ymax*100.0, - nx, ny, - &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], - filename); + DETECTOR_OUT_2D ("Detector_pn", "X position [cm]", "Y position [cm]", xmin * 100.0, xmax * 100.0, ymin * 100.0, ymax * 100.0, nx, ny, &PSD_N[0][0], + &PSD_p[0][0], &PSD_p2[0][0], filename); %} FINALLY %{ - Table_Free(&T); - destroy_darr2d(PSD_N); - destroy_darr2d(PSD_p); - destroy_darr2d(PSD_p2); + Table_Free (&T); + destroy_darr2d (PSD_N); + destroy_darr2d (PSD_p); + destroy_darr2d (PSD_p2); %} MCDISPLAY %{ - box(0,0,0,xwidth,yheight,zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/Laue_crystal_BC.comp b/mcxtrace-comps/contrib/Laue_crystal_BC.comp index bb8ec779b7..51627db2b9 100644 --- a/mcxtrace-comps/contrib/Laue_crystal_BC.comp +++ b/mcxtrace-comps/contrib/Laue_crystal_BC.comp @@ -100,233 +100,222 @@ DECLARE double f_rel; double f_nt; t_Table m_t; - t_Table f0_t; + t_Table f0_t; %} INITIALIZE %{ - int status; - if (material){ - if ((status=Table_Read(&(m_t),material,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,material); - exit(-1); - } - char **header_parsed; - header_parsed=Table_ParseHeader(m_t.header,"Z","A[r]","rho","Z/A","sigma[a]",NULL); - if(header_parsed[2]){rho=strtod(header_parsed[2],NULL);} - if(header_parsed[0]){Z=strtod(header_parsed[0],NULL);} - if(header_parsed[1]){At=strtod(header_parsed[1],NULL);} - }else{ - fprintf(stderr,"Error(%s): No material file specified\n",NAME_CURRENT_COMP); + int status; + if (material) { + if ((status = Table_Read (&(m_t), material, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, material); + exit (-1); } - if(form_factors){ - if ((status=Table_Read(&(f0_t),form_factors,0))==-1){ - fprintf(stderr,"Error(%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,form_factors); - exit(-1); - } + char** header_parsed; + header_parsed = Table_ParseHeader (m_t.header, "Z", "A[r]", "rho", "Z/A", "sigma[a]", NULL); + if (header_parsed[2]) { + rho = strtod (header_parsed[2], NULL); + } + if (header_parsed[0]) { + Z = strtod (header_parsed[0], NULL); + } + if (header_parsed[1]) { + At = strtod (header_parsed[1], NULL); } - - printf("INFO: (%s): initialized\n",NAME_CURRENT_COMP); + } else { + fprintf (stderr, "Error(%s): No material file specified\n", NAME_CURRENT_COMP); + } + if (form_factors) { + if ((status = Table_Read (&(f0_t), form_factors, 0)) == -1) { + fprintf (stderr, "Error(%s): Could not parse file \"%s\"\n", NAME_CURRENT_COMP, form_factors); + exit (-1); + } + } + + printf ("INFO: (%s): initialized\n", NAME_CURRENT_COMP); %} TRACE %{ - double E; // (keV) x-ray energy - double K; // length of k-vector - double tin; // 'time' of intersection of ray with y=0 plane (which include the crystal surface) - double x_int,y_int,z_int; // intersection with the y=0 plane - double dist; // distance from position at t=0 to the y=0 plane - double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). - - /* get the photon's kvector and energy */ - K=sqrt(kx*kx+ky*ky+kz*kz); - E = K2E*K; /* use built-in constants for consistency */ - /* make unit vector in the direction of k :*/ - double k0hat[3] = {kx/K, ky/K, kz/K}; // unit vector in the direction of k-vector. - double kk0[3]={kx,ky,kz}; - -#ifdef MCDEBUG - fflush(NULL); // make sure all error messages are in order! -#endif - - /*intersection calculation*/ - tin = (-y-thickness/2)/k0hat[1]; - if (tin > 1e-8){ - /* check whether our intersection lies within the boundaries of the crystal*/ - x_int=x+k0hat[0]*tin; - y_int=y+k0hat[1]*tin; - z_int=z+k0hat[2]*tin; + double E; // (keV) x-ray energy + double K; // length of k-vector + double tin; // 'time' of intersection of ray with y=0 plane (which include the crystal surface) + double x_int, y_int, z_int; // intersection with the y=0 plane + double dist; // distance from position at t=0 to the y=0 plane + double f00, f0h, fp, fpp; // atomic form factors for Q=0 is (f00 + fp + i*fpp) and for Q= ha*+kb*+lc* it is (f0h + fp + i*fpp). + + /* get the photon's kvector and energy */ + K = sqrt (kx * kx + ky * ky + kz * kz); + E = K2E * K; /* use built-in constants for consistency */ + /* make unit vector in the direction of k :*/ + double k0hat[3] = { kx / K, ky / K, kz / K }; // unit vector in the direction of k-vector. + double kk0[3] = { kx, ky, kz }; - if (fabs(x_int)<=width/2 && fabs(z_int)<=length/2){ - dist=sqrt(SQR(x-x_int)+SQR(y-y_int)+SQR(z-z_int)); - PROP_DL(dist); /* now the photon is on the crystal surface, ready to be reflected... */ - SCATTER; - double d=cbrt(V)/(sqrt(h*h+k*k+l*l));/*this is valid only for cubic structures*/ - f00 = Z; - f0h = Table_Value(f0_t,1/(2*d),Z); - fp = Table_Value(m_t,E,1)-Z; - fpp = Table_Value(m_t,E,2); + #ifdef MCDEBUG + fflush (NULL); // make sure all error messages are in order! + #endif - double nhat[3]={0,1,0}; // crystal surface normal is yhat, and is set up by mcxtrace to be into the crystal + /*intersection calculation*/ + tin = (-y - thickness / 2) / k0hat[1]; + if (tin > 1e-8) { + /* check whether our intersection lies within the boundaries of the crystal*/ + x_int = x + k0hat[0] * tin; + y_int = y + k0hat[1] * tin; + z_int = z + k0hat[2] * tin; - double alpha[3]={alphax, alphay, alphaz}; - int fail; + if (fabs (x_int) <= width / 2 && fabs (z_int) <= length / 2) { + dist = sqrt (SQR (x - x_int) + SQR (y - y_int) + SQR (z - z_int)); + PROP_DL (dist); /* now the photon is on the crystal surface, ready to be reflected... */ + SCATTER; + double d = cbrt (V) / (sqrt (h * h + k * k + l * l)); /*this is valid only for cubic structures*/ + f00 = Z; + f0h = Table_Value (f0_t, 1 / (2 * d), Z); + fp = Table_Value (m_t, E, 1) - Z; + fpp = Table_Value (m_t, E, 2); - double complex chi0, chih; - double k0mag, hscale, thetaB; - - Mx_CubicCrystalChi(&chi0, &chih, &k0mag, &hscale, &thetaB, - f00, f0h, fp, fpp, V, h, k, l, - debye_waller_B, E, - crystal_type,structure_factor_scale_r,structure_factor_scale_i); + double nhat[3] = { 0, 1, 0 }; // crystal surface normal is yhat, and is set up by mcxtrace to be into the crystal - double selsum, transmitlim, borrmannlim; - selsum=transmission_sampling+forward_diffraction_sampling+laue_sampling; - transmitlim=transmission_sampling/selsum; - borrmannlim=(transmission_sampling+forward_diffraction_sampling)/selsum; - double sampling_selector=rand01(); + double alpha[3] = { alphax, alphay, alphaz }; + int fail; - double Rsig, Rpi, Tsig, Tpi, Asig, Api; - double kh[3], sig_axis[3], pi_axis[3]; - fail=Mx_LaueReflectivityBC(&Rsig, &Rpi, &Tsig, &Tpi, &Asig, &Api, kh, // these are the return values - k0hat, nhat, alpha, - chi0, chih, chih, k0mag, hscale, thetaB, thickness - ); + double complex chi0, chih; + double k0mag, hscale, thetaB; - cross(sig_axis, k0hat, kh, 1); /* kin x kout is sigma direction for diffracted ray */ - /* pi is a vector perpendicular to k_in and sig i.e. the direction of pi polarization incoming */ - cross(pi_axis, k0hat, sig_axis, 1); - - /* resolve incoming polarization into sig and pi bits, and scale by sqrt(reflectivity) which is amplitude scale */ - double Esig=(Ex*sig_axis[0]+Ey*sig_axis[1]+Ez*sig_axis[2]); - double Epi= (Ex* pi_axis[0]+Ey* pi_axis[1]+Ez* pi_axis[2]); - if(Esig==0 && Epi==0) { /* someone didn't set the polarization direction; set it now to a random value and it will propagate */ - double psi=rand01()*PI/2; - Esig=cos(psi); Epi=sin(psi); - } - - /*, now hack out Borrmann (forward diffraction) estimate. Not really right! - Total forward beam is exact, just partitioning into split beam is estimated */ - double sigpart=2*Tsig/(Tsig+Asig); // should be 1 far away from anything interesting, 2 for enhanced - Asig=(1-((1-sigpart)*(1-sigpart)))*Tsig; // when sigpart=1, all is normal atten, otherwise share - Tsig=Tsig-Asig; - double pipart=2*Tpi/(Tpi+Api); // should be 1 far away from anything interesting, 2 for enhanced - Api=(1-((1-pipart)*(1-pipart)))*Tpi; // when sigpart=1, all is normal atten, otherwise share - Tpi=Tpi-Api; + Mx_CubicCrystalChi (&chi0, &chih, &k0mag, &hscale, &thetaB, f00, f0h, fp, fpp, V, h, k, l, debye_waller_B, E, crystal_type, structure_factor_scale_r, + structure_factor_scale_i); -#ifdef MCDEBUG - fprintf(stderr, "Laue %s, Rsig=%.3e, Rpi=%.3e, sigpart=%.3f, pipart=%.3f, Tsig=%.3e, Tpi=%.3e, Asig=%.3e, Api=%.3e\n", - NAME_CURRENT_COMP, Rsig, Rpi, sigpart, pipart, Tsig, Tpi, Asig, Api); -#endif + double selsum, transmitlim, borrmannlim; + selsum = transmission_sampling + forward_diffraction_sampling + laue_sampling; + transmitlim = transmission_sampling / selsum; + borrmannlim = (transmission_sampling + forward_diffraction_sampling) / selsum; + double sampling_selector = rand01 (); - /* R is the Laue reflected intensity, - T is the transmitted Bormann intensity, - A is the undiffracted beam */ + double Rsig, Rpi, Tsig, Tpi, Asig, Api; + double kh[3], sig_axis[3], pi_axis[3]; + fail = Mx_LaueReflectivityBC (&Rsig, &Rpi, &Tsig, &Tpi, &Asig, &Api, kh, // these are the return values + k0hat, nhat, alpha, chi0, chih, chih, k0mag, hscale, thetaB, thickness); - /* In all cases, the emergence from the crystal is only exactly right in the symmetric case, for now. - The assumption is the rays move along the crystal planes. - The real answer requires getting the normal to the dispersion surface, which I have not implemented. - For usual, thin crystals, this makes only a tiny difference.*/ - -#ifdef MCDEBUG - fprintf(stderr, "LAUE component %s: Rsig=%.3e Rpi=%.3e Tsig=%.3f Tpi=%.3e\n", - NAME_CURRENT_COMP, Rsig, Rpi, Tsig, Tpi); -#endif + cross (sig_axis, k0hat, kh, 1); /* kin x kout is sigma direction for diffracted ray */ + /* pi is a vector perpendicular to k_in and sig i.e. the direction of pi polarization incoming */ + cross (pi_axis, k0hat, sig_axis, 1); - if(transmission_sampling && ( sampling_selector<= transmitlim )) { - double EsigT=Esig*sqrt(Asig); /* Tsig & Tpi are the un-diffracted beam intensity */ - double EpiT=Epi*sqrt(Api); - double RT=EsigT*EsigT+EpiT*EpiT; /* projected transmitted, squared back to intensity */ - Ex=EpiT*pi_axis[0]+EsigT*sig_axis[0]; - Ey=EpiT*pi_axis[1]+EsigT*sig_axis[1]; - Ez=EpiT*pi_axis[2]+EsigT*sig_axis[2]; - p*=RT/(transmission_sampling/selsum); - PROP_DL(thickness/k0hat[1]); // continue through crystal to back side in original direction - } else { // doing either forward diffraction (Borrman anomalous transmission) or Laue diffraction. - // for both types of diffraction, the beam propagates inside the crystal - // along the lattice plane (perpendicular to diffraction plane and alpha) (B&C figure 2, again). - double propdir[3]; - cross(propdir, sig_axis, alpha, 1); // beam appears to move this direction - // now bump beam to back side of crystal - x+=thickness*propdir[0]/propdir[1]; - y+=thickness; - z+=thickness*propdir[2]/propdir[1]; - // now, decide between forward (Borrman) diffraction and Laue - if(forward_diffraction_sampling && (sampling_selector < borrmannlim)) { - // transmitted direction is same as incoming, no k vector change - double EsigT=Esig*sqrt(Tsig); /* Fsig and Fpi are the Borrmann intensity */ - double EpiT=Epi*sqrt(Tpi); - double RT=EsigT*EsigT+EpiT*EpiT; /* projected Borrmann, squared back to intensity */ - Ex=EpiT*pi_axis[0]+EsigT*sig_axis[0]; - Ey=EpiT*pi_axis[1]+EsigT*sig_axis[1]; - Ez=EpiT*pi_axis[2]+EsigT*sig_axis[2]; - p*=RT/(forward_diffraction_sampling/selsum); - } else if(laue_sampling) { - /* update outgoing direction vector kx, ky, kz = kh, fully scaled diffracted k vector */ - kx=kh[0]; ky=kh[1]; kz=kh[2]; - double pi1_axis[3]; - cross(pi1_axis, kh, sig_axis,1); - /* pi1 is now a vector perpendicular to k_out and sig i.e. the direction of pi polarization outgoing */ - double EsigR=Esig*sqrt(Rsig); /* Rsig and Rpi are the Laue intensity */ - double EpiR=Epi*sqrt(Rpi); - double RR=EsigR*EsigR+EpiR*EpiR; /* projected diffracted reflectivity, squared back to intensity */ - Ex=EpiR*pi1_axis[0]+EsigR*sig_axis[0]; - Ey=EpiR*pi1_axis[1]+EsigR*sig_axis[1]; - Ez=EpiR*pi1_axis[2]+EsigR*sig_axis[2]; - p*=RR/(laue_sampling/selsum); - } else { - // no diffraction sampled! illegal. - fprintf(stderr, "LAUE component %s: all sampling choices set to zero. cannot run.\n", - NAME_CURRENT_COMP); - exit(1); - } - NORM(Ex, Ey, Ez); - if (p==0) ABSORB; - } + /* resolve incoming polarization into sig and pi bits, and scale by sqrt(reflectivity) which is amplitude scale */ + double Esig = (Ex * sig_axis[0] + Ey * sig_axis[1] + Ez * sig_axis[2]); + double Epi = (Ex * pi_axis[0] + Ey * pi_axis[1] + Ez * pi_axis[2]); + if (Esig == 0 && Epi == 0) { /* someone didn't set the polarization direction; set it now to a random value and it will propagate */ + double psi = rand01 () * PI / 2; + Esig = cos (psi); + Epi = sin (psi); + } + + /*, now hack out Borrmann (forward diffraction) estimate. Not really right! + Total forward beam is exact, just partitioning into split beam is estimated */ + double sigpart = 2 * Tsig / (Tsig + Asig); // should be 1 far away from anything interesting, 2 for enhanced + Asig = (1 - ((1 - sigpart) * (1 - sigpart))) * Tsig; // when sigpart=1, all is normal atten, otherwise share + Tsig = Tsig - Asig; + double pipart = 2 * Tpi / (Tpi + Api); // should be 1 far away from anything interesting, 2 for enhanced + Api = (1 - ((1 - pipart) * (1 - pipart))) * Tpi; // when sigpart=1, all is normal atten, otherwise share + Tpi = Tpi - Api; + + #ifdef MCDEBUG + fprintf (stderr, "Laue %s, Rsig=%.3e, Rpi=%.3e, sigpart=%.3f, pipart=%.3f, Tsig=%.3e, Tpi=%.3e, Asig=%.3e, Api=%.3e\n", NAME_CURRENT_COMP, Rsig, Rpi, + sigpart, pipart, Tsig, Tpi, Asig, Api); + #endif + + /* R is the Laue reflected intensity, + T is the transmitted Bormann intensity, + A is the undiffracted beam */ + + /* In all cases, the emergence from the crystal is only exactly right in the symmetric case, for now. + The assumption is the rays move along the crystal planes. + The real answer requires getting the normal to the dispersion surface, which I have not implemented. + For usual, thin crystals, this makes only a tiny difference.*/ + + #ifdef MCDEBUG + fprintf (stderr, "LAUE component %s: Rsig=%.3e Rpi=%.3e Tsig=%.3f Tpi=%.3e\n", NAME_CURRENT_COMP, Rsig, Rpi, Tsig, Tpi); + #endif + + if (transmission_sampling && (sampling_selector <= transmitlim)) { + double EsigT = Esig * sqrt (Asig); /* Tsig & Tpi are the un-diffracted beam intensity */ + double EpiT = Epi * sqrt (Api); + double RT = EsigT * EsigT + EpiT * EpiT; /* projected transmitted, squared back to intensity */ + Ex = EpiT * pi_axis[0] + EsigT * sig_axis[0]; + Ey = EpiT * pi_axis[1] + EsigT * sig_axis[1]; + Ez = EpiT * pi_axis[2] + EsigT * sig_axis[2]; + p *= RT / (transmission_sampling / selsum); + PROP_DL (thickness / k0hat[1]); // continue through crystal to back side in original direction + } else { // doing either forward diffraction (Borrman anomalous transmission) or Laue diffraction. + // for both types of diffraction, the beam propagates inside the crystal + // along the lattice plane (perpendicular to diffraction plane and alpha) (B&C figure 2, again). + double propdir[3]; + cross (propdir, sig_axis, alpha, 1); // beam appears to move this direction + // now bump beam to back side of crystal + x += thickness * propdir[0] / propdir[1]; + y += thickness; + z += thickness * propdir[2] / propdir[1]; + // now, decide between forward (Borrman) diffraction and Laue + if (forward_diffraction_sampling && (sampling_selector < borrmannlim)) { + // transmitted direction is same as incoming, no k vector change + double EsigT = Esig * sqrt (Tsig); /* Fsig and Fpi are the Borrmann intensity */ + double EpiT = Epi * sqrt (Tpi); + double RT = EsigT * EsigT + EpiT * EpiT; /* projected Borrmann, squared back to intensity */ + Ex = EpiT * pi_axis[0] + EsigT * sig_axis[0]; + Ey = EpiT * pi_axis[1] + EsigT * sig_axis[1]; + Ez = EpiT * pi_axis[2] + EsigT * sig_axis[2]; + p *= RT / (forward_diffraction_sampling / selsum); + } else if (laue_sampling) { + /* update outgoing direction vector kx, ky, kz = kh, fully scaled diffracted k vector */ + kx = kh[0]; + ky = kh[1]; + kz = kh[2]; + double pi1_axis[3]; + cross (pi1_axis, kh, sig_axis, 1); + /* pi1 is now a vector perpendicular to k_out and sig i.e. the direction of pi polarization outgoing */ + double EsigR = Esig * sqrt (Rsig); /* Rsig and Rpi are the Laue intensity */ + double EpiR = Epi * sqrt (Rpi); + double RR = EsigR * EsigR + EpiR * EpiR; /* projected diffracted reflectivity, squared back to intensity */ + Ex = EpiR * pi1_axis[0] + EsigR * sig_axis[0]; + Ey = EpiR * pi1_axis[1] + EsigR * sig_axis[1]; + Ez = EpiR * pi1_axis[2] + EsigR * sig_axis[2]; + p *= RR / (laue_sampling / selsum); } else { - RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); + // no diffraction sampled! illegal. + fprintf (stderr, "LAUE component %s: all sampling choices set to zero. cannot run.\n", NAME_CURRENT_COMP); + exit (1); } + NORM (Ex, Ey, Ez); + if (p == 0) + ABSORB; + } + } else { + RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); } - + } %} MCDISPLAY %{ - double w2=width/2., l2=length/2., t2=thickness/2., wl2=(width+length)/2.0, arrow=wl2/10.0; - magnify(""); - multiline(9, - -w2, -t2, -l2, -w2, t2, -l2, w2, t2, -l2, w2, -t2, -l2, // front face - w2, -t2, l2, // connector - -w2, -t2, l2, -w2, t2, l2, w2, t2, l2, w2, -t2, l2 // back face - ); - line(-w2, -t2, -l2, -w2, -t2, l2); // extra edges - line(-w2, t2, -l2, -w2, t2, l2); // extra edges - line( w2, t2, -l2, w2, t2, l2); // extra edges - - multiline(8, - 0.,0.,0., - alphax*wl2, alphay*wl2, alphaz*wl2, - alphax*wl2, alphay*(wl2-arrow)+alphaz*arrow, alphaz*(wl2-arrow)+alphay*arrow, - alphax*wl2, alphay*(wl2-arrow)-alphaz*arrow, alphaz*(wl2-arrow)+alphay*arrow, - alphax*wl2, alphay*wl2, alphaz*wl2, - alphax*(wl2-arrow)+alphaz*arrow, alphay*wl2, alphaz*(wl2-arrow)+alphax*arrow, - alphax*(wl2-arrow)-alphaz*arrow, alphay*wl2, alphaz*(wl2-arrow)+alphax*arrow, - alphax*wl2, alphay*wl2, alphaz*wl2 - ); // display alpha vector assuming it is mostly in the 'z' direction + double w2 = width / 2., l2 = length / 2., t2 = thickness / 2., wl2 = (width + length) / 2.0, arrow = wl2 / 10.0; + magnify (""); + multiline (9, -w2, -t2, -l2, -w2, t2, -l2, w2, t2, -l2, w2, -t2, -l2, // front face + w2, -t2, l2, // connector + -w2, -t2, l2, -w2, t2, l2, w2, t2, l2, w2, -t2, l2 // back face + ); + line (-w2, -t2, -l2, -w2, -t2, l2); // extra edges + line (-w2, t2, -l2, -w2, t2, l2); // extra edges + line (w2, t2, -l2, w2, t2, l2); // extra edges - double nx=0., ny=1.0/2.0, nz=0.; // short surface normal - multiline(8, - 0.,0.,0., - nx*wl2, ny*wl2, nz*wl2, - nx*wl2, ny*(wl2-arrow)+nz*arrow, nz*(wl2-arrow)+ny*arrow, - nx*wl2, ny*(wl2-arrow)+nz*arrow, nz*(wl2-arrow)-ny*arrow, - nx*wl2, ny*wl2, nz*wl2, - nx*(wl2-arrow)+ny*arrow, ny*(wl2-arrow)+nx*arrow, nz*wl2, - nx*(wl2-arrow)-ny*arrow, ny*(wl2-arrow)+nx*arrow, nz*wl2, - nx*wl2, ny*wl2, nz*wl2 - ); // display surface normal vector, pointing in forward beam direction, - // half the size of alpha vector to make it easy to tell them apart + multiline (8, 0., 0., 0., alphax* wl2, alphay* wl2, alphaz* wl2, alphax* wl2, alphay*(wl2 - arrow) + alphaz * arrow, alphaz*(wl2 - arrow) + alphay * arrow, + alphax* wl2, alphay*(wl2 - arrow) - alphaz * arrow, alphaz*(wl2 - arrow) + alphay * arrow, alphax* wl2, alphay* wl2, alphaz* wl2, + alphax*(wl2 - arrow) + alphaz * arrow, alphay* wl2, alphaz*(wl2 - arrow) + alphax * arrow, alphax*(wl2 - arrow) - alphaz * arrow, alphay* wl2, + alphaz*(wl2 - arrow) + alphax * arrow, alphax* wl2, alphay* wl2, alphaz* wl2); // display alpha vector assuming it is mostly in the 'z' direction + double nx = 0., ny = 1.0 / 2.0, nz = 0.; // short surface normal + multiline (8, 0., 0., 0., nx* wl2, ny* wl2, nz* wl2, nx* wl2, ny*(wl2 - arrow) + nz * arrow, nz*(wl2 - arrow) + ny * arrow, nx* wl2, + ny*(wl2 - arrow) + nz * arrow, nz*(wl2 - arrow) - ny * arrow, nx* wl2, ny* wl2, nz* wl2, nx*(wl2 - arrow) + ny * arrow, + ny*(wl2 - arrow) + nx * arrow, nz* wl2, nx*(wl2 - arrow) - ny * arrow, ny*(wl2 - arrow) + nx * arrow, nz* wl2, nx* wl2, ny* wl2, + nz* wl2); // display surface normal vector, pointing in forward beam direction, + // half the size of alpha vector to make it easy to tell them apart %} END diff --git a/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp b/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp index dae6d665f8..42e4d142e1 100644 --- a/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp +++ b/mcxtrace-comps/contrib/Mirror_toroid_pothole.comp @@ -45,130 +45,153 @@ SETTING PARAMETERS (zdepth=0.1, xwidth=0.01, radius, radius_o,R0=0,string coatin SHARE %{ #include - %include "read_table-lib" + %include "read_table-lib" %include "reflectivity-lib" struct potholestruct { - double e_min,e_max,e_step,theta_min,theta_max,theta_step; - int use_reflec_table; + double e_min, e_max, e_step, theta_min, theta_max, theta_step; + int use_reflec_table; }; %} DECLARE %{ - struct potholestruct prms; - t_Table reflec_table; - t_Reflec re; + struct potholestruct prms; + t_Table reflec_table; + t_Reflec re; %} INITIALIZE %{ - if (coating && strlen(coating)){ - char **header_parsed; - t_Table *tp=&reflec_table; - /* read 1st block data from file into tp */ - if (Table_Read(tp, coating, 1) <= 0) - { - exit(fprintf(stderr,"Error: %s: cannot read file %s\n",NAME_CURRENT_COMP, coating)); - } - header_parsed = Table_ParseHeader(tp->header, - "e_min=","e_max=","e_step=","theta_min=","theta_max=","theta_step=",NULL); - if (header_parsed[0] && header_parsed[1] && header_parsed[2] && - header_parsed[3] && header_parsed[4] && header_parsed[5]) - { - prms.e_min=strtod(header_parsed[0],NULL); - prms.e_max=strtod(header_parsed[1],NULL); - prms.e_step=strtod(header_parsed[2],NULL); - prms.theta_min=strtod(header_parsed[3],NULL); - prms.theta_max=strtod(header_parsed[4],NULL); - prms.theta_step=strtod(header_parsed[5],NULL); - } else { - exit(fprintf(stderr,"Error: %s: wrong/missing header line(s) in file %s\n", NAME_CURRENT_COMP, coating)); - } - if (((tp->rows-2)*prms.e_step) > (prms.e_max-prms.e_min)) - { - exit(fprintf(stderr,"Error: %s: e_step does not match e_min and e_max in file %s (<)\n",NAME_CURRENT_COMP, coating)); - } - if ((prms.e_max-prms.e_min) > ((tp->rows)*prms.e_step)) - { - exit(fprintf(stderr,"Error: %s: e_step does not match e_min and e_max in file %s (>)\n",NAME_CURRENT_COMP, coating)); - } - if (((tp->columns-2)*prms.theta_step) > (prms.theta_max-prms.theta_min)) - { - exit(fprintf(stderr,"Error: %s: theta_step does not match theta_min and theta_max in file %s (<)\n",NAME_CURRENT_COMP,coating)); - } - if ((prms.theta_max-prms.theta_min) > ((tp->columns)*prms.theta_step)) - { - exit(fprintf(stderr,"Error: %s: theta_step does not match theta_min and theta_max in file %s (>)\n",NAME_CURRENT_COMP,coating)); - } - prms.use_reflec_table=1; - }else{ - prms.use_reflec_table=0; + if (coating && strlen (coating)) { + char** header_parsed; + t_Table* tp = &reflec_table; + /* read 1st block data from file into tp */ + if (Table_Read (tp, coating, 1) <= 0) { + exit (fprintf (stderr, "Error: %s: cannot read file %s\n", NAME_CURRENT_COMP, coating)); } + header_parsed = Table_ParseHeader (tp->header, "e_min=", "e_max=", "e_step=", "theta_min=", "theta_max=", "theta_step=", NULL); + if (header_parsed[0] && header_parsed[1] && header_parsed[2] && header_parsed[3] && header_parsed[4] && header_parsed[5]) { + prms.e_min = strtod (header_parsed[0], NULL); + prms.e_max = strtod (header_parsed[1], NULL); + prms.e_step = strtod (header_parsed[2], NULL); + prms.theta_min = strtod (header_parsed[3], NULL); + prms.theta_max = strtod (header_parsed[4], NULL); + prms.theta_step = strtod (header_parsed[5], NULL); + } else { + exit (fprintf (stderr, "Error: %s: wrong/missing header line(s) in file %s\n", NAME_CURRENT_COMP, coating)); + } + if (((tp->rows - 2) * prms.e_step) > (prms.e_max - prms.e_min)) { + exit (fprintf (stderr, "Error: %s: e_step does not match e_min and e_max in file %s (<)\n", NAME_CURRENT_COMP, coating)); + } + if ((prms.e_max - prms.e_min) > ((tp->rows) * prms.e_step)) { + exit (fprintf (stderr, "Error: %s: e_step does not match e_min and e_max in file %s (>)\n", NAME_CURRENT_COMP, coating)); + } + if (((tp->columns - 2) * prms.theta_step) > (prms.theta_max - prms.theta_min)) { + exit (fprintf (stderr, "Error: %s: theta_step does not match theta_min and theta_max in file %s (<)\n", NAME_CURRENT_COMP, coating)); + } + if ((prms.theta_max - prms.theta_min) > ((tp->columns) * prms.theta_step)) { + exit (fprintf (stderr, "Error: %s: theta_step does not match theta_min and theta_max in file %s (>)\n", NAME_CURRENT_COMP, coating)); + } + prms.use_reflec_table = 1; + } else { + prms.use_reflec_table = 0; + } %} TRACE %{ - int status; - double l0,l1,l2,l3,tx,ty,tz; - + int status; + double l0, l1, l2, l3, tx, ty, tz; + + do { + /*TODO check the permutation of coordinates and the actual position of the cylinder.*/ + status = cylinder_intersect (&l0, &l1, x, z, y - radius, kx, kz, ky, radius, zdepth); + if (!status) + break; /* no interaction with the cylinder*/ + // if(status & (02|04)) break; /*we exit the top/bottom of the cylinder*/ + if (status & (24)) + break; + // if(status & (8|16)) break; /*we exit the top/bottom of the cylinder*/ + /*if the first intersection is behind the particle - this means the ray is on the wrong side of the mirror*/ + if (l1 < 0) + break; do { - /*TODO check the permutation of coordinates and the actual position of the cylinder.*/ - status= cylinder_intersect(&l0,&l1,x,z,y-radius,kx,kz,ky, radius, zdepth); - if (!status) break; /* no interaction with the cylinder*/ - //if(status & (02|04)) break; /*we exit the top/bottom of the cylinder*/ - if(status & (24)) break; - //if(status & (8|16)) break; /*we exit the top/bottom of the cylinder*/ - /*if the first intersection is behind the particle - this means the ray is on the wrong side of the mirror*/ - if(l1<0) break; - do { - /*this to do a test propagation*/ - double op[12]; - op[0]=x; op[1]=y; op[2]=z; op[3]=kx; op[4]=ky; op[5]=kz; op[6]=phi; op[7]=t; op[8]=Ex; op[9]=Ey; op[10]=Ez; op[12]=p; - mcPROP_DL(l1); - (tx)=x;(ty)=y; (tz)=z; - x=op[0]; y=op[1]; z=op[2]; kx=op[3]; ky=op[4]; kz=op[5]; phi=op[6]; t=op[7]; Ex=op[8]; Ey=op[9]; Ez=op[10]; p=op[12]; - }while(0); - /*check mirror limits of intersectio point.*/ - if(tz<-zdepth/2.0 || tz>zdepth/2.0) break; - /*check if the width is OK*/ - double xmax=acos(xwidth/(2.0*radius))*radius; - if(tx<-xmax || tx>xmax) break; - - status=ellipsoid_intersect(&l2,&l3,x,y-radius,z,kx,ky,kz,radius,radius,radius_o+radius,NULL); - if (!status) break; - if(l3<0) break; /*This shouldn't be possible*/ - /*the mirror is indeed hit*/ - PROP_DL(l3); - SCATTER; - - double nx,ny,nz; - nx=2*x/(radius*radius); - ny=2*(y-radius)/(radius*radius);/*ellipsoid is displaced to put Origin on the surface*/ - nz=2*z/((radius+radius_o)*(radius+radius_o)); - NORM(nx,ny,nz); - - /*By definition the normal vector points out from the ellipsoid*/ - double s=scalar_prod(kx,ky,kz,nx,ny,nz); - double k=sqrt(scalar_prod(kx,ky,kz,kx,ky,kz)); - - kx=kx-2*s*nx; - ky=ky-2*s*ny; - kz=kz-2*s*nz; - - /*find energy, and glancing angle*/ - if (prms.use_reflec_table){ - /*the get ref. by call to Table_value2d*/ - double R; - double theta=RAD2DEG*(acos(s/k)-M_PI_2); - double e=K2E*k; - R=Table_Value2d(reflec_table,fabs((e-prms.e_min)/prms.e_step), fabs(((theta-prms.theta_min)/prms.theta_step))); - p*=R; - /*update phase - as an approximation turn by 180 deg.*/ - phi+=M_PI; - }else{ - p*=R0; - } - }while(0); + /*this to do a test propagation*/ + double op[12]; + op[0] = x; + op[1] = y; + op[2] = z; + op[3] = kx; + op[4] = ky; + op[5] = kz; + op[6] = phi; + op[7] = t; + op[8] = Ex; + op[9] = Ey; + op[10] = Ez; + op[12] = p; + mcPROP_DL (l1); + (tx) = x; + (ty) = y; + (tz) = z; + x = op[0]; + y = op[1]; + z = op[2]; + kx = op[3]; + ky = op[4]; + kz = op[5]; + phi = op[6]; + t = op[7]; + Ex = op[8]; + Ey = op[9]; + Ez = op[10]; + p = op[12]; + } while (0); + /*check mirror limits of intersectio point.*/ + if (tz < -zdepth / 2.0 || tz > zdepth / 2.0) + break; + /*check if the width is OK*/ + double xmax = acos (xwidth / (2.0 * radius)) * radius; + if (tx < -xmax || tx > xmax) + break; + + status = ellipsoid_intersect (&l2, &l3, x, y - radius, z, kx, ky, kz, radius, radius, radius_o + radius, NULL); + if (!status) + break; + if (l3 < 0) + break; /*This shouldn't be possible*/ + /*the mirror is indeed hit*/ + PROP_DL (l3); + SCATTER; + + double nx, ny, nz; + nx = 2 * x / (radius * radius); + ny = 2 * (y - radius) / (radius * radius); /*ellipsoid is displaced to put Origin on the surface*/ + nz = 2 * z / ((radius + radius_o) * (radius + radius_o)); + NORM (nx, ny, nz); + + /*By definition the normal vector points out from the ellipsoid*/ + double s = scalar_prod (kx, ky, kz, nx, ny, nz); + double k = sqrt (scalar_prod (kx, ky, kz, kx, ky, kz)); + + kx = kx - 2 * s * nx; + ky = ky - 2 * s * ny; + kz = kz - 2 * s * nz; + + /*find energy, and glancing angle*/ + if (prms.use_reflec_table) { + /*the get ref. by call to Table_value2d*/ + double R; + double theta = RAD2DEG * (acos (s / k) - M_PI_2); + double e = K2E * k; + R = Table_Value2d (reflec_table, fabs ((e - prms.e_min) / prms.e_step), fabs (((theta - prms.theta_min) / prms.theta_step))); + p *= R; + /*update phase - as an approximation turn by 180 deg.*/ + phi += M_PI; + } else { + p *= R0; + } + } while (0); %} @@ -176,71 +199,69 @@ TRACE MCDISPLAY %{ - //See parametric equation used for the 3D view here: https://en.wikipedia.org/wiki/Torus#Geometry or here:https://web.maths.unsw.edu.au/~rsw/Torus/index.php - //Because the code above doesn't actually reproduce a torus exactly but uses a cylinder and an ellipsoid to construct the toroidal mirror (to avoid solving quartic equations), - //there's a slight mismatch between the 3D view (exact torus) and the actual toroidal mirror coded above. - - double x0,y0,z0,x1,y1,z1,theta,phi,phi_total,theta_total,theta_beginning,phi_beginning,theta_end,phi_end; - //int N=50; too many points, makes the 3D matlab viewer bug - int N=10; - - phi_total = RAD2DEG*(zdepth/(radius+radius_o));//degrees - theta_total = RAD2DEG*(xwidth/radius); - - //fix theta=0 and phi=180 degrees respectively - theta=0; - phi=180; - - //then start at the beginning of the theta or phi range in order to map the flat surface (xwidth, zdepth) on the torus. - theta_beginning=theta-(theta_total/2); - phi_beginning=phi-(phi_total/2); - - theta_end = theta+(theta_total/2); - phi_end = phi+(phi_total/2); - + // See parametric equation used for the 3D view here: https://en.wikipedia.org/wiki/Torus#Geometry or here:https://web.maths.unsw.edu.au/~rsw/Torus/index.php + // Because the code above doesn't actually reproduce a torus exactly but uses a cylinder and an ellipsoid to construct the toroidal mirror (to avoid solving + // quartic equations), there's a slight mismatch between the 3D view (exact torus) and the actual toroidal mirror coded above. + + double x0, y0, z0, x1, y1, z1, theta, phi, phi_total, theta_total, theta_beginning, phi_beginning, theta_end, phi_end; + // int N=50; too many points, makes the 3D matlab viewer bug + int N = 10; + + phi_total = RAD2DEG * (zdepth / (radius + radius_o)); // degrees + theta_total = RAD2DEG * (xwidth / radius); + + // fix theta=0 and phi=180 degrees respectively + theta = 0; + phi = 180; + + // then start at the beginning of the theta or phi range in order to map the flat surface (xwidth, zdepth) on the torus. + theta_beginning = theta - (theta_total / 2); + phi_beginning = phi - (phi_total / 2); + + theta_end = theta + (theta_total / 2); + phi_end = phi + (phi_total / 2); + theta = theta_beginning; - phi = phi_beginning; - - while (theta<=theta_end){ - phi=180-(phi_total/2); - x0=radius*sin(DEG2RAD*theta); - y0=(radius_o+radius*cos(DEG2RAD*theta))*cos(DEG2RAD*phi)+(radius_o+radius); - z0=(radius_o+radius*cos(DEG2RAD*theta))*sin(DEG2RAD*phi); - while (phi<=phi_end){ - x1=radius*sin(DEG2RAD*theta); - y1=(radius_o+radius*cos(DEG2RAD*theta))*cos(DEG2RAD*phi)+(radius_o+radius); - z1=(radius_o+radius*cos(DEG2RAD*theta))*sin(DEG2RAD*phi); - line(x0,y0,z0,x1,y1,z1); - x0=x1; - y0=y1; - z0=z1; - phi+=phi_total/N; - - } - theta+=theta_total/N; - } - + phi = phi_beginning; + + while (theta <= theta_end) { + phi = 180 - (phi_total / 2); + x0 = radius * sin (DEG2RAD * theta); + y0 = (radius_o + radius * cos (DEG2RAD * theta)) * cos (DEG2RAD * phi) + (radius_o + radius); + z0 = (radius_o + radius * cos (DEG2RAD * theta)) * sin (DEG2RAD * phi); + while (phi <= phi_end) { + x1 = radius * sin (DEG2RAD * theta); + y1 = (radius_o + radius * cos (DEG2RAD * theta)) * cos (DEG2RAD * phi) + (radius_o + radius); + z1 = (radius_o + radius * cos (DEG2RAD * theta)) * sin (DEG2RAD * phi); + line (x0, y0, z0, x1, y1, z1); + x0 = x1; + y0 = y1; + z0 = z1; + phi += phi_total / N; + } + theta += theta_total / N; + } + theta = theta_beginning; - phi = phi_beginning; - - while (phi<=phi_end){ - theta=-(theta_total/2); - x0=radius*sin(DEG2RAD*theta); - y0=(radius_o+radius*cos(DEG2RAD*theta))*cos(DEG2RAD*phi)+(radius_o+radius); - z0=(radius_o+radius*cos(DEG2RAD*theta))*sin(DEG2RAD*phi); - while (theta<=theta_end){ - x1=radius*sin(DEG2RAD*theta); - y1=(radius_o+radius*cos(DEG2RAD*theta))*cos(DEG2RAD*phi)+(radius_o+radius); - z1=(radius_o+radius*cos(DEG2RAD*theta))*sin(DEG2RAD*phi); - line(x0,y0,z0,x1,y1,z1); - x0=x1; - y0=y1; - z0=z1; - theta+=theta_total/N; - } - phi+=phi_total/N; + phi = phi_beginning; + + while (phi <= phi_end) { + theta = -(theta_total / 2); + x0 = radius * sin (DEG2RAD * theta); + y0 = (radius_o + radius * cos (DEG2RAD * theta)) * cos (DEG2RAD * phi) + (radius_o + radius); + z0 = (radius_o + radius * cos (DEG2RAD * theta)) * sin (DEG2RAD * phi); + while (theta <= theta_end) { + x1 = radius * sin (DEG2RAD * theta); + y1 = (radius_o + radius * cos (DEG2RAD * theta)) * cos (DEG2RAD * phi) + (radius_o + radius); + z1 = (radius_o + radius * cos (DEG2RAD * theta)) * sin (DEG2RAD * phi); + line (x0, y0, z0, x1, y1, z1); + x0 = x1; + y0 = y1; + z0 = z1; + theta += theta_total / N; + } + phi += phi_total / N; } - %} diff --git a/mcxtrace-comps/contrib/PSD_monitor_rad.comp b/mcxtrace-comps/contrib/PSD_monitor_rad.comp index b5175430d4..b4aea6980e 100644 --- a/mcxtrace-comps/contrib/PSD_monitor_rad.comp +++ b/mcxtrace-comps/contrib/PSD_monitor_rad.comp @@ -61,24 +61,25 @@ INITIALIZE %{ int i; - PSDr_N = create_darr1d(nr); - PSDr_p = create_darr1d(nr); - PSDr_p2 = create_darr1d(nr); - PSDr_av_p = create_darr1d(nr); - PSDr_av_p2 = create_darr1d(nr); - - for (i=0; i QQ[NumberOfDatapoints - 1])) { @@ -222,7 +222,7 @@ TRACE Intensity = (Slope * q + Offset) / ForwardScattering; - p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * Intensity * exp(- Absorption * (l + l1)); + p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * Intensity * exp (-Absorption * (l + l1)); SCATTER; } @@ -230,7 +230,7 @@ TRACE MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSCylinders.comp b/mcxtrace-comps/contrib/SAXSCylinders.comp index 0a28e52b8a..f9a4fe3b51 100644 --- a/mcxtrace-comps/contrib/SAXSCylinders.comp +++ b/mcxtrace-comps/contrib/SAXSCylinders.comp @@ -58,109 +58,109 @@ DECLARE INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Prefactor = NumberDensity * pow(PI * Height * pow(R, 2), 2) * pow(DeltaRho, 2); + Prefactor = NumberDensity * pow (PI * Height * pow (R, 2), 2) * pow (DeltaRho, 2); - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Formfactor1; - double Formfactor2; - double Intensity; - double SolidAngle; - double qx; - double qy; - double qz; - double q; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - - // variables needed for integration over alpha - int i; - const int NumberOfSteps = 30; - double Alpha; - const double AlphaMin = 0; - const double AlphaMax = PI / 2.0; - const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfSteps); - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Compute scattering - Intensity = 0.0; - - for (i = 0; i < NumberOfSteps; ++i) { - Alpha = (i + 0.5) * AlphaStep; - - Formfactor1 = gsl_sf_bessel_J1(q * R * sin(Alpha)) / (q * R * sin(Alpha)); - Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); - - Intensity += sin(Alpha) * Prefactor * pow(2 * Formfactor1 * Formfactor2, 2) * AlphaStep; - } - - p *= l_full * SolidAngle / (4.0 * PI) * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Formfactor1; + double Formfactor2; + double Intensity; + double SolidAngle; + double qx; + double qy; + double qz; + double q; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + + // variables needed for integration over alpha + int i; + const int NumberOfSteps = 30; + double Alpha; + const double AlphaMin = 0; + const double AlphaMax = PI / 2.0; + const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfSteps); + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Compute scattering + Intensity = 0.0; + + for (i = 0; i < NumberOfSteps; ++i) { + Alpha = (i + 0.5) * AlphaStep; + + Formfactor1 = gsl_sf_bessel_J1 (q * R * sin (Alpha)) / (q * R * sin (Alpha)); + Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); + + Intensity += sin (Alpha) * Prefactor * pow (2 * Formfactor1 * Formfactor2, 2) * AlphaStep; + } + + p *= l_full * SolidAngle / (4.0 * PI) * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSEllipticCylinders.comp b/mcxtrace-comps/contrib/SAXSEllipticCylinders.comp index e990387c14..a21f97e048 100644 --- a/mcxtrace-comps/contrib/SAXSEllipticCylinders.comp +++ b/mcxtrace-comps/contrib/SAXSEllipticCylinders.comp @@ -60,121 +60,121 @@ DECLARE INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Prefactor = NumberDensity * pow(PI * Height * R1 * R2, 2) * pow(DeltaRho, 2); + Prefactor = NumberDensity * pow (PI * Height * R1 * R2, 2) * pow (DeltaRho, 2); - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; %} TRACE %{ - double l0; - double l1; - double l_full; - double l; - double l_1; - double Formfactor1; - double Formfactor2; - double Intensity; - double SolidAngle; - double qx; - double qy; - double qz; - double q; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - double ProjectedRadius; - char Intersect = 0; - - /* variables needed for integration over alpha */ - int i; - const int NumberOfStepsInAlpha = 30; - double Alpha; - const double AlphaMin = 0.0; - const double AlphaMax = PI / 2.0; - const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfStepsInAlpha); - - /* Variables needed in integration over beta */ - int j; - const int NumberOfStepsInBeta = 30; - double Beta; - const double BetaMin = 0.0; - const double BetaMax = PI / 2.0; - const double BetaStep = (BetaMax - BetaMin) / (1.0 * NumberOfStepsInBeta); - - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - /* Generate new direction of photon */ - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - /* Compute q */ - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - /* Compute scattering */ - Intensity = 0.0; - - for (i = 0; i < NumberOfStepsInAlpha; ++i) { - Alpha = (i + 0.5) * AlphaStep; - - for (j = 0; j < NumberOfStepsInBeta; ++j) { - Beta = (j + 0.5) * BetaStep; - ProjectedRadius = sqrt(pow(R1 * sin(Beta), 2) + pow(R2 * cos(Beta), 2)); - - Formfactor1 = gsl_sf_bessel_J1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); - Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); - - Intensity += 2 / PI * sin(Alpha) * Prefactor * pow(2 * Formfactor1 * Formfactor2, 2) * AlphaStep * BetaStep; - } - } - - p *= l_full * SolidAngle / (4.0 * PI) * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + double l0; + double l1; + double l_full; + double l; + double l_1; + double Formfactor1; + double Formfactor2; + double Intensity; + double SolidAngle; + double qx; + double qy; + double qz; + double q; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + double ProjectedRadius; + char Intersect = 0; + + /* variables needed for integration over alpha */ + int i; + const int NumberOfStepsInAlpha = 30; + double Alpha; + const double AlphaMin = 0.0; + const double AlphaMax = PI / 2.0; + const double AlphaStep = (AlphaMax - AlphaMin) / (1.0 * NumberOfStepsInAlpha); + + /* Variables needed in integration over beta */ + int j; + const int NumberOfStepsInBeta = 30; + double Beta; + const double BetaMin = 0.0; + const double BetaMax = PI / 2.0; + const double BetaStep = (BetaMax - BetaMin) / (1.0 * NumberOfStepsInBeta); + + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + /* Generate new direction of photon */ + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + /* Compute q */ + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + /* Compute scattering */ + Intensity = 0.0; + + for (i = 0; i < NumberOfStepsInAlpha; ++i) { + Alpha = (i + 0.5) * AlphaStep; + + for (j = 0; j < NumberOfStepsInBeta; ++j) { + Beta = (j + 0.5) * BetaStep; + ProjectedRadius = sqrt (pow (R1 * sin (Beta), 2) + pow (R2 * cos (Beta), 2)); + + Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha)); + Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); + + Intensity += 2 / PI * sin (Alpha) * Prefactor * pow (2 * Formfactor1 * Formfactor2, 2) * AlphaStep * BetaStep; + } + } + + p *= l_full * SolidAngle / (4.0 * PI) * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSLiposomes.comp b/mcxtrace-comps/contrib/SAXSLiposomes.comp index 86e541d9d6..357d34b2ab 100644 --- a/mcxtrace-comps/contrib/SAXSLiposomes.comp +++ b/mcxtrace-comps/contrib/SAXSLiposomes.comp @@ -57,50 +57,51 @@ SETTING PARAMETERS (Radius = 800.0, Thickness = 38.89, SigmaRadius = 0.20, nRadi SHARE %{ // Functions used for compution the intensity from a given liposome -#pragma acc routine - double FormfactorSphere(double q, double R) - { - return 3 * (sin(q * R) - q * R * cos(q * R)) / pow(q * R, 3); + #pragma acc routine + double + FormfactorSphere (double q, double R) { + return 3 * (sin (q * R) - q * R * cos (q * R)) / pow (q * R, 3); } -#pragma acc routine - double IntensityOfLiposome(double q, double R, double ThicknessHead, double ThicknessTail, double ThicknessCH3, double DeltaRhoHead, double DeltaRhoCH2, double DeltaRhoCH3) - { - double RHeadOut,RTailOut,RCH3Out; - double RCH3In,RTailIn,RHeadIn; - double VolumeHeadOut, VolumeTailOut, VolumeCH3,VolumeTailIn,VolumeHeadIn; - double AmplitudeHeadOut,AmplitudeTailOut,AmplitudeCH3,AmplitudeTailIn,AmplitudeHeadIn; + #pragma acc routine + double + IntensityOfLiposome (double q, double R, double ThicknessHead, double ThicknessTail, double ThicknessCH3, double DeltaRhoHead, double DeltaRhoCH2, + double DeltaRhoCH3) { + double RHeadOut, RTailOut, RCH3Out; + double RCH3In, RTailIn, RHeadIn; + double VolumeHeadOut, VolumeTailOut, VolumeCH3, VolumeTailIn, VolumeHeadIn; + double AmplitudeHeadOut, AmplitudeTailOut, AmplitudeCH3, AmplitudeTailIn, AmplitudeHeadIn; double Intensity; RHeadOut = R + ThicknessHead + ThicknessTail + ThicknessCH3; RTailOut = R + ThicknessTail + ThicknessCH3; - RCH3Out = R + ThicknessCH3; - RCH3In = R - ThicknessCH3; - RTailIn = R - ThicknessTail - ThicknessCH3; - RHeadIn = R - ThicknessHead - ThicknessTail - ThicknessCH3; + RCH3Out = R + ThicknessCH3; + RCH3In = R - ThicknessCH3; + RTailIn = R - ThicknessTail - ThicknessCH3; + RHeadIn = R - ThicknessHead - ThicknessTail - ThicknessCH3; - VolumeHeadOut = 4.0 / 3.0 * PI * (pow(RHeadOut, 3) - pow(RTailOut, 3)); - VolumeTailOut = 4.0 / 3.0 * PI * (pow(RTailOut, 3) - pow(RCH3Out, 3)); - VolumeCH3 = 4.0 / 3.0 * PI * (pow(RCH3Out, 3) - pow(RCH3In, 3)); - VolumeTailIn = 4.0 / 3.0 * PI * (pow(RCH3In, 3) - pow(RTailIn, 3)); - VolumeHeadIn = 4.0 / 3.0 * PI * (pow(RTailIn, 3) - pow(RHeadIn, 3)); + VolumeHeadOut = 4.0 / 3.0 * PI * (pow (RHeadOut, 3) - pow (RTailOut, 3)); + VolumeTailOut = 4.0 / 3.0 * PI * (pow (RTailOut, 3) - pow (RCH3Out, 3)); + VolumeCH3 = 4.0 / 3.0 * PI * (pow (RCH3Out, 3) - pow (RCH3In, 3)); + VolumeTailIn = 4.0 / 3.0 * PI * (pow (RCH3In, 3) - pow (RTailIn, 3)); + VolumeHeadIn = 4.0 / 3.0 * PI * (pow (RTailIn, 3) - pow (RHeadIn, 3)); - AmplitudeHeadOut = DeltaRhoHead * VolumeHeadOut * - (pow(RHeadOut, 3) * FormfactorSphere(q, RHeadOut) - pow(RTailOut, 3) * FormfactorSphere(q, RTailOut)) / (pow(RHeadOut, 3) - pow(RTailOut, 3)); + AmplitudeHeadOut = DeltaRhoHead * VolumeHeadOut * (pow (RHeadOut, 3) * FormfactorSphere (q, RHeadOut) - pow (RTailOut, 3) * FormfactorSphere (q, RTailOut)) + / (pow (RHeadOut, 3) - pow (RTailOut, 3)); - AmplitudeTailOut = DeltaRhoCH2 * VolumeTailOut * - (pow(RTailOut, 3) * FormfactorSphere(q, RTailOut) - pow(RCH3Out, 3) * FormfactorSphere(q, RCH3Out)) / (pow(RTailOut, 3) - pow(RCH3Out, 3)); + AmplitudeTailOut = DeltaRhoCH2 * VolumeTailOut * (pow (RTailOut, 3) * FormfactorSphere (q, RTailOut) - pow (RCH3Out, 3) * FormfactorSphere (q, RCH3Out)) + / (pow (RTailOut, 3) - pow (RCH3Out, 3)); - AmplitudeCH3 = DeltaRhoCH3 * VolumeCH3 * - (pow(RCH3Out, 3) * FormfactorSphere(q, RCH3Out) - pow(RCH3In, 3) * FormfactorSphere(q, RCH3In)) / (pow(RCH3Out, 3) - pow(RCH3In, 3)); + AmplitudeCH3 = DeltaRhoCH3 * VolumeCH3 * (pow (RCH3Out, 3) * FormfactorSphere (q, RCH3Out) - pow (RCH3In, 3) * FormfactorSphere (q, RCH3In)) + / (pow (RCH3Out, 3) - pow (RCH3In, 3)); - AmplitudeTailIn = DeltaRhoCH2 * VolumeTailIn * - (pow(RCH3In, 3) * FormfactorSphere(q, RCH3In) - pow(RTailIn, 3) * FormfactorSphere(q, RTailIn)) / (pow(RCH3In, 3) - pow(RTailIn, 3)); + AmplitudeTailIn = DeltaRhoCH2 * VolumeTailIn * (pow (RCH3In, 3) * FormfactorSphere (q, RCH3In) - pow (RTailIn, 3) * FormfactorSphere (q, RTailIn)) + / (pow (RCH3In, 3) - pow (RTailIn, 3)); - AmplitudeHeadIn = DeltaRhoHead * VolumeHeadIn * - (pow(RTailIn, 3) * FormfactorSphere(q, RTailIn) - pow(RHeadIn, 3) * FormfactorSphere(q, RHeadIn)) / (pow(RTailIn, 3) - pow(RHeadIn, 3)); + AmplitudeHeadIn = DeltaRhoHead * VolumeHeadIn * (pow (RTailIn, 3) * FormfactorSphere (q, RTailIn) - pow (RHeadIn, 3) * FormfactorSphere (q, RHeadIn)) + / (pow (RTailIn, 3) - pow (RHeadIn, 3)); - Intensity = pow(AmplitudeHeadOut + AmplitudeTailOut + AmplitudeCH3 + AmplitudeTailIn + AmplitudeHeadIn, 2); + Intensity = pow (AmplitudeHeadOut + AmplitudeTailOut + AmplitudeCH3 + AmplitudeTailIn + AmplitudeHeadIn, 2); return Intensity; } @@ -108,147 +109,146 @@ SHARE DECLARE %{ - double Absorption; - double NumberDensity; - - double RMin; - double RMax; - double RStep; - - // Scattering lengths - double DeltaRhoHead; - double DeltaRhoCH2Tail; - double DeltaRhoCH3Tail; - - // Thickness - double ThicknessOfHead; - double ThicknessOfCH2Tail; - double ThicknessOfCH3Tail; + double Absorption; + double NumberDensity; + + double RMin; + double RMax; + double RStep; + + // Scattering lengths + double DeltaRhoHead; + double DeltaRhoCH2Tail; + double DeltaRhoCH3Tail; + + // Thickness + double ThicknessOfHead; + double ThicknessOfCH2Tail; + double ThicknessOfCH3Tail; %} INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } - - Absorption = AbsorptionCrosssection; + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } + Absorption = AbsorptionCrosssection; - RMin = Radius - 3.0 * SigmaRadius * Radius; + RMin = Radius - 3.0 * SigmaRadius * Radius; - if (RMin < Thickness / 2.0) { - RMin = Thickness / 2.0; - } + if (RMin < Thickness / 2.0) { + RMin = Thickness / 2.0; + } - RMax = Radius + 3.0 * SigmaRadius * Radius; + RMax = Radius + 3.0 * SigmaRadius * Radius; - RStep = (RMax - RMin) / (1.0f * nRadius); + RStep = (RMax - RMin) / (1.0f * nRadius); - // Molecular properties of liposomes - const double ScatteringLengthOfWater = 2.82E-12; - const double VolumeOfWater = 30.0; + // Molecular properties of liposomes + const double ScatteringLengthOfWater = 2.82E-12; + const double VolumeOfWater = 30.0; - double RhoWater = ScatteringLengthOfWater / VolumeOfWater; - double RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; - double RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; - double RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; + double RhoWater = ScatteringLengthOfWater / VolumeOfWater; + double RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; + double RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; + double RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; - DeltaRhoHead = RhoHead - RhoWater; - DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; - DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; + DeltaRhoHead = RhoHead - RhoWater; + DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; + DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; - ThicknessOfHead = Thickness * VolumeOfHeadgroup / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); - ThicknessOfCH2Tail = Thickness * VolumeOfCH2Tail / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); - ThicknessOfCH3Tail = Thickness * VolumeOfCH3Tail / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); + ThicknessOfHead = Thickness * VolumeOfHeadgroup / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); + ThicknessOfCH2Tail = Thickness * VolumeOfCH2Tail / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); + ThicknessOfCH3Tail = Thickness * VolumeOfCH3Tail / (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail); %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Intensity; - double Weight1; - double Weight2; - double IntensityPart; - double SolidAngle; - double qx; - double qy; - double qz; - double k; - double q; - double dl; - double kx_i; - double ky_i; - double kz_i; - int Intersect = 0; - double R; - int i; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q=sqrt(qx*qx+qy*qy+qz*qz); - // Compute scattering - Intensity = 0.0; - Weight1 = 1.0 / (SigmaRadius * Radius * sqrt(2.0 * PI)); - for (i = 0; i < nRadius; i++) { - R = RMin + RStep* (i + 0.5); - IntensityPart = IntensityOfLiposome(q, R, ThicknessOfHead, ThicknessOfCH2Tail, ThicknessOfCH3Tail, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); - Weight2 = exp(- pow((R - Radius) / (sqrt(2.0) * SigmaRadius * Radius), 2)); - - Intensity += Weight1 * Weight2 * IntensityPart * RStep; - } - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Intensity; + double Weight1; + double Weight2; + double IntensityPart; + double SolidAngle; + double qx; + double qy; + double qz; + double k; + double q; + double dl; + double kx_i; + double ky_i; + double kz_i; + int Intersect = 0; + double R; + int i; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (qx * qx + qy * qy + qz * qz); + // Compute scattering + Intensity = 0.0; + Weight1 = 1.0 / (SigmaRadius * Radius * sqrt (2.0 * PI)); + for (i = 0; i < nRadius; i++) { + R = RMin + RStep * (i + 0.5); + IntensityPart = IntensityOfLiposome (q, R, ThicknessOfHead, ThicknessOfCH2Tail, ThicknessOfCH3Tail, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); + Weight2 = exp (-pow ((R - Radius) / (sqrt (2.0) * SigmaRadius * Radius), 2)); + + Intensity += Weight1 * Weight2 * IntensityPart * RStep; + } + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSNanodiscs.comp b/mcxtrace-comps/contrib/SAXSNanodiscs.comp index b4f8b539ee..a6b7c89eac 100644 --- a/mcxtrace-comps/contrib/SAXSNanodiscs.comp +++ b/mcxtrace-comps/contrib/SAXSNanodiscs.comp @@ -60,23 +60,23 @@ NOACC /*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/ SHARE %{ -#include + #include /* Functions used to compute the formfactor of a nanodisc */ - double ND_FormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) - { - double ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2)); + double + ND_FormfactorCylinder (double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) { + double ProjectedRadius = sqrt (pow (MajorSemiAxis * sin (Beta), 2) + pow (MinorSemiAxis * cos (Beta), 2)); - double Formfactor1 = gsl_sf_bessel_J1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); - double Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); + double Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha)); + double Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); return 2 * Formfactor1 * Formfactor2; } - double ND_IntensityOfEmptyNanodiscs(double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, - double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3, - double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail) - { + double + ND_IntensityOfEmptyNanodiscs (double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, double HeightOfBelt, double HeightOfLipids, + double HeightOfTails, double HeightOfCH3, double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, + double DeltaRhoCH3Tail) { double Intensity; double IntensityPart; @@ -93,9 +93,9 @@ SHARE double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt; double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt; - double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; - double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; - double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; + double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; + double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; + double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; // Variables needed for integration over alpha @@ -121,183 +121,187 @@ SHARE Alpha = (i + 0.5) * AlphaStep; for (j = 0; j < NumberOfStepsInBeta; ++j) { - Beta = (j + 0.5) * BetaStep; + Beta = (j + 0.5) * BetaStep; - // Compute formfactors - FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * ND_FormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) - - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) / - (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); + // Compute formfactors + FormfactorOfBelt + = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis + * ND_FormfactorCylinder (q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) + - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) + / (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); - FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) - - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) / - (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); + FormfactorOfHeads + = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) + - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta)) + / (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) - - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3 , Alpha, Beta)) / - (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); + FormfactorOfCH2Tail + = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) + - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta)) + / (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH3Tail = ND_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); + FormfactorOfCH3Tail = ND_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - // Compute amplitudes - AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; - AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; - AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; - AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; + // Compute amplitudes + AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; + AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; + AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; + AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; - // Perform integration - IntensityPart = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); + // Perform integration + IntensityPart = pow (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); - Intensity += 2 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep; + Intensity += 2 / PI * sin (Alpha) * IntensityPart * AlphaStep * BetaStep; } } return Intensity; } - %} DECLARE %{ - // Declarations - double Absorption; - double NumberDensity; - - // Scattering lengths - double DeltaRhoHead; - double DeltaRhoBelt; - double DeltaRhoCH2Tail; - double DeltaRhoCH3Tail; - - // Geometric properties - double MajorSemiAxis; - double MinorSemiAxis; - double ThicknessOfBelt; - - double HeightOfLipids; - double HeightOfTails; - double HeightOfCH3; + // Declarations + double Absorption; + double NumberDensity; + + // Scattering lengths + double DeltaRhoHead; + double DeltaRhoBelt; + double DeltaRhoCH2Tail; + double DeltaRhoCH3Tail; + + // Geometric properties + double MajorSemiAxis; + double MinorSemiAxis; + double ThicknessOfBelt; + + double HeightOfLipids; + double HeightOfTails; + double HeightOfCH3; %} INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; - // Scattering properties of different components - const double ScatteringLengthOfWater = 2.82E-12; - const double VolumeOfWater = 30.0; - // Scattering lengths - double RhoWater; - double RhoBelt; - double RhoHead; - double RhoCH2Tail; - double RhoCH3Tail; + // Scattering properties of different components + const double ScatteringLengthOfWater = 2.82E-12; + const double VolumeOfWater = 30.0; + // Scattering lengths + double RhoWater; + double RhoBelt; + double RhoHead; + double RhoCH2Tail; + double RhoCH3Tail; - RhoWater = ScatteringLengthOfWater / VolumeOfWater; - RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; - RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; - RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; - RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; + RhoWater = ScatteringLengthOfWater / VolumeOfWater; + RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; + RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; + RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; + RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; - DeltaRhoBelt = RhoBelt - RhoWater; - DeltaRhoHead = RhoHead - RhoWater; - DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; - DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; + DeltaRhoBelt = RhoBelt - RhoWater; + DeltaRhoHead = RhoHead - RhoWater; + DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; + DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; - // Geometric properties of different components - const double AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; + // Geometric properties of different components + const double AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; - MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio)); - MajorSemiAxis = MinorSemiAxis * AxisRatio; + MinorSemiAxis = sqrt (AreaOfLipids / (PI * AxisRatio)); + MajorSemiAxis = MinorSemiAxis * AxisRatio; - HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; + HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; - ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; + ThicknessOfBelt = sqrt (pow (MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double q; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Compute scattering - Intensity = exp(- pow(q * Roughness, 2))* ND_IntensityOfEmptyNanodiscs(q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, - HeightOfTails, HeightOfCH3, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double q; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Compute scattering + Intensity = exp (-pow (q * Roughness, 2)) + * ND_IntensityOfEmptyNanodiscs (q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3, + DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSNanodiscsFast.comp b/mcxtrace-comps/contrib/SAXSNanodiscsFast.comp index 7b8266b40f..c022c3a7f5 100644 --- a/mcxtrace-comps/contrib/SAXSNanodiscsFast.comp +++ b/mcxtrace-comps/contrib/SAXSNanodiscsFast.comp @@ -65,23 +65,23 @@ NOACC /*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/ SHARE %{ -#include + #include /* Functions used to compute the formfactor of a nanodisc */ - double NDF_FormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) - { - double ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2)); + double + NDF_FormfactorCylinder (double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) { + double ProjectedRadius = sqrt (pow (MajorSemiAxis * sin (Beta), 2) + pow (MinorSemiAxis * cos (Beta), 2)); - double Formfactor1 = gsl_sf_bessel_J1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); - double Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); + double Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha)); + double Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); return 2 * Formfactor1 * Formfactor2; } - double NDF_IntensityOfEmptyNanodiscs(double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, - double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3, - double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail) - { + double + NDF_IntensityOfEmptyNanodiscs (double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, double HeightOfBelt, double HeightOfLipids, + double HeightOfTails, double HeightOfCH3, double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, + double DeltaRhoCH3Tail) { // Declaration double Intensity; double IntensityPart; @@ -99,9 +99,9 @@ SHARE const double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt; const double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt; - const double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; - const double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; - const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; const double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; // Variables needed for integration over alpha @@ -127,219 +127,222 @@ SHARE Alpha = (i + 0.5) * AlphaStep; for (j = 0; j < NumberOfStepsInBeta; ++j) { - Beta = (j + 0.5) * BetaStep; + Beta = (j + 0.5) * BetaStep; - // Compute formfactors - FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * NDF_FormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) - - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) / - (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); + // Compute formfactors + FormfactorOfBelt + = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis + * NDF_FormfactorCylinder (q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) + - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) + / (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); - FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) - - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) / - (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); + FormfactorOfHeads + = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) + - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta)) + / (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) - - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3 , Alpha, Beta)) / - (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); + FormfactorOfCH2Tail + = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) + - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta)) + / (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH3Tail = NDF_FormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - - // Compute amplitudes - AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; - AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; - AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; - AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; + FormfactorOfCH3Tail = NDF_FormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - // Perform integration - IntensityPart = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); + // Compute amplitudes + AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; + AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; + AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; + AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; - Intensity += 2.0 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep; + // Perform integration + IntensityPart = pow (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); + + Intensity += 2.0 / PI * sin (Alpha) * IntensityPart * AlphaStep * BetaStep; } } return Intensity; } - %} DECLARE %{ - // Declarations - double RhoWater; - double RhoBelt; - double RhoHead; - double RhoCH2Tail; - double RhoCH3Tail; - - double DeltaRhoHead; - double DeltaRhoBelt; - double DeltaRhoCH2Tail; - double DeltaRhoCH3Tail; - - double MajorSemiAxis; - double MinorSemiAxis; - double ThicknessOfBelt; - - double HeightOfBelt; - double HeightOfLipids; - double HeightOfTails; - double HeightOfCH3; - - double AreaOfLipids; - double *qArray; - double *IArray; - - double NumberDensity; - double Absorption; + // Declarations + double RhoWater; + double RhoBelt; + double RhoHead; + double RhoCH2Tail; + double RhoCH3Tail; + + double DeltaRhoHead; + double DeltaRhoBelt; + double DeltaRhoCH2Tail; + double DeltaRhoCH3Tail; + + double MajorSemiAxis; + double MinorSemiAxis; + double ThicknessOfBelt; + + double HeightOfBelt; + double HeightOfLipids; + double HeightOfTails; + double HeightOfCH3; + + double AreaOfLipids; + double* qArray; + double* IArray; + + double NumberDensity; + double Absorption; %} INITIALIZE %{ - const double ScatteringLengthOfWater = 2.82E-12; - const double VolumeOfWater = 30.0; + const double ScatteringLengthOfWater = 2.82E-12; + const double VolumeOfWater = 30.0; - const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); - int i; - double qDummy; + const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); + int i; + double qDummy; - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; - // Scattering properties of different components - RhoWater = ScatteringLengthOfWater / VolumeOfWater; - RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; - RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; - RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; - RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; + // Scattering properties of different components + RhoWater = ScatteringLengthOfWater / VolumeOfWater; + RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; + RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; + RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; + RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; - DeltaRhoBelt = RhoBelt - RhoWater; - DeltaRhoHead = RhoHead - RhoWater; - DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; - DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; + DeltaRhoBelt = RhoBelt - RhoWater; + DeltaRhoHead = RhoHead - RhoWater; + DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; + DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; - // Geometric properties of different components - AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; + // Geometric properties of different components + AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; - MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio)); - MajorSemiAxis = MinorSemiAxis * AxisRatio; + MinorSemiAxis = sqrt (AreaOfLipids / (PI * AxisRatio)); + MajorSemiAxis = MinorSemiAxis * AxisRatio; - HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; + HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; - ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; + ThicknessOfBelt = sqrt (pow (MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; - // Compute scattering from nanodiscs in predecided points - qArray = (double *) calloc(NumberOfQBins, sizeof(double)); - IArray = (double *) calloc(NumberOfQBins, sizeof(double)); + // Compute scattering from nanodiscs in predecided points + qArray = (double*)calloc (NumberOfQBins, sizeof (double)); + IArray = (double*)calloc (NumberOfQBins, sizeof (double)); - for (i = 0; i < NumberOfQBins; ++i) { - qDummy = qMin + (i + 0.5) * qStep; + for (i = 0; i < NumberOfQBins; ++i) { + qDummy = qMin + (i + 0.5) * qStep; - qArray[i] = qDummy; - IArray[i] = NDF_IntensityOfEmptyNanodiscs(qDummy, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, - HeightOfTails, HeightOfCH3, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); - } + qArray[i] = qDummy; + IArray[i] = NDF_IntensityOfEmptyNanodiscs (qDummy, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3, + DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail); + } %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double q; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - int i; - double Slope; - double Offset; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Discard photon, if q is out of range - if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { - ABSORB; - } - - // Find the first value of q in the curve larger than that of the photon - i = 1; - - while (q > qArray[i]) { - ++i; - } - - // Do a linear interpolation - Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); - Offset = IArray[i] - Slope * qArray[i]; - - Intensity = (Slope * q + Offset) * exp(- pow(q * Roughness, 2)); - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double q; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + int i; + double Slope; + double Offset; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Discard photon, if q is out of range + if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { + ABSORB; + } + + // Find the first value of q in the curve larger than that of the photon + i = 1; + + while (q > qArray[i]) { + ++i; + } + + // Do a linear interpolation + Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); + Offset = IArray[i] - Slope * qArray[i]; + + Intensity = (Slope * q + Offset) * exp (-pow (q * Roughness, 2)); + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSNanodiscsWithTags.comp b/mcxtrace-comps/contrib/SAXSNanodiscsWithTags.comp index 961fa47782..63f6606387 100644 --- a/mcxtrace-comps/contrib/SAXSNanodiscsWithTags.comp +++ b/mcxtrace-comps/contrib/SAXSNanodiscsWithTags.comp @@ -65,22 +65,22 @@ NOACC /*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/ SHARE %{ -#include + #include /* Functions used for compution the intensity from a given liposome */ - double NDTFormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) - { - double ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2)); + double + NDTFormfactorCylinder (double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) { + double ProjectedRadius = sqrt (pow (MajorSemiAxis * sin (Beta), 2) + pow (MinorSemiAxis * cos (Beta), 2)); - double Formfactor1 = gsl_sf_bessel_J1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); - double Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); + double Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha)); + double Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); return 2 * Formfactor1 * Formfactor2; } // Functions used to compute the characteristics of the histidine tags - double NDTDebye(double q, double r) - { + double + NDTDebye (double q, double r) { // Variables needed in functions double ReturnValue; @@ -90,41 +90,40 @@ SHARE // Computation Dummy1 = (q * q) * (r * r); - Dummy2 = 2 * ((exp(-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1)); + Dummy2 = 2 * ((exp (-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1)); ReturnValue = Dummy2; return ReturnValue; } - - double NDTXiRodWithoutEndcaps(double q, double r, double l, double Alpha) - { + double + NDTXiRodWithoutEndcaps (double q, double r, double l, double Alpha) { // Variables needed in function double ReturnValue; // Computation - ReturnValue = gsl_sf_bessel_J0(q * r * sin(Alpha)) * sin(q * l * cos(Alpha) / 2.f) / (q * l * cos(Alpha) / 2.f); + ReturnValue = gsl_sf_bessel_J0 (q * r * sin (Alpha)) * sin (q * l * cos (Alpha) / 2.f) / (q * l * cos (Alpha) / 2.f); return ReturnValue; } - double NDTPsiHammouda(double x) - { + double + NDTPsiHammouda (double x) { // Variables used in function double ReturnValue; // Computation - ReturnValue = (1.0 - exp(-pow(x, 2))) / pow(x, 2); + ReturnValue = (1.0 - exp (-pow (x, 2))) / pow (x, 2); return ReturnValue; } /* Computation of the intensity from a given nanodisc */ - double NDT_IntensityOfEmptyNanodiscsWithTags(double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, - double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3, double RadiusOfGyrationForHisTag, double VolumeOfOneHisTag, - double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail, double DeltaRhoTag) - { + double + NDT_IntensityOfEmptyNanodiscsWithTags (double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, double HeightOfBelt, double HeightOfLipids, + double HeightOfTails, double HeightOfCH3, double RadiusOfGyrationForHisTag, double VolumeOfOneHisTag, + double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail, double DeltaRhoTag) { double Intensity; double IntensityPart; double IntensityDisc; @@ -141,11 +140,11 @@ SHARE double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt; double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt; - double AverageRadiusOfBelt = sqrt(pow(OuterMajorSemiAxis, 2) + pow(OuterMinorSemiAxis, 2)); + double AverageRadiusOfBelt = sqrt (pow (OuterMajorSemiAxis, 2) + pow (OuterMinorSemiAxis, 2)); - double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; - double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; - double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; + double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; + double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; + double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; double Autocorrelation; @@ -176,46 +175,50 @@ SHARE Alpha = (i + 0.5) * AlphaStep; for (j = 0; j < NumberOfStepsInBeta; ++j) { - Beta = (j + 0.5) * BetaStep; + Beta = (j + 0.5) * BetaStep; - // Compute formfactors - FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * NDTFormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) - - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) / - (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); + // Compute formfactors + FormfactorOfBelt + = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis + * NDTFormfactorCylinder (q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) + - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) + / (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); - FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) - - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) / - (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); + FormfactorOfHeads + = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) + - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta)) + / (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) - - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3 , Alpha, Beta)) / - (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); + FormfactorOfCH2Tail + = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) + - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta)) + / (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH3Tail = NDTFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); + FormfactorOfCH3Tail = NDTFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - // Compute amplitudes - AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; - AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; - AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; - AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; + // Compute amplitudes + AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; + AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; + AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; + AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; - // Perform integration - IntensityDisc = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); + // Perform integration + IntensityDisc = pow (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); - // Add the histidine tags - FormfactorOfTags = NDTXiRodWithoutEndcaps(q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha); + // Add the histidine tags + FormfactorOfTags = NDTXiRodWithoutEndcaps (q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha); - Crosscorrelation = 2.0 * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * pow(FormfactorOfTags, 2) * pow(NDTPsiHammouda(q * RadiusOfGyrationForHisTag), 2); + Crosscorrelation = 2.0 * pow (DeltaRhoTag * VolumeOfOneHisTag, 2) * pow (FormfactorOfTags, 2) * pow (NDTPsiHammouda (q * RadiusOfGyrationForHisTag), 2); - Autocorrelation = 2.0 * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * NDTDebye(q, RadiusOfGyrationForHisTag); + Autocorrelation = 2.0 * pow (DeltaRhoTag * VolumeOfOneHisTag, 2) * NDTDebye (q, RadiusOfGyrationForHisTag); - Disccorrelation = 4.0 * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) * - DeltaRhoTag * VolumeOfOneHisTag * FormfactorOfTags * NDTPsiHammouda(q * RadiusOfGyrationForHisTag); + Disccorrelation = 4.0 * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) * DeltaRhoTag * VolumeOfOneHisTag + * FormfactorOfTags * NDTPsiHammouda (q * RadiusOfGyrationForHisTag); - IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation; + IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation; - // Return the final intensity - Intensity += 2 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep; + // Return the final intensity + Intensity += 2 / PI * sin (Alpha) * IntensityPart * AlphaStep * BetaStep; } } @@ -225,149 +228,151 @@ SHARE DECLARE %{ - // Declarations - double RhoWater; - double RhoBelt; - double RhoHead; - double RhoCH2Tail; - double RhoCH3Tail; - double RhoTag; - - double DeltaRhoHead; - double DeltaRhoBelt; - double DeltaRhoCH2Tail; - double DeltaRhoCH3Tail; - double DeltaRhoTag; - - double MajorSemiAxis; - double MinorSemiAxis; - double ThicknessOfBelt; - - double HeightOfBelt; - double HeightOfLipids; - double HeightOfTails; - double HeightOfCH3; - - double NumberDensity; - double Absorption; + // Declarations + double RhoWater; + double RhoBelt; + double RhoHead; + double RhoCH2Tail; + double RhoCH3Tail; + double RhoTag; + + double DeltaRhoHead; + double DeltaRhoBelt; + double DeltaRhoCH2Tail; + double DeltaRhoCH3Tail; + double DeltaRhoTag; + + double MajorSemiAxis; + double MinorSemiAxis; + double ThicknessOfBelt; + + double HeightOfBelt; + double HeightOfLipids; + double HeightOfTails; + double HeightOfCH3; + + double NumberDensity; + double Absorption; %} INITIALIZE %{ - const double ScatteringLengthOfWater = 2.82E-12; - const double VolumeOfWater = 30.0; + const double ScatteringLengthOfWater = 2.82E-12; + const double VolumeOfWater = 30.0; - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; - // Scattering properties of different components - RhoWater = ScatteringLengthOfWater / VolumeOfWater; - RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; - RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; - RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; - RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; - RhoTag = ScatteringLengthOfOneHisTag / VolumeOfOneHisTag; + // Scattering properties of different components + RhoWater = ScatteringLengthOfWater / VolumeOfWater; + RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; + RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; + RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; + RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; + RhoTag = ScatteringLengthOfOneHisTag / VolumeOfOneHisTag; - DeltaRhoBelt = RhoBelt - RhoWater; - DeltaRhoHead = RhoHead - RhoWater; - DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; - DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; - DeltaRhoTag = RhoTag - RhoWater; + DeltaRhoBelt = RhoBelt - RhoWater; + DeltaRhoHead = RhoHead - RhoWater; + DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; + DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; + DeltaRhoTag = RhoTag - RhoWater; - // Geometric properties of different components - const double AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; + // Geometric properties of different components + const double AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; - MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio)); - MajorSemiAxis = MinorSemiAxis * AxisRatio; + MinorSemiAxis = sqrt (AreaOfLipids / (PI * AxisRatio)); + MajorSemiAxis = MinorSemiAxis * AxisRatio; - HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; + HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; - ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; + ThicknessOfBelt = sqrt (pow (MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double q; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - /* Compute properties of photon */ - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - /* Compute scattering */ - Intensity = exp(- pow(q * Roughness, 2)) * NDT_IntensityOfEmptyNanodiscsWithTags(q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3, - RadiusOfGyrationForHisTag, VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail, DeltaRhoTag); - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double q; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + /* Compute properties of photon */ + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + /* Compute scattering */ + Intensity = exp (-pow (q * Roughness, 2)) + * NDT_IntensityOfEmptyNanodiscsWithTags (q, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, + HeightOfCH3, RadiusOfGyrationForHisTag, VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, + DeltaRhoCH3Tail, DeltaRhoTag); + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSNanodiscsWithTagsFast.comp b/mcxtrace-comps/contrib/SAXSNanodiscsWithTagsFast.comp index 262644d4ee..daca7f78e0 100644 --- a/mcxtrace-comps/contrib/SAXSNanodiscsWithTagsFast.comp +++ b/mcxtrace-comps/contrib/SAXSNanodiscsWithTagsFast.comp @@ -69,22 +69,22 @@ NOACC /*X-ray Parameters (x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p)*/ SHARE %{ -#include + #include /* Functions used for compution the intensity from a given liposome */ - double NDTFFormfactorCylinder(double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) - { - double const ProjectedRadius = sqrt(pow(MajorSemiAxis * sin(Beta), 2) + pow(MinorSemiAxis * cos(Beta), 2)); + double + NDTFFormfactorCylinder (double q, double MajorSemiAxis, double MinorSemiAxis, double Height, double Alpha, double Beta) { + double const ProjectedRadius = sqrt (pow (MajorSemiAxis * sin (Beta), 2) + pow (MinorSemiAxis * cos (Beta), 2)); - double const Formfactor1 = gsl_sf_bessel_J1(q * ProjectedRadius * sin(Alpha)) / (q * ProjectedRadius * sin(Alpha)); - double const Formfactor2 = sin(q * Height * cos(Alpha) / 2.0) / (q * Height * cos(Alpha) / 2.0); + double const Formfactor1 = gsl_sf_bessel_J1 (q * ProjectedRadius * sin (Alpha)) / (q * ProjectedRadius * sin (Alpha)); + double const Formfactor2 = sin (q * Height * cos (Alpha) / 2.0) / (q * Height * cos (Alpha) / 2.0); return 2 * Formfactor1 * Formfactor2; } // Functions used to compute the characteristics of the histidine tags - double NDTFDebye(double q, double r) - { + double + NDTFDebye (double q, double r) { // Variables needed in functions double ReturnValue; @@ -94,41 +94,41 @@ SHARE // Computation Dummy1 = (q * q) * (r * r); - Dummy2 = 2 * ((exp(-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1)); + Dummy2 = 2 * ((exp (-Dummy1) + Dummy1 - 1) / (Dummy1 * Dummy1)); ReturnValue = Dummy2; return ReturnValue; } - - double NDTFXiRodWithoutEndcaps(double q, double r, double l, double Alpha) - { + double + NDTFXiRodWithoutEndcaps (double q, double r, double l, double Alpha) { // Variables needed in function double ReturnValue; // Computation - ReturnValue = gsl_sf_bessel_J0(q * r * sin(Alpha)) * sin(q * l * cos(Alpha) / 2.f) / (q * l * cos(Alpha) / 2.f); + ReturnValue = gsl_sf_bessel_J0 (q * r * sin (Alpha)) * sin (q * l * cos (Alpha) / 2.f) / (q * l * cos (Alpha) / 2.f); return ReturnValue; } - double NDTFPsiHammouda(double x) - { + double + NDTFPsiHammouda (double x) { // Variables used in function double ReturnValue; // Computation - ReturnValue = (1.0 - exp(-pow(x, 2))) / pow(x, 2); + ReturnValue = (1.0 - exp (-pow (x, 2))) / pow (x, 2); return ReturnValue; } /* Computation of the intensity from a given nanodisc */ - double NDTF_IntensityOfEmptyNanodiscsWithTags(double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, - double HeightOfBelt, double HeightOfLipids, double HeightOfTails, double HeightOfCH3, double RadiusOfGyrationForHisTag, double VolumeOfOneHisTag, - double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail, double DeltaRhoTag) - { + double + NDTF_IntensityOfEmptyNanodiscsWithTags (double q, double MajorSemiAxis, double MinorSemiAxis, double ThicknessOfBelt, double HeightOfBelt, + double HeightOfLipids, double HeightOfTails, double HeightOfCH3, double RadiusOfGyrationForHisTag, + double VolumeOfOneHisTag, double DeltaRhoBelt, double DeltaRhoHead, double DeltaRhoCH2Tail, double DeltaRhoCH3Tail, + double DeltaRhoTag) { double Intensity; double IntensityPart; double IntensityDisc; @@ -145,11 +145,11 @@ SHARE const double OuterMajorSemiAxis = MajorSemiAxis + ThicknessOfBelt; const double OuterMinorSemiAxis = MinorSemiAxis + ThicknessOfBelt; - const double AverageRadiusOfBelt = sqrt(pow(OuterMajorSemiAxis, 2) + pow(OuterMinorSemiAxis, 2)); + const double AverageRadiusOfBelt = sqrt (pow (OuterMajorSemiAxis, 2) + pow (OuterMinorSemiAxis, 2)); - const double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; - const double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; - const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfBelt = PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfHeads = PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis; + const double VolumeOfCH2Tail = PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; const double VolumeOfCH3Tail = PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis; double Autocorrelation; @@ -180,46 +180,50 @@ SHARE Alpha = (i + 0.5) * AlphaStep; for (j = 0; j < NumberOfStepsInBeta; ++j) { - Beta = (j + 0.5) * BetaStep; + Beta = (j + 0.5) * BetaStep; + + // Compute formfactors + FormfactorOfBelt + = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis + * NDTFFormfactorCylinder (q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) + - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) + / (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); - // Compute formfactors - FormfactorOfBelt = (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis * NDTFFormfactorCylinder(q, OuterMajorSemiAxis, OuterMinorSemiAxis, HeightOfBelt, Alpha, Beta) - - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfBelt, Alpha, Beta)) / - (PI * HeightOfBelt * OuterMajorSemiAxis * OuterMinorSemiAxis - PI * HeightOfBelt * MajorSemiAxis * MinorSemiAxis); + FormfactorOfHeads + = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) + - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta)) + / (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); - FormfactorOfHeads = (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfLipids, Alpha, Beta) - - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails , Alpha, Beta)) / - (PI * HeightOfLipids * MajorSemiAxis * MinorSemiAxis - PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis); + FormfactorOfCH2Tail + = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) + - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta)) + / (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); - FormfactorOfCH2Tail = (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfTails, Alpha, Beta) - - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis * NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3 , Alpha, Beta)) / - (PI * HeightOfTails * MajorSemiAxis * MinorSemiAxis - PI * HeightOfCH3 * MajorSemiAxis * MinorSemiAxis); + FormfactorOfCH3Tail = NDTFFormfactorCylinder (q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - FormfactorOfCH3Tail = NDTFFormfactorCylinder(q, MajorSemiAxis, MinorSemiAxis, HeightOfCH3, Alpha, Beta); - - // Compute amplitudes - AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; - AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; - AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; - AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; + // Compute amplitudes + AmplitudeOfBelt = DeltaRhoBelt * VolumeOfBelt * FormfactorOfBelt; + AmplitudeOfHeads = DeltaRhoHead * VolumeOfHeads * FormfactorOfHeads; + AmplitudeOfCH2Tail = DeltaRhoCH2Tail * VolumeOfCH2Tail * FormfactorOfCH2Tail; + AmplitudeOfCH3Tail = DeltaRhoCH3Tail * VolumeOfCH3Tail * FormfactorOfCH3Tail; - // Perform integration - IntensityDisc = pow(AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); + // Perform integration + IntensityDisc = pow (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail, 2); - // Add the histidine tags - FormfactorOfTags = NDTFXiRodWithoutEndcaps(q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha); + // Add the histidine tags + FormfactorOfTags = NDTFXiRodWithoutEndcaps (q, AverageRadiusOfBelt + RadiusOfGyrationForHisTag, HeightOfBelt, Alpha); - Crosscorrelation = 2.0 * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * pow(FormfactorOfTags, 2) * pow(NDTFPsiHammouda(q * RadiusOfGyrationForHisTag), 2); + Crosscorrelation = 2.0 * pow (DeltaRhoTag * VolumeOfOneHisTag, 2) * pow (FormfactorOfTags, 2) * pow (NDTFPsiHammouda (q * RadiusOfGyrationForHisTag), 2); - Autocorrelation = 2.0 * pow(DeltaRhoTag * VolumeOfOneHisTag, 2) * NDTFDebye(q, RadiusOfGyrationForHisTag); + Autocorrelation = 2.0 * pow (DeltaRhoTag * VolumeOfOneHisTag, 2) * NDTFDebye (q, RadiusOfGyrationForHisTag); - Disccorrelation = 4.0 * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) * - DeltaRhoTag * VolumeOfOneHisTag * FormfactorOfTags * NDTFPsiHammouda(q * RadiusOfGyrationForHisTag); + Disccorrelation = 4.0 * (AmplitudeOfBelt + AmplitudeOfHeads + AmplitudeOfCH2Tail + AmplitudeOfCH3Tail) * DeltaRhoTag * VolumeOfOneHisTag + * FormfactorOfTags * NDTFPsiHammouda (q * RadiusOfGyrationForHisTag); - IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation; + IntensityPart = IntensityDisc + Crosscorrelation + Autocorrelation + Disccorrelation; - // Return the final intensity - Intensity += 2 / PI * sin(Alpha) * IntensityPart * AlphaStep * BetaStep; + // Return the final intensity + Intensity += 2 / PI * sin (Alpha) * IntensityPart * AlphaStep * BetaStep; } } @@ -229,186 +233,187 @@ SHARE DECLARE %{ - // Declarations - double RhoWater; - double RhoBelt; - double RhoHead; - double RhoCH2Tail; - double RhoCH3Tail; - double RhoTag; - - double DeltaRhoHead; - double DeltaRhoBelt; - double DeltaRhoCH2Tail; - double DeltaRhoCH3Tail; - double DeltaRhoTag; - - double MajorSemiAxis; - double MinorSemiAxis; - double ThicknessOfBelt; - - double HeightOfBelt; - double HeightOfLipids; - double HeightOfTails; - double HeightOfCH3; - - double AreaOfLipids; - double *qArray; - double *IArray; - - double NumberDensity; - double Absorption; + // Declarations + double RhoWater; + double RhoBelt; + double RhoHead; + double RhoCH2Tail; + double RhoCH3Tail; + double RhoTag; + + double DeltaRhoHead; + double DeltaRhoBelt; + double DeltaRhoCH2Tail; + double DeltaRhoCH3Tail; + double DeltaRhoTag; + + double MajorSemiAxis; + double MinorSemiAxis; + double ThicknessOfBelt; + + double HeightOfBelt; + double HeightOfLipids; + double HeightOfTails; + double HeightOfCH3; + + double AreaOfLipids; + double* qArray; + double* IArray; + + double NumberDensity; + double Absorption; %} INITIALIZE %{ - const double ScatteringLengthOfWater = 2.82E-12; - const double VolumeOfWater = 30.0; + const double ScatteringLengthOfWater = 2.82E-12; + const double VolumeOfWater = 30.0; - const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); - int i; - double qDummy; + const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); + int i; + double qDummy; - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; - // Scattering properties of different components - RhoWater = ScatteringLengthOfWater / VolumeOfWater; - RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; - RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; - RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; - RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; - RhoTag = ScatteringLengthOfOneHisTag / VolumeOfOneHisTag; + // Scattering properties of different components + RhoWater = ScatteringLengthOfWater / VolumeOfWater; + RhoBelt = ScatteringLengthOfOneMSP / VolumeOfOneMSP; + RhoHead = ScatteringLengthOfHeadgroup / VolumeOfHeadgroup; + RhoCH2Tail = ScatteringLengthOfCH2Tail / VolumeOfCH2Tail; + RhoCH3Tail = ScatteringLengthOfCH3Tail / VolumeOfCH3Tail; + RhoTag = ScatteringLengthOfOneHisTag / VolumeOfOneHisTag; - DeltaRhoBelt = RhoBelt - RhoWater; - DeltaRhoHead = RhoHead - RhoWater; - DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; - DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; - DeltaRhoTag = RhoTag - RhoWater; + DeltaRhoBelt = RhoBelt - RhoWater; + DeltaRhoHead = RhoHead - RhoWater; + DeltaRhoCH2Tail = RhoCH2Tail - RhoWater; + DeltaRhoCH3Tail = RhoCH3Tail - RhoWater; + DeltaRhoTag = RhoTag - RhoWater; - // Geometric properties of different components - AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; + // Geometric properties of different components + AreaOfLipids = NumberOfLipids * AreaPerLipidHeadgroup / 2.0; - MinorSemiAxis = sqrt(AreaOfLipids / (PI * AxisRatio)); - MajorSemiAxis = MinorSemiAxis * AxisRatio; + MinorSemiAxis = sqrt (AreaOfLipids / (PI * AxisRatio)); + MajorSemiAxis = MinorSemiAxis * AxisRatio; - HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; - HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; + HeightOfLipids = 2.0 * (VolumeOfHeadgroup + VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfTails = 2.0 * (VolumeOfCH2Tail + VolumeOfCH3Tail) / AreaPerLipidHeadgroup; + HeightOfCH3 = 2.0 * VolumeOfCH3Tail / AreaPerLipidHeadgroup; - ThicknessOfBelt = sqrt(pow(MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; + ThicknessOfBelt = sqrt (pow (MinorSemiAxis + MajorSemiAxis, 2) / 4.0 + 2.0 * VolumeOfOneMSP / (PI * HeightOfMSP)) - (MajorSemiAxis + MinorSemiAxis) / 2.0; - // Compute scattering from nanodiscs in predecided points - qArray = (double *) calloc(NumberOfQBins, sizeof(double)); - IArray = (double *) calloc(NumberOfQBins, sizeof(double)); + // Compute scattering from nanodiscs in predecided points + qArray = (double*)calloc (NumberOfQBins, sizeof (double)); + IArray = (double*)calloc (NumberOfQBins, sizeof (double)); - for (i = 0; i < NumberOfQBins; ++i) { - qDummy = qMin + (i + 0.5) * qStep; + for (i = 0; i < NumberOfQBins; ++i) { + qDummy = qMin + (i + 0.5) * qStep; - qArray[i] = qDummy; - IArray[i] = NDTF_IntensityOfEmptyNanodiscsWithTags(qDummy, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, HeightOfCH3, RadiusOfGyrationForHisTag, - VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, DeltaRhoCH3Tail, DeltaRhoTag); - } + qArray[i] = qDummy; + IArray[i] = NDTF_IntensityOfEmptyNanodiscsWithTags (qDummy, MajorSemiAxis, MinorSemiAxis, ThicknessOfBelt, HeightOfMSP, HeightOfLipids, HeightOfTails, + HeightOfCH3, RadiusOfGyrationForHisTag, VolumeOfOneHisTag, DeltaRhoBelt, DeltaRhoHead, DeltaRhoCH2Tail, + DeltaRhoCH3Tail, DeltaRhoTag); + } %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double q; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - int i; - double Slope; - double Offset; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - /* Compute properties of photon */ - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Discard photon, if q is out of range - if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { - ABSORB; - } - - // Find the first value of q in the curve larger than that of the photon - i = 1; - - while (q > qArray[i]) { - ++i; - } - - // Do a linear interpolation - Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); - Offset = IArray[i] - Slope * qArray[i]; - - Intensity = (Slope * q + Offset) * exp(- pow(q * Roughness, 2)); - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double q; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + int i; + double Slope; + double Offset; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + /* Compute properties of photon */ + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Discard photon, if q is out of range + if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { + ABSORB; + } + + // Find the first value of q in the curve larger than that of the photon + i = 1; + + while (q > qArray[i]) { + ++i; + } + + // Do a linear interpolation + Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); + Offset = IArray[i] - Slope * qArray[i]; + + Intensity = (Slope * q + Offset) * exp (-pow (q * Roughness, 2)); + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSPDB.comp b/mcxtrace-comps/contrib/SAXSPDB.comp index 9b0eff2034..4f11ab0461 100644 --- a/mcxtrace-comps/contrib/SAXSPDB.comp +++ b/mcxtrace-comps/contrib/SAXSPDB.comp @@ -65,12 +65,12 @@ EXTEND %{ DECLARE %{ - double Absorption; - double NumberDensity; + double Absorption; + double NumberDensity; - // Protein properties - // cdouble Matrix[SAXSPDBOrderOfHarmonics+1][SAXSPDBOrderOfHarmonics+1]; - ProteinStruct Protein; + // Protein properties + // cdouble Matrix[SAXSPDBOrderOfHarmonics+1][SAXSPDBOrderOfHarmonics+1]; + ProteinStruct Protein; %} INITIALIZE @@ -79,154 +79,148 @@ INITIALIZE NumberDensity = Concentration * 6.02214129e19; // Standard sample handling - if (!xwidth || !yheight || !zdepth) { - exit(fprintf(stderr, "SAXSPDB: %s: ERROR: Sample has no volume - check parameters.\n", NAME_CURRENT_COMP)); - } + if (!xwidth || !yheight || !zdepth) { + exit (fprintf (stderr, "SAXSPDB: %s: ERROR: Sample has no volume - check parameters.\n", NAME_CURRENT_COMP)); + } // count the number of residues - Absorption = AbsorptionCrosssection; - Protein.NumberOfResidues = CountResidues(PDBFilepath); - Protein.Beads = calloc(Protein.NumberOfResidues,sizeof(BeadStruct)); - if (Protein.Beads == NULL) - exit(fprintf(stderr, "SAXSPDB: %s: ERROR: memory allocation\n", NAME_CURRENT_COMP)); + Absorption = AbsorptionCrosssection; + Protein.NumberOfResidues = CountResidues (PDBFilepath); + Protein.Beads = calloc (Protein.NumberOfResidues, sizeof (BeadStruct)); + if (Protein.Beads == NULL) + exit (fprintf (stderr, "SAXSPDB: %s: ERROR: memory allocation\n", NAME_CURRENT_COMP)); // initialize the protein from the PDB - ReadAminoPDB(PDBFilepath, &Protein); - MPI_MASTER( - printf("SAXSPDB: %s: Scattering from %s with %d residues\n", - NAME_CURRENT_COMP, PDBFilepath, Protein.NumberOfResidues); - ); - + ReadAminoPDB (PDBFilepath, &Protein); + MPI_MASTER (printf ("SAXSPDB: %s: Scattering from %s with %d residues\n", NAME_CURRENT_COMP, PDBFilepath, Protein.NumberOfResidues);); %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double q; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - double Slope; - double Offset; - int i,j,ResidueID; - cdouble** Matrix=(cdouble**)matrix2d_new(sizeof(cdouble),SAXSPDBOrderOfHarmonics+1,SAXSPDBOrderOfHarmonics+1); - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "SAXSPDB %s: Photon already inside sample - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Compute scattering for a given q-value - // ResetMatrix(Matrix, SAXSPDBOrderOfHarmonics); - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) - for (j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) - Matrix[i][j] = cplx(0,0); - - for (ResidueID = 0; ResidueID < Protein.NumberOfResidues; ++ResidueID) { - // ExpandStructure(Matrix, &Protein, ResidueID, qDummy, RhoSolvent); - double Legendre[SAXSPDBOrderOfHarmonics + 1]; - double Bessel[SAXSPDBOrderOfHarmonics + 1]; - - // Residue information - BeadStruct Residue = Protein.Beads[ResidueID]; - const double Volume = Residue.Volume; - const double DeltaRhoProtein = Residue.ScatteringLength - Volume * RhoSolvent; - - double X,Y,Z; - Protein_Beads_get_coords(Residue, &X, &Y, &Z); - X = (X * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].xv) / DeltaRhoProtein; - - Y = (Y * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].yv) / DeltaRhoProtein; - - Z = (Z * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].zv) / DeltaRhoProtein; - - // Convert bead position to spherical coordinates - const double Radius = sqrt(pow(X, 2) + pow(Y, 2) + pow(Z, 2)); - const double Theta = acos(Z / Radius); - const double C = acos(X / (Radius * sin(Theta))) * Sign(Y); - - // Expand protein structure on harmonics - gsl_sf_bessel_jl_array(SAXSPDBOrderOfHarmonics, q * Radius, Bessel); - - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { - gsl_sf_legendre_sphPlm_array(SAXSPDBOrderOfHarmonics, i, cos(Theta), &Legendre[i]); - - for(j = i; j <= SAXSPDBOrderOfHarmonics; ++j) { - Matrix[j][i] = cadd(Matrix[j][i],rmul(sqrt(4.0 * PI) * DeltaRhoProtein * Bessel[j] * Legendre[j] , cmul(cpow(cplx(0,1), cplx(j,0)), Polar(1.0, -i * C)))); - } - } - - } // for ResidueID - - Intensity = 0; - // ComputeIntensity(Matrix, SAXSPDBOrderOfHarmonics); - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { - for (j = 0; j <= i; ++j) { - Intensity += ((j > 0) + 1.0) * pow(cabs(Matrix[i][j]), 2); - } - } - - // Compute new weight - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double q; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + double Slope; + double Offset; + int i, j, ResidueID; + cdouble** Matrix = (cdouble**)matrix2d_new (sizeof (cdouble), SAXSPDBOrderOfHarmonics + 1, SAXSPDBOrderOfHarmonics + 1); + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "SAXSPDB %s: Photon already inside sample - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Compute scattering for a given q-value + // ResetMatrix(Matrix, SAXSPDBOrderOfHarmonics); + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) + for (j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) + Matrix[i][j] = cplx (0, 0); + + for (ResidueID = 0; ResidueID < Protein.NumberOfResidues; ++ResidueID) { + // ExpandStructure(Matrix, &Protein, ResidueID, qDummy, RhoSolvent); + double Legendre[SAXSPDBOrderOfHarmonics + 1]; + double Bessel[SAXSPDBOrderOfHarmonics + 1]; + + // Residue information + BeadStruct Residue = Protein.Beads[ResidueID]; + const double Volume = Residue.Volume; + const double DeltaRhoProtein = Residue.ScatteringLength - Volume * RhoSolvent; + + double X, Y, Z; + Protein_Beads_get_coords (Residue, &X, &Y, &Z); + X = (X * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].xv) / DeltaRhoProtein; + + Y = (Y * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].yv) / DeltaRhoProtein; + + Z = (Z * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].zv) / DeltaRhoProtein; + + // Convert bead position to spherical coordinates + const double Radius = sqrt (pow (X, 2) + pow (Y, 2) + pow (Z, 2)); + const double Theta = acos (Z / Radius); + const double C = acos (X / (Radius * sin (Theta))) * Sign (Y); + + // Expand protein structure on harmonics + gsl_sf_bessel_jl_array (SAXSPDBOrderOfHarmonics, q * Radius, Bessel); + + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { + gsl_sf_legendre_sphPlm_array (SAXSPDBOrderOfHarmonics, i, cos (Theta), &Legendre[i]); + + for (j = i; j <= SAXSPDBOrderOfHarmonics; ++j) { + Matrix[j][i] = cadd (Matrix[j][i], + rmul (sqrt (4.0 * PI) * DeltaRhoProtein * Bessel[j] * Legendre[j], cmul (cpow (cplx (0, 1), cplx (j, 0)), Polar (1.0, -i * C)))); + } + } + + } // for ResidueID + + Intensity = 0; + // ComputeIntensity(Matrix, SAXSPDBOrderOfHarmonics); + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { + for (j = 0; j <= i; ++j) { + Intensity += ((j > 0) + 1.0) * pow (cabs (Matrix[i][j]), 2); + } + } + + // Compute new weight + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSPDBFast.comp b/mcxtrace-comps/contrib/SAXSPDBFast.comp index a9a9d7b5cb..bd38019003 100644 --- a/mcxtrace-comps/contrib/SAXSPDBFast.comp +++ b/mcxtrace-comps/contrib/SAXSPDBFast.comp @@ -58,491 +58,484 @@ NOACC SHARE %{ %include "read_table-lib"; // for Open_File - - #include - #include - %include "mccode-complex-lib" -#ifndef SAXSPDB -#define SAXSPDB - #define SAXSPDBOrderOfHarmonics 21 - -void** matrix2d_new( size_t esize, size_t idim, size_t jdim ) { - size_t const ptrs_size = sizeof(void*) * idim; - size_t const row_size = esize * jdim; - void **const rows = malloc( ptrs_size ); - for ( size_t i = 0; i < idim; ++i ) - rows[i] = malloc( row_size ); - return rows; -} - // Simple mathematical functions - int Sign(double x) { - int Sign; - - if (x > 0) { - Sign = 1; - } else if (x < 0) { - Sign = -1; - } else { - Sign = 0; - } - - return Sign; - } - - void complex_print_matrix(cdouble **Matrix, int N, int M) - { - int i,j; - for (i = 0; i < N; ++i) - { - for (j = 0; j < M; ++j) - { - cdouble z = Matrix[i][j]; - fprintf(stderr, - "(%.12e,%.12e)%s", - creal(z), - cimag(z), - (j < M - 1) ? " " : "\n"); - } + + #include + #include + %include "mccode-complex-lib" + #ifndef SAXSPDB + #define SAXSPDB + #define SAXSPDBOrderOfHarmonics 21 + + void** + matrix2d_new (size_t esize, size_t idim, size_t jdim) { + size_t const ptrs_size = sizeof (void*) * idim; + size_t const row_size = esize * jdim; + void** const rows = malloc (ptrs_size); + for (size_t i = 0; i < idim; ++i) + rows[i] = malloc (row_size); + return rows; + } + // Simple mathematical functions + int + Sign (double x) { + int Sign; + + if (x > 0) { + Sign = 1; + } else if (x < 0) { + Sign = -1; + } else { + Sign = 0; } - } - - cdouble Polar(double R, double Concentration) { - cdouble Polar; - - Polar = cplx(R * cos(Concentration) , R* sin(Concentration)); - - return Polar; - } - - // Protein structs - struct Bead - { - double x; - double y; - double z; - - double xv; - double yv; - double zv; - - double Volume; - double ScatteringLength; - char Atom; - }; - typedef struct Bead BeadStruct; - - struct Protein - { - BeadStruct *Beads; - int NumberOfResidues; - }; - typedef struct Protein ProteinStruct; - - // functions for the INITIALIZE ---------------------------------------------- - - // Function used to determine the number of residues in the .pdb-file - int CountResidues(char *PDBFilepath) - { - // Declarations - double Dummy1; - double Dummy2; - double Dummy3; - char Line[65535]; - char DummyChar; - char Atom; - int NumberOfResidues = 0; - int ResidueID; - int PreviousResidueID = 0; - FILE *PDBFile; - - // I/O - PDBFile = Open_File(PDBFilepath, "r",NULL); - if (PDBFile == NULL) { - exit(fprintf(stderr, "SAXSPDBFast: %s: ERROR: Cannot open %s... \n", __FILE__, PDBFilepath)); - } - - while (fgets(Line, sizeof(Line), PDBFile) != NULL) { - ResidueID = 0; - if (strncmp(Line, "ATOM", 4)) continue; - if (sscanf(Line, "ATOM%*18c%d%*4c%lf%lf%lf", &ResidueID, &Dummy1, &Dummy2, &Dummy3) == 4) { - if (ResidueID != PreviousResidueID && ResidueID != 0) ++NumberOfResidues; - PreviousResidueID = ResidueID; - } - } - fclose(PDBFile); - return NumberOfResidues; - } // CountResidues - - // Function used to read .pdb-file - int ReadAminoPDB(char *PDBFilename, ProteinStruct *Protein) - { - // Declarations and input - int NumberOfResidues = Protein->NumberOfResidues; - BeadStruct *Residue = Protein->Beads; - FILE *PDBFile; - - int i = 0; - int PreviousResidueID = 0; - int ResidueID = 0; - - double Weight = 0.0; - double W = 0.0; - - double Aweight = 0.0; - double A = 0.0; - - double x; - double y; - double z; - - double X = 0.0; - double Y = 0.0; - double Z = 0.0; - - double XA = 0.0; - double YA = 0.0; - double ZA = 0.0; - - char Atom; - - char Buffer[65535]; - char DummyChar; - - // Atomic weighing factors - const double WH = 5.15; - const double WC = 16.44; - const double WN = 2.49; - const double WO = 9.13; - const double WS = 19.86; - const double WP = 5.73; - - // Scattering lengths - const double AH = 1 * 2.82e-13; - const double AD = 1 * 2.82e-13; - const double AC = 6 * 2.82e-13; - const double AN = 7 * 2.82e-13; - const double AO = 8 * 2.82e-13; - const double AP = 15 * 2.82e-13; - const double AS = 16 * 2.82e-13; - - // Program - if (NumberOfResidues <= 0 || (PDBFile = Open_File(PDBFilename, "r",NULL)) == 0) { - exit(printf("ERROR: Cannot open file: %s. \n", PDBFilename)); - } - - while (fgets(Buffer, sizeof(Buffer), PDBFile) != NULL) { - // a typical line is: - // ATOM 8726 N VAL B 576 76.450 47.214 58.026 1.00111.85 N - Atom = 0; - ResidueID = 0; - if (strncmp(Buffer, "ATOM", 4)) continue; - if (!sscanf(Buffer,"ATOM%*9c%c%*8c%d%*4c%lf%lf%lf%*23c%c", &DummyChar, &ResidueID, &x, &y, &z, &Atom)) { - fprintf(stderr, "SAXSPDBFast: %s: ReadAminoPDB: [%i] invalid PDB line %s\n", __FILE__, i, Buffer); - continue; - } - if (ResidueID != PreviousResidueID && ResidueID != 0) { + return Sign; + } - if (PreviousResidueID != 0 && Aweight && Weight) { + void + complex_print_matrix (cdouble** Matrix, int N, int M) { + int i, j; + for (i = 0; i < N; ++i) { + for (j = 0; j < M; ++j) { + cdouble z = Matrix[i][j]; + fprintf (stderr, "(%.12e,%.12e)%s", creal (z), cimag (z), (j < M - 1) ? " " : "\n"); + } + } + } - // Assign center of scattering - Residue[i].xv = X / Weight; - Residue[i].yv = Y / Weight; - Residue[i].zv = Z / Weight; + cdouble + Polar (double R, double Concentration) { + cdouble Polar; - // Assign center of mass - Residue[i].x = XA / Aweight; - Residue[i].y = YA / Aweight; - Residue[i].z = ZA / Aweight; + Polar = cplx (R * cos (Concentration), R * sin (Concentration)); - // Other residue attributes - Residue[i].Volume = Weight; - Residue[i].ScatteringLength = Aweight; - Residue[i].Atom = Atom; + return Polar; + } - X = Y = Z = Weight = 0.0; - XA = YA = ZA = Aweight = 0.0; + // Protein structs + struct Bead { + double x; + double y; + double z; + + double xv; + double yv; + double zv; + + double Volume; + double ScatteringLength; + char Atom; + }; + typedef struct Bead BeadStruct; + + struct Protein { + BeadStruct* Beads; + int NumberOfResidues; + }; + typedef struct Protein ProteinStruct; + + // functions for the INITIALIZE ---------------------------------------------- + + // Function used to determine the number of residues in the .pdb-file + int + CountResidues (char* PDBFilepath) { + // Declarations + double Dummy1; + double Dummy2; + double Dummy3; + char Line[65535]; + char DummyChar; + char Atom; + int NumberOfResidues = 0; + int ResidueID; + int PreviousResidueID = 0; + FILE* PDBFile; + + // I/O + PDBFile = Open_File (PDBFilepath, "r", NULL); + if (PDBFile == NULL) { + exit (fprintf (stderr, "SAXSPDBFast: %s: ERROR: Cannot open %s... \n", __FILE__, PDBFilepath)); + } - ++i; + while (fgets (Line, sizeof (Line), PDBFile) != NULL) { + ResidueID = 0; + if (strncmp (Line, "ATOM", 4)) + continue; + if (sscanf (Line, "ATOM%*18c%d%*4c%lf%lf%lf", &ResidueID, &Dummy1, &Dummy2, &Dummy3) == 4) { + if (ResidueID != PreviousResidueID && ResidueID != 0) + ++NumberOfResidues; + PreviousResidueID = ResidueID; + } + } + fclose (PDBFile); + return NumberOfResidues; + } // CountResidues + + // Function used to read .pdb-file + int + ReadAminoPDB (char* PDBFilename, ProteinStruct* Protein) { + // Declarations and input + int NumberOfResidues = Protein->NumberOfResidues; + BeadStruct* Residue = Protein->Beads; + FILE* PDBFile; + + int i = 0; + int PreviousResidueID = 0; + int ResidueID = 0; + + double Weight = 0.0; + double W = 0.0; + + double Aweight = 0.0; + double A = 0.0; + + double x; + double y; + double z; + + double X = 0.0; + double Y = 0.0; + double Z = 0.0; + + double XA = 0.0; + double YA = 0.0; + double ZA = 0.0; + + char Atom; + + char Buffer[65535]; + char DummyChar; + + // Atomic weighing factors + const double WH = 5.15; + const double WC = 16.44; + const double WN = 2.49; + const double WO = 9.13; + const double WS = 19.86; + const double WP = 5.73; + + // Scattering lengths + const double AH = 1 * 2.82e-13; + const double AD = 1 * 2.82e-13; + const double AC = 6 * 2.82e-13; + const double AN = 7 * 2.82e-13; + const double AO = 8 * 2.82e-13; + const double AP = 15 * 2.82e-13; + const double AS = 16 * 2.82e-13; + + // Program + if (NumberOfResidues <= 0 || (PDBFile = Open_File (PDBFilename, "r", NULL)) == 0) { + exit (printf ("ERROR: Cannot open file: %s. \n", PDBFilename)); + } - } + while (fgets (Buffer, sizeof (Buffer), PDBFile) != NULL) { + // a typical line is: + // ATOM 8726 N VAL B 576 76.450 47.214 58.026 1.00111.85 N + Atom = 0; + ResidueID = 0; + if (strncmp (Buffer, "ATOM", 4)) + continue; + if (!sscanf (Buffer, "ATOM%*9c%c%*8c%d%*4c%lf%lf%lf%*23c%c", &DummyChar, &ResidueID, &x, &y, &z, &Atom)) { + fprintf (stderr, "SAXSPDBFast: %s: ReadAminoPDB: [%i] invalid PDB line %s\n", __FILE__, i, Buffer); + continue; + } + + if (ResidueID != PreviousResidueID && ResidueID != 0) { + + if (PreviousResidueID != 0 && Aweight && Weight) { + + // Assign center of scattering + Residue[i].xv = X / Weight; + Residue[i].yv = Y / Weight; + Residue[i].zv = Z / Weight; + + // Assign center of mass + Residue[i].x = XA / Aweight; + Residue[i].y = YA / Aweight; + Residue[i].z = ZA / Aweight; + + // Other residue attributes + Residue[i].Volume = Weight; + Residue[i].ScatteringLength = Aweight; + Residue[i].Atom = Atom; + + X = Y = Z = Weight = 0.0; + XA = YA = ZA = Aweight = 0.0; + + ++i; + } - PreviousResidueID = ResidueID; - } + PreviousResidueID = ResidueID; + } + + // Finish the final amino acid + if (i == NumberOfResidues - 1 && Aweight && Weight) { + Residue[i].xv = X / Weight; + Residue[i].yv = Y / Weight; + Residue[i].zv = Z / Weight; + + // Assign center of mass + Residue[i].x = XA / Aweight; + Residue[i].y = YA / Aweight; + Residue[i].z = ZA / Aweight; + + // Other residue attributes + Residue[i].Volume = Weight; + Residue[i].ScatteringLength = Aweight; + Residue[i].Atom = 'X'; + } + + switch (Atom) { + case 'C': + A = AC; + W = WC; + break; + + case 'N': + A = AN; + W = WN; + break; + + case 'O': + A = AO; + W = WO; + break; + + case 'S': + A = AS; + W = WS; + break; + + case 'H': + A = AH; + W = WH; + break; + + case 'P': + A = AP; + W = WP; + break; + + default: + A = 0.0; + W = 0.0; + } + + Weight += W; + Aweight += A; + + X += W * x; + Y += W * y; + Z += W * z; + + XA += A * x; + YA += A * y; + ZA += A * z; + } - // Finish the final amino acid - if (i == NumberOfResidues - 1 && Aweight && Weight) { - Residue[i].xv = X / Weight; - Residue[i].yv = Y / Weight; - Residue[i].zv = Z / Weight; + fclose (PDBFile); - // Assign center of mass - Residue[i].x = XA / Aweight; - Residue[i].y = YA / Aweight; - Residue[i].z = ZA / Aweight; + return (NumberOfResidues); + } // ReadAminoPDB\ - // Other residue attributes - Residue[i].Volume = Weight; - Residue[i].ScatteringLength = Aweight; - Residue[i].Atom = 'X'; - } - switch(Atom) { - case 'C': - A = AC; - W = WC; - break; - - case 'N': - A = AN; - W = WN; - break; - - case 'O': - A = AO; - W = WO; - break; - - case 'S': - A = AS; - W = WS; - break; - - case 'H': - A = AH; - W = WH; - break; - - case 'P': - A = AP; - W = WP; - break; - - default: - A = 0.0; - W = 0.0; - } - - Weight += W; - Aweight += A; - - X += W * x; - Y += W * y; - Z += W * z; - - XA += A * x; - YA += A * y; - ZA += A * z; - } - - fclose(PDBFile); - - return(NumberOfResidues); - } // ReadAminoPDB\ - - -#endif /*SAXSPDB*/ + #endif /*SAXSPDB*/ %} DECLARE %{ - double Absorption; - double NumberDensity; + double Absorption; + double NumberDensity; - // Arrays for storing q and I(q) - DArray1d qArray; - DArray1d IArray; + // Arrays for storing q and I(q) + DArray1d qArray; + DArray1d IArray; %} INITIALIZE %{ - // Protein properties - ProteinStruct Protein; - int qbin; - - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; - - // Standard sample handling - if (!xwidth || !yheight || !zdepth) { - exit(fprintf(stderr, "SAXSPDBFast: %s: ERROR: Sample has no volume - check parameters.\n", NAME_CURRENT_COMP)); - } + // Protein properties + ProteinStruct Protein; + int qbin; + + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; + + // Standard sample handling + if (!xwidth || !yheight || !zdepth) { + exit (fprintf (stderr, "SAXSPDBFast: %s: ERROR: Sample has no volume - check parameters.\n", NAME_CURRENT_COMP)); + } // count the number of residues - Absorption = AbsorptionCrosssection; - Protein.NumberOfResidues = CountResidues(PDBFilepath); - Protein.Beads = calloc(Protein.NumberOfResidues,sizeof(BeadStruct)); - if (Protein.Beads == NULL) - exit(fprintf(stderr, "SAXSPDB: %s: ERROR: memory allocation\n", NAME_CURRENT_COMP)); + Absorption = AbsorptionCrosssection; + Protein.NumberOfResidues = CountResidues (PDBFilepath); + Protein.Beads = calloc (Protein.NumberOfResidues, sizeof (BeadStruct)); + if (Protein.Beads == NULL) + exit (fprintf (stderr, "SAXSPDB: %s: ERROR: memory allocation\n", NAME_CURRENT_COMP)); - qArray = create_darr1d(NumberOfQBins); - IArray = create_darr1d(NumberOfQBins); + qArray = create_darr1d (NumberOfQBins); + IArray = create_darr1d (NumberOfQBins); // initialize the protein from the PDB - ReadAminoPDB(PDBFilepath, &Protein); - MPI_MASTER( - printf("SAXSPDBFast: %s: Initializing scattering from %s with %d residues on %d Q-values\n", - NAME_CURRENT_COMP, PDBFilepath, Protein.NumberOfResidues, NumberOfQBins); - ); - - // Computing scattering profile I(q) - for (qbin = 0; qbin < NumberOfQBins; ++qbin) { - int i,j,ResidueID; - double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); - double q = qMin + qStep * (qbin + 0.5); - cdouble** Matrix = (cdouble**)matrix2d_new(sizeof(cdouble),SAXSPDBOrderOfHarmonics+1,SAXSPDBOrderOfHarmonics+1); - // init Matrix = 0 - // ResetMatrix(Matrix, SAXSPDBOrderOfHarmonics); - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) - for (j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) - Matrix[i][j] = cplx(0,0); - - - for (ResidueID = 0; ResidueID < Protein.NumberOfResidues; ++ResidueID) { - // ExpandStructure(Matrix, &Protein, ResidueID, qDummy, RhoSolvent); - double Legendre[SAXSPDBOrderOfHarmonics + 1]; - double Bessel[SAXSPDBOrderOfHarmonics + 1]; - - // Residue information - const double Volume = Protein.Beads[ResidueID].Volume; - const double DeltaRhoProtein = Protein.Beads[ResidueID].ScatteringLength - Volume * RhoSolvent; - - const double x = (Protein.Beads[ResidueID].x * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].xv) / DeltaRhoProtein; - - const double y = (Protein.Beads[ResidueID].y * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].yv) / DeltaRhoProtein; - - const double z = (Protein.Beads[ResidueID].z * Protein.Beads[ResidueID].ScatteringLength - - RhoSolvent * Volume * Protein.Beads[ResidueID].zv) / DeltaRhoProtein; - - // Convert bead position to spherical coordinates - const double Radius = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - const double Theta = acos(z / Radius); - const double C = acos(x / (Radius * sin(Theta))) * Sign(y); - - // Expand protein structure on harmonics - gsl_sf_bessel_jl_array(SAXSPDBOrderOfHarmonics, q * Radius, Bessel); - - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { - gsl_sf_legendre_sphPlm_array(SAXSPDBOrderOfHarmonics, i, cos(Theta), &Legendre[i]); - - for(j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) { - if (j < i) Matrix[j][i] = cplx(0,0); - else - Matrix[j][i] = cadd(Matrix[j][i],rmul(sqrt(4.0 * PI) * DeltaRhoProtein * Bessel[j] * Legendre[j] , cmul(cpow(cplx(0,1), cplx(j,0)), Polar(1.0, -i * C)))); - } - } - - } // for ResidueID - - qArray[qbin] = q; - // IArray[qbin] = ComputeIntensity(Matrix, SAXSPDBOrderOfHarmonics); - IArray[qbin] = 0; - for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { - for (j = 0; j <= i; ++j) { - IArray[qbin] += ((j > 0) + 1.0) * pow(cabs(Matrix[i][j]), 2); - } - } - // printf("I(q=%g) = %g\n", q, IArray[qbin]); - } // for qbin - MPI_MASTER( - printf("SAXSPDBFast: %s: %s initialization I(q) done\n", NAME_CURRENT_COMP, PDBFilepath); - ); + ReadAminoPDB (PDBFilepath, &Protein); + MPI_MASTER (printf ("SAXSPDBFast: %s: Initializing scattering from %s with %d residues on %d Q-values\n", NAME_CURRENT_COMP, PDBFilepath, + Protein.NumberOfResidues, NumberOfQBins);); + + // Computing scattering profile I(q) + for (qbin = 0; qbin < NumberOfQBins; ++qbin) { + int i, j, ResidueID; + double qStep = (qMax - qMin) / (1.0 * NumberOfQBins); + double q = qMin + qStep * (qbin + 0.5); + cdouble** Matrix = (cdouble**)matrix2d_new (sizeof (cdouble), SAXSPDBOrderOfHarmonics + 1, SAXSPDBOrderOfHarmonics + 1); + // init Matrix = 0 + // ResetMatrix(Matrix, SAXSPDBOrderOfHarmonics); + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) + for (j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) + Matrix[i][j] = cplx (0, 0); + + for (ResidueID = 0; ResidueID < Protein.NumberOfResidues; ++ResidueID) { + // ExpandStructure(Matrix, &Protein, ResidueID, qDummy, RhoSolvent); + double Legendre[SAXSPDBOrderOfHarmonics + 1]; + double Bessel[SAXSPDBOrderOfHarmonics + 1]; + + // Residue information + const double Volume = Protein.Beads[ResidueID].Volume; + const double DeltaRhoProtein = Protein.Beads[ResidueID].ScatteringLength - Volume * RhoSolvent; + + const double x + = (Protein.Beads[ResidueID].x * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].xv) / DeltaRhoProtein; + + const double y + = (Protein.Beads[ResidueID].y * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].yv) / DeltaRhoProtein; + + const double z + = (Protein.Beads[ResidueID].z * Protein.Beads[ResidueID].ScatteringLength - RhoSolvent * Volume * Protein.Beads[ResidueID].zv) / DeltaRhoProtein; + + // Convert bead position to spherical coordinates + const double Radius = sqrt (pow (x, 2) + pow (y, 2) + pow (z, 2)); + const double Theta = acos (z / Radius); + const double C = acos (x / (Radius * sin (Theta))) * Sign (y); + + // Expand protein structure on harmonics + gsl_sf_bessel_jl_array (SAXSPDBOrderOfHarmonics, q * Radius, Bessel); + + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { + gsl_sf_legendre_sphPlm_array (SAXSPDBOrderOfHarmonics, i, cos (Theta), &Legendre[i]); + + for (j = 0; j <= SAXSPDBOrderOfHarmonics; ++j) { + if (j < i) + Matrix[j][i] = cplx (0, 0); + else + Matrix[j][i] = cadd (Matrix[j][i], + rmul (sqrt (4.0 * PI) * DeltaRhoProtein * Bessel[j] * Legendre[j], cmul (cpow (cplx (0, 1), cplx (j, 0)), Polar (1.0, -i * C)))); + } + } + } // for ResidueID + + qArray[qbin] = q; + // IArray[qbin] = ComputeIntensity(Matrix, SAXSPDBOrderOfHarmonics); + IArray[qbin] = 0; + for (i = 0; i <= SAXSPDBOrderOfHarmonics; ++i) { + for (j = 0; j <= i; ++j) { + IArray[qbin] += ((j > 0) + 1.0) * pow (cabs (Matrix[i][j]), 2); + } + } + // printf("I(q=%g) = %g\n", q, IArray[qbin]); + } // for qbin + MPI_MASTER (printf ("SAXSPDBFast: %s: %s initialization I(q) done\n", NAME_CURRENT_COMP, PDBFilepath);); %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double q; - double Intensity; - double Weight; - double IntensityPart; - double SolidAngle; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - double Slope; - double Offset; - int i; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "SAXSPDBFast: %s: Photon already inside sample - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Discard photon, if q is out of range - if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { - ABSORB; - } - - // Find the first value of q in the curve larger than that of the photon - i = 1; - - while (q > qArray[i]) { - ++i; - } - - // Do a linear interpolation - Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); - Offset = IArray[i] - Slope * qArray[i]; - - Intensity = (Slope * q + Offset); - - p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1)); - - SCATTER; + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double q; + double Intensity; + double Weight; + double IntensityPart; + double SolidAngle; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + double Slope; + double Offset; + int i; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "SAXSPDBFast: %s: Photon already inside sample - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Discard photon, if q is out of range + if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) { + ABSORB; + } + + // Find the first value of q in the curve larger than that of the photon + i = 1; + + while (q > qArray[i]) { + ++i; + } + + // Do a linear interpolation + Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]); + Offset = IArray[i] - Slope * qArray[i]; + + Intensity = (Slope * q + Offset); + + p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp (-Absorption * (l + l1)); + + SCATTER; } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSQMonitor.comp b/mcxtrace-comps/contrib/SAXSQMonitor.comp index 937f8686e9..2fd72eac79 100644 --- a/mcxtrace-comps/contrib/SAXSQMonitor.comp +++ b/mcxtrace-comps/contrib/SAXSQMonitor.comp @@ -50,170 +50,149 @@ SETTING PARAMETERS (string RFilename="RDetector", string qFilename="QDetector", DECLARE %{ - DArray1d Nofq; - DArray1d Iofq; - DArray1d IofqSquared; + DArray1d Nofq; + DArray1d Iofq; + DArray1d IofqSquared; - DArray1d NofR; - DArray1d IofR; - DArray1d IofRSquared; + DArray1d NofR; + DArray1d IofR; + DArray1d IofRSquared; %} INITIALIZE %{ - Nofq=create_darr1d(NumberOfBins); - Iofq=create_darr1d(NumberOfBins); - IofqSquared=create_darr1d(NumberOfBins); - NofR=create_darr1d(NumberOfBins); - IofR=create_darr1d(NumberOfBins); - IofRSquared=create_darr1d(NumberOfBins); - - if (DistanceFromSample <= 0) - exit(fprintf(stderr,"%s: ERROR: you must set the DistanceFromSample > 0.\n", NAME_CURRENT_COMP)); - - if (qMax<=0 && LambdaMin) { - double TwoThetaMax = atan(RadiusDetector / DistanceFromSample); - qMax = 4 * PI * sin(TwoThetaMax / 2.0) / LambdaMin; - } - - // Use instance name for monitor output if no input was given - if (!strcmp(RFilename,"\0")) sprintf(RFilename,"%s_R", NAME_CURRENT_COMP); - if (!strcmp(qFilename,"\0")) sprintf(qFilename,"%s_Q", NAME_CURRENT_COMP); + Nofq = create_darr1d (NumberOfBins); + Iofq = create_darr1d (NumberOfBins); + IofqSquared = create_darr1d (NumberOfBins); + NofR = create_darr1d (NumberOfBins); + IofR = create_darr1d (NumberOfBins); + IofRSquared = create_darr1d (NumberOfBins); + + if (DistanceFromSample <= 0) + exit (fprintf (stderr, "%s: ERROR: you must set the DistanceFromSample > 0.\n", NAME_CURRENT_COMP)); + + if (qMax <= 0 && LambdaMin) { + double TwoThetaMax = atan (RadiusDetector / DistanceFromSample); + qMax = 4 * PI * sin (TwoThetaMax / 2.0) / LambdaMin; + } + + // Use instance name for monitor output if no input was given + if (!strcmp (RFilename, "\0")) + sprintf (RFilename, "%s_R", NAME_CURRENT_COMP); + if (!strcmp (qFilename, "\0")) + sprintf (qFilename, "%s_Q", NAME_CURRENT_COMP); %} TRACE %{ - int i; - double TwoTheta; - double Lambda; - - double R; - double RLow; - double RHigh; - - double q; - double qLow; - double qHigh; - - double TwoThetaLow; - double TwoThetaHigh; - double AreaOfSlice; - - PROP_Z0; - - R = sqrt(x*x+y*y); - - // Computation of q - if (Lambda0 <= 0.0) { - Lambda = 2.0 * PI / sqrt(kx*kx+ky*ky+kz*kz); - } else { - Lambda = Lambda0; - } - - TwoTheta = atan(R / DistanceFromSample); - q = 4.0 * PI * sin(TwoTheta / 2.0) / Lambda; - - // Put photon in the correct r-bin - if (R < RadiusDetector) { - i = floor(NumberOfBins * R / RadiusDetector); - - RLow = RadiusDetector / NumberOfBins * i; - RHigh = RadiusDetector / NumberOfBins * (i + 1); - - TwoThetaLow = atan(RLow / DistanceFromSample); - TwoThetaHigh = atan(RHigh / DistanceFromSample); - - AreaOfSlice = fabs((cos(2.0 * TwoThetaLow) - cos(2.0 * TwoThetaHigh)) * 2.0 * PI); - -#pragma acc atomic - NofR[i] = NofR[i] + 1; - double p_A=p/AreaOfSlice; -#pragma acc atomic - IofR[i] += p_A; - double p2_A2= p*p/(AreaOfSlice*AreaOfSlice); -#pragma acc atomic - IofRSquared[i] += p2_A2; - } - - // Put photon in the correct q-bin - if (q < qMax) { - i = floor(NumberOfBins * q / qMax); - - qLow = qMax / NumberOfBins * i; - qHigh = qMax / NumberOfBins * (i + 1); - - TwoThetaLow = asin(qLow * Lambda / (4.0 * PI)); - TwoThetaHigh = asin(qHigh * Lambda / (4.0 * PI)); - - AreaOfSlice = fabs((cos(2.0 * TwoThetaLow) - cos(2.0 * TwoThetaHigh)) * 2.0 * PI); - -#pragma acc atomic - Nofq[i] = Nofq[i] + 1; - double p_A=p/AreaOfSlice; -#pragma acc atomic - Iofq[i] += p_A; - double p2_A2= p*p/(AreaOfSlice*AreaOfSlice); -#pragma acc atomic - IofqSquared[i] += p2_A2; - - SCATTER; - } - - // Restore xray if requested - if (restore_xray) { - RESTORE_XRAY(INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); - } + int i; + double TwoTheta; + double Lambda; + + double R; + double RLow; + double RHigh; + + double q; + double qLow; + double qHigh; + + double TwoThetaLow; + double TwoThetaHigh; + double AreaOfSlice; + + PROP_Z0; + + R = sqrt (x * x + y * y); + + // Computation of q + if (Lambda0 <= 0.0) { + Lambda = 2.0 * PI / sqrt (kx * kx + ky * ky + kz * kz); + } else { + Lambda = Lambda0; + } + + TwoTheta = atan (R / DistanceFromSample); + q = 4.0 * PI * sin (TwoTheta / 2.0) / Lambda; + + // Put photon in the correct r-bin + if (R < RadiusDetector) { + i = floor (NumberOfBins * R / RadiusDetector); + + RLow = RadiusDetector / NumberOfBins * i; + RHigh = RadiusDetector / NumberOfBins * (i + 1); + + TwoThetaLow = atan (RLow / DistanceFromSample); + TwoThetaHigh = atan (RHigh / DistanceFromSample); + + AreaOfSlice = fabs ((cos (2.0 * TwoThetaLow) - cos (2.0 * TwoThetaHigh)) * 2.0 * PI); + + #pragma acc atomic + NofR[i] = NofR[i] + 1; + double p_A = p / AreaOfSlice; + #pragma acc atomic + IofR[i] += p_A; + double p2_A2 = p * p / (AreaOfSlice * AreaOfSlice); + #pragma acc atomic + IofRSquared[i] += p2_A2; + } + + // Put photon in the correct q-bin + if (q < qMax) { + i = floor (NumberOfBins * q / qMax); + + qLow = qMax / NumberOfBins * i; + qHigh = qMax / NumberOfBins * (i + 1); + + TwoThetaLow = asin (qLow * Lambda / (4.0 * PI)); + TwoThetaHigh = asin (qHigh * Lambda / (4.0 * PI)); + + AreaOfSlice = fabs ((cos (2.0 * TwoThetaLow) - cos (2.0 * TwoThetaHigh)) * 2.0 * PI); + + #pragma acc atomic + Nofq[i] = Nofq[i] + 1; + double p_A = p / AreaOfSlice; + #pragma acc atomic + Iofq[i] += p_A; + double p2_A2 = p * p / (AreaOfSlice * AreaOfSlice); + #pragma acc atomic + IofqSquared[i] += p2_A2; + + SCATTER; + } + + // Restore xray if requested + if (restore_xray) { + RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); + } %} SAVE %{ // Output I(r) - DETECTOR_OUT_1D( - "QMonitor - Radially averaged distribution", - "Radius [m]", - "I(r)", - "r", - 0.0, - RadiusDetector, - NumberOfBins, - NofR, - IofR, - IofRSquared, - RFilename - ); - - // Output I(q) - char fname1[256]; - snprintf(fname1, 256, "%s_q", NAME_CURRENT_COMP); - // we use mcdetector_out_nD in order to update the [comp_name]_I symbol for the 2nd output - mcdetector_out_1D( - "QMonitor - Distribution in q (Radially averaged)", - "q [1 / AA]", - "I(q)", - "q", - 0.0, - qMax, - NumberOfBins, - Nofq, - Iofq, - IofqSquared, - qFilename, fname1,POS_A_CURRENT_COMP,ROT_A_CURRENT_COMP,INDEX_CURRENT_COMP - ); + DETECTOR_OUT_1D ("QMonitor - Radially averaged distribution", "Radius [m]", "I(r)", "r", 0.0, RadiusDetector, NumberOfBins, NofR, IofR, IofRSquared, RFilename); + + // Output I(q) + char fname1[256]; + snprintf (fname1, 256, "%s_q", NAME_CURRENT_COMP); + // we use mcdetector_out_nD in order to update the [comp_name]_I symbol for the 2nd output + mcdetector_out_1D ("QMonitor - Distribution in q (Radially averaged)", "q [1 / AA]", "I(q)", "q", 0.0, qMax, NumberOfBins, Nofq, Iofq, IofqSquared, qFilename, + fname1, POS_A_CURRENT_COMP, ROT_A_CURRENT_COMP, INDEX_CURRENT_COMP); %} FINALLY %{ - destroy_darr1d(Nofq); - destroy_darr1d(Iofq); - destroy_darr1d(IofqSquared); - destroy_darr1d(NofR); - destroy_darr1d(IofR); - destroy_darr1d(IofRSquared); + destroy_darr1d (Nofq); + destroy_darr1d (Iofq); + destroy_darr1d (IofqSquared); + destroy_darr1d (NofR); + destroy_darr1d (IofR); + destroy_darr1d (IofRSquared); %} MCDISPLAY %{ - circle("xy", 0, 0, 0, RadiusDetector); + circle ("xy", 0, 0, 0, RadiusDetector); %} END diff --git a/mcxtrace-comps/contrib/SAXSShells.comp b/mcxtrace-comps/contrib/SAXSShells.comp index 83cac3bac5..c7987cdd3e 100644 --- a/mcxtrace-comps/contrib/SAXSShells.comp +++ b/mcxtrace-comps/contrib/SAXSShells.comp @@ -48,119 +48,119 @@ SETTING PARAMETERS (R = 100.0, Thickness = 5.0, Concentration = 0.01, DeltaRho = DECLARE %{ - double Prefactor; - double Absorption; - double q; - double NumberDensity; + double Prefactor; + double Absorption; + double q; + double NumberDensity; - double RBig; - double RSmall; + double RBig; + double RSmall; - double VolumeBigSphere; - double VolumeSmallSphere; - double Volume; + double VolumeBigSphere; + double VolumeSmallSphere; + double Volume; %} INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - if (Thickness >= R) { - printf("%s: Thickness of shell larger than radius of shell!\n", NAME_CURRENT_COMP); - } + if (Thickness >= R) { + printf ("%s: Thickness of shell larger than radius of shell!\n", NAME_CURRENT_COMP); + } - RBig = R + Thickness / 2.0; - RSmall = R - Thickness / 2.0; + RBig = R + Thickness / 2.0; + RSmall = R - Thickness / 2.0; - VolumeBigSphere = 4.0 / 3.0 * PI * pow(RBig, 3); - VolumeSmallSphere = 4.0 / 3.0 * PI * pow(RSmall, 3); + VolumeBigSphere = 4.0 / 3.0 * PI * pow (RBig, 3); + VolumeSmallSphere = 4.0 / 3.0 * PI * pow (RSmall, 3); - Volume = VolumeBigSphere - VolumeSmallSphere; + Volume = VolumeBigSphere - VolumeSmallSphere; - Prefactor = NumberDensity * pow(Volume, 2) * pow(DeltaRho, 2); + Prefactor = NumberDensity * pow (Volume, 2) * pow (DeltaRho, 2); - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double FormfactorBigSphere; - double FormfactorSmallSphere; - double Formfactor; - double SolidAngle; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Compute scattering - FormfactorBigSphere = 3.0 * (sin(q * RBig) - q * RBig * cos(q * RBig)) / pow(q * RBig, 3); - FormfactorSmallSphere = 3.0 * (sin(q * RSmall) - q * RSmall * cos(q * RSmall)) / pow(q * RSmall, 3); - Formfactor = (FormfactorBigSphere * VolumeBigSphere - FormfactorSmallSphere * VolumeSmallSphere) / (VolumeBigSphere - VolumeSmallSphere); - - p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * pow(Formfactor, 2) * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double FormfactorBigSphere; + double FormfactorSmallSphere; + double Formfactor; + double SolidAngle; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Compute scattering + FormfactorBigSphere = 3.0 * (sin (q * RBig) - q * RBig * cos (q * RBig)) / pow (q * RBig, 3); + FormfactorSmallSphere = 3.0 * (sin (q * RSmall) - q * RSmall * cos (q * RSmall)) / pow (q * RSmall, 3); + Formfactor = (FormfactorBigSphere * VolumeBigSphere - FormfactorSmallSphere * VolumeSmallSphere) / (VolumeBigSphere - VolumeSmallSphere); + + p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * pow (Formfactor, 2) * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END diff --git a/mcxtrace-comps/contrib/SAXSSpheres.comp b/mcxtrace-comps/contrib/SAXSSpheres.comp index 142743f499..c96f9bfdd2 100644 --- a/mcxtrace-comps/contrib/SAXSSpheres.comp +++ b/mcxtrace-comps/contrib/SAXSSpheres.comp @@ -47,97 +47,97 @@ SETTING PARAMETERS (R = 100.0, Concentration = 0.01, DeltaRho = 1.0e-14, Absorpt DECLARE %{ - double Prefactor; - double Absorption; - double q; - double NumberDensity; + double Prefactor; + double Absorption; + double q; + double NumberDensity; %} INITIALIZE %{ - // Rescale concentration into number of aggregates per m^3 times 10^-4 - NumberDensity = Concentration * 6.02214129e19; + // Rescale concentration into number of aggregates per m^3 times 10^-4 + NumberDensity = Concentration * 6.02214129e19; - // Computations - if (!xwidth || !yheight || !zdepth) { - printf("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); - } + // Computations + if (!xwidth || !yheight || !zdepth) { + printf ("%s: Sample has no volume, check parameters!\n", NAME_CURRENT_COMP); + } - Prefactor = NumberDensity * pow(4.0 / 3.0 * PI * pow(R, 3), 2) * pow(DeltaRho, 2); + Prefactor = NumberDensity * pow (4.0 / 3.0 * PI * pow (R, 3), 2) * pow (DeltaRho, 2); - Absorption = AbsorptionCrosssection; + Absorption = AbsorptionCrosssection; %} TRACE %{ - // Declarations - double l0; - double l1; - double l_full; - double l; - double l_1; - double Formfactor; - double SolidAngle; - double qx; - double qy; - double qz; - double k; - double dl; - double kx_i; - double ky_i; - double kz_i; - char Intersect = 0; - - // Computation - Intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); - - if (Intersect) { - - if (l0 < 0.0) { - fprintf(stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); - ABSORB; - } - - // Compute properties of photon - k = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2)); - l_full = l1 - l0; - dl = rand01() * (l1 - l0) + l0; - PROP_DL(dl); - l = dl - l0; - - // Store properties of incoming photon - kx_i = kx; - ky_i = ky; - kz_i = kz; - - // Generate new direction of photon - randvec_target_circle(&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); - - NORM(kx, ky, kz); - - kx *= k; - ky *= k; - kz *= k; - - // Compute q - qx = kx_i - kx; - qy = ky_i - ky; - qz = kz_i - kz; - - q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2)); - - // Compute scattering - Formfactor = 3.0 * (sin(q * R) - q * R * cos(q * R)) / pow(q * R, 3); - - p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * pow(Formfactor, 2) * exp(- Absorption * (l + l1)); - - SCATTER; - } + // Declarations + double l0; + double l1; + double l_full; + double l; + double l_1; + double Formfactor; + double SolidAngle; + double qx; + double qy; + double qz; + double k; + double dl; + double kx_i; + double ky_i; + double kz_i; + char Intersect = 0; + + // Computation + Intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); + + if (Intersect) { + + if (l0 < 0.0) { + fprintf (stderr, "Photon already inside sample %s - absorbing...\n", NAME_CURRENT_COMP); + ABSORB; + } + + // Compute properties of photon + k = sqrt (pow (kx, 2) + pow (ky, 2) + pow (kz, 2)); + l_full = l1 - l0; + dl = rand01 () * (l1 - l0) + l0; + PROP_DL (dl); + l = dl - l0; + + // Store properties of incoming photon + kx_i = kx; + ky_i = ky; + kz_i = kz; + + // Generate new direction of photon + randvec_target_circle (&kx, &ky, &kz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius); + + NORM (kx, ky, kz); + + kx *= k; + ky *= k; + kz *= k; + + // Compute q + qx = kx_i - kx; + qy = ky_i - ky; + qz = kz_i - kz; + + q = sqrt (pow (qx, 2) + pow (qy, 2) + pow (qz, 2)); + + // Compute scattering + Formfactor = 3.0 * (sin (q * R) - q * R * cos (q * R)) / pow (q * R, 3); + + p *= l_full * SolidAngle / (4.0 * PI) * Prefactor * pow (Formfactor, 2) * exp (-Absorption * (l + l1)); + + SCATTER; + } %} MCDISPLAY %{ - box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0); + box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); %} END