Changeset 994d77f in sasmodels for sasmodels/models/core_shell_cylinder.c
- Timestamp:
- Oct 30, 2014 10:33:53 AM (9 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:
- ef2861b
- Parents:
- d087487b
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_cylinder.c
r5d4777d r994d77f 1 real form_volume(real radius, real thickness, reallength);2 real Iq(real q, real core_sld, real shell_sld, realsolvent_sld,3 real radius, real thickness, reallength);4 real Iqxy(real qx, real qy, real core_sld, real shell_sld, realsolvent_sld,5 real radius, real thickness, real length, real theta, realphi);1 double form_volume(double radius, double thickness, double length); 2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 3 double radius, double thickness, double length); 4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 5 double radius, double thickness, double length, double theta, double phi); 6 6 7 7 // twovd = 2 * volume * delta_rho 8 8 // besarg = q * R * sin(alpha) 9 9 // siarg = q * L/2 * cos(alpha) 10 real _cyl(real twovd, real besarg, realsiarg);11 real _cyl(real twovd, real besarg, realsiarg)10 double _cyl(double twovd, double besarg, double siarg); 11 double _cyl(double twovd, double besarg, double siarg) 12 12 { 13 const real bj = (besarg == REAL(0.0) ? REAL(0.5): J1(besarg)/besarg);14 const real si = (siarg == REAL(0.0) ? REAL(1.0): sin(siarg)/siarg);13 const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 14 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 15 15 return twovd*si*bj; 16 16 } 17 17 18 real form_volume(real radius, real thickness, reallength)18 double form_volume(double radius, double thickness, double length) 19 19 { 20 20 return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness); 21 21 } 22 22 23 real Iq(realq,24 realcore_sld,25 realshell_sld,26 realsolvent_sld,27 realradius,28 realthickness,29 reallength)23 double Iq(double q, 24 double core_sld, 25 double shell_sld, 26 double solvent_sld, 27 double radius, 28 double thickness, 29 double length) 30 30 { 31 31 // precalculate constants 32 const realcore_qr = q*radius;33 const real core_qh = q*REAL(0.5)*length;34 const real core_twovd = REAL(2.0)* form_volume(radius,0,length)32 const double core_qr = q*radius; 33 const double core_qh = q*0.5*length; 34 const double core_twovd = 2.0 * form_volume(radius,0,length) 35 35 * (core_sld-shell_sld); 36 const realshell_qr = q*(radius + thickness);37 const real shell_qh = q*(REAL(0.5)*length + thickness);38 const real shell_twovd = REAL(2.0)* form_volume(radius,thickness,length)36 const double shell_qr = q*(radius + thickness); 37 const double shell_qh = q*(0.5*length + thickness); 38 const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 39 39 * (shell_sld-solvent_sld); 40 real total = REAL(0.0);41 // reallower=0, upper=M_PI_2;40 double total = 0.0; 41 // double lower=0, upper=M_PI_2; 42 42 for (int i=0; i<76 ;i++) { 43 43 // translate a point in [-1,1] to a point in [lower,upper] 44 //const realalpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;45 realsn, cn;46 const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);44 //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 45 double sn, cn; 46 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 47 47 SINCOS(alpha, sn, cn); 48 const realfq = _cyl(core_twovd, core_qr*sn, core_qh*cn)48 const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 49 49 + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 50 50 total += Gauss76Wt[i] * fq * fq * sn; 51 51 } 52 52 // translate dx in [-1,1] to dx in [lower,upper] 53 //const realform = (upper-lower)/2.0*total;54 return REAL(1.0e-4)* total * M_PI_4;53 //const double form = (upper-lower)/2.0*total; 54 return 1.0e-4 * total * M_PI_4; 55 55 } 56 56 57 57 58 real Iqxy(real qx, realqy,59 realcore_sld,60 realshell_sld,61 realsolvent_sld,62 realradius,63 realthickness,64 reallength,65 realtheta,66 realphi)58 double Iqxy(double qx, double qy, 59 double core_sld, 60 double shell_sld, 61 double solvent_sld, 62 double radius, 63 double thickness, 64 double length, 65 double theta, 66 double phi) 67 67 { 68 realsn, cn; // slots to hold sincos function output68 double sn, cn; // slots to hold sincos function output 69 69 70 70 // Compute angle alpha between q and the cylinder axis … … 72 72 // # The following correction factor exists in sasview, but it can't be 73 73 // # right, so we are leaving it out for now. 74 // const realcorrection = fabs(cn)*M_PI_2;75 const realq = sqrt(qx*qx+qy*qy);76 const realcos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);77 const realalpha = acos(cos_val);74 // const double correction = fabs(cn)*M_PI_2; 75 const double q = sqrt(qx*qx+qy*qy); 76 const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 77 const double alpha = acos(cos_val); 78 78 79 const realcore_qr = q*radius;80 const real core_qh = q*REAL(0.5)*length;81 const real core_twovd = REAL(2.0)* form_volume(radius,0,length)79 const double core_qr = q*radius; 80 const double core_qh = q*0.5*length; 81 const double core_twovd = 2.0 * form_volume(radius,0,length) 82 82 * (core_sld-shell_sld); 83 const realshell_qr = q*(radius + thickness);84 const real shell_qh = q*(REAL(0.5)*length + thickness);85 const real shell_twovd = REAL(2.0)* form_volume(radius,thickness,length)83 const double shell_qr = q*(radius + thickness); 84 const double shell_qh = q*(0.5*length + thickness); 85 const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 86 86 * (shell_sld-solvent_sld); 87 87 88 88 SINCOS(alpha, sn, cn); 89 const realfq = _cyl(core_twovd, core_qr*sn, core_qh*cn)89 const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 90 90 + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 91 return REAL(1.0e-4)* fq * fq;91 return 1.0e-4 * fq * fq; 92 92 }
Note: See TracChangeset
for help on using the changeset viewer.