Changeset 34c2649 in sasview for sansmodels/src/sans/models/libigor/libCylinder.h
- Timestamp:
- Nov 21, 2011 4:52:25 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:
- 1810613
- Parents:
- afea8fe
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/libigor/libCylinder.h
r3a084629 r34c2649 67 67 68 68 /////////functions for WRC implementation of flexible cylinders 69 static double Sk_WR(double q, double L, double b); 70 static double w_WR(double x); 71 static double u1(double q, double L, double b); 72 static double u_WR(double q, double L, double b); 73 static double Rgsquarezero(double q, double L, double b); 74 static double Rgsquareshort(double q, double L, double b); 75 static double Rgsquare(double q, double L, double b); 76 static double AlphaSquare(double x); 77 static double miu(double x); 78 static double Sdebye(double q, double L, double b); 79 static double Sdebye1(double q, double L, double b); 80 static double Sexv(double q, double L, double b); 81 static double Sexvnew(double q, double L, double b); 82 static double a2short(double q, double L, double b, double p1short, double p2short, double q0); 83 static double a1short(double q, double L, double b, double p1short, double p2short, double q0); 84 static double a2long(double q, double L, double b, double p1, double p2, double q0); 85 static double sech_WR(double x); 86 static double a1long(double q, double L, double b, double p1, double p2, double q0); 87 static double gammaln(double xx); 69 static double 70 gammaln(double xx) 71 { 72 double x,y,tmp,ser; 73 static double cof[6]={76.18009172947146,-86.50532032941677, 74 24.01409824083091,-1.231739572450155, 75 0.1208650973866179e-2,-0.5395239384953e-5}; 76 int j; 77 78 y=x=xx; 79 tmp=x+5.5; 80 tmp -= (x+0.5)*log(tmp); 81 ser=1.000000000190015; 82 for (j=0;j<=5;j++) ser += cof[j]/++y; 83 return -tmp+log(2.5066282746310005*ser/x); 84 } 85 86 // 87 static double 88 AlphaSquare(double x) 89 { 90 double yy; 91 yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 92 93 return (yy); 94 } 95 96 // 97 static double 98 Rgsquarezero(double q, double L, double b) 99 { 100 double yy; 101 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)))); 102 103 return (yy); 104 } 105 106 // 107 static double 108 Rgsquareshort(double q, double L, double b) 109 { 110 double yy; 111 yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 112 113 return (yy); 114 } 115 116 // 117 static double 118 Rgsquare(double q, double L, double b) 119 { 120 double yy; 121 yy = AlphaSquare(L/b)*L*b/6.0; 122 123 return (yy); 124 } 125 126 // ?? funciton is not used - but should the log actually be log10??? 127 static double 128 miu(double x) 129 { 130 double yy; 131 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))); 132 133 return (yy); 134 } 135 136 //WR named this w (too generic) 137 static double 138 w_WR(double x) 139 { 140 double yy; 141 yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 142 143 return (yy); 144 } 145 146 // 147 static double 148 u1(double q, double L, double b) 149 { 150 double yy; 151 152 yy = Rgsquareshort(q,L,b)*q*q; 153 154 return (yy); 155 } 156 157 // was named u 158 static double 159 u_WR(double q, double L, double b) 160 { 161 double yy; 162 yy = Rgsquare(q,L,b)*q*q; 163 return (yy); 164 } 165 166 167 // 168 static double 169 Sdebye1(double q, double L, double b) 170 { 171 double yy; 172 yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 173 174 return (yy); 175 } 176 177 178 // 179 static double 180 Sdebye(double q, double L, double b) 181 { 182 double yy; 183 yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 184 185 return (yy); 186 } 187 188 // 189 static double 190 Sexv(double q, double L, double b) 191 { 192 double yy,C1,C2,C3,miu,Rg2; 193 C1=1.22; 194 C2=0.4288; 195 C3=-1.651; 196 miu = 0.585; 197 198 Rg2 = Rgsquare(q,L,b); 199 200 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))); 201 202 return (yy); 203 } 204 205 // this must be WR modified version 206 static double 207 Sexvnew(double q, double L, double b) 208 { 209 double yy,C1,C2,C3,miu; 210 double del=1.05,C_star2,Rg2; 211 212 C1=1.22; 213 C2=0.4288; 214 C3=-1.651; 215 miu = 0.585; 216 217 //calculating the derivative to decide on the corection (cutoff) term? 218 // I have modified this from WRs original code 219 220 if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 221 C_star2 = 0.0; 222 } else { 223 C_star2 = 1.0; 224 } 225 226 Rg2 = Rgsquare(q,L,b); 227 228 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))); 229 230 return (yy); 231 } 232 233 // these are the messy ones 234 static double 235 a2short(double q, double L, double b, double p1short, double p2short, double q0) 236 { 237 double yy,Rg2_sh; 238 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 239 double pi; 240 241 E = 2.718281828459045091; 242 pi = 4.0*atan(1.0); 243 Rg2_sh = Rgsquareshort(q,L,b); 244 Rg2_sh2 = Rg2_sh*Rg2_sh; 245 t1 = ((q0*q0*Rg2_sh)/(b*b)); 246 Et1 = pow(E,t1); 247 Emt1 =pow(E,-t1); 248 q02 = q0*q0; 249 q0p = pow(q0,(-4.0 + p2short) ); 250 251 //E is the number e 252 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))))))); 253 254 return (yy); 255 } 256 257 // 258 static double 259 a1short(double q, double L, double b, double p1short, double p2short, double q0) 260 { 261 double yy,Rg2_sh; 262 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 263 double pi; 264 265 E = 2.718281828459045091; 266 pi = 4.0*atan(1.0); 267 Rg2_sh = Rgsquareshort(q,L,b); 268 Rg2_sh2 = Rg2_sh*Rg2_sh; 269 b3 = b*b*b; 270 t1 = ((q0*q0*Rg2_sh)/(b*b)); 271 Et1 = pow(E,t1); 272 Emt1 =pow(E,-t1); 273 q02 = q0*q0; 274 q0p = pow(q0,(-4.0 + p1short) ); 275 276 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)))))); 277 278 return(yy); 279 } 280 281 282 //need to define this on my own 283 static double 284 sech_WR(double x) 285 { 286 return(1/cosh(x)); 287 } 288 289 // this one will be lots of trouble 290 static double 291 a2long(double q, double L, double b, double p1, double p2, double q0) 292 { 293 double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 294 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 295 double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 296 297 pi = 4.0*atan(1.0); 298 E = 2.718281828459045091; 299 if( L/b > 10.0) { 300 C = 3.06/pow((L/b),0.44); 301 } else { 302 C = 1.0; 303 } 304 305 C1 = 1.22; 306 C2 = 0.4288; 307 C3 = -1.651; 308 C4 = 1.523; 309 C5 = 0.1477; 310 miu = 0.585; 311 312 Rg2 = Rgsquare(q,L,b); 313 Rg22 = Rg2*Rg2; 314 b2 = b*b; 315 b3 = b*b*b; 316 b4 = b3*b; 317 q02 = q0*q0; 318 q03 = q0*q0*q0; 319 q04 = q03*q0; 320 q05 = q04*q0; 321 322 t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 323 324 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; 325 326 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); 327 328 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); 329 330 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); 331 332 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); 333 334 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)); 335 336 t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 337 338 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; 339 340 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); 341 342 343 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)))))))))); 344 345 return (yy); 346 } 347 // 348 static double 349 a1long(double q, double L, double b, double p1, double p2, double q0) 350 { 351 double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 352 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 353 double E,pi; 354 double b2,b3,b4,q02,q03,q04,q05,Rg22; 355 356 pi = 4.0*atan(1.0); 357 E = 2.718281828459045091; 358 359 if( L/b > 10.0) { 360 C = 3.06/pow((L/b),0.44); 361 } else { 362 C = 1.0; 363 } 364 365 C1 = 1.22; 366 C2 = 0.4288; 367 C3 = -1.651; 368 C4 = 1.523; 369 C5 = 0.1477; 370 miu = 0.585; 371 372 Rg2 = Rgsquare(q,L,b); 373 Rg22 = Rg2*Rg2; 374 b2 = b*b; 375 b3 = b*b*b; 376 b4 = b3*b; 377 q02 = q0*q0; 378 q03 = q0*q0*q0; 379 q04 = q03*q0; 380 q05 = q04*q0; 381 382 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)))); 383 384 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)))))); 385 386 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)))); 387 388 t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 389 390 t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 391 392 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))); 393 394 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)); 395 396 t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 397 398 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)))))); 399 400 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)))))); 401 402 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)); 403 404 t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 405 406 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)))); 407 408 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)))))); 409 410 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)))); 411 412 413 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))))))))))); 414 415 return (yy); 416 } 417 418 419 static double 420 Sk_WR(double q, double L, double b) 421 { 422 // 423 double p1,p2,p1short,p2short,q0; 424 double C,ans,q0short,Sexvmodify,pi; 425 426 pi = 4.0*atan(1.0); 427 428 p1 = 4.12; 429 p2 = 4.42; 430 p1short = 5.36; 431 p2short = 5.62; 432 q0 = 3.1; 433 // 434 q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 435 436 // 437 if(L/b > 10.0) { 438 C = 3.06/pow((L/b),0.44); 439 } else { 440 C = 1.0; 441 } 442 // 443 444 if( L > 4*b ) { // Longer Chains 445 if (q*b <= 3.1) { //Modified by Yun on Oct. 15, 446 Sexvmodify = Sexvnew(q, L, b); 447 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); 448 } else { //q(i)*b > 3.1 449 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); 450 } 451 } else { //L <= 4*b Shorter Chains 452 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 453 if (q*b<=0.01) { 454 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 455 } else { 456 ans = Sdebye1(q,L,b); 457 } 458 } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 459 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); 460 } 461 } 462 463 return(ans); 464 //return(a2long(q, L, b, p1, p2, q0)); 465 } 466
Note: See TracChangeset
for help on using the changeset viewer.