Changeset ff85cc5 in sasview for sansmodels/src
- Timestamp:
- Aug 28, 2009 12:24:34 PM (15 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:
- a50ef5e
- Parents:
- e2afadf
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/libigor/libStructureFactor.c
rd2e8e5d rff85cc5 1 /* SimpleFit.c1 /* SimpleFit.c 2 2 3 3 A simplified project designed to act as a template for your curve fitting function. … … 7 7 #include "StandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h 8 8 #include "libStructureFactor.h" 9 10 9 11 //Hard Sphere Structure Factor 10 12 // … … 15 17 double calp,cbeta,cgam,prefac,c,vstruc; 16 18 double r,phi,struc; 17 19 18 20 r = dp[0]; 19 21 phi = dp[1]; … … 26 28 // 27 29 // calculate the structure factor 28 // 30 // 29 31 a = 2.0*q*r; 30 32 asq = a*a; … … 40 42 vstruc = 1.0/(1.0-c); 41 43 struc = vstruc; 42 44 43 45 return(struc); 44 46 } … … 54 56 double lam,lam2,test,mu,alpha,beta; 55 57 double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 56 58 57 59 qv= q; 58 60 rad = dp[0]; … … 60 62 eps = dp[2]; 61 63 tau = dp[3]; 62 64 63 65 eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps); 64 66 65 67 sig = 2.0 * rad; 66 68 aa = sig/(1.0 - eps); … … 114 116 // 115 117 sq = 1.0/(aq*aq +bq*bq); 116 118 117 119 return(sq); 118 120 } … … 132 134 double sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamm; 133 135 double x,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck; 134 136 135 137 x= q; 136 138 137 139 req = dp[0]; 138 140 phis = dp[1]; 139 141 edibkb = dp[2]; 140 142 lambda = dp[3]; 141 143 142 144 sigma = req*2.; 143 145 eta = phis; 144 146 eta2 = eta*eta; 145 147 eta3 = eta*eta2; 146 eta4 = eta*eta3; 147 etam1 = 1. - eta; 148 eta4 = eta*eta3; 149 etam1 = 1. - eta; 148 150 etam14 = etam1*etam1*etam1*etam1; 149 151 alpha = ( pow((1. + 2.*eta),2) + eta3*( eta-4.0 ) )/etam14; … … 152 154 // 153 155 // calculate the structure factor 154 156 155 157 sk = x*sigma; 156 158 sk2 = sk*sk; … … 165 167 ck = -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3; 166 168 struc = 1./(1.-ck); 167 169 168 170 return(struc); 169 171 } 170 172 171 // Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions 173 // Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions 172 174 // 173 175 double … … 182 184 double pi; 183 185 int ierr; 184 186 185 187 pi = 4.0*atan(1.); 186 188 QQ= q; 187 188 diam=dp[0] *2; //in dp[0] coming from python is radius : cahnged on Mar. 30, 2009189 190 diam=dp[0]; //in � 189 191 zz = dp[1]; //# of charges 190 VolFrac=dp[2]; 192 VolFrac=dp[2]; 191 193 Temp=dp[3]; //in degrees Kelvin 192 194 csalt=dp[4]; //in molarity … … 202 204 Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0); //in m^3 203 205 cs=csalt*6.022E23*1E3; //# salt molecules/m^3 204 206 205 207 // Compute the derived values of : 206 // Ionic strength IonSt (in C^2/m^3) 208 // Ionic strength IonSt (in C^2/m^3) 207 209 // Kappa (Debye-Huckel screening length in m) 208 210 // and gamma Exp(-k) … … 211 213 // Kappa=2/SIdiam // Use to compare with HP paper 212 214 gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2)); 213 214 // Finally set up dimensionless parameters 215 216 // Finally set up dimensionless parameters 215 217 Qdiam=QQ*diam; 216 218 gMSAWave[6] = Kappa*SIdiam; 217 219 gMSAWave[4] = VolFrac; 218 220 219 221 //Function sqhpa(qq) {this is where Hayter-Penfold program began} 220 222 221 223 // FIRST CALCULATE COUPLING 222 224 223 225 ss=pow(gMSAWave[4],(1.0/3.0)); 224 226 gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss); 225 227 226 228 // CALCULATE COEFFICIENTS, CHECK ALL IS WELL 227 229 // AND IF SO CALCULATE S(Q*SIG) 228 230 229 231 ierr=0; 230 232 ierr=sqcoef(ierr); … … 237 239 // print "Please report HPMSA problem with above error code" 238 240 } 239 241 240 242 return(SofQ); 241 243 } … … 267 269 int 268 270 sqcoef(int ir) 269 { 271 { 270 272 int itm=40,ix,ig,ii; 271 273 double acc=5.0E-6,del,e1,e2,f1,f2; 272 274 273 275 // WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" 274 276 f1=0; //these were never properly initialized... 275 277 f2=0; 276 278 277 279 ig=1; 278 280 if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) { … … 290 292 gMSAWave[10]=fmin(gMSAWave[4],0.20); 291 293 if ((ig!=1) || ( gMSAWave[9]>=0.15)) { 292 ii=0; 294 ii=0; 293 295 do { 294 296 ii=ii+1; 295 297 if(ii>itm) { 296 298 ir=-1; 297 return ir; 299 return ir; 298 300 } 299 301 if (gMSAWave[10]<=0.0) { … … 347 349 int 348 350 sqfun(int ix, int ir) 349 { 351 { 350 352 double acc=1.0e-6; 351 353 double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek; … … 365 367 p2=0; 366 368 p3=0; 367 369 368 370 // CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981) 369 370 reta = gMSAWave[16]; 371 372 reta = gMSAWave[16]; 371 373 eta2 = reta*reta; 372 374 eta3 = eta2*reta; … … 379 381 ibig=1; 380 382 } 381 383 382 384 gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]); 383 385 rgek = gMSAWave[11]; … … 389 391 d = 1.0-reta; 390 392 d2 = d*d; 391 dak = d/rak; 393 dak = d/rak; 392 394 dd2 = 1.0/d2; 393 395 dd4 = dd2*dd2; … … 400 402 eta21 = 2.0*reta+1.0; 401 403 eta22 = eta21*eta21; 402 404 403 405 // ALPHA(I) 404 406 405 407 al1 = -eta21*dak; 406 408 al2 = (14.0*eta2-4.0*reta-1.0)*dak2; 407 409 al3 = 36.0*eta2*dak4; 408 410 409 411 // BETA(I) 410 412 411 413 be1 = -(eta2+7.0*reta+1.0)*dak; 412 414 be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2; 413 415 be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4; 414 416 415 417 // NU(I) 416 418 417 419 vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak; 418 420 vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2; … … 420 422 vu4 = vu1+e24*rak*vu3; 421 423 vu5 = eta6d*(vu2+4.0*vu3); 422 424 423 425 // PHI(I) 424 426 425 427 ph1 = eta6d/rak; 426 428 ph2 = d-e12*dak2; 427 429 428 430 // TAU(I) 429 431 430 432 ta1 = (reta+5.0)/(5.0*rak); 431 433 ta2 = eta2d*dak2; … … 433 435 ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2); 434 436 ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2; 435 437 436 438 // double PRECISION SINH(K), COSH(K) 437 439 438 440 ex1 = exp(rak); 439 441 ex2 = 0.0; … … 445 447 ckma = ck-1.0-rak*sk; 446 448 skma = sk-rak*ck; 447 449 448 450 // a(I) 449 451 450 452 a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4; 451 453 if (ibig==0) { … … 453 455 a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4; 454 456 } 455 457 456 458 // b(I) 457 459 458 460 b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4; 459 461 if (ibig==0) { … … 461 463 b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4; 462 464 } 463 465 464 466 // V(I) 465 467 466 468 v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45; 467 469 if (ibig==0) { … … 469 471 v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45; 470 472 } 471 472 473 474 473 475 // P(I) 474 476 475 477 pp1 = ph1*ph1; 476 478 pp2 = ph2*ph2; … … 482 484 p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2; 483 485 } 484 486 485 487 // T(I) 486 488 487 489 t1 = ta3+ta4*a1+ta5*b1; 488 490 if (ibig!=0) { 489 491 490 492 // VERY LARGE SCREENING: ASYMPTOTIC SOLUTION 491 493 492 494 v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45; 493 495 t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0; … … 516 518 } 517 519 gMSAWave[10] = gMSAWave[16]; 518 520 519 521 } else { 520 522 521 523 t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk); 522 524 t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0; 523 525 524 526 // MU(i) 525 527 526 528 um1 = t2*a2-e12*v2*v2; 527 529 um2 = t1*a2+t2*a1-e24*v1*v2; … … 530 532 um5 = t1*a3+t3*a1-e24*v1*v3; 531 533 um6 = t3*a3-e12*v3*v3; 532 534 533 535 // GILLAN CONDITION ? 534 536 // … … 538 540 // 539 541 if ((ix==1) || (ix==3)) { 540 542 541 543 // NO - CALCULATE REMAINING COEFFICIENTS. 542 544 543 545 // LAMBDA(I) 544 546 545 547 al1 = e12*p2*p2; 546 548 al2 = e24*p1*p2-b2-b2; … … 549 551 al5 = e24*p1*p3-b3-b3-ak2; 550 552 al6 = e12*p3*p3; 551 553 552 554 // OMEGA(I) 553 555 554 556 w16 = um1*al6-al1*um6; 555 557 w15 = um1*al5-al1*um5; … … 557 559 w13 = um1*al3-al1*um3; 558 560 w12 = um1*al2-al1*um2; 559 561 560 562 w26 = um2*al6-al2*um6; 561 563 w25 = um2*al5-al2*um5; 562 564 w24 = um2*al4-al2*um4; 563 565 564 566 w36 = um3*al6-al3*um6; 565 567 w35 = um3*al5-al3*um5; … … 571 573 w3526 = w35+w26; 572 574 w3425 = w34+w25; 573 575 574 576 // QUARTIC COEFFICIENTS 575 577 576 578 w4 = w16*w16-w13*w36; 577 579 w3 = 2.0*w16*w15-w13*w3526-w12*w36; … … 579 581 w1 = 2.0*w15*w14-w13*w24-w12*w3425; 580 582 w0 = w14*w14-w12*w24; 581 583 582 584 // ESTIMATE THE STARTING VALUE OF f 583 585 584 586 if (ix==1) { 585 587 // LARGE K … … 600 602 fap = -(pg+p2*ca)/p3; 601 603 } 602 604 603 605 // AND REFINE IT ACCORDING TO NEWTON 604 606 ii=0; … … 616 618 del=fabs((fap-fa)/fa); 617 619 } while (del>acc); 618 620 619 621 ir = ir+ii; 620 622 fa = fap; … … 648 650 } 649 651 650 // called as DiamCyl(hcyl,rcyl) 651 //modified from Igor NIST package XOP 652 // called as DiamCyl(hcyl,rcyl) 652 653 double 653 DiamCyl(double dp[], double q)654 { 655 double hcyl, rcyl; 654 DiamCyl(double hcyl, double rcyl) 655 { 656 656 657 double diam,a,b,t1,t2,ddd; 657 658 double pi; 658 659 rcyl = dp[0]; 660 hcyl = dp[1]; 659 661 660 pi = 4.0*atan(1.0); 662 if (rcyl == 0 || hcyl == 0) { 663 return 0.0; 664 } 661 665 662 a = rcyl; 666 663 b = hcyl/2.0; … … 669 666 ddd = 3.0*t1*t2; 670 667 diam = pow(ddd,(1.0/3.0)); 671 672 return(diam /2); //return radius668 669 return(diam); 673 670 } 674 671 … … 680 677 //returns DIAMETER 681 678 // called as DiamEllip(aa,bb) 682 683 //modified from Igor NIST package XOP684 679 double 685 DiamEllip(double dp[], double q)686 { 687 double aa, bb; 680 DiamEllip(double aa, double bb) 681 { 682 688 683 double ee,e1,bd,b1,bL,b2,del,ddd,diam; 689 690 aa = dp[0]; 691 bb = dp[1]; 692 if (aa == 0 || bb == 0) { 693 return 0.0; 694 } 695 if (aa == bb) { 696 return aa; 697 } 684 698 685 if(aa>bb) { 699 686 ee = (aa*aa - bb*bb)/(aa*aa); … … 701 688 ee = (bb*bb - aa*aa)/(bb*bb); 702 689 } 703 690 704 691 bd = 1.0-ee; 705 692 e1 = sqrt(ee); … … 708 695 b2 = 1.0 + bd/2/e1*log(bL); 709 696 del = 0.75*b1*b2; 710 697 711 698 ddd = 2.0*(del+1.0)*aa*bb*bb; //volume is always calculated correctly 712 699 diam = pow(ddd,(1.0/3.0)); 713 714 return(diam /2); //return radius700 701 return(diam); 715 702 } 716 703 717 704 double 718 705 sqhcal(double qq) 719 { 720 double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; 706 { 707 double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; 721 708 // WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" 722 709 723 710 etaz = gMSAWave[10]; 724 711 akz = gMSAWave[12]; … … 733 720 sk = 0.5*(x1-x2); 734 721 ak2 = akz*akz; 735 722 736 723 if (qq<=0.0) { 737 724 SofQ = -1.0/gMSAWave[0];
Note: See TracChangeset
for help on using the changeset viewer.