/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2008, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SasView_binary_hard_sphere * * %Identification * Written by: Jose Robledo * Based on sasmodels from SasView * Origin: FZJ / DTU / ESS DMSC * * * SasView binary_hard_sphere model component as sample description. * * %Description * * SasView_binary_hard_sphere component, generated from binary_hard_sphere.c in sasmodels. * * Example: * SasView_binary_hard_sphere(radius_lg, radius_sm, volfraction_lg, volfraction_sm, sld_lg, sld_sm, sld_solvent, * model_scale=1.0, model_abs=0.0, xwidth=0.01, yheight=0.01, zdepth=0.005, R=0, * int target_index=1, target_x=0, target_y=0, target_z=1, * focus_xw=0.5, focus_yh=0.5, focus_aw=0, focus_ah=0, focus_r=0, * pd_radius_lg=0.0, pd_radius_sm=0.0) * * %Parameters * INPUT PARAMETERS: * radius_lg: [AA] ([0, inf]) radius of large particle. * radius_sm: [AA] ([0, inf]) radius of small particle. * volfraction_lg: [] ([0, 1]) volume fraction of large particle. * volfraction_sm: [] ([0, 1]) volume fraction of small particle. * sld_lg: [1e-6/AA^2] ([-inf, inf]) scattering length density of large particle. * sld_sm: [1e-6/AA^2] ([-inf, inf]) scattering length density of small particle. * sld_solvent: [1e-6/AA^2] ([-inf, inf]) Solvent scattering length density. * Optional parameters: * model_abs: [ ] Absorption cross section density at 2200 m/s. * model_scale: [ ] Global scale factor for scattering kernel. For systems without inter-particle interference, the form factors can be related to the scattering intensity by the particle volume fraction. * xwidth: [m] ([-inf, inf]) Horiz. dimension of sample, as a width. * yheight: [m] ([-inf, inf]) vert . dimension of sample, as a height for cylinder/box * zdepth: [m] ([-inf, inf]) depth of sample * R: [m] Outer radius of sample in (x,z) plane for cylinder/sphere. * target_x: [m] relative focus target position. * target_y: [m] relative focus target position. * target_z: [m] relative focus target position. * target_index: [ ] Relative index of component to focus at, e.g. next is +1. * focus_xw: [m] horiz. dimension of a rectangular area. * focus_yh: [m], vert. dimension of a rectangular area. * focus_aw: [deg], horiz. angular dimension of a rectangular area. * focus_ah: [deg], vert. angular dimension of a rectangular area. * focus_r: [m] case of circular focusing, focusing radius. * pd_radius_lg: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_radius_sm: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable * * %Link * %End *******************************************************************************/ DEFINE COMPONENT SasView_binary_hard_sphere SETTING PARAMETERS ( radius_lg=100, radius_sm=25, volfraction_lg=0.1, volfraction_sm=0.2, sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36, model_scale=1.0, model_abs=0.0, xwidth=0.01, yheight=0.01, zdepth=0.005, R=0, target_x=0, target_y=0, target_z=1, int target_index=1, focus_xw=0.5, focus_yh=0.5, focus_aw=0, focus_ah=0, focus_r=0, pd_radius_lg=0.0, pd_radius_sm=0.0) SHARE %{ %include "sas_kernel_header.c" /* BEGIN Required header for SASmodel binary_hard_sphere */ #define HAS_Iq #define FORM_VOL #ifndef SAS_HAVE_sas_3j1x_x #define SAS_HAVE_sas_3j1x_x #line 1 "sas_3j1x_x" /** * Spherical Bessel function 3*j1(x)/x * * Used for low q to avoid cancellation error. * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but * in this case it is likely cancellation errors in the original expression * using double precision that are the source. */ double sas_3j1x_x(double q); // The choice of the number of terms in the series and the cutoff value for // switching between series and direct calculation depends on the numeric // precision. // // Point where direct calculation reaches machine precision: // // single machine precision eps 3e-8 at qr=1.1 ** // double machine precision eps 4e-16 at qr=1.1 // // Point where Taylor series reaches machine precision (eps), where taylor // series matches direct calculation (cross) and the error at that point: // // prec n eps cross error // single 3 0.28 0.4 6.2e-7 // single 4 0.68 0.7 2.3e-7 // single 5 1.18 1.2 7.5e-8 // double 3 0.01 0.03 2.3e-13 // double 4 0.06 0.1 3.1e-14 // double 5 0.16 0.2 5.0e-15 // // ** Note: relative error on single precision starts increase on the direct // method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3. This should be // safe for the sans range, with objects of 100 nm supported to a q of 0.1 // while maintaining 5 digits of precision. For usans/sesans, the objects // are larger but the q is smaller, so again it should be fine. // // See explore/sph_j1c.py for code to explore these ranges. // Use 4th order series #if FLOAT_SIZE>4 #define SPH_J1C_CUTOFF 0.1 #else #define SPH_J1C_CUTOFF 0.7 #endif #pragma acc routine seq double sas_3j1x_x(double q) { // 2017-05-18 PAK - support negative q if (fabs(q) < SPH_J1C_CUTOFF) { const double q2 = q*q; return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); } else { double sin_q, cos_q; SINCOS(q, sin_q, cos_q); return 3.0*(sin_q/q - cos_q)/(q*q); } } #endif // SAS_HAVE_sas_3j1x_x #ifndef SAS_HAVE_binary_hard_sphere #define SAS_HAVE_binary_hard_sphere #line 1 "binary_hard_sphere" double form_volume_binary_hard_sphere(void); double Iq_binary_hard_sphere(double q, double lg_radius, double sm_radius, double lg_vol_frac, double sm_vol_frac, double lg_sld, double sm_sld, double solvent_sld ); void calculate_psfs(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12 ); double form_volume_binary_hard_sphere(void) { return 1.0; } double Iq_binary_hard_sphere(double q, double lg_radius, double sm_radius, double lg_vol_frac, double sm_vol_frac, double lg_sld, double sm_sld, double solvent_sld) { double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten; //my local names double psf11,psf12,psf22; double phi1,phi2,phr,a3; double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; r2 = lg_radius; r1 = sm_radius; phi2 = lg_vol_frac; phi1 = sm_vol_frac; rho2 = lg_sld; rho1 = sm_sld; rhos = solvent_sld; phi = phi1 + phi2; aa = r1/r2; //calculate the number fraction of larger spheres (eqn 2 in reference) a3=aa*aa*aa; phr=phi2/phi; nf2 = phr*a3/(1.0-phr+phr*a3); // calculate the PSF's here calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12); // /* do form factor calculations */ v1 = M_4PI_3*r1*r1*r1; v2 = M_4PI_3*r2*r2*r2; n1 = phi1/v1; n2 = phi2/v2; qr1 = r1*q; qr2 = r2*q; sc1 = sas_3j1x_x(qr1); sc2 = sas_3j1x_x(qr2); b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; inten = n1*b1*b1*psf11; inten += sqrt(n1*n2)*2.0*b1*b2*psf12; inten += n2*b2*b2*psf22; ///* convert I(1/A) to (1/cm) */ inten *= 1.0e8; ///*convert rho^2 in 10^-6A to A*/ inten *= 1.0e-12; return(inten); } void calculate_psfs(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12) { // variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 // calculate constant terms double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; double a1,a2i,a2,b1,b2,b12,gm1,gm12; double yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; s2 = 2.0*r2; // s1 = aa*s2; why is this never used? check original paper? v = phi; a3 = aa*aa*aa; v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v; v2=(nf2/(nf2+(1.-nf2)*a3))*v; g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v); g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v); g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v); wmv = 1/(1.-v); wmv3 = wmv*wmv*wmv; wmv4 = wmv*wmv3; a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3; a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3; a2=a2i/a3; b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12); b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12); b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12; gm1=(v1*a1+a3*v2*a2)*.5; gm12=2.*gm1*(1.-aa)/aa; //c //c calculate the direct correlation functions and print results //c // do 20 j=1,npts yy=qval*s2; //c calculate direct correlation functions //c ----c11 ay=aa*yy; ay2 = ay*ay; ay3 = ay*ay*ay; t1=a1*(sin(ay)-ay*cos(ay)); t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay; t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3; f11=24.*v1*(t1+t2+t3)/ay3; //c ------c22 y2=yy*yy; y3=yy*y2; tt1=a2*(sin(yy)-yy*cos(yy)); tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy; tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3; f22=24.*v2*(tt1+tt2+tt3)/y3; //c -----c12 yl=.5*yy*(1.-aa); yl3=yl*yl*yl; wma3 = (1.-aa)*(1.-aa)*(1.-aa); y1=aa*yy; y13 = y1*y1*y1; ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3); t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1)); t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1; t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1); t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.); t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1; t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1); t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1); t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1); ttt2=sin(yl)*(t21+t22+t23)/(y13*y1); ttt3=cos(yl)*(t31+t32+t33)/(y13*y1); ttt4=a1*(t41+t42)/y1; f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3); c11=f11; c22=f22; c12=f12; *s11=1./(1.+c11-(c12)*c12/(1.+c22)); *s22=1./(1.+c22-(c12)*c12/(1.+c11)); *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12)); return; } #endif // SAS_HAVE_binary_hard_sphere /* END Required header for SASmodel binary_hard_sphere */ %} DECLARE %{ double shape; double my_a_v; %} INITIALIZE %{ shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere */ if (xwidth && yheight && zdepth) shape = 1; else if (R > 0 && yheight) shape = 0; else if (R > 0 && !yheight) shape = 2; if (shape < 0) exit (fprintf (stderr, "SasView_model: %s: sample has invalid dimensions.\n" "ERROR Please check parameter values.\n", NAME_CURRENT_COMP)); /* now compute target coords if a component index is supplied */ if (!target_index && !target_x && !target_y && !target_z) target_index = 1; if (target_index) { Coords ToTarget; ToTarget = coords_sub (POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP); ToTarget = rot_apply (ROT_A_CURRENT_COMP, ToTarget); coords_get (ToTarget, &target_x, &target_y, &target_z); } if (!(target_x || target_y || target_z)) { printf ("SasView_model: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); target_z = 1; } my_a_v = model_abs * 2200 * 100; /* Is not yet divided by v. 100: Convert barns -> fm^2 */ %} TRACE %{ double t0, t1, v, l_full, l, l_1, dt, d_phi, my_s; double aim_x = 0, aim_y = 0, aim_z = 1, axis_x, axis_y, axis_z; double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z; double f, solid_angle, vx_i, vy_i, vz_i, q, qx, qy, qz; char intersect = 0; /* Intersection neutron trajectory / sample (sample surface) */ if (shape == 0) { intersect = cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, R, yheight); } else if (shape == 1) { intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); } else if (shape == 2) { intersect = sphere_intersect (&t0, &t1, x, y, z, vx, vy, vz, R); } if (intersect) { if (t0 < 0) ABSORB; /* Neutron enters at t=t0. */ v = sqrt (vx * vx + vy * vy + vz * vz); l_full = v * (t1 - t0); /* Length of full path through sample */ dt = rand01 () * (t1 - t0) + t0; /* Time of scattering */ PROP_DT (dt); /* Point of scattering */ l = v * (dt - t0); /* Penetration in sample */ vx_i = vx; vy_i = vy; vz_i = vz; if ((target_x || target_y || target_z)) { aim_x = target_x - x; /* Vector pointing at target (anal./det.) */ aim_y = target_y - y; aim_z = target_z - z; } if (focus_aw && focus_ah) { randvec_target_rect_angular (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); } else if (focus_xw && focus_yh) { randvec_target_rect (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); } else { randvec_target_circle (&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r); } NORM (vx, vy, vz); vx *= v; vy *= v; vz *= v; qx = V2K * (vx_i - vx); qy = V2K * (vy_i - vy); qz = V2K * (vz_i - vz); q = sqrt (qx * qx + qy * qy + qz * qz); double trace_radius_lg = radius_lg; double trace_radius_sm = radius_sm; if (pd_radius_lg != 0.0 || pd_radius_sm != 0.0) { trace_radius_lg = (randnorm () * pd_radius_lg + 1.0) * radius_lg; trace_radius_sm = (randnorm () * pd_radius_sm + 1.0) * radius_sm; } // Sample dependent. Retrieved from SasView.///////////////////// float Iq_out; Iq_out = 1; Iq_out = Iq_binary_hard_sphere (q, trace_radius_lg, trace_radius_sm, volfraction_lg, volfraction_sm, sld_lg, sld_sm, sld_solvent); float vol; vol = 1; // Scale by 1.0E2 [SasView: 1/cm -> McStas: 1/m] Iq_out = model_scale * Iq_out / vol * 1.0E2; l_1 = v * t1; p *= l_full * solid_angle / (4 * PI) * Iq_out * exp (-my_a_v * (l + l_1) / v); SCATTER; } %} MCDISPLAY %{ if (shape == 0) { /* cylinder */ cylinder (0, 0, 0, R, yheight, 0, 0, 1, 0); } else if (shape == 1) { /* box */ box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); } else if (shape == 2) { /* sphere */ sphere (0, 0, 0, R); } %} END