Changeset 5f0dcab in sasview
- Timestamp:
- Dec 7, 2011 3:39:05 PM (13 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- feadd6f
- Parents:
- 2605da22
- Location:
- sansmodels/src/sans/models/libigor
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/libigor/libCylinder.c
r34c2649 r5f0dcab 8 8 #include "GaussWeights.h" 9 9 #include "libCylinder.h" 10 11 /////////functions for WRC implementation of flexible cylinders 12 static double 13 gammaln(double xx) 14 { 15 double x,y,tmp,ser; 16 static double cof[6]={76.18009172947146,-86.50532032941677, 17 24.01409824083091,-1.231739572450155, 18 0.1208650973866179e-2,-0.5395239384953e-5}; 19 int j; 20 21 y=x=xx; 22 tmp=x+5.5; 23 tmp -= (x+0.5)*log(tmp); 24 ser=1.000000000190015; 25 for (j=0;j<=5;j++) ser += cof[j]/++y; 26 return -tmp+log(2.5066282746310005*ser/x); 27 } 28 /////////functions for WRC implementation of flexible cylinders 29 30 // 31 static double 32 AlphaSquare(double x) 33 { 34 double yy; 35 yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 36 37 return (yy); 38 } 39 40 // 41 static double 42 Rgsquarezero(double q, double L, double b) 43 { 44 double yy; 45 yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); 46 47 return (yy); 48 } 49 50 // 51 static double 52 Rgsquareshort(double q, double L, double b) 53 { 54 double yy; 55 yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 56 57 return (yy); 58 } 59 60 // 61 static double 62 Rgsquare(double q, double L, double b) 63 { 64 double yy; 65 yy = AlphaSquare(L/b)*L*b/6.0; 66 67 return (yy); 68 } 69 70 // ?? funciton is not used - but should the log actually be log10??? 71 static double 72 miu(double x) 73 { 74 double yy; 75 yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); 76 77 return (yy); 78 } 79 80 //WR named this w (too generic) 81 static double 82 w_WR(double x) 83 { 84 double yy; 85 yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 86 87 return (yy); 88 } 89 90 // 91 static double 92 u1(double q, double L, double b) 93 { 94 double yy; 95 96 yy = Rgsquareshort(q,L,b)*q*q; 97 98 return (yy); 99 } 100 101 // was named u 102 static double 103 u_WR(double q, double L, double b) 104 { 105 double yy; 106 yy = Rgsquare(q,L,b)*q*q; 107 return (yy); 108 } 109 110 111 // 112 static double 113 Sdebye1(double q, double L, double b) 114 { 115 double yy; 116 yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 117 118 return (yy); 119 } 120 121 122 // 123 static double 124 Sdebye(double q, double L, double b) 125 { 126 double yy; 127 yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 128 129 return (yy); 130 } 131 132 // 133 static double 134 Sexv(double q, double L, double b) 135 { 136 double yy,C1,C2,C3,miu,Rg2; 137 C1=1.22; 138 C2=0.4288; 139 C3=-1.651; 140 miu = 0.585; 141 142 Rg2 = Rgsquare(q,L,b); 143 144 yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 145 146 return (yy); 147 } 148 149 // this must be WR modified version 150 static double 151 Sexvnew(double q, double L, double b) 152 { 153 double yy,C1,C2,C3,miu; 154 double del=1.05,C_star2,Rg2; 155 156 C1=1.22; 157 C2=0.4288; 158 C3=-1.651; 159 miu = 0.585; 160 161 //calculating the derivative to decide on the corection (cutoff) term? 162 // I have modified this from WRs original code 163 164 if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 165 C_star2 = 0.0; 166 } else { 167 C_star2 = 1.0; 168 } 169 170 Rg2 = Rgsquare(q,L,b); 171 172 yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 173 174 return (yy); 175 } 176 177 // these are the messy ones 178 static double 179 a2short(double q, double L, double b, double p1short, double p2short, double q0) 180 { 181 double yy,Rg2_sh; 182 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 183 double pi; 184 185 E = 2.718281828459045091; 186 pi = 4.0*atan(1.0); 187 Rg2_sh = Rgsquareshort(q,L,b); 188 Rg2_sh2 = Rg2_sh*Rg2_sh; 189 t1 = ((q0*q0*Rg2_sh)/(b*b)); 190 Et1 = pow(E,t1); 191 Emt1 =pow(E,-t1); 192 q02 = q0*q0; 193 q0p = pow(q0,(-4.0 + p2short) ); 194 195 //E is the number e 196 yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); 197 198 return (yy); 199 } 200 201 // 202 static double 203 a1short(double q, double L, double b, double p1short, double p2short, double q0) 204 { 205 double yy,Rg2_sh; 206 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 207 double pi; 208 209 E = 2.718281828459045091; 210 pi = 4.0*atan(1.0); 211 Rg2_sh = Rgsquareshort(q,L,b); 212 Rg2_sh2 = Rg2_sh*Rg2_sh; 213 b3 = b*b*b; 214 t1 = ((q0*q0*Rg2_sh)/(b*b)); 215 Et1 = pow(E,t1); 216 Emt1 =pow(E,-t1); 217 q02 = q0*q0; 218 q0p = pow(q0,(-4.0 + p1short) ); 219 220 yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); 221 222 return(yy); 223 } 224 225 226 //need to define this on my own 227 static double 228 sech_WR(double x) 229 { 230 return(1/cosh(x)); 231 } 232 233 // this one will be lots of trouble 234 static double 235 a2long(double q, double L, double b, double p1, double p2, double q0) 236 { 237 double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 238 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 239 double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 240 241 pi = 4.0*atan(1.0); 242 E = 2.718281828459045091; 243 if( L/b > 10.0) { 244 C = 3.06/pow((L/b),0.44); 245 } else { 246 C = 1.0; 247 } 248 249 C1 = 1.22; 250 C2 = 0.4288; 251 C3 = -1.651; 252 C4 = 1.523; 253 C5 = 0.1477; 254 miu = 0.585; 255 256 Rg2 = Rgsquare(q,L,b); 257 Rg22 = Rg2*Rg2; 258 b2 = b*b; 259 b3 = b*b*b; 260 b4 = b3*b; 261 q02 = q0*q0; 262 q03 = q0*q0*q0; 263 q04 = q03*q0; 264 q05 = q04*q0; 265 266 t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 267 268 t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 269 270 t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5); 271 272 t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 273 274 t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 275 276 t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 277 278 t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 279 280 t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 281 282 t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 283 284 t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 285 286 287 yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 288 289 return (yy); 290 } 291 // 292 static double 293 a1long(double q, double L, double b, double p1, double p2, double q0) 294 { 295 double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 296 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 297 double E,pi; 298 double b2,b3,b4,q02,q03,q04,q05,Rg22; 299 300 pi = 4.0*atan(1.0); 301 E = 2.718281828459045091; 302 303 if( L/b > 10.0) { 304 C = 3.06/pow((L/b),0.44); 305 } else { 306 C = 1.0; 307 } 308 309 C1 = 1.22; 310 C2 = 0.4288; 311 C3 = -1.651; 312 C4 = 1.523; 313 C5 = 0.1477; 314 miu = 0.585; 315 316 Rg2 = Rgsquare(q,L,b); 317 Rg22 = Rg2*Rg2; 318 b2 = b*b; 319 b3 = b*b*b; 320 b4 = b3*b; 321 q02 = q0*q0; 322 q03 = q0*q0*q0; 323 q04 = q03*q0; 324 q05 = q04*q0; 325 326 t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 327 328 t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 329 330 t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 331 332 t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 333 334 t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 335 336 t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); 337 338 t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 339 340 t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 341 342 t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 343 344 t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 345 346 t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 347 348 t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 349 350 t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 351 352 t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 353 354 t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 355 356 357 yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); 358 359 return (yy); 360 } 361 362 363 static double 364 Sk_WR(double q, double L, double b) 365 { 366 // 367 double p1,p2,p1short,p2short,q0; 368 double C,ans,q0short,Sexvmodify,pi; 369 370 pi = 4.0*atan(1.0); 371 372 p1 = 4.12; 373 p2 = 4.42; 374 p1short = 5.36; 375 p2short = 5.62; 376 q0 = 3.1; 377 // 378 q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 379 380 // 381 if(L/b > 10.0) { 382 C = 3.06/pow((L/b),0.44); 383 } else { 384 C = 1.0; 385 } 386 // 387 388 if( L > 4*b ) { // Longer Chains 389 if (q*b <= 3.1) { //Modified by Yun on Oct. 15, 390 Sexvmodify = Sexvnew(q, L, b); 391 ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 392 } else { //q(i)*b > 3.1 393 ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); 394 } 395 } else { //L <= 4*b Shorter Chains 396 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 397 if (q*b<=0.01) { 398 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 399 } else { 400 ans = Sdebye1(q,L,b); 401 } 402 } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 403 ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); 404 } 405 } 406 407 return(ans); 408 //return(a2long(q, L, b, p1, p2, q0)); 409 } 10 410 11 411 /* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x … … 20 420 CylinderForm(double dp[], double q) 21 421 { 422 22 423 int i; 23 424 double Pi; -
sansmodels/src/sans/models/libigor/libCylinder.h
r6a2c931 r5f0dcab 68 68 double CSPPKernel(double dp[], double mu, double uu); 69 69 70 /////////functions for WRC implementation of flexible cylinders71 static double72 gammaln(double xx)73 {74 double x,y,tmp,ser;75 static double cof[6]={76.18009172947146,-86.50532032941677,76 24.01409824083091,-1.231739572450155,77 0.1208650973866179e-2,-0.5395239384953e-5};78 int j;79 70 80 y=x=xx;81 tmp=x+5.5;82 tmp -= (x+0.5)*log(tmp);83 ser=1.000000000190015;84 for (j=0;j<=5;j++) ser += cof[j]/++y;85 return -tmp+log(2.5066282746310005*ser/x);86 }87 88 //89 static double90 AlphaSquare(double x)91 {92 double yy;93 yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );94 95 return (yy);96 }97 98 //99 static double100 Rgsquarezero(double q, double L, double b)101 {102 double yy;103 yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));104 105 return (yy);106 }107 108 //109 static double110 Rgsquareshort(double q, double L, double b)111 {112 double yy;113 yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);114 115 return (yy);116 }117 118 //119 static double120 Rgsquare(double q, double L, double b)121 {122 double yy;123 yy = AlphaSquare(L/b)*L*b/6.0;124 125 return (yy);126 }127 128 // ?? funciton is not used - but should the log actually be log10???129 static double130 miu(double x)131 {132 double yy;133 yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x)));134 135 return (yy);136 }137 138 //WR named this w (too generic)139 static double140 w_WR(double x)141 {142 double yy;143 yy = 0.5*(1 + tanh((x - 1.523)/0.1477));144 145 return (yy);146 }147 148 //149 static double150 u1(double q, double L, double b)151 {152 double yy;153 154 yy = Rgsquareshort(q,L,b)*q*q;155 156 return (yy);157 }158 159 // was named u160 static double161 u_WR(double q, double L, double b)162 {163 double yy;164 yy = Rgsquare(q,L,b)*q*q;165 return (yy);166 }167 168 169 //170 static double171 Sdebye1(double q, double L, double b)172 {173 double yy;174 yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );175 176 return (yy);177 }178 179 180 //181 static double182 Sdebye(double q, double L, double b)183 {184 double yy;185 yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));186 187 return (yy);188 }189 190 //191 static double192 Sexv(double q, double L, double b)193 {194 double yy,C1,C2,C3,miu,Rg2;195 C1=1.22;196 C2=0.4288;197 C3=-1.651;198 miu = 0.585;199 200 Rg2 = Rgsquare(q,L,b);201 202 yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));203 204 return (yy);205 }206 207 // this must be WR modified version208 static double209 Sexvnew(double q, double L, double b)210 {211 double yy,C1,C2,C3,miu;212 double del=1.05,C_star2,Rg2;213 214 C1=1.22;215 C2=0.4288;216 C3=-1.651;217 miu = 0.585;218 219 //calculating the derivative to decide on the corection (cutoff) term?220 // I have modified this from WRs original code221 222 if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {223 C_star2 = 0.0;224 } else {225 C_star2 = 1.0;226 }227 228 Rg2 = Rgsquare(q,L,b);229 230 yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));231 232 return (yy);233 }234 235 // these are the messy ones236 static double237 a2short(double q, double L, double b, double p1short, double p2short, double q0)238 {239 double yy,Rg2_sh;240 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;241 double pi;242 243 E = 2.718281828459045091;244 pi = 4.0*atan(1.0);245 Rg2_sh = Rgsquareshort(q,L,b);246 Rg2_sh2 = Rg2_sh*Rg2_sh;247 t1 = ((q0*q0*Rg2_sh)/(b*b));248 Et1 = pow(E,t1);249 Emt1 =pow(E,-t1);250 q02 = q0*q0;251 q0p = pow(q0,(-4.0 + p2short) );252 253 //E is the number e254 yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2)))))));255 256 return (yy);257 }258 259 //260 static double261 a1short(double q, double L, double b, double p1short, double p2short, double q0)262 {263 double yy,Rg2_sh;264 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;265 double pi;266 267 E = 2.718281828459045091;268 pi = 4.0*atan(1.0);269 Rg2_sh = Rgsquareshort(q,L,b);270 Rg2_sh2 = Rg2_sh*Rg2_sh;271 b3 = b*b*b;272 t1 = ((q0*q0*Rg2_sh)/(b*b));273 Et1 = pow(E,t1);274 Emt1 =pow(E,-t1);275 q02 = q0*q0;276 q0p = pow(q0,(-4.0 + p1short) );277 278 yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2))))));279 280 return(yy);281 }282 283 284 //need to define this on my own285 static double286 sech_WR(double x)287 {288 return(1/cosh(x));289 }290 291 // this one will be lots of trouble292 static double293 a2long(double q, double L, double b, double p1, double p2, double q0)294 {295 double yy,C1,C2,C3,C4,C5,miu,C,Rg2;296 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;297 double E,b2,b3,b4,q02,q03,q04,q05,Rg22;298 299 pi = 4.0*atan(1.0);300 E = 2.718281828459045091;301 if( L/b > 10.0) {302 C = 3.06/pow((L/b),0.44);303 } else {304 C = 1.0;305 }306 307 C1 = 1.22;308 C2 = 0.4288;309 C3 = -1.651;310 C4 = 1.523;311 C5 = 0.1477;312 miu = 0.585;313 314 Rg2 = Rgsquare(q,L,b);315 Rg22 = Rg2*Rg2;316 b2 = b*b;317 b3 = b*b*b;318 b4 = b3*b;319 q02 = q0*q0;320 q03 = q0*q0*q0;321 q04 = q03*q0;322 q05 = q04*q0;323 324 t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));325 326 t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;327 328 t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);329 330 t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);331 332 t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);333 334 t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);335 336 t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));337 338 t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));339 340 t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;341 342 t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);343 344 345 yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));346 347 return (yy);348 }349 //350 static double351 a1long(double q, double L, double b, double p1, double p2, double q0)352 {353 double yy,C,C1,C2,C3,C4,C5,miu,Rg2;354 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;355 double E,pi;356 double b2,b3,b4,q02,q03,q04,q05,Rg22;357 358 pi = 4.0*atan(1.0);359 E = 2.718281828459045091;360 361 if( L/b > 10.0) {362 C = 3.06/pow((L/b),0.44);363 } else {364 C = 1.0;365 }366 367 C1 = 1.22;368 C2 = 0.4288;369 C3 = -1.651;370 C4 = 1.523;371 C5 = 0.1477;372 miu = 0.585;373 374 Rg2 = Rgsquare(q,L,b);375 Rg22 = Rg2*Rg2;376 b2 = b*b;377 b3 = b*b*b;378 b4 = b3*b;379 q02 = q0*q0;380 q03 = q0*q0*q0;381 q04 = q03*q0;382 q05 = q04*q0;383 384 t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));385 386 t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));387 388 t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));389 390 t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));391 392 t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));393 394 t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));395 396 t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));397 398 t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));399 400 t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));401 402 t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));403 404 t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));405 406 t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));407 408 t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));409 410 t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));411 412 t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));413 414 415 yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));416 417 return (yy);418 }419 420 421 static double422 Sk_WR(double q, double L, double b)423 {424 //425 double p1,p2,p1short,p2short,q0;426 double C,ans,q0short,Sexvmodify,pi;427 428 pi = 4.0*atan(1.0);429 430 p1 = 4.12;431 p2 = 4.42;432 p1short = 5.36;433 p2short = 5.62;434 q0 = 3.1;435 //436 q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);437 438 //439 if(L/b > 10.0) {440 C = 3.06/pow((L/b),0.44);441 } else {442 C = 1.0;443 }444 //445 446 if( L > 4*b ) { // Longer Chains447 if (q*b <= 3.1) { //Modified by Yun on Oct. 15,448 Sexvmodify = Sexvnew(q, L, b);449 ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);450 } else { //q(i)*b > 3.1451 ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);452 }453 } else { //L <= 4*b Shorter Chains454 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {455 if (q*b<=0.01) {456 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;457 } else {458 ans = Sdebye1(q,L,b);459 }460 } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)461 ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L);462 }463 }464 465 return(ans);466 //return(a2long(q, L, b, p1, p2, q0));467 }468 -
sansmodels/src/sans/models/libigor/libSphere.c
r34c2649 r5f0dcab 8 8 #include "GaussWeights.h" 9 9 #include "libSphere.h" 10 11 12 static double 13 gammln(double xx) { 14 double x,y,tmp,ser; 15 static double cof[6]={76.18009172947146,-86.50532032941677, 16 24.01409824083091,-1.231739572450155, 17 0.1208650973866179e-2,-0.5395239384953e-5}; 18 int j; 19 20 y=x=xx; 21 tmp=x+5.5; 22 tmp -= (x+0.5)*log(tmp); 23 ser=1.000000000190015; 24 for (j=0;j<=5;j++) ser += cof[j]/++y; 25 return -tmp+log(2.5066282746310005*ser/x); 26 } 27 28 static double 29 LogNormal_distr(double sig, double mu, double pt) 30 { 31 double retval,pi; 32 33 pi = 4.0*atan(1.0); 34 retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 35 return(retval); 36 } 37 38 static double 39 Gauss_distr(double sig, double avg, double pt) 40 { 41 double retval,Pi; 42 43 Pi = 4.0*atan(1.0); 44 retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0); 45 return(retval); 46 } 47 48 static double SchulzPoint(double x, double avg, double zz) { 49 double dr; 50 dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 51 return (exp(dr)); 52 }; 53 10 54 11 55 // scattering from a sphere - hardly needs to be an XOP... -
sansmodels/src/sans/models/libigor/libSphere.h
r6a2c931 r5f0dcab 49 49 double SC_Integrand(double w[], double qq, double xx, double yy); 50 50 double SCeval(double Theta, double Phi, double temp3, double temp4, double temp5); 51 52 53 54 static double55 gammln(double xx) {56 double x,y,tmp,ser;57 static double cof[6]={76.18009172947146,-86.50532032941677,58 24.01409824083091,-1.231739572450155,59 0.1208650973866179e-2,-0.5395239384953e-5};60 int j;61 62 y=x=xx;63 tmp=x+5.5;64 tmp -= (x+0.5)*log(tmp);65 ser=1.000000000190015;66 for (j=0;j<=5;j++) ser += cof[j]/++y;67 return -tmp+log(2.5066282746310005*ser/x);68 }69 70 static double71 LogNormal_distr(double sig, double mu, double pt)72 {73 double retval,pi;74 75 pi = 4.0*atan(1.0);76 retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );77 return(retval);78 }79 80 static double81 Gauss_distr(double sig, double avg, double pt)82 {83 double retval,Pi;84 85 Pi = 4.0*atan(1.0);86 retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);87 return(retval);88 }89 90 static double SchulzPoint(double x, double avg, double zz) {91 double dr;92 dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0));93 return (exp(dr));94 };95
Note: See TracChangeset
for help on using the changeset viewer.