Changeset c5ac2b2 in sasmodels for sasmodels/models
- Timestamp:
- Jul 18, 2016 2:05:12 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:
- 4a21b121
- Parents:
- 256dfe1
- Location:
- sasmodels/models
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_multi_shell.c
rabdd01c rc5ac2b2 1 FourShell(double dp[], double q)2 {3 // variables are:4 //[0] scale factor5 //[1] radius of core [ï¿œ]6 //[2] SLD of the core [ï¿œ-2]7 //[3] thickness of shell 1 [ï¿œ]8 //[4] SLD of shell 19 //[5] thickness of shell 2 [ï¿œ]10 //[6] SLD of shell 211 //[7] thickness of shell 312 //[8] SLD of shell 313 //[9] thickness of shell 314 //[10] SLD of shell 315 //[11] SLD of solvent16 //[12] background [cm-1]17 18 double x,pi;19 double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names20 double bes,f,vol,qr,contr,f2;21 double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;22 23 pi = 4.0*atan(1.0);24 x=q;25 26 scale = dp[0];27 rcore = dp[1];28 rhocore = dp[2];29 thick1 = dp[3];30 rhoshel1 = dp[4];31 thick2 = dp[5];32 rhoshel2 = dp[6];33 thick3 = dp[7];34 rhoshel3 = dp[8];35 thick4 = dp[9];36 rhoshel4 = dp[10];37 rhosolv = dp[11];38 bkg = dp[12];39 40 // core first, then add in shells41 qr=x*rcore;42 contr = rhocore-rhoshel1;43 if(qr == 0){44 bes = 1.0;45 }else{46 bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);47 }48 vol = 4.0*pi/3.0*rcore*rcore*rcore;49 f = vol*bes*contr;50 //now the shell (1)51 qr=x*(rcore+thick1);52 contr = rhoshel1-rhoshel2;53 if(qr == 0){54 bes = 1.0;55 }else{56 bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);57 }58 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);59 f += vol*bes*contr;60 //now the shell (2)61 qr=x*(rcore+thick1+thick2);62 contr = rhoshel2-rhoshel3;63 if(qr == 0){64 bes = 1.0;65 }else{66 bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);67 }68 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);69 f += vol*bes*contr;70 //now the shell (3)71 qr=x*(rcore+thick1+thick2+thick3);72 contr = rhoshel3-rhoshel4;73 if(qr == 0){74 bes = 1.0;75 }else{76 bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);77 }78 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);79 f += vol*bes*contr;80 //now the shell (4)81 qr=x*(rcore+thick1+thick2+thick3+thick4);82 contr = rhoshel4-rhosolv;83 if(qr == 0){84 bes = 1.0;85 }else{86 bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);87 }88 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);89 f += vol*bes*contr;90 91 92 // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]93 f2 = f*f/vol*1.0e8;94 95 //scale if desired96 f2 *= scale;97 // then add in the background98 f2 += bkg;99 100 return(f2);101 }102 103 // above is cut and past with no changes from libigor four shell104 static double105 f_exp(double q, double r, double sld_in, double sld_out,106 double thickness, double A)107 {108 const double vol = 4.0/3.0 * M_PI * r * r * r;109 const double qr = q * r;110 const double alpha = A * r/thickness;111 const double bes = sph_j1c(qr);112 const double B = (sld_out - sld_in)/expm1(A);113 const double C = sld_in - B;114 double fun;115 if (qr == 0.0) {116 fun = 1.0;117 } else if (fabs(A) > 0.0) {118 const double qrsq = qr*qr;119 const double alphasq = alpha*alpha;120 const double sumsq = alphasq + qrsq;121 double sinqr, cosqr;122 SINCOS(qr, sinqr, cosqr);123 fun = -3.0(124 ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq125 - (alpha*sinqr/qr - cosqr)126 ) / sumsq;127 } else {128 fun = bes;129 }130 return vol * (B*fun + C*bes);131 }132 133 static double134 f_linear(double q, double r, double sld, double slope)135 {136 const double vol = 4.0/3.0 * M_PI * r * r * r;137 const double qr = q * r;138 const double bes = sph_j1c(qr);139 double fun = 0.0;140 if (qr > 0.0) {141 const double qrsq = qr*qr;142 double sinqr, cosqr;143 SINCOS(qr, sinqr, cosqr);144 // Jae-He's code seems to simplify to this145 // fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);146 // Rederiving the math, we get the following instead:147 fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);148 }149 return vol * (sld*bes + fun);150 }151 1 152 2 static double … … 154 4 { 155 5 const double bes = sph_j1c(q * r); 156 const double vol = 4.0/3.0 * M_PI * r * r * r;6 const double vol = M_4PI_3 * cube(r); 157 7 return sld * vol * bes; 158 8 } 159 9 160 static double 10 double 11 form_volume(double core_radius, double n, double thickness[]); 12 double 161 13 form_volume(double core_radius, double n, double thickness[]) 162 14 { 163 int i;164 15 double r = core_radius; 165 for (i =0; i < n; i++) {16 for (int i=0; i < n; i++) { 166 17 r += thickness[i]; 167 18 } 168 return 4.0/3.0 * M_PI * r * r * r;19 return M_4PI_3 * cube(r); 169 20 } 170 21 171 static double 172 Iq(double q, double core_sld, double core_radius, double solvent_sld, 173 double n, double in_sld[], double out_sld[], double thickness[], 174 double A[]) 22 double 23 Iq(double q, double core_sld, double core_radius, 24 double solvent_sld, double num_shells, double sld[], double thickness[]); 25 double 26 Iq(double q, double core_sld, double core_radius, 27 double solvent_sld, double num_shells, double sld[], double thickness[]) 175 28 { 176 int i; 29 const int n = (int)ceil(num_shells); 30 double f, r, last_sld; 177 31 r = core_radius; 178 f = f_constant(q, r, core_sld); 179 for (i=0; i<n; i++){ 180 const double r0 = r; 32 last_sld = core_sld; 33 f = 0.; 34 for (int i=0; i<n; i++) { 35 f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r); 36 last_sld = sld[i]; 181 37 r += thickness[i]; 182 if (r == r0) {183 // no thickness, so nothing to add184 } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) {185 f -= f_constant(q, r0, sld_in[i]);186 f += f_constant(q, r, sld_in[i]);187 } else if (fabs(A[i]) < 1.0e-4) {188 const double slope = (sld_out[i] - sld_in[i])/thickness[i];189 f -= f_linear(q, r0, sld_in[i], slope);190 f += f_linear(q, r, sld_out[i], slope);191 } else {192 f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);193 f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);194 }195 38 } 196 f -= f_constant(q, r, solvent_sld); 197 f2 = f * f * 1.0e-4; 198 199 return f2; 39 f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r); 40 return f * f * 1.0e-4; 200 41 } -
sasmodels/models/core_multi_shell.py
r46ed760 rc5ac2b2 87 87 88 88 # ["name", "units", default, [lower, upper], "type","description"], 89 parameters = [["volfraction", "", 0.05, [0,1],"", 90 "volume fraction of spheres"], 91 ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 89 parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 92 90 "Core scattering length density"], 93 91 ["radius", "Ang", 200., [0, inf], "volume", … … 95 93 ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 96 94 "Solvent scattering length density"], 97 ["n", "", 1, [0, 4], "volume",95 ["n", "", 1, [0, 10], "volume", 98 96 "number of shells"], 99 97 ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", … … 103 101 ] 104 102 105 #source = ["lib/sph_j1c.c", "onion.c"]103 source = ["lib/sph_j1c.c", "core_multi_shell.c"] 106 104 107 def Iq(q, *args, **kw): 108 return q 109 110 def Iqxy(qx, *args, **kw): 111 return qx 112 113 114 def profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness): 105 def profile(sld_core, radius, sld_solvent, n, sld, thickness): 115 106 """ 116 107 SLD profile … … 120 111 121 112 """ 122 # r = []123 # beta = []124 # # for core at r=0125 # r.append(0)126 # beta.append(self.params['sld_core0'])127 # # for core at r=rad_core128 # r.append(self.params['rad_core0'])129 # beta.append(self.params['sld_core0'])130 #131 # # for shells132 # for n in range(1, self.n_shells+1):133 # # Left side of each shells134 # r0 = r[len(r)-1]135 # r.append(r0)136 # exec "beta.append(self.params['sld_shell%s'% str(n)])"137 #138 # # Right side of each shells139 # exec "r0 += self.params['thick_shell%s'% str(n)]"140 # r.append(r0)141 # exec "beta.append(self.params['sld_shell%s'% str(n)])"142 #143 # # for solvent144 # r0 = r[len(r)-1]145 # r.append(r0)146 # beta.append(self.params['sld_solv'])147 # r_solv = 5*r0/4148 # r.append(r_solv)149 # beta.append(self.params['sld_solv'])150 #151 # return r, beta152 # above is directly from old code -- below is alotered from Woitek's first153 # cut an the onion. Seems like we should be consistant?154 155 113 total_radius = 1.25*(sum(thickness[:n]) + radius + 1) 156 114 … … 172 130 r.append(r0 + thickness[k]) 173 131 beta.append(sld[k]) 174 # add in the solvent132 # add in the solvent 175 133 r.append(r[-1]) 176 134 beta.append(sld_solvent) … … 187 145 return 1.0, 1.0 188 146 189 demo = dict(volfraction = 1.0, 190 sld_core = 6.4, 147 demo = dict(sld_core = 6.4, 191 148 radius = 60, 192 149 sld_solvent = 6.4, 193 n = 1, 194 sld = [2.0], 195 thickness = [10]) 150 n = 2, 151 sld = [2.0, 3.0], 152 thickness = 20, 153 thickness1_pd = 0.3, 154 thickness2_pd = 0.3, 155 thickness1_pd_n = 10, 156 thickness2_pd_n = 10, 157 )
Note: See TracChangeset
for help on using the changeset viewer.