Changeset 2a0b2b1 in sasmodels for sasmodels/models/core_shell_bicelle_elliptical.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_bicelle_elliptical.c
rdedcf34 r2a0b2b1 17 17 double thick_face, 18 18 double length, 19 double rhoc,20 double rhoh,21 double rhor,22 double rhosolv)19 double sld_core, 20 double sld_face, 21 double sld_rim, 22 double sld_solvent) 23 23 { 24 double si1,si2,be1,be2;25 24 // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 26 25 // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 27 // const double uplim = M_PI_4;28 26 const double halfheight = 0.5*length; 29 //const double va = 0.0;30 //const double vb = 1.0;31 // inner integral limits32 //const double vaj=0.0;33 //const double vbj=M_PI;34 35 27 const double r_major = r_minor * x_core; 36 28 const double r2A = 0.5*(square(r_major) + square(r_minor)); 37 29 const double r2B = 0.5*(square(r_major) - square(r_minor)); 38 const double dr1 = (rhoc-rhoh) *M_PI*r_minor*r_major*(2.0*halfheight);;39 const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);40 const double dr3 = (rhoh-rhor) *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);41 //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);42 //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);43 //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);30 const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 31 const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 32 const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 33 const double dr1 = vol1*(sld_core-sld_face); 34 const double dr2 = vol2*(sld_rim-sld_solvent); 35 const double dr3 = vol3*(sld_face-sld_rim); 44 36 45 37 //initialize integral … … 47 39 for(int i=0;i<76;i++) { 48 40 //setup inner integral over the ellipsoidal cross-section 49 // since we generate these lots of times, why not store them somewhere? 50 //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 51 const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 52 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 53 double inner_sum=0; 54 double sinarg1 = q*halfheight*cos_alpha; 55 double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 56 si1 = sas_sinx_x(sinarg1); 57 si2 = sas_sinx_x(sinarg2); 41 //const double va = 0.0; 42 //const double vb = 1.0; 43 //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 44 const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 45 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 46 const double qab = q*sin_theta; 47 const double qc = q*cos_theta; 48 const double si1 = sas_sinx_x(halfheight*qc); 49 const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 50 double inner_sum=0.0; 58 51 for(int j=0;j<76;j++) { 59 52 //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 60 //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 61 const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 62 const double rr = sqrt(r2A - r2B*cos(beta)); 63 double besarg1 = q*rr*sin_alpha; 64 double besarg2 = q*(rr+thick_rim)*sin_alpha; 65 be1 = sas_2J1x_x(besarg1); 66 be2 = sas_2J1x_x(besarg2); 67 inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 68 dr2*si2*be2 + 69 dr3*si2*be1); 53 // inner integral limits 54 //const double vaj=0.0; 55 //const double vbj=M_PI; 56 //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 57 const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 58 const double rr = sqrt(r2A - r2B*cos(phi)); 59 const double be1 = sas_2J1x_x(rr*qab); 60 const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 61 const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 62 63 inner_sum += Gauss76Wt[j] * fq * fq; 70 64 } 71 65 //now calculate outer integral … … 83 77 double thick_face, 84 78 double length, 85 double rhoc,86 double rhoh,87 double rhor,88 double rhosolv,79 double sld_core, 80 double sld_face, 81 double sld_rim, 82 double sld_solvent, 89 83 double theta, 90 84 double phi, 91 85 double psi) 92 86 { 93 // THIS NEEDS TESTING94 87 double q, xhat, yhat, zhat; 95 88 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 96 const double dr1 = rhoc-rhoh; 97 const double dr2 = rhor-rhosolv; 98 const double dr3 = rhoh-rhor; 89 const double qa = q*xhat; 90 const double qb = q*yhat; 91 const double qc = q*zhat; 92 93 const double dr1 = sld_core-sld_face; 94 const double dr2 = sld_rim-sld_solvent; 95 const double dr3 = sld_face-sld_rim; 99 96 const double r_major = r_minor*x_core; 100 97 const double halfheight = 0.5*length; … … 104 101 105 102 // Compute effective radius in rotated coordinates 106 const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));107 const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)108 + square((r_minor+thick_rim)* yhat));109 const double be1 = sas_2J1x_x( q *r_hat );110 const double be2 = sas_2J1x_x( q *rshell_hat );111 const double si1 = sas_sinx_x( q*halfheight*zhat);112 const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat);113 const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1);114 return 1.0e-4 * Aq;103 const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 104 const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 105 + square((r_minor+thick_rim)*qb)); 106 const double be1 = sas_2J1x_x( qr_hat ); 107 const double be2 = sas_2J1x_x( qrshell_hat ); 108 const double si1 = sas_sinx_x( halfheight*qc ); 109 const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 110 const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1; 111 return 1.0e-4 * fq*fq; 115 112 } 116
Note: See TracChangeset
for help on using the changeset viewer.