Changeset 994d77f in sasmodels for sasmodels/models/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/cylinder.c
r0496031 r994d77f 1 real form_volume(real radius, reallength);2 real Iq(real q, real sld, real solvent_sld, real radius, reallength);3 real Iqxy(real qx, real qy, real sld, realsolvent_sld,4 real radius, real length, real theta, realphi);1 double form_volume(double radius, double length); 2 double Iq(double q, double sld, double solvent_sld, double radius, double length); 3 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 double radius, double length, double theta, double phi); 5 5 6 6 // twovd = 2 * volume * delta_rho 7 7 // besarg = q * R * sin(alpha) 8 8 // siarg = q * L/2 * cos(alpha) 9 real _cyl(real twovd, real besarg, realsiarg);10 real _cyl(real twovd, real besarg, realsiarg)9 double _cyl(double twovd, double besarg, double siarg); 10 double _cyl(double twovd, double besarg, double siarg) 11 11 { 12 const real bj = (besarg == REAL(0.0) ? REAL(0.5): J1(besarg)/besarg);13 const real si = (siarg == REAL(0.0) ? REAL(1.0): sin(siarg)/siarg);12 const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 13 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 14 14 return twovd*si*bj; 15 15 } 16 16 17 real form_volume(real radius, reallength)17 double form_volume(double radius, double length) 18 18 { 19 19 return M_PI*radius*radius*length; 20 20 } 21 21 22 real Iq(realq,23 realsld,24 realsolvent_sld,25 realradius,26 reallength)22 double Iq(double q, 23 double sld, 24 double solvent_sld, 25 double radius, 26 double length) 27 27 { 28 const realqr = q*radius;29 const real qh = q*REAL(0.5)*length;30 const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);31 real total = REAL(0.0);32 // reallower=0, upper=M_PI_2;28 const double qr = q*radius; 29 const double qh = q*0.5*length; 30 const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length); 31 double total = 0.0; 32 // double lower=0, upper=M_PI_2; 33 33 for (int i=0; i<76 ;i++) { 34 34 // translate a point in [-1,1] to a point in [lower,upper] 35 //const realalpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;36 const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);37 realsn, cn;35 //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 36 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 37 double sn, cn; 38 38 SINCOS(alpha, sn, cn); 39 const realfq = _cyl(twovd, qr*sn, qh*cn);39 const double fq = _cyl(twovd, qr*sn, qh*cn); 40 40 total += Gauss76Wt[i] * fq * fq * sn; 41 41 } 42 42 // translate dx in [-1,1] to dx in [lower,upper] 43 //const realform = (upper-lower)/2.0*total;44 return REAL(1.0e-4)* total * M_PI_4;43 //const double form = (upper-lower)/2.0*total; 44 return 1.0e-4 * total * M_PI_4; 45 45 } 46 46 47 47 48 real Iqxy(real qx, realqy,49 realsld,50 realsolvent_sld,51 realradius,52 reallength,53 realtheta,54 realphi)48 double Iqxy(double qx, double qy, 49 double sld, 50 double solvent_sld, 51 double radius, 52 double length, 53 double theta, 54 double phi) 55 55 { 56 56 // TODO: check that radius<0 and length<0 give zero scattering. 57 57 // This should be the case since the polydispersity weight vector should 58 58 // be zero length, and this function never called. 59 realsn, cn; // slots to hold sincos function output59 double sn, cn; // slots to hold sincos function output 60 60 61 61 // Compute angle alpha between q and the cylinder axis 62 62 SINCOS(theta*M_PI_180, sn, cn); 63 const realq = sqrt(qx*qx+qy*qy);64 const realcos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);65 const realalpha = acos(cos_val);63 const double q = sqrt(qx*qx+qy*qy); 64 const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 65 const double alpha = acos(cos_val); 66 66 67 const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);67 const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length); 68 68 SINCOS(alpha, sn, cn); 69 const real fq = _cyl(twovd, q*radius*sn, q*REAL(0.5)*length*cn);70 return REAL(1.0e-4)* fq * fq;69 const double fq = _cyl(twovd, q*radius*sn, q*0.5*length*cn); 70 return 1.0e-4 * fq * fq; 71 71 }
Note: See TracChangeset
for help on using the changeset viewer.