# Changeset 7e0b281 in sasmodels

Ignore:
Timestamp:
Apr 17, 2017 6:34:52 PM (15 months ago)
Branches:
master, ESS_GUI, beta_approx, cuda-test, ticket-1084, ticket-1102-pinhole, ticket-1105_mass_surface_fractal, ticket-1112, ticket-608-user-defined-weights
Children:
64ca163
Parents:
cb038a2
Message:

adjust paracrystal equations to better match terminology in Matsuoka paper

Files:
5 edited

Unmodified
Removed
• ## explore/bccpy.py

 rcb038a2 SLD_SOLVENT = 6.3 # Note: using Matsuoka 1990; this is different from what # is in the sasmodels/models code  (see bcc vs bcc_old). # The difference is that the sign of phi and theta seem to be # negative in the old vs. the new, yielding a pattern that is # swapped left to right and top to bottom. def sc(qa, qb, qc): return qa, qb, qc def bcc(qa, qb, qc): a1 = (+qa + qb + qc)/2 a2 = (-qa - qb + qc)/2 a3 = (-qa + qb - qc)/2 return a1, a2, a3 def bcc_old(qa, qb, qc): a1 = (+qa + qb - qc)/2.0 a2 = (+qa - qb + qc)/2.0 a3 = (-qa + qb + qc)/2.0 return a1, a2, a3 def fcc(qa, qb, qc): a1 = ( 0. + qb + qc)/2 a2 = (-qa + 0. + qc)/2 a3 = (-qa + qb + 0.)/2 return a1, a2, a3 def fcc_old(qa, qb, qc): a1 = ( qa + qb + 0.)/2 a2 = ( qa + 0. + qc)/2 a3 = ( 0. + qb + qc)/2 return a1, a2, a3 KERNEL = bcc def kernel(q, dnn, d_factor, theta, phi): """ qc = q*cos(theta) if 0: # sc a1, a2, a3 = qa, qb, qc dcos = dnn if 1: # bcc a1 = +qa - qc + qb a2 = +qa + qc - qb a3 = -qa + qc + qb dcos = dnn/2 if 0: # fcc a1 = qb + qa a2 = qa + qc a3 = qb + qc dcos = dnn/2 arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2) exp_arg = exp(-arg) den = [((exp_arg - 2*cos(dcos*a))*exp_arg + 1.0) for a in (a1, a2, a3)] Sq = -expm1(-2*arg)**3/np.prod(den, axis=0) return Sq a1, a2, a3 = KERNEL(qa, qb, qc) # Note: paper says that different directions can have different distortion # factors.  Easy enough to add to the code.  This would definitely break # 8-fold symmetry. arg = -0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2) exp_arg = exp(arg) den = [((exp_arg - 2*cos(dnn*a))*exp_arg + 1.0) for a in (a1, a2, a3)] Zq = -expm1(2*arg)**3/np.prod(den, axis=0) return Zq def integrand(theta, phi): evals[0] += 1 Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) return Sq*sin(theta) Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) return Zq*sin(theta) ans = dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi) print("dblquad evals =", evals[0]) Aw = w[None, :] * w[:, None] sin_theta = np.fmax(abs(sin(Atheta)), 1e-6) Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) print("gauss %d evals ="%n, n**2) return np.sum(Sq*Aw*sin_theta)*8/(4*pi) return np.sum(Zq*Aw*sin_theta)*8/(4*pi) phi = np.linspace(0, pi/2, n) Atheta, Aphi = np.meshgrid(theta, phi) Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) Sq *= abs(sin(Atheta)) Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) Zq *= abs(sin(Atheta)) dx, dy = theta[1]-theta[0], phi[1]-phi[0] print("rect", n, np.sum(Sq)*dx*dy*8/(4*pi)) print("trapz", n, np.trapz(np.trapz(Sq, dx=dx), dx=dy)*8/(4*pi)) print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi)) print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi)) print("rect", n, np.sum(Zq)*dx*dy*8/(4*pi)) print("trapz", n, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*8/(4*pi)) print("simpson", n, simps(simps(Zq, dx=dx), dx=dy)*8/(4*pi)) print("romb", n, romb(romb(Zq, dx=dx), dx=dy)*8/(4*pi)) print("gridded %d evals ="%n, n**2) #phi = np.linspace(0, pi/2, n) Atheta, Aphi = np.meshgrid(theta, phi) Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) Sq *= abs(sin(Atheta)) pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6))) Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) Zq *= abs(sin(Atheta)) pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6))) pylab.axis('tight') pylab.title("BCC S(q) for q=%g, dnn=%g d_factor=%g" % (q, dnn, d_factor)) pylab.title("%s Z(q) for q=%g, dnn=%g d_factor=%g" % (KERNEL.__name__, q, dnn, d_factor)) pylab.xlabel("theta (degrees)") pylab.ylabel("phi (degrees)")
• ## explore/sc.c

 rfdd56a1 static double _sq_sc(double qa, double qb, double qc, double dnn, double d_factor) sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) { // Rewriting equations for efficiency, accuracy and readability, and so const double a3 = qc; const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); // Numerator: (1 - exp(a)^2)^3 //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 const double exp_arg = exp(-arg); const double Sq = -cube(expm1(-2.0*arg)) const double exp_arg = exp(arg); const double Zq = -cube(expm1(2.0*arg)) / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); return Sq; return Zq; } // occupied volume fraction calculated from lattice symmetry and sphere radius static double _sc_volume_fraction(double radius, double dnn) sc_volume_fraction(double radius, double dnn) { return sphere_volume(radius/dnn); const double qa = qab*cos_phi; const double qb = qab*sin_phi; const double fq = _sq_sc(qa, qb, qc, dnn, d_factor); inner_sum += fq; const double form = sc_Zq(qa, qb, qc, dnn, d_factor); inner_sum += form; } inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx outer_sum *= theta_m/(n*n); #endif double Sq; double Zq; if (sym > 0.) { Sq = outer_sum/M_PI_2; Zq = outer_sum/M_PI_2; } else { Sq = outer_sum/(4.0*M_PI); Zq = outer_sum/(4.0*M_PI); } return Zq; const double Pq = sphere_form(q, radius, sld, solvent_sld); return _sc_volume_fraction(radius, dnn) * Pq * Sq; return sc_volume_fraction(radius, dnn) * Pq * Zq; } q = sqrt(qa*qa + qb*qb + qc*qc); const double Pq = sphere_form(q, radius, sld, solvent_sld); const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor); return _sc_volume_fraction(radius, dnn) * Pq * Sq; const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); return sc_volume_fraction(radius, dnn) * Pq * Zq; }
• ## sasmodels/models/bcc_paracrystal.c

 r2a0b2b1 static double _sq_bcc(double qa, double qb, double qc, double dnn, double d_factor) bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) { // Rewriting equations for efficiency, accuracy and readability, and so // code is reusable between 1D and 2D models. const double a1 = +qa - qc + qb; const double a2 = +qa + qc - qb; const double a3 = -qa + qc + qb; const double half_dnn = 0.5*dnn; const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); #if 0  // Equations as written in Matsuoka const double a1 = (+qa + qb + qc)/2.0; const double a2 = (-qa - qb + qc)/2.0; const double a3 = (-qa + qb - qc)/2.0; #else const double a1 = (+qa + qb - qc)/2.0; const double a2 = (+qa - qb + qc)/2.0; const double a3 = (-qa + qb + qc)/2.0; #endif #if 1 //         => (-(exp(2a) - 1))^3 //         => -expm1(2a)^3 // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 const double exp_arg = exp(-arg); const double Sq = -cube(expm1(-2.0*arg)) / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0)); // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); const double exp_arg = exp(arg); const double Zq = -cube(expm1(2.0*arg)) / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); #else // Alternate form, which perhaps is more approachable const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); const double sinh_qd = sinh(arg); const double cosh_qd = cosh(arg); const double Sq = sinh_qd/(cosh_qd - cos(half_dnn*a1)) * sinh_qd/(cosh_qd - cos(half_dnn*a2)) * sinh_qd/(cosh_qd - cos(half_dnn*a3)); const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) * sinh_qd/(cosh_qd - cos(dnn*a2)) * sinh_qd/(cosh_qd - cos(dnn*a3)); #endif return Sq; return Zq; } // occupied volume fraction calculated from lattice symmetry and sphere radius static double _bcc_volume_fraction(double radius, double dnn) bcc_volume_fraction(double radius, double dnn) { return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); const double qa = qab*cos_phi; const double qb = qab*sin_phi; const double fq = _sq_bcc(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * fq; const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * form; } inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx } outer_sum *= theta_m; const double Sq = outer_sum/(4.0*M_PI); const double Zq = outer_sum/(4.0*M_PI); const double Pq = sphere_form(q, radius, sld, solvent_sld); return _bcc_volume_fraction(radius, dnn) * Pq * Sq; return bcc_volume_fraction(radius, dnn) * Pq * Zq; } q = sqrt(qa*qa + qb*qb + qc*qc); const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); const double Pq = sphere_form(q, radius, sld, solvent_sld); const double Sq = _sq_bcc(qa, qb, qc, dnn, d_factor); return _bcc_volume_fraction(radius, dnn) * Pq * Sq; return bcc_volume_fraction(radius, dnn) * Pq * Zq; }
• ## sasmodels/models/fcc_paracrystal.c

 r2a0b2b1 static double _sq_fcc(double qa, double qb, double qc, double dnn, double d_factor) fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) { // Rewriting equations for efficiency, accuracy and readability, and so // code is reusable between 1D and 2D models. const double a1 = qb + qa; const double a2 = qa + qc; const double a3 = qb + qc; const double half_dnn = 0.5*dnn; const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); #if 0  // Equations as written in Matsuoka const double a1 = ( qa + qb)/2.0; const double a2 = (-qa + qc)/2.0; const double a3 = (-qa + qb)/2.0; #else const double a1 = ( qa + qb)/2.0; const double a2 = ( qa + qc)/2.0; const double a3 = ( qb + qc)/2.0; #endif // Numerator: (1 - exp(a)^2)^3 //         => (-(exp(2a) - 1))^3 //         => -expm1(2a)^3 // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 const double exp_arg = exp(-arg); const double Sq = -cube(expm1(-2.0*arg)) / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0)); // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); const double exp_arg = exp(arg); const double Zq = -cube(expm1(2.0*arg)) / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); return Sq; return Zq; } // occupied volume fraction calculated from lattice symmetry and sphere radius static double _fcc_volume_fraction(double radius, double dnn) fcc_volume_fraction(double radius, double dnn) { return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); const double qa = qab*cos_phi; const double qb = qab*sin_phi; const double fq = _sq_fcc(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * fq; const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * form; } inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx } outer_sum *= theta_m; const double Sq = outer_sum/(4.0*M_PI); const double Zq = outer_sum/(4.0*M_PI); const double Pq = sphere_form(q, radius, sld, solvent_sld); return _fcc_volume_fraction(radius, dnn) * Pq * Sq; return fcc_volume_fraction(radius, dnn) * Pq * Zq; } q = sqrt(qa*qa + qb*qb + qc*qc); const double Pq = sphere_form(q, radius, sld, solvent_sld); const double Sq = _sq_fcc(qa, qb, qc, dnn, d_factor); return _fcc_volume_fraction(radius, dnn) * Pq * Sq; const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); return fcc_volume_fraction(radius, dnn) * Pq * Zq; }
• ## sasmodels/models/sc_paracrystal.c

 r2a0b2b1 static double _sq_sc(double qa, double qb, double qc, double dnn, double d_factor) sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) { // Rewriting equations for efficiency, accuracy and readability, and so // code is reusable between 1D and 2D models. const double a1 = qa; const double a2 = qb; const double a3 = qc; const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); // Numerator: (1 - exp(a)^2)^3 //         => (-(exp(2a) - 1))^3 //         => -expm1(2a)^3 // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 const double exp_arg = exp(-arg); const double Sq = -cube(expm1(-2.0*arg)) // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); const double exp_arg = exp(arg); const double Zq = -cube(expm1(2.0*arg)) / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); return Sq; return Zq; } // occupied volume fraction calculated from lattice symmetry and sphere radius static double _sc_volume_fraction(double radius, double dnn) sc_volume_fraction(double radius, double dnn) { return sphere_volume(radius/dnn); const double qa = qab*cos_phi; const double qb = qab*sin_phi; const double fq = _sq_sc(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * fq; const double form = sc_Zq(qa, qb, qc, dnn, d_factor); inner_sum += Gauss150Wt[j] * form; } inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx } outer_sum *= theta_m; const double Sq = outer_sum/M_PI_2; const double Zq = outer_sum/M_PI_2; const double Pq = sphere_form(q, radius, sld, solvent_sld); return _sc_volume_fraction(radius, dnn) * Pq * Sq; return sc_volume_fraction(radius, dnn) * Pq * Zq; } q = sqrt(qa*qa + qb*qb + qc*qc); const double Pq = sphere_form(q, radius, sld, solvent_sld); const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor); return _sc_volume_fraction(radius, dnn) * Pq * Sq; const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); return sc_volume_fraction(radius, dnn) * Pq * Zq; }
Note: See TracChangeset for help on using the changeset viewer.