Changeset 994d77f in sasmodels for sasmodels/models
- Timestamp:
- Oct 30, 2014 12:33:53 PM (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:
- ef2861b
- Parents:
- d087487b
- Location:
- sasmodels/models
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/capped_cylinder.c
rf4cf580 r994d77f 1 real form_volume(real radius, real cap_radius, reallength);2 real Iq(real q, real sld, realsolvent_sld,3 real radius, real cap_radius, reallength);4 real Iqxy(real qx, real qy, real sld, realsolvent_sld,5 real radius, real cap_radius, real length, real theta, realphi);1 double form_volume(double radius, double cap_radius, double length); 2 double Iq(double q, double sld, double solvent_sld, 3 double radius, double cap_radius, double length); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 5 double radius, double cap_radius, double length, double theta, double phi); 6 6 7 7 // Integral over a convex lens kernel for t in [h/R,1]. See the docs for … … 13 13 // length is the cylinder length, or the separation between the lens halves 14 14 // alpha is the angle of the cylinder wrt q. 15 real _cap_kernel(real q, real h, real cap_radius, reallength,16 real sin_alpha, realcos_alpha);17 real _cap_kernel(real q, real h, real cap_radius, reallength,18 real sin_alpha, realcos_alpha)15 double _cap_kernel(double q, double h, double cap_radius, double length, 16 double sin_alpha, double cos_alpha); 17 double _cap_kernel(double q, double h, double cap_radius, double length, 18 double sin_alpha, double cos_alpha) 19 19 { 20 20 // For speed, we are pre-calculating terms which are constant over the 21 21 // kernel. 22 const real upper = REAL(1.0);23 const reallower = h/cap_radius; // integral lower bound22 const double upper = 1.0; 23 const double lower = h/cap_radius; // integral lower bound 24 24 // cos term in integral is: 25 25 // cos (q (R t - h + L/2) cos(alpha)) … … 29 29 // m = q R cos(alpha) 30 30 // b = q(L/2-h) cos(alpha) 31 const realm = q*cap_radius*cos_alpha; // cos argument slope32 const real b = q*(REAL(0.5)*length-h)*cos_alpha; // cos argument intercept33 const realqrst = q*cap_radius*sin_alpha; // Q*R*sin(theta)34 real total = REAL(0.0);31 const double m = q*cap_radius*cos_alpha; // cos argument slope 32 const double b = q*(0.5*length-h)*cos_alpha; // cos argument intercept 33 const double qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta) 34 double total = 0.0; 35 35 for (int i=0; i<76 ;i++) { 36 36 // translate a point in [-1,1] to a point in [lower,upper] 37 //const realt = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;38 const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower);39 const real radical = REAL(1.0)- t*t;40 const realarg = qrst*sqrt(radical); // cap bessel function arg41 const real be = (arg == REAL(0.0) ? REAL(0.5): J1(arg)/arg);42 const realFq = cos(m*t + b) * radical * be;37 //const double t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 38 const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower); 39 const double radical = 1.0 - t*t; 40 const double arg = qrst*sqrt(radical); // cap bessel function arg 41 const double be = (arg == 0.0 ? 0.5 : J1(arg)/arg); 42 const double Fq = cos(m*t + b) * radical * be; 43 43 total += Gauss76Wt[i] * Fq; 44 44 } 45 45 // translate dx in [-1,1] to dx in [lower,upper] 46 //const realform = (upper-lower)/2.0*total;47 const real integral = REAL(0.5)*(upper-lower)*total;48 return REAL(4.0)*M_PI*cap_radius*cap_radius*cap_radius*integral;46 //const double form = (upper-lower)/2.0*total; 47 const double integral = 0.5*(upper-lower)*total; 48 return 4.0*M_PI*cap_radius*cap_radius*cap_radius*integral; 49 49 } 50 50 51 real form_volume(real radius, real cap_radius, reallength)51 double form_volume(double radius, double cap_radius, double length) 52 52 { 53 53 // cap radius should never be less than radius when this is called … … 60 60 // = pi r^2 L + pi hc (r^2 + hc^2/3) 61 61 // = pi * (r^2 (L+hc) + hc^3/3) 62 const realhc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius);63 return M_PI*(radius*radius*(length+hc) + REAL(0.333333333333333)*hc*hc*hc);62 const double hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius); 63 return M_PI*(radius*radius*(length+hc) + 0.333333333333333*hc*hc*hc); 64 64 } 65 65 66 real Iq(realq,67 realsld,68 realsolvent_sld,69 realradius,70 realcap_radius,71 reallength)66 double Iq(double q, 67 double sld, 68 double solvent_sld, 69 double radius, 70 double cap_radius, 71 double length) 72 72 { 73 realsn, cn; // slots to hold sincos function output73 double sn, cn; // slots to hold sincos function output 74 74 75 75 // Exclude invalid inputs. 76 if (cap_radius < radius) return REAL(-1.0);76 if (cap_radius < radius) return -1.0; 77 77 78 const real lower = REAL(0.0);79 const realupper = M_PI_2;80 const realh = sqrt(cap_radius*cap_radius - radius*radius); // negative h81 real total = REAL(0.0);78 const double lower = 0.0; 79 const double upper = M_PI_2; 80 const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 81 double total = 0.0; 82 82 for (int i=0; i<76 ;i++) { 83 83 // translate a point in [-1,1] to a point in [lower,upper] 84 const real alpha= REAL(0.5)*(Gauss76Z[i]*(upper-lower) + upper + lower);84 const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower); 85 85 SINCOS(alpha, sn, cn); 86 86 87 const realcap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);87 const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 88 88 89 89 // The following is CylKernel() / sin(alpha), but we are doing it in place 90 90 // to avoid sin(alpha)/sin(alpha) for alpha = 0. It is also a teensy bit 91 91 // faster since we don't multiply and divide sin(alpha). 92 const realbesarg = q*radius*sn;93 const real siarg = q*REAL(0.5)*length*cn;92 const double besarg = q*radius*sn; 93 const double siarg = q*0.5*length*cn; 94 94 // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 95 const real bj = (besarg == REAL(0.0) ? REAL(0.5): J1(besarg)/besarg);96 const real si = (siarg == REAL(0.0) ? REAL(1.0): sin(siarg)/siarg);97 const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si;95 const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 96 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 97 const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 98 98 99 99 // Volume weighted average F(q) 100 const realAq = cyl_Fq + cap_Fq;100 const double Aq = cyl_Fq + cap_Fq; 101 101 total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration 102 102 } 103 103 // translate dx in [-1,1] to dx in [lower,upper] 104 const real form = total * REAL(0.5)*(upper-lower);104 const double form = total * 0.5*(upper-lower); 105 105 106 106 // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 107 107 // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 108 108 // The additional volume factor is for polydisperse volume normalization. 109 const reals = (sld - solvent_sld);110 return REAL(1.0e-4)* form * s * s; // form_volume(radius, cap_radius, length);109 const double s = (sld - solvent_sld); 110 return 1.0e-4 * form * s * s; // form_volume(radius, cap_radius, length); 111 111 } 112 112 113 113 114 real Iqxy(real qx, realqy,115 realsld,116 realsolvent_sld,117 realradius,118 realcap_radius,119 reallength,120 realtheta,121 realphi)114 double Iqxy(double qx, double qy, 115 double sld, 116 double solvent_sld, 117 double radius, 118 double cap_radius, 119 double length, 120 double theta, 121 double phi) 122 122 { 123 realsn, cn; // slots to hold sincos function output123 double sn, cn; // slots to hold sincos function output 124 124 125 125 // Exclude invalid inputs. 126 if (cap_radius < radius) return REAL(-1.0);126 if (cap_radius < radius) return -1.0; 127 127 128 128 // Compute angle alpha between q and the cylinder axis … … 130 130 // # The following correction factor exists in sasview, but it can't be 131 131 // # right, so we are leaving it out for now. 132 const realq = sqrt(qx*qx+qy*qy);133 const realcos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);134 const realalpha = acos(cos_val); // rod angle relative to q132 const double q = sqrt(qx*qx+qy*qy); 133 const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 134 const double alpha = acos(cos_val); // rod angle relative to q 135 135 SINCOS(alpha, sn, cn); 136 136 137 const realh = sqrt(cap_radius*cap_radius - radius*radius); // negative h138 const realcap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);137 const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 138 const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 139 139 140 const realbesarg = q*radius*sn;141 const real siarg = q*REAL(0.5)*length*cn;140 const double besarg = q*radius*sn; 141 const double siarg = q*0.5*length*cn; 142 142 // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 143 const real bj = (besarg == REAL(0.0) ? REAL(0.5): J1(besarg)/besarg);144 const real si = (siarg == REAL(0.0) ? REAL(1.0): sin(siarg)/siarg);145 const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si;143 const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 144 const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 145 const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 146 146 147 147 // Volume weighted average F(q) 148 const realAq = cyl_Fq + cap_Fq;148 const double Aq = cyl_Fq + cap_Fq; 149 149 150 150 // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 151 const reals = (sld - solvent_sld);152 return REAL(1.0e-4)* Aq * Aq * s * s; // form_volume(radius, cap_radius, length);151 const double s = (sld - solvent_sld); 152 return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); 153 153 } -
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 } -
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 } -
sasmodels/models/cylinder_clone.c
r5d4777d 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, real solvent_sld, 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, double radius, double length, double theta, double phi); 4 4 5 5 … … 7 7 // besarg = q * R * sin(alpha) 8 8 // siarg = q * L/2 * cos(alpha) 9 real _cyl(real twovd, real besarg, real siarg, realalpha);10 real _cyl(real twovd, real besarg, real siarg, realalpha)9 double _cyl(double twovd, double besarg, double siarg, double alpha); 10 double _cyl(double twovd, double besarg, double siarg, double alpha) 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 real Iq(realq,22 realsldCyl,23 realsldSolv,24 realradius,25 reallength)21 double Iq(double q, 22 double sldCyl, 23 double sldSolv, 24 double radius, 25 double length) 26 26 { 27 const realqr = 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 // reallower=0, upper=M_PI_2;27 const double qr = q*radius; 28 const double qh = q*0.5*length; 29 const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length); 30 double total = 0.0; 31 // double lower=0, upper=M_PI_2; 32 32 for (int i=0; i<76 ;i++) { 33 33 // translate a point in [-1,1] to a point in [lower,upper] 34 //const realalpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;35 const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);36 realsn, cn;34 //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 35 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 36 double sn, cn; 37 37 SINCOS(alpha, sn, cn); 38 const realfq = _cyl(twovd, qr*sn, qh*cn, alpha);38 const double fq = _cyl(twovd, qr*sn, qh*cn, alpha); 39 39 total += Gauss76Wt[i] * fq * fq * sn; 40 40 } 41 41 // translate dx in [-1,1] to dx in [lower,upper] 42 //const realform = (upper-lower)/2.0*total;43 return REAL(1.0e8)* total * M_PI_4;42 //const double form = (upper-lower)/2.0*total; 43 return 1.0e8 * total * M_PI_4; 44 44 } 45 45 46 real Iqxy(real qx, realqy,47 realsldCyl,48 realsldSolv,49 realradius,50 reallength,51 realcyl_theta,52 realcyl_phi)46 double Iqxy(double qx, double qy, 47 double sldCyl, 48 double sldSolv, 49 double radius, 50 double length, 51 double cyl_theta, 52 double cyl_phi) 53 53 { 54 realsn, cn; // slots to hold sincos function output54 double sn, cn; // slots to hold sincos function output 55 55 56 56 // Compute angle alpha between q and the cylinder axis … … 58 58 // # The following correction factor exists in sasview, but it can't be 59 59 // # right, so we are leaving it out for now. 60 const realspherical_integration = fabs(cn)*M_PI_2;61 const realq = sqrt(qx*qx+qy*qy);62 const realcos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q);63 const realalpha = acos(cos_val);60 const double spherical_integration = fabs(cn)*M_PI_2; 61 const double q = sqrt(qx*qx+qy*qy); 62 const double cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q); 63 const double alpha = acos(cos_val); 64 64 65 const realqr = q*radius;66 const real qh = q*REAL(0.5)*length;67 const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length);65 const double qr = q*radius; 66 const double qh = q*0.5*length; 67 const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length); 68 68 SINCOS(alpha, sn, cn); 69 const realfq = _cyl(twovd, qr*sn, qh*cn, alpha);70 return REAL(1.0e8)* fq * fq * spherical_integration;69 const double fq = _cyl(twovd, qr*sn, qh*cn, alpha); 70 return 1.0e8 * fq * fq * spherical_integration; 71 71 } -
sasmodels/models/ellipsoid.c
r5d4777d r994d77f 1 real form_volume(real rpolar, realrequatorial);2 real Iq(real q, real sld, real solvent_sld, real rpolar, realrequatorial);3 real Iqxy(real qx, real qy, real sld, realsolvent_sld,4 real rpolar, real requatorial, real theta, realphi);1 double form_volume(double rpolar, double requatorial); 2 double Iq(double q, double sld, double solvent_sld, double rpolar, double requatorial); 3 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 double rpolar, double requatorial, double theta, double phi); 5 5 6 real _ellipsoid_kernel(real q, real rpolar, real requatorial, realcos_alpha);7 real _ellipsoid_kernel(real q, real rpolar, real requatorial, realcos_alpha)6 double _ellipsoid_kernel(double q, double rpolar, double requatorial, double cos_alpha); 7 double _ellipsoid_kernel(double q, double rpolar, double requatorial, double cos_alpha) 8 8 { 9 realsn, cn;10 realratio = rpolar/requatorial;11 const real u = q*requatorial*sqrt(REAL(1.0)12 + cos_alpha*cos_alpha*(ratio*ratio - REAL(1.0)));9 double sn, cn; 10 double ratio = rpolar/requatorial; 11 const double u = q*requatorial*sqrt(1.0 12 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 13 13 SINCOS(u, sn, cn); 14 const real f = ( u==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-u*cn)/(u*u*u) );14 const double f = ( u==0.0 ? 1.0 : 3.0*(sn-u*cn)/(u*u*u) ); 15 15 return f*f; 16 16 } 17 17 18 real form_volume(real rpolar, realrequatorial)18 double form_volume(double rpolar, double requatorial) 19 19 { 20 return REAL(1.333333333333333)*M_PI*rpolar*requatorial*requatorial;20 return 1.333333333333333*M_PI*rpolar*requatorial*requatorial; 21 21 } 22 22 23 real Iq(realq,24 realsld,25 realsolvent_sld,26 realrpolar,27 realrequatorial)23 double Iq(double q, 24 double sld, 25 double solvent_sld, 26 double rpolar, 27 double requatorial) 28 28 { 29 //const real lower = REAL(0.0);30 //const real upper = REAL(1.0);31 real total = REAL(0.0);29 //const double lower = 0.0; 30 //const double upper = 1.0; 31 double total = 0.0; 32 32 for (int i=0;i<76;i++) { 33 //const realcos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;34 const real cos_alpha = REAL(0.5)*(Gauss76Z[i] + REAL(1.0));33 //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 34 const double cos_alpha = 0.5*(Gauss76Z[i] + 1.0); 35 35 total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 36 36 } 37 //const realform = (upper-lower)/2*total;38 const real form = REAL(0.5)*total;39 const reals = (sld - solvent_sld) * form_volume(rpolar, requatorial);40 return REAL(1.0e-4)* form * s * s;37 //const double form = (upper-lower)/2*total; 38 const double form = 0.5*total; 39 const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 40 return 1.0e-4 * form * s * s; 41 41 } 42 42 43 real Iqxy(real qx, realqy,44 realsld,45 realsolvent_sld,46 realrpolar,47 realrequatorial,48 realtheta,49 realphi)43 double Iqxy(double qx, double qy, 44 double sld, 45 double solvent_sld, 46 double rpolar, 47 double requatorial, 48 double theta, 49 double phi) 50 50 { 51 realsn, cn;51 double sn, cn; 52 52 53 const realq = sqrt(qx*qx + qy*qy);53 const double q = sqrt(qx*qx + qy*qy); 54 54 SINCOS(theta*M_PI_180, sn, cn); 55 const realcos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);56 const realform = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha);57 const reals = (sld - solvent_sld) * form_volume(rpolar, requatorial);55 const double cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 56 const double form = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 57 const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 58 58 59 return REAL(1.0e-4)* form * s * s;59 return 1.0e-4 * form * s * s; 60 60 } 61 61 -
sasmodels/models/lamellar.py
r19dcb933 r994d77f 77 77 # This should perhaps be volume normalized? 78 78 form_volume = """ 79 return REAL(1.0);79 return 1.0; 80 80 """ 81 81 82 82 Iq = """ 83 const realsub = sld - solvent_sld;84 const realqsq = q*q;85 return REAL(4.0e-4)*M_PI*sub*sub/qsq * (REAL(1.0)-cos(q*thickness))83 const double sub = sld - solvent_sld; 84 const double qsq = q*q; 85 return 4.0e-4*M_PI*sub*sub/qsq * (1.0-cos(q*thickness)) 86 86 / (thickness*qsq); 87 87 """ … … 89 89 Iqxy = """ 90 90 // never called since no orientation or magnetic parameters. 91 return REAL(-1.0);91 return -1.0; 92 92 """ 93 93 -
sasmodels/models/lib/J1.c
r14de349 r994d77f 1 real J1(realx);2 real J1(realx)1 double J1(double x); 2 double J1(double x) 3 3 { 4 const realax = fabs(x);5 if (ax < REAL(8.0)) {6 const realy = x*x;7 const real ans1 = x*(REAL(72362614232.0)8 + y*( REAL(-7895059235.0)9 + y*( REAL(242396853.1)10 + y*( REAL(-2972611.439)11 + y*( REAL(15704.48260)12 + y*( REAL(-30.16036606)))))));13 const real ans2 = REAL(144725228442.0)14 + y*( REAL(2300535178.0)15 + y*( REAL(18583304.74)16 + y*( REAL(99447.43394)17 + y*( REAL(376.9991397)4 const double ax = fabs(x); 5 if (ax < 8.0) { 6 const double y = x*x; 7 const double ans1 = x*(72362614232.0 8 + y*(-7895059235.0 9 + y*(242396853.1 10 + y*(-2972611.439 11 + y*(15704.48260 12 + y*(-30.16036606)))))); 13 const double ans2 = 144725228442.0 14 + y*(2300535178.0 15 + y*(18583304.74 16 + y*(99447.43394 17 + y*(376.9991397 18 18 + y)))); 19 19 return ans1/ans2; 20 20 } else { 21 const real y = REAL(64.0)/(ax*ax);22 const real xx = ax - REAL(2.356194491);23 const real ans1 = REAL(1.0)24 + y*( REAL(0.183105e-2)25 + y*( REAL(-0.3516396496e-4)26 + y*( REAL(0.2457520174e-5)27 + y* REAL(-0.240337019e-6))));28 const real ans2 = REAL(0.04687499995)29 + y*( REAL(-0.2002690873e-3)30 + y*( REAL(0.8449199096e-5)31 + y*( REAL(-0.88228987e-6)32 + y* REAL(0.105787412e-6))));33 realsn,cn;21 const double y = 64.0/(ax*ax); 22 const double xx = ax - 2.356194491; 23 const double ans1 = 1.0 24 + y*(0.183105e-2 25 + y*(-0.3516396496e-4 26 + y*(0.2457520174e-5 27 + y*-0.240337019e-6))); 28 const double ans2 = 0.04687499995 29 + y*(-0.2002690873e-3 30 + y*(0.8449199096e-5 31 + y*(-0.88228987e-6 32 + y*0.105787412e-6))); 33 double sn,cn; 34 34 SINCOS(xx, sn, cn); 35 const real ans = sqrt(REAL(0.636619772)/ax) * (cn*ans1 - (REAL(8.0)/ax)*sn*ans2);36 return (x < REAL(0.0)) ? -ans : ans;35 const double ans = sqrt(0.636619772/ax) * (cn*ans1 - (8.0/ax)*sn*ans2); 36 return (x < 0.0) ? -ans : ans; 37 37 } 38 38 } -
sasmodels/models/lib/gauss150.c
r14de349 r994d77f 9 9 10 10 // Gaussians 11 constant realGauss150Z[150]={12 REAL(-0.9998723404457334),13 REAL(-0.9993274305065947),14 REAL(-0.9983473449340834),15 REAL(-0.9969322929775997),16 REAL(-0.9950828645255290),17 REAL(-0.9927998590434373),18 REAL(-0.9900842691660192),19 REAL(-0.9869372772712794),20 REAL(-0.9833602541697529),21 REAL(-0.9793547582425894),22 REAL(-0.9749225346595943),23 REAL(-0.9700655145738374),24 REAL(-0.9647858142586956),25 REAL(-0.9590857341746905),26 REAL(-0.9529677579610971),27 REAL(-0.9464345513503147),28 REAL(-0.9394889610042837),29 REAL(-0.9321340132728527),30 REAL(-0.9243729128743136),31 REAL(-0.9162090414984952),32 REAL(-0.9076459563329236),33 REAL(-0.8986873885126239),34 REAL(-0.8893372414942055),35 REAL(-0.8795995893549102),36 REAL(-0.8694786750173527),37 REAL(-0.8589789084007133),38 REAL(-0.8481048644991847),39 REAL(-0.8368612813885015),40 REAL(-0.8252530581614230),41 REAL(-0.8132852527930605),42 REAL(-0.8009630799369827),43 REAL(-0.7882919086530552),44 REAL(-0.7752772600680049),45 REAL(-0.7619248049697269),46 REAL(-0.7482403613363824),47 REAL(-0.7342298918013638),48 REAL(-0.7198995010552305),49 REAL(-0.7052554331857488),50 REAL(-0.6903040689571928),51 REAL(-0.6750519230300931),52 REAL(-0.6595056411226444),53 REAL(-0.6436719971150083),54 REAL(-0.6275578900977726),55 REAL(-0.6111703413658551),56 REAL(-0.5945164913591590),57 REAL(-0.5776035965513142),58 REAL(-0.5604390262878617),59 REAL(-0.5430302595752546),60 REAL(-0.5253848818220803),61 REAL(-0.5075105815339176),62 REAL(-0.4894151469632753),63 REAL(-0.4711064627160663),64 REAL(-0.4525925063160997),65 REAL(-0.4338813447290861),66 REAL(-0.4149811308476706),67 REAL(-0.3959000999390257),68 REAL(-0.3766465660565522),69 REAL(-0.3572289184172501),70 REAL(-0.3376556177463400),71 REAL(-0.3179351925907259),72 REAL(-0.2980762356029071),73 REAL(-0.2780873997969574),74 REAL(-0.2579773947782034),75 REAL(-0.2377549829482451),76 REAL(-0.2174289756869712),77 REAL(-0.1970082295132342),78 REAL(-0.1765016422258567),79 REAL(-0.1559181490266516),80 REAL(-0.1352667186271445),81 REAL(-0.1145563493406956),82 REAL(-0.0937960651617229),83 REAL(-0.0729949118337358),84 REAL(-0.0521619529078925),85 REAL(-0.0313062657937972),86 REAL(-0.0104369378042598),87 REAL(0.0104369378042598),88 REAL(0.0313062657937972),89 REAL(0.0521619529078925),90 REAL(0.0729949118337358),91 REAL(0.0937960651617229),92 REAL(0.1145563493406956),93 REAL(0.1352667186271445),94 REAL(0.1559181490266516),95 REAL(0.1765016422258567),96 REAL(0.1970082295132342),97 REAL(0.2174289756869712),98 REAL(0.2377549829482451),99 REAL(0.2579773947782034),100 REAL(0.2780873997969574),101 REAL(0.2980762356029071),102 REAL(0.3179351925907259),103 REAL(0.3376556177463400),104 REAL(0.3572289184172501),105 REAL(0.3766465660565522),106 REAL(0.3959000999390257),107 REAL(0.4149811308476706),108 REAL(0.4338813447290861),109 REAL(0.4525925063160997),110 REAL(0.4711064627160663),111 REAL(0.4894151469632753),112 REAL(0.5075105815339176),113 REAL(0.5253848818220803),114 REAL(0.5430302595752546),115 REAL(0.5604390262878617),116 REAL(0.5776035965513142),117 REAL(0.5945164913591590),118 REAL(0.6111703413658551),119 REAL(0.6275578900977726),120 REAL(0.6436719971150083),121 REAL(0.6595056411226444),122 REAL(0.6750519230300931),123 REAL(0.6903040689571928),124 REAL(0.7052554331857488),125 REAL(0.7198995010552305),126 REAL(0.7342298918013638),127 REAL(0.7482403613363824),128 REAL(0.7619248049697269),129 REAL(0.7752772600680049),130 REAL(0.7882919086530552),131 REAL(0.8009630799369827),132 REAL(0.8132852527930605),133 REAL(0.8252530581614230),134 REAL(0.8368612813885015),135 REAL(0.8481048644991847),136 REAL(0.8589789084007133),137 REAL(0.8694786750173527),138 REAL(0.8795995893549102),139 REAL(0.8893372414942055),140 REAL(0.8986873885126239),141 REAL(0.9076459563329236),142 REAL(0.9162090414984952),143 REAL(0.9243729128743136),144 REAL(0.9321340132728527),145 REAL(0.9394889610042837),146 REAL(0.9464345513503147),147 REAL(0.9529677579610971),148 REAL(0.9590857341746905),149 REAL(0.9647858142586956),150 REAL(0.9700655145738374),151 REAL(0.9749225346595943),152 REAL(0.9793547582425894),153 REAL(0.9833602541697529),154 REAL(0.9869372772712794),155 REAL(0.9900842691660192),156 REAL(0.9927998590434373),157 REAL(0.9950828645255290),158 REAL(0.9969322929775997),159 REAL(0.9983473449340834),160 REAL(0.9993274305065947),161 REAL(0.9998723404457334)11 constant double Gauss150Z[150]={ 12 -0.9998723404457334, 13 -0.9993274305065947, 14 -0.9983473449340834, 15 -0.9969322929775997, 16 -0.9950828645255290, 17 -0.9927998590434373, 18 -0.9900842691660192, 19 -0.9869372772712794, 20 -0.9833602541697529, 21 -0.9793547582425894, 22 -0.9749225346595943, 23 -0.9700655145738374, 24 -0.9647858142586956, 25 -0.9590857341746905, 26 -0.9529677579610971, 27 -0.9464345513503147, 28 -0.9394889610042837, 29 -0.9321340132728527, 30 -0.9243729128743136, 31 -0.9162090414984952, 32 -0.9076459563329236, 33 -0.8986873885126239, 34 -0.8893372414942055, 35 -0.8795995893549102, 36 -0.8694786750173527, 37 -0.8589789084007133, 38 -0.8481048644991847, 39 -0.8368612813885015, 40 -0.8252530581614230, 41 -0.8132852527930605, 42 -0.8009630799369827, 43 -0.7882919086530552, 44 -0.7752772600680049, 45 -0.7619248049697269, 46 -0.7482403613363824, 47 -0.7342298918013638, 48 -0.7198995010552305, 49 -0.7052554331857488, 50 -0.6903040689571928, 51 -0.6750519230300931, 52 -0.6595056411226444, 53 -0.6436719971150083, 54 -0.6275578900977726, 55 -0.6111703413658551, 56 -0.5945164913591590, 57 -0.5776035965513142, 58 -0.5604390262878617, 59 -0.5430302595752546, 60 -0.5253848818220803, 61 -0.5075105815339176, 62 -0.4894151469632753, 63 -0.4711064627160663, 64 -0.4525925063160997, 65 -0.4338813447290861, 66 -0.4149811308476706, 67 -0.3959000999390257, 68 -0.3766465660565522, 69 -0.3572289184172501, 70 -0.3376556177463400, 71 -0.3179351925907259, 72 -0.2980762356029071, 73 -0.2780873997969574, 74 -0.2579773947782034, 75 -0.2377549829482451, 76 -0.2174289756869712, 77 -0.1970082295132342, 78 -0.1765016422258567, 79 -0.1559181490266516, 80 -0.1352667186271445, 81 -0.1145563493406956, 82 -0.0937960651617229, 83 -0.0729949118337358, 84 -0.0521619529078925, 85 -0.0313062657937972, 86 -0.0104369378042598, 87 0.0104369378042598, 88 0.0313062657937972, 89 0.0521619529078925, 90 0.0729949118337358, 91 0.0937960651617229, 92 0.1145563493406956, 93 0.1352667186271445, 94 0.1559181490266516, 95 0.1765016422258567, 96 0.1970082295132342, 97 0.2174289756869712, 98 0.2377549829482451, 99 0.2579773947782034, 100 0.2780873997969574, 101 0.2980762356029071, 102 0.3179351925907259, 103 0.3376556177463400, 104 0.3572289184172501, 105 0.3766465660565522, 106 0.3959000999390257, 107 0.4149811308476706, 108 0.4338813447290861, 109 0.4525925063160997, 110 0.4711064627160663, 111 0.4894151469632753, 112 0.5075105815339176, 113 0.5253848818220803, 114 0.5430302595752546, 115 0.5604390262878617, 116 0.5776035965513142, 117 0.5945164913591590, 118 0.6111703413658551, 119 0.6275578900977726, 120 0.6436719971150083, 121 0.6595056411226444, 122 0.6750519230300931, 123 0.6903040689571928, 124 0.7052554331857488, 125 0.7198995010552305, 126 0.7342298918013638, 127 0.7482403613363824, 128 0.7619248049697269, 129 0.7752772600680049, 130 0.7882919086530552, 131 0.8009630799369827, 132 0.8132852527930605, 133 0.8252530581614230, 134 0.8368612813885015, 135 0.8481048644991847, 136 0.8589789084007133, 137 0.8694786750173527, 138 0.8795995893549102, 139 0.8893372414942055, 140 0.8986873885126239, 141 0.9076459563329236, 142 0.9162090414984952, 143 0.9243729128743136, 144 0.9321340132728527, 145 0.9394889610042837, 146 0.9464345513503147, 147 0.9529677579610971, 148 0.9590857341746905, 149 0.9647858142586956, 150 0.9700655145738374, 151 0.9749225346595943, 152 0.9793547582425894, 153 0.9833602541697529, 154 0.9869372772712794, 155 0.9900842691660192, 156 0.9927998590434373, 157 0.9950828645255290, 158 0.9969322929775997, 159 0.9983473449340834, 160 0.9993274305065947, 161 0.9998723404457334 162 162 }; 163 163 164 constant realGauss150Wt[150]={165 REAL(0.0003276086705538),166 REAL(0.0007624720924706),167 REAL(0.0011976474864367),168 REAL(0.0016323569986067),169 REAL(0.0020663664924131),170 REAL(0.0024994789888943),171 REAL(0.0029315036836558),172 REAL(0.0033622516236779),173 REAL(0.0037915348363451),174 REAL(0.0042191661429919),175 REAL(0.0046449591497966),176 REAL(0.0050687282939456),177 REAL(0.0054902889094487),178 REAL(0.0059094573005900),179 REAL(0.0063260508184704),180 REAL(0.0067398879387430),181 REAL(0.0071507883396855),182 REAL(0.0075585729801782),183 REAL(0.0079630641773633),184 REAL(0.0083640856838475),185 REAL(0.0087614627643580),186 REAL(0.0091550222717888),187 REAL(0.0095445927225849),188 REAL(0.0099300043714212),189 REAL(0.0103110892851360),190 REAL(0.0106876814158841),191 REAL(0.0110596166734735),192 REAL(0.0114267329968529),193 REAL(0.0117888704247183),194 REAL(0.0121458711652067),195 REAL(0.0124975796646449),196 REAL(0.0128438426753249),197 REAL(0.0131845093222756),198 REAL(0.0135194311690004),199 REAL(0.0138484622795371),200 REAL(0.0141714592928592),201 REAL(0.0144882814685445),202 REAL(0.0147987907597169),203 REAL(0.0151028518701744),204 REAL(0.0154003323133401),205 REAL(0.0156911024699895),206 REAL(0.0159750356447283),207 REAL(0.0162520081211971),208 REAL(0.0165218992159766),209 REAL(0.0167845913311726),210 REAL(0.0170399700056559),211 REAL(0.0172879239649355),212 REAL(0.0175283451696437),213 REAL(0.0177611288626114),214 REAL(0.0179861736145128),215 REAL(0.0182033813680609),216 REAL(0.0184126574807331),217 REAL(0.0186139107660094),218 REAL(0.0188070535331042),219 REAL(0.0189920016251754),220 REAL(0.0191686744559934),221 REAL(0.0193369950450545),222 REAL(0.0194968900511231),223 REAL(0.0196482898041878),224 REAL(0.0197911283358190),225 REAL(0.0199253434079123),226 REAL(0.0200508765398072),227 REAL(0.0201676730337687),228 REAL(0.0202756819988200),229 REAL(0.0203748563729175),230 REAL(0.0204651529434560),231 REAL(0.0205465323660984),232 REAL(0.0206189591819181),233 REAL(0.0206824018328499),234 REAL(0.0207368326754401),235 REAL(0.0207822279928917),236 REAL(0.0208185680053983),237 REAL(0.0208458368787627),238 REAL(0.0208640227312962),239 REAL(0.0208731176389954),240 REAL(0.0208731176389954),241 REAL(0.0208640227312962),242 REAL(0.0208458368787627),243 REAL(0.0208185680053983),244 REAL(0.0207822279928917),245 REAL(0.0207368326754401),246 REAL(0.0206824018328499),247 REAL(0.0206189591819181),248 REAL(0.0205465323660984),249 REAL(0.0204651529434560),250 REAL(0.0203748563729175),251 REAL(0.0202756819988200),252 REAL(0.0201676730337687),253 REAL(0.0200508765398072),254 REAL(0.0199253434079123),255 REAL(0.0197911283358190),256 REAL(0.0196482898041878),257 REAL(0.0194968900511231),258 REAL(0.0193369950450545),259 REAL(0.0191686744559934),260 REAL(0.0189920016251754),261 REAL(0.0188070535331042),262 REAL(0.0186139107660094),263 REAL(0.0184126574807331),264 REAL(0.0182033813680609),265 REAL(0.0179861736145128),266 REAL(0.0177611288626114),267 REAL(0.0175283451696437),268 REAL(0.0172879239649355),269 REAL(0.0170399700056559),270 REAL(0.0167845913311726),271 REAL(0.0165218992159766),272 REAL(0.0162520081211971),273 REAL(0.0159750356447283),274 REAL(0.0156911024699895),275 REAL(0.0154003323133401),276 REAL(0.0151028518701744),277 REAL(0.0147987907597169),278 REAL(0.0144882814685445),279 REAL(0.0141714592928592),280 REAL(0.0138484622795371),281 REAL(0.0135194311690004),282 REAL(0.0131845093222756),283 REAL(0.0128438426753249),284 REAL(0.0124975796646449),285 REAL(0.0121458711652067),286 REAL(0.0117888704247183),287 REAL(0.0114267329968529),288 REAL(0.0110596166734735),289 REAL(0.0106876814158841),290 REAL(0.0103110892851360),291 REAL(0.0099300043714212),292 REAL(0.0095445927225849),293 REAL(0.0091550222717888),294 REAL(0.0087614627643580),295 REAL(0.0083640856838475),296 REAL(0.0079630641773633),297 REAL(0.0075585729801782),298 REAL(0.0071507883396855),299 REAL(0.0067398879387430),300 REAL(0.0063260508184704),301 REAL(0.0059094573005900),302 REAL(0.0054902889094487),303 REAL(0.0050687282939456),304 REAL(0.0046449591497966),305 REAL(0.0042191661429919),306 REAL(0.0037915348363451),307 REAL(0.0033622516236779),308 REAL(0.0029315036836558),309 REAL(0.0024994789888943),310 REAL(0.0020663664924131),311 REAL(0.0016323569986067),312 REAL(0.0011976474864367),313 REAL(0.0007624720924706),314 REAL(0.0003276086705538)164 constant double Gauss150Wt[150]={ 165 0.0003276086705538, 166 0.0007624720924706, 167 0.0011976474864367, 168 0.0016323569986067, 169 0.0020663664924131, 170 0.0024994789888943, 171 0.0029315036836558, 172 0.0033622516236779, 173 0.0037915348363451, 174 0.0042191661429919, 175 0.0046449591497966, 176 0.0050687282939456, 177 0.0054902889094487, 178 0.0059094573005900, 179 0.0063260508184704, 180 0.0067398879387430, 181 0.0071507883396855, 182 0.0075585729801782, 183 0.0079630641773633, 184 0.0083640856838475, 185 0.0087614627643580, 186 0.0091550222717888, 187 0.0095445927225849, 188 0.0099300043714212, 189 0.0103110892851360, 190 0.0106876814158841, 191 0.0110596166734735, 192 0.0114267329968529, 193 0.0117888704247183, 194 0.0121458711652067, 195 0.0124975796646449, 196 0.0128438426753249, 197 0.0131845093222756, 198 0.0135194311690004, 199 0.0138484622795371, 200 0.0141714592928592, 201 0.0144882814685445, 202 0.0147987907597169, 203 0.0151028518701744, 204 0.0154003323133401, 205 0.0156911024699895, 206 0.0159750356447283, 207 0.0162520081211971, 208 0.0165218992159766, 209 0.0167845913311726, 210 0.0170399700056559, 211 0.0172879239649355, 212 0.0175283451696437, 213 0.0177611288626114, 214 0.0179861736145128, 215 0.0182033813680609, 216 0.0184126574807331, 217 0.0186139107660094, 218 0.0188070535331042, 219 0.0189920016251754, 220 0.0191686744559934, 221 0.0193369950450545, 222 0.0194968900511231, 223 0.0196482898041878, 224 0.0197911283358190, 225 0.0199253434079123, 226 0.0200508765398072, 227 0.0201676730337687, 228 0.0202756819988200, 229 0.0203748563729175, 230 0.0204651529434560, 231 0.0205465323660984, 232 0.0206189591819181, 233 0.0206824018328499, 234 0.0207368326754401, 235 0.0207822279928917, 236 0.0208185680053983, 237 0.0208458368787627, 238 0.0208640227312962, 239 0.0208731176389954, 240 0.0208731176389954, 241 0.0208640227312962, 242 0.0208458368787627, 243 0.0208185680053983, 244 0.0207822279928917, 245 0.0207368326754401, 246 0.0206824018328499, 247 0.0206189591819181, 248 0.0205465323660984, 249 0.0204651529434560, 250 0.0203748563729175, 251 0.0202756819988200, 252 0.0201676730337687, 253 0.0200508765398072, 254 0.0199253434079123, 255 0.0197911283358190, 256 0.0196482898041878, 257 0.0194968900511231, 258 0.0193369950450545, 259 0.0191686744559934, 260 0.0189920016251754, 261 0.0188070535331042, 262 0.0186139107660094, 263 0.0184126574807331, 264 0.0182033813680609, 265 0.0179861736145128, 266 0.0177611288626114, 267 0.0175283451696437, 268 0.0172879239649355, 269 0.0170399700056559, 270 0.0167845913311726, 271 0.0165218992159766, 272 0.0162520081211971, 273 0.0159750356447283, 274 0.0156911024699895, 275 0.0154003323133401, 276 0.0151028518701744, 277 0.0147987907597169, 278 0.0144882814685445, 279 0.0141714592928592, 280 0.0138484622795371, 281 0.0135194311690004, 282 0.0131845093222756, 283 0.0128438426753249, 284 0.0124975796646449, 285 0.0121458711652067, 286 0.0117888704247183, 287 0.0114267329968529, 288 0.0110596166734735, 289 0.0106876814158841, 290 0.0103110892851360, 291 0.0099300043714212, 292 0.0095445927225849, 293 0.0091550222717888, 294 0.0087614627643580, 295 0.0083640856838475, 296 0.0079630641773633, 297 0.0075585729801782, 298 0.0071507883396855, 299 0.0067398879387430, 300 0.0063260508184704, 301 0.0059094573005900, 302 0.0054902889094487, 303 0.0050687282939456, 304 0.0046449591497966, 305 0.0042191661429919, 306 0.0037915348363451, 307 0.0033622516236779, 308 0.0029315036836558, 309 0.0024994789888943, 310 0.0020663664924131, 311 0.0016323569986067, 312 0.0011976474864367, 313 0.0007624720924706, 314 0.0003276086705538 315 315 }; -
sasmodels/models/lib/gauss20.c
r14de349 r994d77f 9 9 10 10 // Gaussians 11 constant realGauss20Wt[20]={12 REAL(.0176140071391521),13 REAL(.0406014298003869),14 REAL(.0626720483341091),15 REAL(.0832767415767047),16 REAL(.10193011981724),17 REAL(.118194531961518),18 REAL(.131688638449177),19 REAL(.142096109318382),20 REAL(.149172986472604),21 REAL(.152753387130726),22 REAL(.152753387130726),23 REAL(.149172986472604),24 REAL(.142096109318382),25 REAL(.131688638449177),26 REAL(.118194531961518),27 REAL(.10193011981724),28 REAL(.0832767415767047),29 REAL(.0626720483341091),30 REAL(.0406014298003869),31 REAL(.0176140071391521)11 constant double Gauss20Wt[20]={ 12 .0176140071391521, 13 .0406014298003869, 14 .0626720483341091, 15 .0832767415767047, 16 .10193011981724, 17 .118194531961518, 18 .131688638449177, 19 .142096109318382, 20 .149172986472604, 21 .152753387130726, 22 .152753387130726, 23 .149172986472604, 24 .142096109318382, 25 .131688638449177, 26 .118194531961518, 27 .10193011981724, 28 .0832767415767047, 29 .0626720483341091, 30 .0406014298003869, 31 .0176140071391521 32 32 }; 33 33 34 constant realGauss20Z[20]={35 REAL(-.993128599185095),36 REAL(-.963971927277914),37 REAL(-.912234428251326),38 REAL(-.839116971822219),39 REAL(-.746331906460151),40 REAL(-.636053680726515),41 REAL(-.510867001950827),42 REAL(-.37370608871542),43 REAL(-.227785851141645),44 REAL(-.076526521133497),45 REAL(.0765265211334973),46 REAL(.227785851141645),47 REAL(.37370608871542),48 REAL(.510867001950827),49 REAL(.636053680726515),50 REAL(.746331906460151),51 REAL(.839116971822219),52 REAL(.912234428251326),53 REAL(.963971927277914),54 REAL(.993128599185095)34 constant double Gauss20Z[20]={ 35 -.993128599185095, 36 -.963971927277914, 37 -.912234428251326, 38 -.839116971822219, 39 -.746331906460151, 40 -.636053680726515, 41 -.510867001950827, 42 -.37370608871542, 43 -.227785851141645, 44 -.076526521133497, 45 .0765265211334973, 46 .227785851141645, 47 .37370608871542, 48 .510867001950827, 49 .636053680726515, 50 .746331906460151, 51 .839116971822219, 52 .912234428251326, 53 .963971927277914, 54 .993128599185095 55 55 }; -
sasmodels/models/lib/gauss76.c
r14de349 r994d77f 9 9 10 10 // Gaussians 11 constant realGauss76Wt[76]={12 REAL(.00126779163408536), //013 REAL(.00294910295364247),14 REAL(.00462793522803742),15 REAL(.00629918049732845),16 REAL(.00795984747723973),17 REAL(.00960710541471375),18 REAL(.0112381685696677),19 REAL(.0128502838475101),20 REAL(.0144407317482767),21 REAL(.0160068299122486),22 REAL(.0175459372914742), //1023 REAL(.0190554584671906),24 REAL(.020532847967908),25 REAL(.0219756145344162),26 REAL(.0233813253070112),27 REAL(.0247476099206597),28 REAL(.026072164497986),29 REAL(.0273527555318275),30 REAL(.028587223650054),31 REAL(.029773487255905),32 REAL(.0309095460374916), //2033 REAL(.0319934843404216),34 REAL(.0330234743977917),35 REAL(.0339977794120564),36 REAL(.0349147564835508),37 REAL(.0357728593807139),38 REAL(.0365706411473296),39 REAL(.0373067565423816),40 REAL(.0379799643084053),41 REAL(.0385891292645067),42 REAL(.0391332242205184), //3043 REAL(.0396113317090621),44 REAL(.0400226455325968),45 REAL(.040366472122844),46 REAL(.0406422317102947),47 REAL(.0408494593018285),48 REAL(.040987805464794),49 REAL(.0410570369162294),50 REAL(.0410570369162294),51 REAL(.040987805464794),52 REAL(.0408494593018285), //4053 REAL(.0406422317102947),54 REAL(.040366472122844),55 REAL(.0400226455325968),56 REAL(.0396113317090621),57 REAL(.0391332242205184),58 REAL(.0385891292645067),59 REAL(.0379799643084053),60 REAL(.0373067565423816),61 REAL(.0365706411473296),62 REAL(.0357728593807139), //5063 REAL(.0349147564835508),64 REAL(.0339977794120564),65 REAL(.0330234743977917),66 REAL(.0319934843404216),67 REAL(.0309095460374916),68 REAL(.029773487255905),69 REAL(.028587223650054),70 REAL(.0273527555318275),71 REAL(.026072164497986),72 REAL(.0247476099206597), //6073 REAL(.0233813253070112),74 REAL(.0219756145344162),75 REAL(.020532847967908),76 REAL(.0190554584671906),77 REAL(.0175459372914742),78 REAL(.0160068299122486),79 REAL(.0144407317482767),80 REAL(.0128502838475101),81 REAL(.0112381685696677),82 REAL(.00960710541471375), //7083 REAL(.00795984747723973),84 REAL(.00629918049732845),85 REAL(.00462793522803742),86 REAL(.00294910295364247),87 REAL(.00126779163408536)//75 (indexed from 0)11 constant double Gauss76Wt[76]={ 12 .00126779163408536, //0 13 .00294910295364247, 14 .00462793522803742, 15 .00629918049732845, 16 .00795984747723973, 17 .00960710541471375, 18 .0112381685696677, 19 .0128502838475101, 20 .0144407317482767, 21 .0160068299122486, 22 .0175459372914742, //10 23 .0190554584671906, 24 .020532847967908, 25 .0219756145344162, 26 .0233813253070112, 27 .0247476099206597, 28 .026072164497986, 29 .0273527555318275, 30 .028587223650054, 31 .029773487255905, 32 .0309095460374916, //20 33 .0319934843404216, 34 .0330234743977917, 35 .0339977794120564, 36 .0349147564835508, 37 .0357728593807139, 38 .0365706411473296, 39 .0373067565423816, 40 .0379799643084053, 41 .0385891292645067, 42 .0391332242205184, //30 43 .0396113317090621, 44 .0400226455325968, 45 .040366472122844, 46 .0406422317102947, 47 .0408494593018285, 48 .040987805464794, 49 .0410570369162294, 50 .0410570369162294, 51 .040987805464794, 52 .0408494593018285, //40 53 .0406422317102947, 54 .040366472122844, 55 .0400226455325968, 56 .0396113317090621, 57 .0391332242205184, 58 .0385891292645067, 59 .0379799643084053, 60 .0373067565423816, 61 .0365706411473296, 62 .0357728593807139, //50 63 .0349147564835508, 64 .0339977794120564, 65 .0330234743977917, 66 .0319934843404216, 67 .0309095460374916, 68 .029773487255905, 69 .028587223650054, 70 .0273527555318275, 71 .026072164497986, 72 .0247476099206597, //60 73 .0233813253070112, 74 .0219756145344162, 75 .020532847967908, 76 .0190554584671906, 77 .0175459372914742, 78 .0160068299122486, 79 .0144407317482767, 80 .0128502838475101, 81 .0112381685696677, 82 .00960710541471375, //70 83 .00795984747723973, 84 .00629918049732845, 85 .00462793522803742, 86 .00294910295364247, 87 .00126779163408536 //75 (indexed from 0) 88 88 }; 89 89 90 constant realGauss76Z[76]={91 REAL(-.999505948362153), //092 REAL(-.997397786355355),93 REAL(-.993608772723527),94 REAL(-.988144453359837),95 REAL(-.981013938975656),96 REAL(-.972229228520377),97 REAL(-.961805126758768),98 REAL(-.949759207710896),99 REAL(-.936111781934811),100 REAL(-.92088586125215),101 REAL(-.904107119545567), //10102 REAL(-.885803849292083),103 REAL(-.866006913771982),104 REAL(-.844749694983342),105 REAL(-.822068037328975),106 REAL(-.7980001871612),107 REAL(-.77258672828181),108 REAL(-.74587051350361),109 REAL(-.717896592387704),110 REAL(-.688712135277641),111 REAL(-.658366353758143), //20112 REAL(-.626910417672267),113 REAL(-.594397368836793),114 REAL(-.560882031601237),115 REAL(-.526420920401243),116 REAL(-.491072144462194),117 REAL(-.454895309813726),118 REAL(-.417951418780327),119 REAL(-.380302767117504),120 REAL(-.342012838966962),121 REAL(-.303146199807908), //30122 REAL(-.263768387584994),123 REAL(-.223945802196474),124 REAL(-.183745593528914),125 REAL(-.143235548227268),126 REAL(-.102483975391227),127 REAL(-.0615595913906112),128 REAL(-.0205314039939986),129 REAL(.0205314039939986),130 REAL(.0615595913906112),131 REAL(.102483975391227), //40132 REAL(.143235548227268),133 REAL(.183745593528914),134 REAL(.223945802196474),135 REAL(.263768387584994),136 REAL(.303146199807908),137 REAL(.342012838966962),138 REAL(.380302767117504),139 REAL(.417951418780327),140 REAL(.454895309813726),141 REAL(.491072144462194), //50142 REAL(.526420920401243),143 REAL(.560882031601237),144 REAL(.594397368836793),145 REAL(.626910417672267),146 REAL(.658366353758143),147 REAL(.688712135277641),148 REAL(.717896592387704),149 REAL(.74587051350361),150 REAL(.77258672828181),151 REAL(.7980001871612), //60152 REAL(.822068037328975),153 REAL(.844749694983342),154 REAL(.866006913771982),155 REAL(.885803849292083),156 REAL(.904107119545567),157 REAL(.92088586125215),158 REAL(.936111781934811),159 REAL(.949759207710896),160 REAL(.961805126758768),161 REAL(.972229228520377), //70162 REAL(.981013938975656),163 REAL(.988144453359837),164 REAL(.993608772723527),165 REAL(.997397786355355),166 REAL(.999505948362153)//7590 constant double Gauss76Z[76]={ 91 -.999505948362153, //0 92 -.997397786355355, 93 -.993608772723527, 94 -.988144453359837, 95 -.981013938975656, 96 -.972229228520377, 97 -.961805126758768, 98 -.949759207710896, 99 -.936111781934811, 100 -.92088586125215, 101 -.904107119545567, //10 102 -.885803849292083, 103 -.866006913771982, 104 -.844749694983342, 105 -.822068037328975, 106 -.7980001871612, 107 -.77258672828181, 108 -.74587051350361, 109 -.717896592387704, 110 -.688712135277641, 111 -.658366353758143, //20 112 -.626910417672267, 113 -.594397368836793, 114 -.560882031601237, 115 -.526420920401243, 116 -.491072144462194, 117 -.454895309813726, 118 -.417951418780327, 119 -.380302767117504, 120 -.342012838966962, 121 -.303146199807908, //30 122 -.263768387584994, 123 -.223945802196474, 124 -.183745593528914, 125 -.143235548227268, 126 -.102483975391227, 127 -.0615595913906112, 128 -.0205314039939986, 129 .0205314039939986, 130 .0615595913906112, 131 .102483975391227, //40 132 .143235548227268, 133 .183745593528914, 134 .223945802196474, 135 .263768387584994, 136 .303146199807908, 137 .342012838966962, 138 .380302767117504, 139 .417951418780327, 140 .454895309813726, 141 .491072144462194, //50 142 .526420920401243, 143 .560882031601237, 144 .594397368836793, 145 .626910417672267, 146 .658366353758143, 147 .688712135277641, 148 .717896592387704, 149 .74587051350361, 150 .77258672828181, 151 .7980001871612, //60 152 .822068037328975, 153 .844749694983342, 154 .866006913771982, 155 .885803849292083, 156 .904107119545567, 157 .92088586125215, 158 .936111781934811, 159 .949759207710896, 160 .961805126758768, 161 .972229228520377, //70 162 .981013938975656, 163 .988144453359837, 164 .993608772723527, 165 .997397786355355, 166 .999505948362153 //75 167 167 }; -
sasmodels/models/sphere.py
r19dcb933 r994d77f 86 86 # This should perhaps be volume normalized? 87 87 form_volume = """ 88 return REAL(1.333333333333333)*M_PI*radius*radius*radius;88 return 1.333333333333333*M_PI*radius*radius*radius; 89 89 """ 90 90 91 91 Iq = """ 92 const realqr = q*radius;93 realsn, cn;92 const double qr = q*radius; 93 double sn, cn; 94 94 SINCOS(qr, sn, cn); 95 const real bes = qr==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-qr*cn)/(qr*qr*qr);96 const realfq = bes * (sld - solvent_sld) * form_volume(radius);97 return REAL(1.0e-4)*fq*fq;95 const double bes = qr==0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr); 96 const double fq = bes * (sld - solvent_sld) * form_volume(radius); 97 return 1.0e-4*fq*fq; 98 98 """ 99 99 … … 101 101 Iqxy = """ 102 102 // never called since no orientation or magnetic parameters. 103 return REAL(-1.0); 103 //return -1.0; 104 return Iq(sqrt(qx*qx + qy*qy), sld, solvent_sld, radius); 104 105 """ 105 106 -
sasmodels/models/triaxial_ellipsoid.c
r5d4777d r994d77f 1 real form_volume(real req_minor, real req_major, realrpolar);2 real Iq(real q, real sld, realsolvent_sld,3 real req_minor, real req_major, realrpolar);4 real Iqxy(real qx, real qy, real sld, realsolvent_sld,5 real req_minor, real req_major, real rpolar, real theta, real phi, realpsi);1 double form_volume(double req_minor, double req_major, double rpolar); 2 double Iq(double q, double sld, double solvent_sld, 3 double req_minor, double req_major, double rpolar); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 5 double req_minor, double req_major, double rpolar, double theta, double phi, double psi); 6 6 7 real form_volume(real req_minor, real req_major, realrpolar)7 double form_volume(double req_minor, double req_major, double rpolar) 8 8 { 9 return REAL(1.333333333333333)*M_PI*req_minor*req_major*rpolar;9 return 1.333333333333333*M_PI*req_minor*req_major*rpolar; 10 10 } 11 11 12 real Iq(realq,13 realsld,14 realsolvent_sld,15 realreq_minor,16 realreq_major,17 realrpolar)12 double Iq(double q, 13 double sld, 14 double solvent_sld, 15 double req_minor, 16 double req_major, 17 double rpolar) 18 18 { 19 // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0); // Exclude invalid region19 // if (req_minor > req_major || req_major > rpolar) return -1.0; // Exclude invalid region 20 20 21 realsn, cn;22 realst, ct;23 //const real lower = REAL(0.0);24 //const real upper = REAL(1.0);25 real outer = REAL(0.0);21 double sn, cn; 22 double st, ct; 23 //const double lower = 0.0; 24 //const double upper = 1.0; 25 double outer = 0.0; 26 26 for (int i=0;i<76;i++) { 27 //const realcos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;28 const real x = REAL(0.5)*(Gauss76Z[i] + REAL(1.0));27 //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 28 const double x = 0.5*(Gauss76Z[i] + 1.0); 29 29 SINCOS(M_PI_2*x, sn, cn); 30 const realacosx2 = req_minor*req_minor*cn*cn;31 const realbsinx2 = req_major*req_major*sn*sn;32 const realc2 = rpolar*rpolar;30 const double acosx2 = req_minor*req_minor*cn*cn; 31 const double bsinx2 = req_major*req_major*sn*sn; 32 const double c2 = rpolar*rpolar; 33 33 34 real inner = REAL(0.0);34 double inner = 0.0; 35 35 for (int j=0;j<76;j++) { 36 const real y = REAL(0.5)*(Gauss76Z[j] + REAL(1.0));37 const real t = q*sqrt(acosx2 + bsinx2*(REAL(1.0)-y*y) + c2*y*y);36 const double y = 0.5*(Gauss76Z[j] + 1.0); 37 const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y); 38 38 SINCOS(t, st, ct); 39 const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );39 const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 40 40 inner += Gauss76Wt[j] * fq * fq ; 41 41 } 42 outer += Gauss76Wt[i] * REAL(0.5)* inner;42 outer += Gauss76Wt[i] * 0.5 * inner; 43 43 } 44 //const realfq2 = (upper-lower)/2*outer;45 const real fq2 = REAL(0.5)*outer;46 const reals = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);47 return REAL(1.0e-4)* fq2 * s * s;44 //const double fq2 = (upper-lower)/2*outer; 45 const double fq2 = 0.5*outer; 46 const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 47 return 1.0e-4 * fq2 * s * s; 48 48 } 49 49 50 real Iqxy(real qx, realqy,51 realsld,52 realsolvent_sld,53 realreq_minor,54 realreq_major,55 realrpolar,56 realtheta,57 realphi,58 realpsi)50 double Iqxy(double qx, double qy, 51 double sld, 52 double solvent_sld, 53 double req_minor, 54 double req_major, 55 double rpolar, 56 double theta, 57 double phi, 58 double psi) 59 59 { 60 // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0); // Exclude invalid region60 // if (req_minor > req_major || req_major > rpolar) return -1.0; // Exclude invalid region 61 61 62 realstheta, ctheta;63 realsphi, cphi;64 realspsi, cpsi;65 realst, ct;62 double stheta, ctheta; 63 double sphi, cphi; 64 double spsi, cpsi; 65 double st, ct; 66 66 67 const realq = sqrt(qx*qx + qy*qy);68 const realqxhat = qx/q;69 const realqyhat = qy/q;67 const double q = sqrt(qx*qx + qy*qy); 68 const double qxhat = qx/q; 69 const double qyhat = qy/q; 70 70 SINCOS(theta*M_PI_180, stheta, ctheta); 71 71 SINCOS(phi*M_PI_180, sphi, cphi); 72 72 SINCOS(psi*M_PI_180, spsi, cpsi); 73 const realcalpha = ctheta*cphi*qxhat + stheta*qyhat;74 const realcnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;75 const realcmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;76 const realt = q*sqrt(req_minor*req_minor*cnu*cnu73 const double calpha = ctheta*cphi*qxhat + stheta*qyhat; 74 const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 75 const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 76 const double t = q*sqrt(req_minor*req_minor*cnu*cnu 77 77 + req_major*req_major*cmu*cmu 78 78 + rpolar*rpolar*calpha*calpha); 79 79 SINCOS(t, st, ct); 80 const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );81 const reals = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);80 const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 81 const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 82 82 83 return REAL(1.0e-4)* fq * fq * s * s;83 return 1.0e-4 * fq * fq * s * s; 84 84 } 85 85
Note: See TracChangeset
for help on using the changeset viewer.