Changeset 2a0b2b1 in sasmodels for sasmodels/models/core_shell_bicelle.c
- Timestamp:
- Apr 14, 2017 8: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.c
rb260926 r2a0b2b1 1 double form_volume(double radius, double thick_rim, double thick_face, double length); 2 double Iq(double q, 3 double radius, 4 double thick_rim, 5 double thick_face, 6 double length, 7 double core_sld, 8 double face_sld, 9 double rim_sld, 10 double solvent_sld); 11 12 13 double Iqxy(double qx, double qy, 14 double radius, 15 double thick_rim, 16 double thick_face, 17 double length, 18 double core_sld, 19 double face_sld, 20 double rim_sld, 21 double solvent_sld, 22 double theta, 23 double phi); 24 25 26 double form_volume(double radius, double thick_rim, double thick_face, double length) 1 static double 2 form_volume(double radius, double thick_rim, double thick_face, double length) 27 3 { 28 return M_PI* (radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);4 return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 29 5 } 30 6 31 7 static double 32 bicelle_kernel(double q, 33 double rad, 34 double radthick, 35 double facthick, 36 double halflength, 37 double rhoc, 38 double rhoh, 39 double rhor, 40 double rhosolv, 41 double sin_alpha, 42 double cos_alpha) 8 bicelle_kernel(double qab, 9 double qc, 10 double radius, 11 double thick_radius, 12 double thick_face, 13 double halflength, 14 double sld_core, 15 double sld_face, 16 double sld_rim, 17 double sld_solvent) 43 18 { 44 const double dr1 = rhoc-rhoh;45 const double dr2 = rhor-rhosolv;46 const double dr3 = rhoh-rhor;47 const double vol1 = M_PI*square(rad )*2.0*(halflength);48 const double vol2 = M_PI*square(rad +radthick)*2.0*(halflength+facthick);49 const double vol3 = M_PI*square(rad )*2.0*(halflength+facthick);19 const double dr1 = sld_core-sld_face; 20 const double dr2 = sld_rim-sld_solvent; 21 const double dr3 = sld_face-sld_rim; 22 const double vol1 = M_PI*square(radius)*2.0*(halflength); 23 const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 24 const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 50 25 51 const double be1 = sas_2J1x_x( q*(rad)*sin_alpha);52 const double be2 = sas_2J1x_x( q*(rad+radthick)*sin_alpha);53 const double si1 = sas_sinx_x( q*(halflength)*cos_alpha);54 const double si2 = sas_sinx_x( q*(halflength+facthick)*cos_alpha);26 const double be1 = sas_2J1x_x((radius)*qab); 27 const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 28 const double si1 = sas_sinx_x((halflength)*qc); 29 const double si2 = sas_sinx_x((halflength+thick_face)*qc); 55 30 56 31 const double t = vol1*dr1*si1*be1 + … … 58 33 vol3*dr3*si2*be1; 59 34 60 const double retval = t*t; 61 62 return retval; 63 35 return t; 64 36 } 65 37 66 38 static double 67 bicelle_integration(double q,68 double rad,69 double radthick,70 double facthick,71 72 double rhoc,73 double rhoh,74 double rhor,75 double rhosolv)39 Iq(double q, 40 double radius, 41 double thick_radius, 42 double thick_face, 43 double length, 44 double sld_core, 45 double sld_face, 46 double sld_rim, 47 double sld_solvent) 76 48 { 77 49 // set up the integration end points … … 79 51 const double halflength = 0.5*length; 80 52 81 double summ= 0.0;53 double total = 0.0; 82 54 for(int i=0;i<N_POINTS_76;i++) { 83 double alpha = (Gauss76Z[i] + 1.0)*uplim; 84 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 SINCOS(alpha, sin_alpha, cos_alpha); 86 double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 87 halflength, rhoc, rhoh, rhor, rhosolv, 88 sin_alpha, cos_alpha); 89 summ += yyy*sin_alpha; 55 double theta = (Gauss76Z[i] + 1.0)*uplim; 56 double sin_theta, cos_theta; // slots to hold sincos function output 57 SINCOS(theta, sin_theta, cos_theta); 58 double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 59 halflength, sld_core, sld_face, sld_rim, sld_solvent); 60 total += Gauss76Wt[i]*fq*fq*sin_theta; 90 61 } 91 62 92 63 // calculate value of integral to return 93 double answer = uplim*summ;94 return answer;64 double answer = total*uplim; 65 return 1.0e-4*answer; 95 66 } 96 67 97 68 static double 98 bicelle_kernel_2d(double qx, double qy,99 100 101 102 103 104 105 106 107 108 69 Iqxy(double qx, double qy, 70 double radius, 71 double thick_rim, 72 double thick_face, 73 double length, 74 double core_sld, 75 double face_sld, 76 double rim_sld, 77 double solvent_sld, 78 double theta, 79 double phi) 109 80 { 110 81 double q, sin_alpha, cos_alpha; 111 82 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 83 const double qab = q*sin_alpha; 84 const double qc = q*cos_alpha; 112 85 113 double answer = bicelle_kernel(q, radius, thick_rim, thick_face,86 double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 114 87 0.5*length, core_sld, face_sld, rim_sld, 115 solvent_sld , sin_alpha, cos_alpha);116 return 1.0e-4* answer;88 solvent_sld); 89 return 1.0e-4*fq*fq; 117 90 } 118 119 double Iq(double q,120 double radius,121 double thick_rim,122 double thick_face,123 double length,124 double core_sld,125 double face_sld,126 double rim_sld,127 double solvent_sld)128 {129 double intensity = bicelle_integration(q, radius, thick_rim, thick_face,130 length, core_sld, face_sld, rim_sld, solvent_sld);131 return intensity*1.0e-4;132 }133 134 135 double Iqxy(double qx, double qy,136 double radius,137 double thick_rim,138 double thick_face,139 double length,140 double core_sld,141 double face_sld,142 double rim_sld,143 double solvent_sld,144 double theta,145 double phi)146 {147 double intensity = bicelle_kernel_2d(qx, qy,148 radius,149 thick_rim,150 thick_face,151 length,152 core_sld,153 face_sld,154 rim_sld,155 solvent_sld,156 theta,157 phi);158 159 return intensity;160 }
Note: See TracChangeset
for help on using the changeset viewer.