Changeset 5d4777d in sasmodels for sasmodels/models/cylinder_clone.c
- Timestamp:
- Sep 2, 2014 1:24:38 AM (10 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:
- f4cf580
- Parents:
- ff7119b
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/cylinder_clone.c
rff7119b r5d4777d 2 2 real Iq(real q, real sld, real solvent_sld, real radius, real length); 3 3 real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 4 5 6 // twovd = 2 * volume * delta_rho 7 // besarg = q * R * sin(alpha) 8 // siarg = q * L/2 * cos(alpha) 9 real _cyl(real twovd, real besarg, real siarg, real alpha); 10 real _cyl(real twovd, real besarg, real siarg, real alpha) 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); 14 return twovd*si*bj; 15 } 4 16 5 17 real form_volume(real radius, real length) … … 7 19 return M_PI*radius*radius*length; 8 20 } 9 10 21 real Iq(real q, 11 22 real sldCyl, … … 14 25 real length) 15 26 { 16 const real h = REAL(0.5)*length; 17 real summ = REAL(0.0); 27 const real qr = q*radius; 28 const real qh = q*REAL(0.5)*length; 29 const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 30 real total = REAL(0.0); 31 // real lower=0, upper=M_PI_2; 18 32 for (int i=0; i<76 ;i++) { 19 //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 20 const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 21 summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi); 33 // translate a point in [-1,1] to a point in [lower,upper] 34 //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 35 const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 36 real sn, cn; 37 SINCOS(alpha, sn, cn); 38 const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 39 total += Gauss76Wt[i] * fq * fq * sn; 22 40 } 23 //const real form = (uplim-lolim)/2.0*summ; 24 const real form = summ * M_PI_4; 25 26 // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 27 // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 28 // The additional volume factor is for polydisperse volume normalization. 29 const real s = (sldCyl - sldSolv) * form_volume(radius, length); 30 return REAL(1.0e8) * form * s * s; 41 // translate dx in [-1,1] to dx in [lower,upper] 42 //const real form = (upper-lower)/2.0*total; 43 return REAL(1.0e8) * total * M_PI_4; 31 44 } 32 45 … … 50 63 const real alpha = acos(cos_val); 51 64 52 // The following is CylKernel() / sin(alpha), but we are doing it in place53 // to avoid sin(alpha)/sin(alpha) for alpha = 0. It is also a teensy bit54 // faster since we don't mulitply and divide sin(alpha).65 const real qr = q*radius; 66 const real qh = q*REAL(0.5)*length; 67 const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 55 68 SINCOS(alpha, sn, cn); 56 const real besarg = q*radius*sn; 57 const real siarg = REAL(0.5)*q*length*cn; 58 // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 59 const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 60 const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 61 const real form = REAL(4.0)*bj*bj*si*si; 62 63 // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 64 // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 65 // The additional volume factor is for polydisperse volume normalization. 66 const real s = (sldCyl - sldSolv) * form_volume(radius, length); 67 return REAL(1.0e8) * form * s * s * spherical_integration; 69 const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 70 return REAL(1.0e8) * fq * fq * spherical_integration; 68 71 }
Note: See TracChangeset
for help on using the changeset viewer.