Changeset e7678b2 in sasmodels for sasmodels/models/core_shell_bicelle.c
- Timestamp:
- Feb 29, 2016 6:21:55 AM (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:
- 73860b6
- Parents:
- deac08c
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_bicelle.c
r8007311 re7678b2 29 29 } 30 30 31 inline double sinc(double x) {return x==0 ? 1.0 : sin(x)/x;} 32 31 33 static double 32 34 bicelle_kernel(double qq, … … 41 43 double dum) 42 44 { 43 double dr1,dr2,dr3; 44 double besarg1,besarg2; 45 double vol1,vol2,vol3; 46 double sinarg1,sinarg2; 47 double t1,t2,t3; 48 double retval,si1,si2,be1,be2; 45 double si1,si2,be1,be2; 49 46 50 dr1 = rhoc-rhoh; 51 dr2 = rhor-rhosolv; 52 dr3= rhoh-rhor; 53 vol1 = M_PI*rad*rad*(2.0*length); 54 vol2 = M_PI*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 55 vol3= M_PI*(rad)*(rad)*(2.0*length+2.0*facthick); 56 besarg1 = qq*rad*sin(dum); 57 besarg2 = qq*(rad+radthick)*sin(dum); 58 sinarg1 = qq*length*cos(dum); 59 sinarg2 = qq*(length+facthick)*cos(dum); 47 const double dr1 = rhoc-rhoh; 48 const double dr2 = rhor-rhosolv; 49 const double dr3 = rhoh-rhor; 50 const double vol1 = M_PI*rad*rad*(2.0*length); 51 const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 52 const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 53 double sn,cn; 54 SINCOS(dum, sn, cn); 55 double besarg1 = qq*rad*sn; 56 double besarg2 = qq*(rad+radthick)*sn; 57 double sinarg1 = qq*length*cn; 58 double sinarg2 = qq*(length+facthick)*cn; 60 59 61 if(besarg1 == 0) { 62 be1 = 0.5; 63 } else { 64 be1 = J1(besarg1)/besarg1; 65 } 66 if(besarg2 == 0) { 67 be2 = 0.5; 68 } else { 69 be2 = J1(besarg2)/besarg2; 70 } 71 if(sinarg1 == 0) { 72 si1 = 1.0; 73 } else { 74 si1 = sin(sinarg1)/sinarg1; 75 } 76 if(sinarg2 == 0) { 77 si2 = 1.0; 78 } else { 79 si2 = sin(sinarg2)/sinarg2; 80 } 81 t1 = 2.0*vol1*dr1*si1*be1; 82 t2 = 2.0*vol2*dr2*si2*be2; 83 t3 = 2.0*vol3*dr3*si2*be1; 60 be1 = J1c(besarg1); 61 be2 = J1c(besarg2); 62 si1 = sinc(sinarg1); 63 si2 = sinc(sinarg2); 84 64 85 retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 86 return(retval); 65 const double t = vol1*dr1*si1*be1 + 66 vol2*dr2*si2*be2 + 67 vol3*dr3*si2*be1; 68 69 const double retval = t*t*sn; 70 71 return(retval); 87 72 88 73 } … … 99 84 double rhosolv) 100 85 { 86 // set up the integration end points 87 const double uplim = M_PI/4; 88 const double halfheight = length/2.0; 101 89 90 double summ = 0.0; 91 for(int i=0;i<N_POINTS_76;i++) { 92 double zi = (Gauss76Z[i] + 1.0)*uplim; 93 double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 94 halfheight, rhoc, rhoh, rhor,rhosolv, zi); 95 summ += yyy; 96 } 102 97 103 double answer,halfheight; 104 double lolim,uplim,summ,yyy,zi; 105 int nord,i; 106 107 // set up the integration end points 108 nord = 76; 109 lolim = 0.0; 110 uplim = M_PI/2; 111 halfheight = length/2.0; 112 113 summ = 0.0; 114 i=0; 115 for(i=0;i<nord;i++) { 116 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 117 yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 118 halfheight, rhoc, rhoh, rhor,rhosolv, zi); 119 summ += yyy; 120 } 121 122 // calculate value of integral to return 123 answer = (uplim-lolim)/2.0*summ; 124 return(answer); 98 // calculate value of integral to return 99 double answer = uplim*summ; 100 return(answer); 125 101 } 126 102 … … 138 114 double phi) 139 115 { 140 double cyl_x, cyl_y;141 double alpha, cos_val;142 double answer;143 144 116 //convert angle degree to radian 145 theta *= M_PI /180.0;146 phi *= M_PI /180.0;117 theta *= M_PI_180; 118 phi *= M_PI_180; 147 119 148 120 // Cylinder orientation 149 c yl_x = cos(theta) * cos(phi);150 c yl_y = sin(theta);121 const double cyl_x = cos(theta) * cos(phi); 122 const double cyl_y = sin(theta); 151 123 152 124 // Compute the angle btw vector q and the axis of the cylinder 153 co s_val = cyl_x*q_x + cyl_y*q_y;154 alpha = acos( cos_val );125 const double cos_val = cyl_x*q_x + cyl_y*q_y; 126 const double alpha = acos( cos_val ); 155 127 156 128 // Get the kernel 157 answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,129 double answer = bicelle_kernel(q, radius, rim_thickness, face_thickness, 158 130 length/2.0, core_sld, face_sld, rim_sld, 159 131 solvent_sld, alpha) / fabs(sin(alpha));
Note: See TracChangeset
for help on using the changeset viewer.