Changeset d4db147 in sasmodels for sasmodels/models
- Timestamp:
- Nov 6, 2017 4:11:48 PM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 9248bf7
- Parents:
- ff31782 (diff), d8ac2ad (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Location:
- sasmodels/models
- Files:
-
- 1 deleted
- 27 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_bicelle_elliptical.c
r74768cb rd4db147 92 92 93 93 // Compute effective radius in rotated coordinates 94 const double qr_hat = sqrt(square(r_major*q a) + square(r_minor*qb));95 const double qrshell_hat = sqrt(square((r_major+thick_rim)*q a)96 + square((r_minor+thick_rim)*q b));94 const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 95 const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 96 + square((r_minor+thick_rim)*qa)); 97 97 const double be1 = sas_2J1x_x( qr_hat ); 98 98 const double be2 = sas_2J1x_x( qrshell_hat ); -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
r74768cb rd4db147 91 91 double sigma) 92 92 { 93 // THIS NEEDS TESTING93 // integrated 2d seems to match 1d reasonably well, except perhaps at very high Q 94 94 // Vol1,2,3 and dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 95 95 const double dr1 = -rhor - rhoh + rhoc + rhosolv; … … 103 103 104 104 // Compute effective radius in rotated coordinates 105 const double qr_hat = sqrt(square(r_major*q a) + square(r_minor*qb));105 const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 106 106 // does this need to be changed for the "missing corners" where there there is no "belt" ? 107 const double qrshell_hat = sqrt(square((r_major+thick_rim)*q a)108 + square((r_minor+thick_rim)*q b));107 const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 108 + square((r_minor+thick_rim)*qa)); 109 109 const double be1 = sas_2J1x_x( qr_hat ); 110 110 const double be2 = sas_2J1x_x( qrshell_hat ); -
sasmodels/models/elliptical_cylinder.c
r74768cb rd4db147 60 60 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 61 61 // Given: radius_major = r_ratio * radius_minor 62 const double qr = radius_minor*sqrt(square(r_ratio*q a) + square(qb));62 const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa)); 63 63 const double be = sas_2J1x_x(qr); 64 64 const double si = sas_sinx_x(qc*0.5*length); -
sasmodels/models/barbell.c
rbecded3 r74768cb 23 23 const double qab_r = radius_bell*qab; // Q*R*sin(theta) 24 24 double total = 0.0; 25 for (int i = 0; i < 76; i++){26 const double t = G auss76Z[i]*zm + zb;25 for (int i = 0; i < GAUSS_N; i++){ 26 const double t = GAUSS_Z[i]*zm + zb; 27 27 const double radical = 1.0 - t*t; 28 28 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 29 29 const double Fq = cos(m*t + b) * radical * bj; 30 total += G auss76Wt[i] * Fq;30 total += GAUSS_W[i] * Fq; 31 31 } 32 32 // translate dx in [-1,1] to dx in [lower,upper] … … 73 73 const double zb = M_PI_4; 74 74 double total = 0.0; 75 for (int i = 0; i < 76; i++){76 const double alpha= G auss76Z[i]*zm + zb;75 for (int i = 0; i < GAUSS_N; i++){ 76 const double alpha= GAUSS_Z[i]*zm + zb; 77 77 double sin_alpha, cos_alpha; // slots to hold sincos function output 78 78 SINCOS(alpha, sin_alpha, cos_alpha); 79 79 const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 80 total += G auss76Wt[i] * Aq * Aq * sin_alpha;80 total += GAUSS_W[i] * Aq * Aq * sin_alpha; 81 81 } 82 82 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/bcc_paracrystal.c
rea60e08 r74768cb 81 81 82 82 double outer_sum = 0.0; 83 for(int i=0; i< 150; i++) {83 for(int i=0; i<GAUSS_N; i++) { 84 84 double inner_sum = 0.0; 85 const double theta = G auss150Z[i]*theta_m + theta_b;85 const double theta = GAUSS_Z[i]*theta_m + theta_b; 86 86 double sin_theta, cos_theta; 87 87 SINCOS(theta, sin_theta, cos_theta); 88 88 const double qc = q*cos_theta; 89 89 const double qab = q*sin_theta; 90 for(int j=0;j< 150;j++) {91 const double phi = G auss150Z[j]*phi_m + phi_b;90 for(int j=0;j<GAUSS_N;j++) { 91 const double phi = GAUSS_Z[j]*phi_m + phi_b; 92 92 double sin_phi, cos_phi; 93 93 SINCOS(phi, sin_phi, cos_phi); … … 95 95 const double qb = qab*sin_phi; 96 96 const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 97 inner_sum += G auss150Wt[j] * form;97 inner_sum += GAUSS_W[j] * form; 98 98 } 99 99 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 100 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;100 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 101 101 } 102 102 outer_sum *= theta_m; -
sasmodels/models/capped_cylinder.c
rbecded3 r74768cb 30 30 const double qab_r = radius_cap*qab; // Q*R*sin(theta) 31 31 double total = 0.0; 32 for (int i=0; i< 76 ;i++) {33 const double t = G auss76Z[i]*zm + zb;32 for (int i=0; i<GAUSS_N; i++) { 33 const double t = GAUSS_Z[i]*zm + zb; 34 34 const double radical = 1.0 - t*t; 35 35 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 36 36 const double Fq = cos(m*t + b) * radical * bj; 37 total += G auss76Wt[i] * Fq;37 total += GAUSS_W[i] * Fq; 38 38 } 39 39 // translate dx in [-1,1] to dx in [lower,upper] … … 95 95 const double zb = M_PI_4; 96 96 double total = 0.0; 97 for (int i=0; i< 76;i++) {98 const double theta = G auss76Z[i]*zm + zb;97 for (int i=0; i<GAUSS_N ;i++) { 98 const double theta = GAUSS_Z[i]*zm + zb; 99 99 double sin_theta, cos_theta; // slots to hold sincos function output 100 100 SINCOS(theta, sin_theta, cos_theta); … … 103 103 const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 104 104 // scale by sin_theta for spherical coord integration 105 total += G auss76Wt[i] * Aq * Aq * sin_theta;105 total += GAUSS_W[i] * Aq * Aq * sin_theta; 106 106 } 107 107 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/core_shell_bicelle.c
rbecded3 r74768cb 52 52 53 53 double total = 0.0; 54 for(int i=0;i< N_POINTS_76;i++) {55 double theta = (G auss76Z[i] + 1.0)*uplim;54 for(int i=0;i<GAUSS_N;i++) { 55 double theta = (GAUSS_Z[i] + 1.0)*uplim; 56 56 double sin_theta, cos_theta; // slots to hold sincos function output 57 57 SINCOS(theta, sin_theta, cos_theta); 58 58 double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 59 59 halflength, sld_core, sld_face, sld_rim, sld_solvent); 60 total += G auss76Wt[i]*fq*fq*sin_theta;60 total += GAUSS_W[i]*fq*fq*sin_theta; 61 61 } 62 62 -
sasmodels/models/core_shell_cylinder.c
rbecded3 r74768cb 30 30 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 31 31 double total = 0.0; 32 for (int i=0; i< 76;i++) {32 for (int i=0; i<GAUSS_N ;i++) { 33 33 // translate a point in [-1,1] to a point in [0, pi/2] 34 //const double theta = ( G auss76Z[i]*(upper-lower) + upper + lower )/2.0;34 //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 35 35 double sin_theta, cos_theta; 36 const double theta = G auss76Z[i]*M_PI_4 + M_PI_4;36 const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 37 37 SINCOS(theta, sin_theta, cos_theta); 38 38 const double qab = q*sin_theta; … … 40 40 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 41 41 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 42 total += G auss76Wt[i] * fq * fq * sin_theta;42 total += GAUSS_W[i] * fq * fq * sin_theta; 43 43 } 44 44 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/core_shell_ellipsoid.c
rbecded3 r74768cb 59 59 const double b = 0.5; 60 60 double total = 0.0; //initialize intergral 61 for(int i=0;i< 76;i++) {62 const double cos_theta = G auss76Z[i]*m + b;61 for(int i=0;i<GAUSS_N;i++) { 62 const double cos_theta = GAUSS_Z[i]*m + b; 63 63 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 64 64 double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, … … 66 66 equat_shell, polar_shell, 67 67 sld_core_shell, sld_shell_solvent); 68 total += G auss76Wt[i] * fq * fq;68 total += GAUSS_W[i] * fq * fq; 69 69 } 70 70 total *= m; -
sasmodels/models/core_shell_parallelepiped.c
rfc0b7aa r74768cb 80 80 // outer integral (with gauss points), integration limits = 0, 1 81 81 double outer_sum = 0; //initialize integral 82 for( int i=0; i< 76; i++) {83 double sigma = 0.5 * ( G auss76Z[i] + 1.0 );82 for( int i=0; i<GAUSS_N; i++) { 83 double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 84 84 double mu_proj = mu * sqrt(1.0-sigma*sigma); 85 85 … … 88 88 const double siCt = sas_sinx_x(mu * sigma * tC_over_b); 89 89 double inner_sum = 0.0; 90 for(int j=0; j< 76; j++) {91 const double uu = 0.5 * ( G auss76Z[j] + 1.0 );90 for(int j=0; j<GAUSS_N; j++) { 91 const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 92 92 double sin_uu, cos_uu; 93 93 SINCOS(M_PI_2*uu, sin_uu, cos_uu); … … 109 109 #endif 110 110 111 inner_sum += G auss76Wt[j] * f * f;111 inner_sum += GAUSS_W[j] * f * f; 112 112 } 113 113 inner_sum *= 0.5; 114 114 // now sum up the outer integral 115 outer_sum += G auss76Wt[i] * inner_sum;115 outer_sum += GAUSS_W[i] * inner_sum; 116 116 } 117 117 outer_sum *= 0.5; -
sasmodels/models/cylinder.c
rbecded3 r74768cb 21 21 22 22 double total = 0.0; 23 for (int i=0; i< 76;i++) {24 const double theta = G auss76Z[i]*zm + zb;23 for (int i=0; i<GAUSS_N ;i++) { 24 const double theta = GAUSS_Z[i]*zm + zb; 25 25 double sin_theta, cos_theta; // slots to hold sincos function output 26 26 // theta (theta,phi) the projection of the cylinder on the detector plane 27 27 SINCOS(theta , sin_theta, cos_theta); 28 28 const double form = fq(q*sin_theta, q*cos_theta, radius, length); 29 total += G auss76Wt[i] * form * form * sin_theta;29 total += GAUSS_W[i] * form * form * sin_theta; 30 30 } 31 31 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/ellipsoid.c
rbecded3 r74768cb 22 22 23 23 // translate a point in [-1,1] to a point in [0, 1] 24 // const double u = G auss76Z[i]*(upper-lower)/2 + (upper+lower)/2;24 // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 25 25 const double zm = 0.5; 26 26 const double zb = 0.5; 27 27 double total = 0.0; 28 for (int i=0;i< 76;i++) {29 const double u = G auss76Z[i]*zm + zb;28 for (int i=0;i<GAUSS_N;i++) { 29 const double u = GAUSS_Z[i]*zm + zb; 30 30 const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 31 31 const double f = sas_3j1x_x(q*r); 32 total += G auss76Wt[i] * f * f;32 total += GAUSS_W[i] * f * f; 33 33 } 34 34 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/elliptical_cylinder.py
reda8b30 r74768cb 43 43 P(q) = \text{scale} <F^2> / V 44 44 45 For 2d data the orientation of the particle is required, described using a different set 46 of angles as in the diagrams below, for further details of the calculation and angular 45 For 2d data the orientation of the particle is required, described using a different set 46 of angles as in the diagrams below, for further details of the calculation and angular 47 47 dispersions see :ref:`orientation` . 48 48 … … 121 121 # pylint: enable=bad-whitespace, line-too-long 122 122 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 124 "elliptical_cylinder.c"] 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 125 124 126 125 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, -
sasmodels/models/fcc_paracrystal.c
rf728001 r74768cb 53 53 54 54 double outer_sum = 0.0; 55 for(int i=0; i< 150; i++) {55 for(int i=0; i<GAUSS_N; i++) { 56 56 double inner_sum = 0.0; 57 const double theta = G auss150Z[i]*theta_m + theta_b;57 const double theta = GAUSS_Z[i]*theta_m + theta_b; 58 58 double sin_theta, cos_theta; 59 59 SINCOS(theta, sin_theta, cos_theta); 60 60 const double qc = q*cos_theta; 61 61 const double qab = q*sin_theta; 62 for(int j=0;j< 150;j++) {63 const double phi = G auss150Z[j]*phi_m + phi_b;62 for(int j=0;j<GAUSS_N;j++) { 63 const double phi = GAUSS_Z[j]*phi_m + phi_b; 64 64 double sin_phi, cos_phi; 65 65 SINCOS(phi, sin_phi, cos_phi); … … 67 67 const double qb = qab*sin_phi; 68 68 const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 69 inner_sum += G auss150Wt[j] * form;69 inner_sum += GAUSS_W[j] * form; 70 70 } 71 71 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 72 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;72 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 73 73 } 74 74 outer_sum *= theta_m; -
sasmodels/models/flexible_cylinder_elliptical.c
r592343f r74768cb 17 17 double sum=0.0; 18 18 19 for(int i=0;i< N_POINTS_76;i++) {20 const double zi = ( G auss76Z[i] + 1.0 )*M_PI_4;19 for(int i=0;i<GAUSS_N;i++) { 20 const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 21 21 double sn, cn; 22 22 SINCOS(zi, sn, cn); 23 23 const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 24 24 const double yyy = sas_2J1x_x(arg); 25 sum += G auss76Wt[i] * yyy * yyy;25 sum += GAUSS_W[i] * yyy * yyy; 26 26 } 27 27 sum *= 0.5; -
sasmodels/models/hollow_cylinder.c
rbecded3 r74768cb 38 38 39 39 double summ = 0.0; //initialize intergral 40 for (int i=0;i< 76;i++) {41 const double cos_theta = 0.5*( G auss76Z[i] * (upper-lower) + lower + upper );40 for (int i=0;i<GAUSS_N;i++) { 41 const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 42 42 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 43 43 const double form = _fq(q*sin_theta, q*cos_theta, 44 44 radius, thickness, length); 45 summ += G auss76Wt[i] * form * form;45 summ += GAUSS_W[i] * form * form; 46 46 } 47 47 -
sasmodels/models/hollow_rectangular_prism.c
r8de1477 r74768cb 39 39 40 40 double outer_sum = 0.0; 41 for(int i=0; i< 76; i++) {41 for(int i=0; i<GAUSS_N; i++) { 42 42 43 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );43 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 44 44 double sin_theta, cos_theta; 45 45 SINCOS(theta, sin_theta, cos_theta); … … 49 49 50 50 double inner_sum = 0.0; 51 for(int j=0; j< 76; j++) {51 for(int j=0; j<GAUSS_N; j++) { 52 52 53 const double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );53 const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 54 54 double sin_phi, cos_phi; 55 55 SINCOS(phi, sin_phi, cos_phi); … … 66 66 const double AP2 = vol_core * termA2 * termB2 * termC2; 67 67 68 inner_sum += G auss76Wt[j] * square(AP1-AP2);68 inner_sum += GAUSS_W[j] * square(AP1-AP2); 69 69 } 70 70 inner_sum *= 0.5 * (v2b-v2a); 71 71 72 outer_sum += G auss76Wt[i] * inner_sum * sin(theta);72 outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 73 73 } 74 74 outer_sum *= 0.5*(v1b-v1a); -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
rab2aea8 r74768cb 1 1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 2 double Iq(double q, double sld, double solvent_sld, double length_a, 2 double Iq(double q, double sld, double solvent_sld, double length_a, 3 3 double b2a_ratio, double c2a_ratio); 4 4 … … 29 29 const double v2a = 0.0; 30 30 const double v2b = M_PI_2; //phi integration limits 31 31 32 32 double outer_sum = 0.0; 33 for(int i=0; i< 76; i++) {34 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );33 for(int i=0; i<GAUSS_N; i++) { 34 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 35 35 36 36 double sin_theta, cos_theta; … … 44 44 45 45 double inner_sum = 0.0; 46 for(int j=0; j< 76; j++) {47 const double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );46 for(int j=0; j<GAUSS_N; j++) { 47 const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 48 48 49 49 double sin_phi, cos_phi; … … 62 62 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 63 63 64 inner_sum += G auss76Wt[j] * square(AL+AT);64 inner_sum += GAUSS_W[j] * square(AL+AT); 65 65 } 66 66 67 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += G auss76Wt[i] * inner_sum * sin_theta;68 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 69 69 } 70 70 -
sasmodels/models/lib/gauss150.c
r994d77f r74768cb 7 7 * 8 8 */ 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 150 15 #define GAUSS_Z Gauss150Z 16 #define GAUSS_W Gauss150Wt 17 18 // Note: using array size 152 so that it is a multiple of 4 9 19 10 20 // Gaussians 11 constant double Gauss150Z[15 0]={21 constant double Gauss150Z[152]={ 12 22 -0.9998723404457334, 13 23 -0.9993274305065947, … … 159 169 0.9983473449340834, 160 170 0.9993274305065947, 161 0.9998723404457334 171 0.9998723404457334, 172 0., 173 0. 162 174 }; 163 175 164 constant double Gauss150Wt[15 0]={176 constant double Gauss150Wt[152]={ 165 177 0.0003276086705538, 166 178 0.0007624720924706, … … 312 324 0.0011976474864367, 313 325 0.0007624720924706, 314 0.0003276086705538 326 0.0003276086705538, 327 0., 328 0. 315 329 }; -
sasmodels/models/lib/gauss20.c
r994d77f r74768cb 7 7 * 8 8 */ 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 20 15 #define GAUSS_Z Gauss20Z 16 #define GAUSS_W Gauss20Wt 9 17 10 18 // Gaussians -
sasmodels/models/lib/gauss76.c
r66d119f r74768cb 7 7 * 8 8 */ 9 #define N_POINTS_76 76 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 76 15 #define GAUSS_Z Gauss76Z 16 #define GAUSS_W Gauss76Wt 10 17 11 18 // Gaussians 12 constant double Gauss76Wt[ N_POINTS_76]={19 constant double Gauss76Wt[76]={ 13 20 .00126779163408536, //0 14 21 .00294910295364247, … … 89 96 }; 90 97 91 constant double Gauss76Z[ N_POINTS_76]={98 constant double Gauss76Z[76]={ 92 99 -.999505948362153, //0 93 100 -.997397786355355, -
sasmodels/models/parallelepiped.c
r9b7b23f r74768cb 23 23 double outer_total = 0; //initialize integral 24 24 25 for( int i=0; i< 76; i++) {26 const double sigma = 0.5 * ( G auss76Z[i] + 1.0 );25 for( int i=0; i<GAUSS_N; i++) { 26 const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 27 27 const double mu_proj = mu * sqrt(1.0-sigma*sigma); 28 28 … … 30 30 // corresponding to angles from 0 to pi/2. 31 31 double inner_total = 0.0; 32 for(int j=0; j< 76; j++) {33 const double uu = 0.5 * ( G auss76Z[j] + 1.0 );32 for(int j=0; j<GAUSS_N; j++) { 33 const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 34 34 double sin_uu, cos_uu; 35 35 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 36 36 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 37 37 const double si2 = sas_sinx_x(mu_proj * cos_uu); 38 inner_total += G auss76Wt[j] * square(si1 * si2);38 inner_total += GAUSS_W[j] * square(si1 * si2); 39 39 } 40 40 inner_total *= 0.5; 41 41 42 42 const double si = sas_sinx_x(mu * c_scaled * sigma); 43 outer_total += G auss76Wt[i] * inner_total * si * si;43 outer_total += GAUSS_W[i] * inner_total * si * si; 44 44 } 45 45 outer_total *= 0.5; -
sasmodels/models/pringle.c
r1e7b0db0 r74768cb 29 29 double sumC = 0.0; // initialize integral 30 30 double r; 31 for (int i=0; i < 76; i++) {32 r = G auss76Z[i]*zm + zb;31 for (int i=0; i < GAUSS_N; i++) { 32 r = GAUSS_Z[i]*zm + zb; 33 33 34 34 const double qrs = r*q_sin_psi; 35 35 const double qrrc = r*r*q_cos_psi; 36 36 37 double y = G auss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);37 double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 38 38 double S, C; 39 39 SINCOS(alpha*qrrc, S, C); … … 86 86 87 87 double sum = 0.0; 88 for (int i = 0; i < 76; i++) {89 double psi = G auss76Z[i]*zm + zb;88 for (int i = 0; i < GAUSS_N; i++) { 89 double psi = GAUSS_Z[i]*zm + zb; 90 90 double sin_psi, cos_psi; 91 91 SINCOS(psi, sin_psi, cos_psi); … … 93 93 double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 94 94 double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 95 sum += G auss76Wt[i] * pringle_kernel;95 sum += GAUSS_W[i] * pringle_kernel; 96 96 } 97 97 -
sasmodels/models/rectangular_prism.c
r8de1477 r74768cb 28 28 29 29 double outer_sum = 0.0; 30 for(int i=0; i< 76; i++) {31 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );30 for(int i=0; i<GAUSS_N; i++) { 31 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 32 32 double sin_theta, cos_theta; 33 33 SINCOS(theta, sin_theta, cos_theta); … … 36 36 37 37 double inner_sum = 0.0; 38 for(int j=0; j< 76; j++) {39 double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );38 for(int j=0; j<GAUSS_N; j++) { 39 double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 40 40 double sin_phi, cos_phi; 41 41 SINCOS(phi, sin_phi, cos_phi); … … 45 45 const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 46 46 const double AP = termA * termB * termC; 47 inner_sum += G auss76Wt[j] * AP * AP;47 inner_sum += GAUSS_W[j] * AP * AP; 48 48 } 49 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += G auss76Wt[i] * inner_sum * sin_theta;50 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 51 51 } 52 52 -
sasmodels/models/sc_paracrystal.c
rf728001 r74768cb 54 54 55 55 double outer_sum = 0.0; 56 for(int i=0; i< 150; i++) {56 for(int i=0; i<GAUSS_N; i++) { 57 57 double inner_sum = 0.0; 58 const double theta = G auss150Z[i]*theta_m + theta_b;58 const double theta = GAUSS_Z[i]*theta_m + theta_b; 59 59 double sin_theta, cos_theta; 60 60 SINCOS(theta, sin_theta, cos_theta); 61 61 const double qc = q*cos_theta; 62 62 const double qab = q*sin_theta; 63 for(int j=0;j< 150;j++) {64 const double phi = G auss150Z[j]*phi_m + phi_b;63 for(int j=0;j<GAUSS_N;j++) { 64 const double phi = GAUSS_Z[j]*phi_m + phi_b; 65 65 double sin_phi, cos_phi; 66 66 SINCOS(phi, sin_phi, cos_phi); … … 68 68 const double qb = qab*sin_phi; 69 69 const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 70 inner_sum += G auss150Wt[j] * form;70 inner_sum += GAUSS_W[j] * form; 71 71 } 72 72 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 73 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;73 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 74 74 } 75 75 outer_sum *= theta_m; -
sasmodels/models/stacked_disks.c
rbecded3 r74768cb 81 81 double halfheight = 0.5*thick_core; 82 82 83 for(int i=0; i< N_POINTS_76; i++) {84 double zi = (G auss76Z[i] + 1.0)*M_PI_4;83 for(int i=0; i<GAUSS_N; i++) { 84 double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 85 85 double sin_alpha, cos_alpha; // slots to hold sincos function output 86 86 SINCOS(zi, sin_alpha, cos_alpha); … … 95 95 solvent_sld, 96 96 d); 97 summ += G auss76Wt[i] * yyy * sin_alpha;97 summ += GAUSS_W[i] * yyy * sin_alpha; 98 98 } 99 99 -
sasmodels/models/triaxial_ellipsoid.c
rbecded3 r74768cb 21 21 const double zb = M_PI_4; 22 22 double outer = 0.0; 23 for (int i=0;i< 76;i++) {24 //const double u = G auss76Z[i]*(upper-lower)/2 + (upper + lower)/2;25 const double phi = G auss76Z[i]*zm + zb;23 for (int i=0;i<GAUSS_N;i++) { 24 //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 25 const double phi = GAUSS_Z[i]*zm + zb; 26 26 const double pa_sinsq_phi = pa*square(sin(phi)); 27 27 … … 29 29 const double um = 0.5; 30 30 const double ub = 0.5; 31 for (int j=0;j< 76;j++) {31 for (int j=0;j<GAUSS_N;j++) { 32 32 // translate a point in [-1,1] to a point in [0, 1] 33 const double usq = square(G auss76Z[j]*um + ub);33 const double usq = square(GAUSS_Z[j]*um + ub); 34 34 const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 35 35 const double fq = sas_3j1x_x(q*r); 36 inner += G auss76Wt[j] * fq * fq;36 inner += GAUSS_W[j] * fq * fq; 37 37 } 38 outer += G auss76Wt[i] * inner; // correcting for dx later38 outer += GAUSS_W[i] * inner; // correcting for dx later 39 39 } 40 40 // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
Note: See TracChangeset
for help on using the changeset viewer.