Changeset 2a0b2b1 in sasmodels for sasmodels/models/core_shell_cylinder.c
- Timestamp:
- Apr 14, 2017 6:30:29 AM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- fdd56a1
- Parents:
- 9901384
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_cylinder.c
r592343f r2a0b2b1 1 double form_volume(double radius, double thickness, double length);2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld,3 double radius, double thickness, double length);4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld,5 double radius, double thickness, double length, double theta, double phi);6 7 1 // vd = volume * delta_rho 8 // besarg = q * R * sin(alpha) 9 // siarg = q * L/2 * cos(alpha) 10 double _cyl(double vd, double besarg, double siarg); 11 double _cyl(double vd, double besarg, double siarg) 2 // besarg = q * R * sin(theta) 3 // siarg = q * L/2 * cos(theta) 4 static double _cyl(double vd, double besarg, double siarg) 12 5 { 13 6 return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 14 7 } 15 8 16 double form_volume(double radius, double thickness, double length) 9 static double 10 form_volume(double radius, double thickness, double length) 17 11 { 18 return M_PI* (radius+thickness)*(radius+thickness)*(length+2.0*thickness);12 return M_PI*square(radius+thickness)*(length+2.0*thickness); 19 13 } 20 14 21 double Iq(double q, 15 static double 16 Iq(double q, 22 17 double core_sld, 23 18 double shell_sld, … … 28 23 { 29 24 // precalculate constants 30 const double core_ qr = q*radius;31 const double core_ qh = q*0.5*length;25 const double core_r = radius; 26 const double core_h = 0.5*length; 32 27 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 33 const double shell_ qr = q*(radius + thickness);34 const double shell_ qh = q*(0.5*length + thickness);28 const double shell_r = (radius + thickness); 29 const double shell_h = (0.5*length + thickness); 35 30 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 36 31 double total = 0.0; 37 // double lower=0, upper=M_PI_2;38 32 for (int i=0; i<76 ;i++) { 39 // translate a point in [-1,1] to a point in [lower,upper] 40 //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 41 double sn, cn; 42 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 43 SINCOS(alpha, sn, cn); 44 const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 45 + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 46 total += Gauss76Wt[i] * fq * fq * sn; 33 // translate a point in [-1,1] to a point in [0, pi/2] 34 //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 35 double sin_theta, cos_theta; 36 const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 37 SINCOS(theta, sin_theta, cos_theta); 38 const double qab = q*sin_theta; 39 const double qc = q*cos_theta; 40 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 41 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 42 total += Gauss76Wt[i] * fq * fq * sin_theta; 47 43 } 48 44 // translate dx in [-1,1] to dx in [lower,upper] … … 64 60 double q, sin_alpha, cos_alpha; 65 61 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 62 const double qab = q*sin_alpha; 63 const double qc = q*cos_alpha; 66 64 67 const double core_ qr = q*radius;68 const double core_ qh = q*0.5*length;65 const double core_r = radius; 66 const double core_h = 0.5*length; 69 67 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 70 const double shell_ qr = q*(radius + thickness);71 const double shell_ qh = q*(0.5*length + thickness);68 const double shell_r = (radius + thickness); 69 const double shell_h = (0.5*length + thickness); 72 70 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 73 71 74 const double fq = _cyl(core_vd, core_ qr*sin_alpha, core_qh*cos_alpha)75 + _cyl(shell_vd, shell_ qr*sin_alpha, shell_qh*cos_alpha);72 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 73 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 76 74 return 1.0e-4 * fq * fq; 77 75 }
Note: See TracChangeset
for help on using the changeset viewer.