Changeset 50e1e40 in sasmodels for sasmodels/models/barbell.c
- Timestamp:
- Mar 1, 2016 5:49:00 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- ad90df9
- Parents:
- a4a7308
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/barbell.c
r3c97ff0 r50e1e40 7 7 8 8 //barbell kernel - same as dumbell 9 double _bell_kernel(double q, double h, double bell_radius, 10 double length, double sin_alpha, double cos_alpha); 11 double _bell_kernel(double q, double h, double bell_radius, 12 double length, double sin_alpha, double cos_alpha) 9 static double 10 _bell_kernel(double q, double h, double bell_radius, 11 double half_length, double sin_alpha, double cos_alpha) 13 12 { 13 // translate a point in [-1,1] to a point in [lower,upper] 14 14 const double upper = 1.0; 15 const double lower = -1.0*h/bell_radius; 15 const double lower = h/bell_radius; 16 const double zm = 0.5*(upper-lower); 17 const double zb = 0.5*(upper+lower); 16 18 19 // cos term in integral is: 20 // cos (q (R t - h + L/2) cos(alpha)) 21 // so turn it into: 22 // cos (m t + b) 23 // where: 24 // m = q R cos(alpha) 25 // b = q(L/2-h) cos(alpha) 26 const double m = q*bell_radius*cos_alpha; // cos argument slope 27 const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 28 const double qrst = q*bell_radius*sin_alpha; // Q*R*sin(theta) 17 29 double total = 0.0; 18 30 for (int i = 0; i < 76; i++){ 19 const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower); 20 const double arg1 = q*cos_alpha*(bell_radius*t+h+length*0.5); 21 const double arg2 = q*bell_radius*sin_alpha*sqrt(1.0-t*t); 22 const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2); 23 const double Fq = cos(arg1)*(1.0-t*t)*be; 31 const double t = Gauss76Z[i]*zm + zb; 32 const double radical = 1.0 - t*t; 33 const double bj = J1c(qrst*sqrt(radical)); 34 const double Fq = cos(m*t + b) * radical * bj; 24 35 total += Gauss76Wt[i] * Fq; 25 36 } 26 const double integral = 0.5*(upper-lower)*total; 27 return 4.0*M_PI*bell_radius*bell_radius*bell_radius*integral; 37 // translate dx in [-1,1] to dx in [lower,upper] 38 const double integral = total*zm; 39 const double bell_Fq = 2*M_PI*cube(bell_radius)*integral; 40 return bell_Fq; 28 41 } 29 42 … … 34 47 35 48 // bell radius should never be less than radius when this is called 36 const double hdist = sqrt( bell_radius*bell_radius - radius*radius);37 const double p1 = 2.0 *bell_radius*bell_radius*bell_radius/3.0;38 const double p2 = bell_radius*bell_radius*hdist;39 const double p3 = hdist*hdist*hdist/3.0;49 const double hdist = sqrt(square(bell_radius) - square(radius)); 50 const double p1 = 2.0/3.0*cube(bell_radius); 51 const double p2 = square(bell_radius)*hdist; 52 const double p3 = cube(hdist)/3.0; 40 53 41 return M_PI* radius*radius*length + 2.0*M_PI*(p1+p2-p3);54 return M_PI*square(radius)*length + 2.0*M_PI*(p1+p2-p3); 42 55 } 43 56 44 double Iq(double q, double sld, 45 double solvent_sld, 46 double bell_radius, 47 double radius, 48 double length) 57 double Iq(double q, double sld, double solvent_sld, 58 double bell_radius, double radius, double length) 49 59 { 50 double sn, cn; // slots to hold sincos function output 60 // Exclude invalid inputs. 61 if (bell_radius < radius) return NAN; 62 const double h = -sqrt(bell_radius*bell_radius - radius*radius); 63 const double half_length = 0.5*length; 51 64 52 if (bell_radius < radius) return NAN; 53 54 const double lower = 0.0; 55 const double upper = M_PI_2; 56 const double h = sqrt(bell_radius*bell_radius-radius*radius); 65 // translate a point in [-1,1] to a point in [0, pi/2] 66 const double zm = M_PI_4; 67 const double zb = M_PI_4; 57 68 double total = 0.0; 58 69 for (int i = 0; i < 76; i++){ 59 const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower); 60 SINCOS(alpha, sn, cn); 70 const double alpha= Gauss76Z[i]*zm + zb; 71 double sin_alpha, cos_alpha; // slots to hold sincos function output 72 SINCOS(alpha, sin_alpha, cos_alpha); 61 73 62 const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn); 74 const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha); 75 const double bj = J1c(q*radius*sin_alpha); 76 const double si = sinc(q*half_length*cos_alpha); 77 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 78 const double Aq = bell_Fq + cyl_Fq; 79 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 80 } 81 // translate dx in [-1,1] to dx in [lower,upper] 82 const double form = total*zm; 63 83 64 const double arg1 = q*length*0.5*cn; 65 const double arg2 = q*radius*sn; 66 // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 67 const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2); 68 const double si = (arg1 == 0.0 ? 1.0 :sin(arg1)/arg1); 69 const double cyl_Fq = M_PI*radius*radius*length*si*2.0*be; 70 71 const double Aq = cyl_Fq + bell_Fq; 72 total += Gauss76Wt[i] * Aq * Aq * sn; 73 } 74 75 const double form = total*(upper-lower)*0.5; 76 77 //Contrast and volume normalization 84 //Contrast 78 85 const double s = (sld - solvent_sld); 79 return form*1.0e-4*s*s; //form_volume(bell_radius,radius,length);86 return 1.0e-4 * s * s * form; 80 87 } 81 88 82 89 83 84 90 double Iqxy(double qx, double qy, 85 double sld, 86 double solvent_sld, 87 double bell_radius, 88 double radius, 89 double length, 90 double theta, 91 double phi) 91 double sld, double solvent_sld, 92 double bell_radius, double radius, double length, 93 double theta, double phi) 92 94 { 93 double sn, cn; // slots to hold sincos function output 95 // Compute angle alpha between q and the cylinder axis 96 double sn, cn; // slots to hold sincos function output 97 SINCOS(theta*M_PI_180, sn, cn); 98 const double q = sqrt(qx*qx+qy*qy); 99 const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 100 const double alpha = acos(cos_val); // rod angle relative to q 94 101 95 102 // Exclude invalid inputs. 96 103 if (bell_radius < radius) return NAN; 104 const double h = -sqrt(square(bell_radius) - square(radius)); 105 const double half_length = 0.5*length; 97 106 98 // Compute angle alpha between q and the cylinder axis 99 SINCOS(theta*M_PI_180, sn, cn); 100 // # The following correction factor exists in sasview, but it can't be 101 // # right, so we are leaving it out for now. 102 const double q = sqrt(qx*qx+qy*qy); 103 const double cos_val = cn*cos(phi*M_PI_180)*qx + sn*qy; 104 const double alpha = acos(cos_val); // rod angle relative to q 105 SINCOS(alpha, sn, cn); 106 107 const double h = sqrt(bell_radius*bell_radius - radius*radius); // negative h 108 const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn)/sn; 109 110 const double besarg = q*radius*sn; 111 const double siarg = q*0.5*length*cn; 112 // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 113 const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 114 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 115 const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 116 117 // Volume weighted average F(q) 107 double sin_alpha, cos_alpha; // slots to hold sincos function output 108 SINCOS(alpha, sin_alpha, cos_alpha); 109 const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha); 110 const double bj = J1c(q*radius*sin_alpha); 111 const double si = sinc(q*half_length*cos_alpha); 112 const double cyl_Fq = M_PI*radius*radius*length*bj*si; 118 113 const double Aq = cyl_Fq + bell_Fq; 119 114 120 // Multiply by contrast^2 , normalize by cylinder volumeand convert to cm-1115 // Multiply by contrast^2 and convert to cm-1 121 116 const double s = (sld - solvent_sld); 122 return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length);117 return 1.0e-4 * square(s * Aq); 123 118 }
Note: See TracChangeset
for help on using the changeset viewer.