Changeset 6e93a02 in sasview for sansmodels/src/libigor/libCylinder.c
- Timestamp:
- Mar 30, 2010 7:53:47 PM (14 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:
- f10063e
- Parents:
- 64017a8
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/libigor/libCylinder.c
r222d75c7 r6e93a02 8 8 #include "GaussWeights.h" 9 9 #include "libCylinder.h" 10 10 11 /* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x 11 12 … … 21 22 int i; 22 23 double Pi; 23 double scale,radius,length,delrho,bkg,halfheight ; //local variables of coefficient wave24 double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv; //local variables of coefficient wave 24 25 int nord=76; //order of integration 25 26 double uplim,lolim; //upper and lower integration limits 26 27 double summ,zi,yyy,answer,vcyl; //running tally of integration 27 28 28 29 Pi = 4.0*atan(1.0); 29 lolim = 0 ;30 lolim = 0.0; 30 31 uplim = Pi/2.0; 31 32 32 33 summ = 0.0; //initialize intergral 33 34 34 35 scale = dp[0]; //make local copies in case memory moves 35 36 radius = dp[1]; 36 37 length = dp[2]; 37 delrho = dp[3]; 38 bkg = dp[4]; 38 sldCyl = dp[3]; 39 sldSolv = dp[4]; 40 bkg = dp[5]; 41 42 delrho = sldCyl-sldSolv; 39 43 halfheight = length/2.0; 40 44 for(i=0;i<nord;i++) { … … 43 47 summ += yyy; 44 48 } 45 49 46 50 answer = (uplim-lolim)/2.0*summ; 47 51 // Multiply by contrast^2 … … 57 61 // add in the background 58 62 answer += bkg; 59 63 60 64 return answer; 61 65 } … … 75 79 { 76 80 int i,j; 77 double Pi ;81 double Pi,slde,sld; 78 82 double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave 79 83 int nord=76; //order of integration … … 87 91 vaj=0.0; 88 92 vbj=Pi; //endpoints of inner integral 89 93 90 94 summ = 0.0; //initialize intergral 91 95 92 96 scale = dp[0]; //make local copies in case memory moves 93 97 ra = dp[1]; 94 98 nu = dp[2]; 95 99 length = dp[3]; 96 delrho = dp[4]; 97 bkg = dp[5]; 100 slde = dp[4]; 101 sld = dp[5]; 102 delrho = slde - sld; 103 bkg = dp[6]; 104 98 105 for(i=0;i<nord;i++) { 99 106 //setup inner integral over the ellipsoidal cross-section 100 107 summj=0; 101 108 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy 102 arg = ra*sqrt(1 -zi*zi);109 arg = ra*sqrt(1.0-zi*zi); 103 110 for(j=0;j<nord;j++) { 104 111 //76 gauss points for the inner integral as well … … 111 118 //divide integral by Pi 112 119 answer /=Pi; 113 120 114 121 //now calculate outer integral 115 122 arg = q*length*zi/2.0; … … 135 142 // add in the background 136 143 answer += bkg; 137 144 138 145 return answer; 139 146 } … … 154 161 { 155 162 int i,j; 156 double Pi ;163 double Pi,slde,sld; 157 164 double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave 158 165 int nordi=76; //order of integration … … 167 174 vaj=0.0; 168 175 vbj=Pi; //endpoints of inner integral 169 176 170 177 summ = 0.0; //initialize intergral 171 178 172 179 scale = dp[0]; //make local copies in case memory moves 173 180 ra = dp[1]; 174 181 nu = dp[2]; 175 182 length = dp[3]; 176 delrho = dp[4]; 177 bkg = dp[5]; 183 slde = dp[4]; 184 sld = dp[5]; 185 delrho = slde - sld; 186 bkg = dp[6]; 187 178 188 for(i=0;i<nordi;i++) { 179 189 //setup inner integral over the ellipsoidal cross-section 180 190 summj=0; 181 191 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy 182 arg = ra*sqrt(1 -zi*zi);192 arg = ra*sqrt(1.0-zi*zi); 183 193 for(j=0;j<nordj;j++) { 184 194 //20 gauss points for the inner integral … … 191 201 //divide integral by Pi 192 202 answer /=Pi; 193 203 194 204 //now calculate outer integral 195 arg = q*length*zi/2 ;205 arg = q*length*zi/2.0; 196 206 if (arg == 0.0){ 197 207 si = 1.0; … … 202 212 summ += yyy; 203 213 } 204 214 205 215 answer = (vb-va)/2.0*summ; 206 216 // Multiply by contrast^2 … … 240 250 double va,vb; //upper and lower integration limits 241 251 double summ,zi,yyy,answer; //running tally of integration 242 double summj,vaj,vbj,zij ; //for the inner integration243 252 double summj,vaj,vbj,zij,slde,sld; //for the inner integration 253 244 254 Pi = 4.0*atan(1.0); 245 255 va = 0.0; … … 247 257 vaj = 0.0; 248 258 vbj = 1.0; //endpoints of inner integral 249 259 250 260 summ = 0.0; //initialize intergral 251 261 252 262 scale = dp[0]; //make local copies in case memory moves 253 263 aa = dp[1]; 254 264 bb = dp[2]; 255 265 cc = dp[3]; 256 delrho = dp[4]; 257 bkg = dp[5]; 266 slde = dp[4]; 267 sld = dp[5]; 268 delrho = slde - sld; 269 bkg = dp[6]; 258 270 for(i=0;i<nordi;i++) { 259 271 //setup inner integral over the ellipsoidal cross-section 260 summj=0 ;272 summj=0.0; 261 273 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy 262 274 for(j=0;j<nordj;j++) { … … 268 280 //now calculate the value of the inner integral 269 281 answer = (vbj-vaj)/2.0*summj; 270 282 271 283 //now calculate outer integral 272 284 yyy = Gauss76Wt[i] * answer; 273 285 summ += yyy; 274 286 } //final scaling is done at the end of the function, after the NT_FP64 case 275 287 276 288 answer = (vb-va)/2.0*summ; 277 289 // Multiply by contrast^2 … … 285 297 // add in the background 286 298 answer += bkg; 287 299 288 300 return answer; 289 301 } … … 310 322 double summ,yyy,answer; //running tally of integration 311 323 double summj,vaj,vbj; //for the inner integration 312 double mu,mudum,arg,sigma,uu,vol ;313 314 324 double mu,mudum,arg,sigma,uu,vol,sldp,sld; 325 326 315 327 // Pi = 4.0*atan(1.0); 316 328 va = 0.0; … … 320 332 321 333 summ = 0.0; //initialize intergral 322 334 323 335 scale = dp[0]; //make local copies in case memory moves 324 336 aa = dp[1]; 325 337 bb = dp[2]; 326 338 cc = dp[3]; 327 delrho = dp[4]; 328 bkg = dp[5]; 329 339 sldp = dp[4]; 340 sld = dp[5]; 341 delrho = sldp - sld; 342 bkg = dp[6]; 343 330 344 mu = q*bb; 331 345 vol = aa*bb*cc; … … 333 347 aa = aa/bb; 334 348 cc = cc/bb; 335 349 336 350 for(i=0;i<nordi;i++) { 337 351 //setup inner integral over the ellipsoidal cross-section 338 summj=0 ;352 summj=0.0; 339 353 sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy 340 354 341 355 for(j=0;j<nordj;j++) { 342 356 //76 gauss points for the inner integral 343 357 uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy 344 mudum = mu*sqrt(1 -sigma*sigma);358 mudum = mu*sqrt(1.0-sigma*sigma); 345 359 yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu); 346 360 summj += yyy; … … 355 369 answer *= sin(arg)*sin(arg)/arg/arg; 356 370 } 357 371 358 372 //now sum up the outer integral 359 373 yyy = Gauss76Wt[i] * answer; 360 374 summ += yyy; 361 375 } //final scaling is done at the end of the function, after the NT_FP64 case 362 376 363 377 answer = (vb-va)/2.0*summ; 364 378 // Multiply by contrast^2 … … 372 386 // add in the background 373 387 answer += bkg; 374 388 375 389 return answer; 376 390 } … … 394 408 int nord=76; //order of integration 395 409 double va,vb,zi; //upper and lower integration limits 396 double summ,answer,pi ; //running tally of integration397 410 double summ,answer,pi,sldc,sld; //running tally of integration 411 398 412 pi = 4.0*atan(1.0); 399 413 va = 0.0; … … 401 415 402 416 summ = 0.0; //initialize intergral 403 417 404 418 scale = dp[0]; //make local copies in case memory moves 405 419 rcore = dp[1]; 406 420 rshell = dp[2]; 407 421 length = dp[3]; 408 delrho = dp[4]; 409 bkg = dp[5]; 410 422 sldc = dp[4]; 423 sld = dp[5]; 424 delrho = sldc - sld; 425 bkg = dp[6]; 426 411 427 for(i=0;i<nord;i++) { 412 428 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 413 429 summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi); 414 430 } 415 431 416 432 answer = (vb-va)/2.0*summ; 417 433 // Multiply by contrast^2 … … 425 441 // add in the background 426 442 answer += bkg; 427 443 428 444 return answer; 429 445 } … … 447 463 int nord=76; //order of integration 448 464 double va,vb,zi; //upper and lower integration limits 449 double summ,answer,pi ; //running tally of integration450 465 double summ,answer,pi,slde,sld; //running tally of integration 466 451 467 pi = 4.0*atan(1.0); 452 468 va = 0.0; … … 454 470 455 471 summ = 0.0; //initialize intergral 456 472 457 473 scale = dp[0]; //make local copies in case memory moves 458 474 nua = dp[1]; 459 475 a = dp[2]; 460 delrho = dp[3]; 461 bkg = dp[4]; 462 476 slde = dp[3]; 477 sld = dp[4]; 478 delrho = slde - sld; 479 bkg = dp[5]; 480 463 481 for(i=0;i<nord;i++) { 464 482 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 465 483 summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi); 466 484 } 467 485 468 486 answer = (vb-va)/2.0*summ; 469 487 // Multiply by contrast^2 … … 477 495 // add in the background 478 496 answer += bkg; 479 497 480 498 return answer; 481 499 } … … 494 512 double uplim,lolim; //upper and lower integration limits 495 513 double summ,zi,yyy,answer,Vpoly; //running tally of integration 496 double range,zz,Pi ;497 514 double range,zz,Pi,sldc,sld; 515 498 516 Pi = 4.0*atan(1.0); 499 517 range = 3.4; 500 518 501 519 summ = 0.0; //initialize intergral 502 520 503 521 scale = dp[0]; //make local copies in case memory moves 504 522 radius = dp[1]; 505 523 length = dp[2]; 506 524 pd = dp[3]; 507 delrho = dp[4]; 508 bkg = dp[5]; 509 525 sldc = dp[4]; 526 sld = dp[5]; 527 delrho = sldc - sld; 528 bkg = dp[6]; 529 510 530 zz = (1.0/pd)*(1.0/pd) - 1.0; 511 531 512 532 lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution 513 if(lolim<0 ) {514 lolim = 0 ;533 if(lolim<0.0) { 534 lolim = 0.0; 515 535 } 516 536 if(pd>0.3) { … … 518 538 } 519 539 uplim = radius*(1.0+range*pd); 520 540 521 541 for(i=0;i<nord;i++) { 522 542 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 524 544 summ += yyy; 525 545 } 526 546 527 547 answer = (uplim-lolim)/2.0*summ; 528 548 //normalize by average cylinder volume … … 536 556 // add in the background 537 557 answer += bkg; 538 558 539 559 return answer; 540 560 } … … 552 572 double uplim,lolim; //upper and lower integration limits 553 573 double summ,zi,yyy,answer,Vpoly; //running tally of integration 554 double range,zz,Pi ;555 556 574 double range,zz,Pi,sldc,sld; 575 576 557 577 Pi = 4.0*atan(1.0); 558 578 range = 3.4; 559 579 560 580 summ = 0.0; //initialize intergral 561 581 562 582 scale = dp[0]; //make local copies in case memory moves 563 583 radius = dp[1]; 564 584 length = dp[2]; 565 585 pd = dp[3]; 566 delrho = dp[4]; 567 bkg = dp[5]; 568 586 sldc = dp[4]; 587 sld = dp[5]; 588 delrho = sldc - sld; 589 bkg = dp[6]; 590 569 591 zz = (1.0/pd)*(1.0/pd) - 1.0; 570 592 571 593 lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution 572 if(lolim<0 ) {573 lolim = 0 ;594 if(lolim<0.0) { 595 lolim = 0.0; 574 596 } 575 597 if(pd>0.3) { … … 577 599 } 578 600 uplim = length*(1.0+range*pd); 579 601 580 602 for(i=0;i<nord;i++) { 581 603 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 583 605 summ += yyy; 584 606 } 585 607 586 608 answer = (uplim-lolim)/2.0*summ; 587 609 //normalize by average cylinder volume (first moment) … … 595 617 // add in the background 596 618 answer += bkg; 597 619 598 620 return answer; 599 621 } … … 613 635 double summ,zi,yyy,answer,Vcyl; //running tally of integration 614 636 double Pi; 615 637 616 638 Pi = 4.0*atan(1.0); 617 639 618 640 lolim = 0.0; 619 641 uplim = Pi/2.0; 620 642 621 643 summ = 0.0; //initialize intergral 622 644 623 645 scale = dp[0]; //make local copies in case memory moves 624 646 rcore = dp[1]; … … 629 651 rhosolv = dp[6]; 630 652 bkg = dp[7]; 631 653 632 654 halfheight = length/2.0; 633 655 634 656 for(i=0;i<nord;i++) { 635 657 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 637 659 summ += yyy; 638 660 } 639 661 640 662 answer = (uplim-lolim)/2.0*summ; 641 // length is the total core length 663 // length is the total core length 642 664 Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick); 643 665 answer /= Vcyl; … … 648 670 // add in the background 649 671 answer += bkg; 650 672 651 673 return answer; 652 674 } … … 667 689 double summ,yyy,answer,Vpoly; //running tally of integration 668 690 double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; 669 691 670 692 Pi = 4.0*atan(1.0); 671 693 672 694 summ = 0.0; //initialize intergral 673 695 Rsqrsumm = 0.0; 674 696 675 697 scale = dp[0]; 676 698 radius = dp[1]; … … 683 705 rhosolv = dp[8]; 684 706 bkg = dp[9]; 685 707 686 708 lolim = exp(log(radius)-(4.*sigma)); 687 if (lolim<0 ) {688 lolim=0 ; //to avoid numerical error when va<0 (-ve r value)709 if (lolim<0.0) { 710 lolim=0.0; //to avoid numerical error when va<0 (-ve r value) 689 711 } 690 712 uplim = exp(log(radius)+(4.*sigma)); 691 713 692 714 for(i=0;i<nord;i++) { 693 715 rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 698 720 Rsqrsumm += Rsqryyy; 699 721 } 700 722 701 723 answer = (uplim-lolim)/2.0*summ; 702 724 Rsqr = (uplim-lolim)/2.0*Rsqrsumm; … … 710 732 // add in the background 711 733 answer += bkg; 712 734 713 735 return answer; 714 736 } … … 726 748 double uplim,lolim; //upper and lower integration limits 727 749 double summ,zi,yyy,answer,oblatevol; //running tally of integration 728 double Pi ;729 750 double Pi,sldc,slds,sld; 751 730 752 Pi = 4.0*atan(1.0); 731 753 732 754 lolim = 0.0; 733 755 uplim = 1.0; 734 756 735 757 summ = 0.0; //initialize intergral 736 737 758 759 738 760 scale = dp[0]; //make local copies in case memory moves 739 761 crmaj = dp[1]; … … 741 763 trmaj = dp[3]; 742 764 trmin = dp[4]; 743 delpc = dp[5]; 744 delps = dp[6]; 745 bkg = dp[7]; 765 sldc = dp[5]; 766 slds = dp[6]; 767 sld = dp[7]; 768 delpc = sldc - slds; //core - shell 769 delps = slds - sld; //shell - solvent 770 bkg = dp[8]; 746 771 747 772 for(i=0;i<nord;i++) { … … 750 775 summ += yyy; 751 776 } 752 777 753 778 answer = (uplim-lolim)/2.0*summ; 754 779 // normalize by particle volume 755 oblatevol = 4 *Pi/3*trmaj*trmaj*trmin;780 oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin; 756 781 answer /= oblatevol; 757 782 758 783 //convert to [cm-1] 759 784 answer *= 1.0e8; … … 762 787 // add in the background 763 788 answer += bkg; 764 789 765 790 return answer; 766 791 } … … 778 803 double uplim,lolim; //upper and lower integration limits 779 804 double summ,zi,yyy,answer,prolatevol; //running tally of integration 780 double Pi ;781 805 double Pi,sldc,slds,sld; 806 782 807 Pi = 4.0*atan(1.0); 783 808 784 809 lolim = 0.0; 785 810 uplim = 1.0; 786 811 787 812 summ = 0.0; //initialize intergral 788 813 789 814 scale = dp[0]; //make local copies in case memory moves 790 815 crmaj = dp[1]; … … 792 817 trmaj = dp[3]; 793 818 trmin = dp[4]; 794 delpc = dp[5]; 795 delps = dp[6]; 796 bkg = dp[7]; 819 sldc = dp[5]; 820 slds = dp[6]; 821 sld = dp[7]; 822 delpc = sldc - slds; //core - shell 823 delps = slds - sld; //shell - sovent 824 bkg = dp[8]; 797 825 798 826 for(i=0;i<nord;i++) { … … 801 829 summ += yyy; 802 830 } 803 831 804 832 answer = (uplim-lolim)/2.0*summ; 805 833 // normalize by particle volume 806 834 prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin; 807 835 answer /= prolatevol; 808 836 809 837 //convert to [cm-1] 810 838 answer *= 1.0e8; … … 813 841 // add in the background 814 842 answer += bkg; 815 843 816 844 return answer; 817 845 } … … 830 858 int nord=76; //order of integration 831 859 double Pi; 832 833 860 861 834 862 Pi = 4.0*atan(1.0); 835 863 836 864 va = 0.0; 837 865 vb = Pi/2.0; 838 866 839 867 summ = 0.0; //initialize intergral 840 868 841 869 scale = dp[0]; 842 870 rcore = dp[1]; … … 849 877 gsd = dp[8]; 850 878 bkg = dp[9]; 851 879 852 880 d=2.0*thick+length; 853 881 halfheight = length/2.0; 854 882 855 883 for(i=0;i<nord;i++) { 856 884 zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0; … … 858 886 summ += yyy; 859 887 } 860 888 861 889 answer = (vb-va)/2.0*summ; 862 // length is the total core length 890 // length is the total core length 863 891 vcyl=Pi*rcore*rcore*(2.0*thick+length)*N; 864 892 answer /= vcyl; … … 869 897 // add in the background 870 898 answer += bkg; 871 899 872 900 return answer; 873 901 } … … 882 910 double scale,del,sig,contr,bkg; //local variables of coefficient wave 883 911 double inten, qval,Pq; 884 double Pi ;885 886 912 double Pi,sldb,sld; 913 914 887 915 Pi = 4.0*atan(1.0); 888 916 scale = dp[0]; 889 917 del = dp[1]; 890 918 sig = dp[2]*del; 891 contr = dp[3]; 892 bkg = dp[4]; 893 qval = q; 894 919 sldb = dp[3]; 920 sld = dp[4]; 921 contr = sldb - sld; 922 bkg = dp[5]; 923 qval=q; 924 895 925 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); 896 926 897 927 inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... 898 928 899 929 inten /= del; //normalize by the thickness (in A) 900 930 901 931 inten *= 1.0e8; // 1/A to 1/cm 902 932 903 933 return(inten+bkg); 904 934 } 905 935 906 936 /* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included 907 ------- 908 ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! 909 910 */ 937 --- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the 938 model is smeared just like any other function 939 */ 911 940 double 912 941 LamellarPS(double dp[], double q) 913 942 { 914 943 double scale,dd,del,sig,contr,NN,Cp,bkg; //local variables of coefficient wave 915 double inten, 916 double Pi,Euler,dQDefault,fii ;944 double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ; 945 double Pi,Euler,dQDefault,fii,sldb,sld; 917 946 int ii,NNint; 918 947 // char buf[256]; 948 949 919 950 Euler = 0.5772156649; // Euler's constant 920 dQDefault = 0.0;//0.0025; //[=] 1/A, q-resolution, default value 951 // dQDefault = 0.0025; //[=] 1/A, q-resolution, default value 952 dQDefault = 0.0; 921 953 dQ = dQDefault; 922 954 923 955 Pi = 4.0*atan(1.0); 924 956 qval = q; 925 957 926 958 scale = dp[0]; 927 959 dd = dp[1]; 928 960 del = dp[2]; 929 961 sig = dp[3]*del; 930 contr = dp[4]; 931 NN = trunc(dp[5]); //be sure that NN is an integer 932 Cp = dp[6]; 933 bkg = dp[7]; 934 962 sldb = dp[4]; 963 sld = dp[5]; 964 contr = sldb - sld; 965 NN = trunc(dp[6]); //be sure that NN is an integer 966 Cp = dp[7]; 967 bkg = dp[8]; 968 935 969 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); 936 970 937 971 NNint = (int)NN; //cast to an integer for the loop 972 973 // sprintf(buf, "qval = %g\r", qval); 974 // XOPNotice(buf); 975 938 976 ii=0; 939 977 Sq = 0.0; 940 978 for(ii=1;ii<(NNint-1);ii+=1) { 941 979 942 980 fii = (double)ii; //do I really need to do this? 943 981 944 982 temp = 0.0; 945 alpha = Cp/4.0/Pi/Pi*(log(Pi* ii) + Euler);983 alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler); 946 984 t1 = 2.0*dQ*dQ*dd*dd*alpha; 947 985 t2 = 2.0*qval*qval*dd*dd*alpha; 948 t3 = dQ*dQ*dd*dd* ii*ii;949 950 temp = 1.0- ii/NN;951 temp *= cos(dd*qval* ii/(1.0+t1));986 t3 = dQ*dQ*dd*dd*fii*fii; 987 988 temp = 1.0-fii/NN; 989 temp *= cos(dd*qval*fii/(1.0+t1)); 952 990 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 953 991 temp /= sqrt(1.0+t1); 954 992 955 993 Sq += temp; 956 994 } 957 995 958 996 Sq *= 2.0; 959 997 Sq += 1.0; 960 998 961 999 inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); 962 1000 963 1001 inten *= 1.0e8; // 1/A to 1/cm 964 1002 965 1003 return(inten+bkg); 966 1004 } … … 968 1006 969 1007 /* LamellarPS_HGX : calculates the form factor of a lamellar structure - with S(q) effects included 970 ------- 971 ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! 972 973 */ 1008 --- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the 1009 model is smeared just like any other function 1010 */ 974 1011 double 975 1012 LamellarPS_HG(double dp[], double q) … … 979 1016 double Pi,Euler,dQDefault,fii; 980 1017 int ii,NNint; 981 982 1018 1019 983 1020 Euler = 0.5772156649; // Euler's constant 984 dQDefault = 0.0; //0.0025; //[=] 1/A, q-resolution, default value 1021 // dQDefault = 0.0025; //[=] 1/A, q-resolution, default value 1022 dQDefault = 0.0; 985 1023 dQ = dQDefault; 986 1024 987 1025 Pi = 4.0*atan(1.0); 988 1026 qval= q; 989 1027 990 1028 scale = dp[0]; 991 1029 dd = dp[1]; … … 998 1036 Cp = dp[8]; 999 1037 bkg = dp[9]; 1000 1001 1038 1039 1002 1040 drh = SLD_H - SLD_S; 1003 1041 drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar 1004 1042 1005 1043 Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); 1006 1044 Pq *= Pq; 1007 1045 Pq *= 4.0/(qval*qval); 1008 1046 1009 1047 NNint = (int)NN; //cast to an integer for the loop 1010 1048 ii=0; 1011 1049 Sq = 0.0; 1012 1050 for(ii=1;ii<(NNint-1);ii+=1) { 1013 1051 1014 1052 fii = (double)ii; //do I really need to do this? 1015 1053 1016 1054 temp = 0.0; 1017 1055 alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); … … 1019 1057 t2 = 2.0*qval*qval*dd*dd*alpha; 1020 1058 t3 = dQ*dQ*dd*dd*ii*ii; 1021 1059 1022 1060 temp = 1.0-ii/NN; 1023 1061 temp *= cos(dd*qval*ii/(1.0+t1)); 1024 1062 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 1025 1063 temp /= sqrt(1.0+t1); 1026 1064 1027 1065 Sq += temp; 1028 1066 } 1029 1067 1030 1068 Sq *= 2.0; 1031 1069 Sq += 1.0; 1032 1070 1033 1071 inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); 1034 1072 1035 1073 inten *= 1.0e8; // 1/A to 1/cm 1036 1074 1037 1075 return(inten+bkg); 1038 1076 } … … 1048 1086 double inten, qval,Pq,drh,drt; 1049 1087 double Pi; 1050 1051 1088 1089 1052 1090 Pi = 4.0*atan(1.0); 1053 1091 qval= q; … … 1059 1097 slds = dp[5]; 1060 1098 bkg = dp[6]; 1061 1062 1099 1100 1063 1101 drh = sldh - slds; 1064 1102 drt = sldt - slds; //correction 13FEB06 by L.Porcar 1065 1103 1066 1104 Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); 1067 1105 Pq *= Pq; 1068 1106 Pq *= 4.0/(qval*qval); 1069 1107 1070 1108 inten = 2.0*Pi*scale*Pq/(qval*qval); //dimensionless... 1071 1109 1072 1110 inten /= 2.0*(delT+delH); //normalize by the bilayer thickness 1073 1111 1074 1112 inten *= 1.0e8; // 1/A to 1/cm 1075 1113 1076 1114 return(inten+bkg); 1077 1115 } … … 1084 1122 FlexExclVolCyl(double dp[], double q) 1085 1123 { 1086 double scale,L,B,bkg,rad,qr,cont ;1124 double scale,L,B,bkg,rad,qr,cont,sldc,slds; 1087 1125 double Pi,flex,crossSect,answer; 1088 1089 1126 1127 1090 1128 Pi = 4.0*atan(1.0); 1091 1129 1092 1130 scale = dp[0]; //make local copies in case memory moves 1093 1131 L = dp[1]; 1094 1132 B = dp[2]; 1095 1133 rad = dp[3]; 1096 cont = dp[4]; 1097 bkg = dp[5]; 1098 1099 1134 sldc = dp[4]; 1135 slds = dp[5]; 1136 cont = sldc-slds; 1137 bkg = dp[6]; 1138 1139 1100 1140 qr = q*rad; 1101 1141 1102 1142 flex = Sk_WR(q,L,B); 1103 1143 1104 1144 crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 1105 1145 flex *= crossSect; … … 1108 1148 flex *= 1.0e8; 1109 1149 answer = scale*flex + bkg; 1110 1150 1111 1151 return answer; 1112 1152 } … … 1119 1159 FlexCyl_Ellip(double dp[], double q) 1120 1160 { 1121 double scale,L,B,bkg,rad,qr,cont,ellRatio ;1161 double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc; 1122 1162 double Pi,flex,crossSect,answer; 1123 1124 1163 1164 1125 1165 Pi = 4.0*atan(1.0); 1126 1166 scale = dp[0]; //make local copies in case memory moves … … 1129 1169 rad = dp[3]; 1130 1170 ellRatio = dp[4]; 1131 cont = dp[5]; 1132 bkg = dp[6]; 1133 1171 sldc = dp[5]; 1172 slds = dp[6]; 1173 bkg = dp[7]; 1174 1175 cont = sldc - slds; 1134 1176 qr = q*rad; 1135 1177 1136 1178 flex = Sk_WR(q,L,B); 1137 1179 1138 1180 crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio)); 1139 1181 flex *= crossSect; … … 1142 1184 flex *= 1.0e8; 1143 1185 answer = scale*flex + bkg; 1144 1186 1145 1187 return answer; 1146 1188 } … … 1151 1193 double uplim,lolim,Pi,summ,arg,zi,yyy,answer; 1152 1194 int i,nord=76; 1153 1195 1154 1196 Pi = 4.0*atan(1.0); 1155 1197 lolim=0.0; 1156 1198 uplim=Pi/2.0; 1157 1199 summ=0.0; 1158 1200 1159 1201 for(i=0;i<nord;i++) { 1160 1202 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1167 1209 answer *= 2.0/Pi; 1168 1210 return(answer); 1169 1211 1170 1212 } 1171 1213 /* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x … … 1177 1219 { 1178 1220 int i; 1179 double scale,radius,length,pd, delrho,bkg,lb; //local variables of coefficient wave1221 double scale,radius,length,pd,bkg,lb,delrho,sldc,slds; //local variables of coefficient wave 1180 1222 int nord=20; //order of integration 1181 1223 double uplim,lolim; //upper and lower integration limits 1182 1224 double summ,zi,yyy,answer,Vpoly; //running tally of integration 1183 1225 double range,zz,Pi; 1184 1226 1185 1227 Pi = 4.0*atan(1.0); 1186 1228 range = 3.4; 1187 1229 1188 1230 summ = 0.0; //initialize intergral 1189 1231 scale = dp[0]; //make local copies in case memory moves … … 1192 1234 lb = dp[3]; 1193 1235 radius = dp[4]; 1194 delrho = dp[5]; 1195 bkg = dp[6]; 1196 1236 sldc = dp[5]; 1237 slds = dp[6]; 1238 bkg = dp[7]; 1239 1240 delrho = sldc - slds; 1197 1241 zz = (1.0/pd)*(1.0/pd) - 1.0; 1198 1242 1199 1243 lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution 1200 if(lolim<0 ) {1201 lolim = 0 ;1244 if(lolim<0.0) { 1245 lolim = 0.0; 1202 1246 } 1203 1247 if(pd>0.3) { … … 1205 1249 } 1206 1250 uplim = length*(1.0+range*pd); 1207 1251 1208 1252 for(i=0;i<nord;i++) { 1209 1253 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1211 1255 summ += yyy; 1212 1256 } 1213 1257 1214 1258 answer = (uplim-lolim)/2.0*summ; 1215 1259 //normalize by average cylinder volume (first moment), using the average length 1216 1260 Vpoly=Pi*radius*radius*length; 1217 1261 answer /= Vpoly; 1218 1262 1219 1263 answer *=delrho*delrho; 1220 1264 1221 1265 //convert to [cm-1] 1222 1266 answer *= 1.0e8; … … 1225 1269 // add in the background 1226 1270 answer += bkg; 1227 1271 1228 1272 return answer; 1229 1273 } … … 1237 1281 { 1238 1282 int i; 1239 double scale,radius,length,pd,delrho,bkg,lb ; //local variables of coefficient wave1283 double scale,radius,length,pd,delrho,bkg,lb,sldc,slds; //local variables of coefficient wave 1240 1284 int nord=76; //order of integration 1241 1285 double uplim,lolim; //upper and lower integration limits 1242 1286 double summ,zi,yyy,answer,Vpoly; //running tally of integration 1243 1287 double range,zz,Pi; 1244 1245 1288 1289 1246 1290 Pi = 4.0*atan(1.0); 1247 1291 range = 3.4; 1248 1292 1249 1293 summ = 0.0; //initialize intergral 1250 1294 1251 1295 scale = dp[0]; //make local copies in case memory moves 1252 1296 length = dp[1]; //radius … … 1254 1298 radius = dp[3]; 1255 1299 pd = dp[4]; 1256 delrho = dp[5]; 1257 bkg = dp[6]; 1258 1300 sldc = dp[5]; 1301 slds = dp[6]; 1302 bkg = dp[7]; 1303 1304 delrho = sldc-slds; 1259 1305 zz = (1.0/pd)*(1.0/pd) - 1.0; 1260 1306 1261 1307 lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution 1262 if(lolim<0 ) {1263 lolim = 0 ;1308 if(lolim<0.0) { 1309 lolim = 0.0; 1264 1310 } 1265 1311 if(pd>0.3) { … … 1267 1313 } 1268 1314 uplim = radius*(1.0+range*pd); 1269 1315 1270 1316 for(i=0;i<nord;i++) { 1271 1317 //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1275 1321 summ += yyy; 1276 1322 } 1277 1323 1278 1324 answer = (uplim-lolim)/2.0*summ; 1279 1325 //normalize by average cylinder volume (second moment), using the average radius 1280 1326 Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0); 1281 1327 answer /= Vpoly; 1282 1328 1283 1329 answer *=delrho*delrho; 1284 1330 1285 1331 //convert to [cm-1] 1286 1332 answer *= 1.0e8; … … 1289 1335 // add in the background 1290 1336 answer += bkg; 1291 1337 1292 1338 return answer; 1293 1339 } … … 1300 1346 double p1,p2,p1short,p2short,q0,qconnect; 1301 1347 double C,epsilon,ans,q0short,Sexvmodify,pi; 1302 1348 1303 1349 pi = 4.0*atan(1.0); 1304 1350 1305 1351 p1 = 4.12; 1306 1352 p2 = 4.42; … … 1309 1355 q0 = 3.1; 1310 1356 qconnect = q0/b; 1311 // 1357 // 1312 1358 q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 1313 1359 1314 1360 // 1315 1361 if(L/b > 10.0) { … … 1321 1367 } 1322 1368 // 1323 1369 1324 1370 if( L > 4*b ) { // Longer Chains 1325 1371 if (q*b <= 3.1) { //Modified by Yun on Oct. 15, … … 1328 1374 } else { //q(i)*b > 3.1 1329 1375 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); 1330 } 1376 } 1331 1377 } else { //L <= 4*b Shorter Chains 1332 1378 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { … … 1340 1386 } 1341 1387 } 1342 1388 1343 1389 return(ans); 1344 1390 //return(a2long(q, L, b, p1, p2, q0)); … … 1351 1397 double yy; 1352 1398 yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 1353 1399 1354 1400 return (yy); 1355 1401 } … … 1360 1406 { 1361 1407 double yy; 1362 1408 1363 1409 yy = Rgsquareshort(q,L,b)*q*q; 1364 1410 1365 1411 return (yy); 1366 1412 } … … 1380 1426 static double 1381 1427 Rgsquarezero(double q, double L, double b) 1382 { 1428 { 1383 1429 double yy; 1384 1430 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)))); 1385 1431 1386 1432 return (yy); 1387 1433 } … … 1390 1436 static double 1391 1437 Rgsquareshort(double q, double L, double b) 1392 { 1438 { 1393 1439 double yy; 1394 1440 yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 1395 1441 1396 1442 return (yy); 1397 1443 } … … 1403 1449 double yy; 1404 1450 yy = AlphaSquare(L/b)*L*b/6.0; 1405 1451 1406 1452 return (yy); 1407 1453 } … … 1410 1456 static double 1411 1457 AlphaSquare(double x) 1412 { 1458 { 1413 1459 double yy; 1414 1460 yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 1415 1461 1416 1462 return (yy); 1417 1463 } … … 1420 1466 static double 1421 1467 miu(double x) 1422 { 1468 { 1423 1469 double yy; 1424 1470 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))); 1425 1471 1426 1472 return (yy); 1427 1473 } … … 1430 1476 static double 1431 1477 Sdebye(double q, double L, double b) 1432 { 1478 { 1433 1479 double yy; 1434 1480 yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 1435 1481 1436 1482 return (yy); 1437 1483 } … … 1440 1486 static double 1441 1487 Sdebye1(double q, double L, double b) 1442 { 1488 { 1443 1489 double yy; 1444 1490 yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 1445 1491 1446 1492 return (yy); 1447 1493 } … … 1450 1496 static double 1451 1497 Sexv(double q, double L, double b) 1452 { 1498 { 1453 1499 double yy,C1,C2,C3,miu,Rg2; 1454 1500 C1=1.22; … … 1456 1502 C3=-1.651; 1457 1503 miu = 0.585; 1458 1504 1459 1505 Rg2 = Rgsquare(q,L,b); 1460 1506 1461 1507 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))); 1462 1508 1463 1509 return (yy); 1464 1510 } … … 1467 1513 static double 1468 1514 Sexvnew(double q, double L, double b) 1469 { 1515 { 1470 1516 double yy,C1,C2,C3,miu; 1471 1517 double del=1.05,C_star2,Rg2; 1472 1518 1473 1519 C1=1.22; 1474 1520 C2=0.4288; 1475 1521 C3=-1.651; 1476 1522 miu = 0.585; 1477 1523 1478 1524 //calculating the derivative to decide on the corection (cutoff) term? 1479 1525 // I have modified this from WRs original code 1480 1526 1481 1527 if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 1482 1528 C_star2 = 0.0; … … 1484 1530 C_star2 = 1.0; 1485 1531 } 1486 1532 1487 1533 Rg2 = Rgsquare(q,L,b); 1488 1534 1489 1535 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))); 1490 1536 1491 1537 return (yy); 1492 1538 } … … 1495 1541 static double 1496 1542 a2short(double q, double L, double b, double p1short, double p2short, double q0) 1497 { 1543 { 1498 1544 double yy,Rg2_sh; 1499 1545 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 1500 1546 double pi; 1501 1547 1502 1548 E = 2.718281828459045091; 1503 1549 pi = 4.0*atan(1.0); … … 1506 1552 t1 = ((q0*q0*Rg2_sh)/(b*b)); 1507 1553 Et1 = pow(E,t1); 1508 Emt1 =pow(E,-t1); 1554 Emt1 =pow(E,-t1); 1509 1555 q02 = q0*q0; 1510 1556 q0p = pow(q0,(-4.0 + p2short) ); 1511 1557 1512 1558 //E is the number e 1513 1559 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))))))); 1514 1560 1515 1561 return (yy); 1516 1562 } … … 1519 1565 static double 1520 1566 a1short(double q, double L, double b, double p1short, double p2short, double q0) 1521 { 1567 { 1522 1568 double yy,Rg2_sh; 1523 1569 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 1524 1570 double pi; 1525 1571 1526 1572 E = 2.718281828459045091; 1527 1573 pi = 4.0*atan(1.0); … … 1531 1577 t1 = ((q0*q0*Rg2_sh)/(b*b)); 1532 1578 Et1 = pow(E,t1); 1533 Emt1 =pow(E,-t1); 1579 Emt1 =pow(E,-t1); 1534 1580 q02 = q0*q0; 1535 1581 q0p = pow(q0,(-4.0 + p1short) ); … … 1547 1593 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 1548 1594 double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 1549 1595 1550 1596 pi = 4.0*atan(1.0); 1551 1597 E = 2.718281828459045091; … … 1555 1601 C = 1.0; 1556 1602 } 1557 1603 1558 1604 C1 = 1.22; 1559 1605 C2 = 0.4288; … … 1562 1608 C5 = 0.1477; 1563 1609 miu = 0.585; 1564 1610 1565 1611 Rg2 = Rgsquare(q,L,b); 1566 1612 Rg22 = Rg2*Rg2; … … 1572 1618 q04 = q03*q0; 1573 1619 q05 = q04*q0; 1574 1620 1575 1621 t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 1576 1622 … … 1581 1627 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); 1582 1628 1583 t5 = (2.0*b4*(((2 *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);1629 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); 1584 1630 1585 1631 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); … … 1602 1648 static double 1603 1649 sech_WR(double x) 1604 { 1650 { 1605 1651 return(1/cosh(x)); 1606 1652 } … … 1609 1655 static double 1610 1656 a1long(double q, double L, double b, double p1, double p2, double q0) 1611 { 1657 { 1612 1658 double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 1613 1659 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 1614 1660 double E,pi; 1615 1661 double b2,b3,b4,q02,q03,q04,q05,Rg22; 1616 1662 1617 1663 pi = 4.0*atan(1.0); 1618 1664 E = 2.718281828459045091; 1619 1665 1620 1666 if( L/b > 10.0) { 1621 1667 C = 3.06/pow((L/b),0.44); … … 1623 1669 C = 1.0; 1624 1670 } 1625 1671 1626 1672 C1 = 1.22; 1627 1673 C2 = 0.4288; … … 1630 1676 C5 = 0.1477; 1631 1677 miu = 0.585; 1632 1678 1633 1679 Rg2 = Rgsquare(q,L,b); 1634 1680 Rg22 = Rg2*Rg2; … … 1640 1686 q04 = q03*q0; 1641 1687 q05 = q04*q0; 1642 1688 1643 1689 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)))); 1644 1690 1645 1691 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)))))); 1646 1692 1647 1693 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)))); 1648 1694 1649 1695 t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 1650 1696 1651 1697 t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 1652 1698 1653 1699 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))); 1654 1700 1655 1701 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)); 1656 1702 1657 1703 t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 1658 1704 1659 1705 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)))))); 1660 1706 1661 1707 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)))))); 1662 1708 1663 1709 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)); 1664 1710 1665 1711 t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 1666 1712 1667 1713 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)))); 1668 1714 1669 1715 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)))))); 1670 1716 1671 1717 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)))); 1672 1673 1718 1719 1674 1720 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))))))))))); 1675 1721 1676 1722 return (yy); 1677 1723 } … … 1683 1729 // 1684 1730 // FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN 1685 // BY (53) AND (56,57) IN CHEN AND 1731 // BY (53) AND (56,57) IN CHEN AND 1686 1732 // KOTLARCHYK REFERENCE 1687 1733 // … … 1695 1741 1696 1742 Pi = 4.0*atan(1.0); 1697 1743 1698 1744 pi43=4.0/3.0*Pi; 1699 1745 aa = crmaj; … … 1719 1765 gfn = gfnc+gfnt; 1720 1766 gfn2 = gfn*gfn; 1721 1767 1722 1768 return (gfn2); 1723 1769 } … … 1729 1775 // 1730 1776 // <OBLATE ELLIPSOID> 1731 // function gfn4 for oblate ellipsoids 1777 // function gfn4 for oblate ellipsoids 1732 1778 double 1733 1779 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) … … 1760 1806 tgfn = gfnc+gfnt; 1761 1807 gfn4 = tgfn*tgfn; 1762 1808 1763 1809 return (gfn4); 1764 1810 } … … 1769 1815 double Pq,vcyl,dl; 1770 1816 double Pi,qr; 1771 1817 1772 1818 Pi = 4.0*atan(1.0); 1773 1819 qr = q*radius; 1774 1820 1775 1821 Pq = Sk_WR(q,zi,lb); //does not have cross section term 1776 1822 if (qr !=0){ 1777 1823 Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 1778 } //else Pk *=1;1824 } 1779 1825 vcyl=Pi*radius*radius*zi; 1780 1826 Pq *= vcyl*vcyl; 1781 1827 1782 1828 dl = SchulzPoint_cpr(zi,length,zz); 1783 return (Pq*dl); 1784 1829 return (Pq*dl); 1830 1785 1831 } 1786 1832 … … 1790 1836 double Pq,vcyl,dr; 1791 1837 double Pi,qr; 1792 1838 1793 1839 Pi = 4.0*atan(1.0); 1794 1840 qr = q*zi; 1795 1841 1796 1842 Pq = Sk_WR(q,Lc,Lb); //does not have cross section term 1797 1843 if (qr !=0){ … … 1801 1847 vcyl=Pi*zi*zi*Lc; 1802 1848 Pq *= vcyl*vcyl; 1803 1849 1804 1850 dr = SchulzPoint_cpr(zi,ravg,zz); 1805 return (Pq*dr); 1806 1851 return (Pq*dr); 1852 1807 1853 } 1808 1854 … … 1813 1859 double lolim,uplim,summ,yyy,zi; 1814 1860 int nord,i; 1815 1816 // set up the integration end points 1861 1862 // set up the integration end points 1817 1863 Pi = 4.0*atan(1.0); 1818 1864 nord = 76; 1819 lolim = 0 ;1865 lolim = 0.0; 1820 1866 uplim = Pi/2.0; 1821 1867 halfheight = length/2.0; 1822 1868 1823 1869 summ = 0.0; // initialize integral 1824 1870 i=0; … … 1828 1874 summ += yyy; 1829 1875 } 1830 1876 1831 1877 // calculate value of integral to return 1832 1878 answer = (uplim-lolim)/2.0*summ; … … 1836 1882 double 1837 1883 CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum) 1838 { 1884 { 1839 1885 // qq is the q-value for the calculation (1/A) 1840 1886 // radius is the core radius of the cylinder (A) 1841 1887 // radthick and facthick are the radial and face layer thicknesses 1842 1888 // rho(n) are the respective SLD's 1843 // length is the *Half* CORE-LENGTH of the cylinder 1889 // length is the *Half* CORE-LENGTH of the cylinder 1844 1890 // dum is the dummy variable for the integration (theta) 1845 1891 1846 1892 double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 1847 1893 double Pi; 1848 1849 Pi = 4.0*atan(1.0); 1850 1894 1895 Pi = 4.0*atan(1.0); 1896 1851 1897 dr1 = rhoc-rhos; 1852 1898 dr2 = rhos-rhosolv; 1853 1899 vol1 = Pi*rad*rad*(2.0*length); 1854 1900 vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 1855 1901 1856 1902 besarg1 = qq*rad*sin(dum); 1857 1903 besarg2 = qq*(rad+radthick)*sin(dum); … … 1873 1919 si1 = sin(sinarg1)/sinarg1; 1874 1920 } 1875 if ( besarg2 == 0.0){1921 if (sinarg2 == 0.0){ 1876 1922 si2 = 1.0; 1877 1923 }else{ … … 1884 1930 retval = ((t1+t2)*(t1+t2))*sin(dum); 1885 1931 return (retval); 1886 1932 1887 1933 } 1888 1934 … … 1894 1940 double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 1895 1941 double Pi; 1896 1942 1897 1943 Pi = 4.0*atan(1.0); 1898 1944 1899 1945 dr1 = rhoc-rhos; 1900 1946 dr2 = rhos-rhosolv; 1901 1947 vol1 = Pi*rcore*rcore*(2.0*length); 1902 1948 vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick); 1903 1949 1904 1950 besarg1 = qq*rcore*sin(dum); 1905 1951 besarg2 = qq*(rcore+thick)*sin(dum); … … 1922 1968 si1 = sin(sinarg1)/sinarg1; 1923 1969 } 1924 if ( besarg2 == 0.0){1970 if (sinarg2 == 0.0){ 1925 1971 si2 = 1.0; 1926 1972 }else{ … … 1932 1978 1933 1979 retval = ((t1+t2)*(t1+t2))*sin(dum); 1934 1980 1935 1981 return (retval); 1936 1982 } … … 1939 1985 Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen) 1940 1986 { 1941 1987 1942 1988 double halfheight,uplim,lolim,zi,summ,yyy,Pi; 1943 1989 double answer,dr,Vcyl; 1944 1990 int i,nord; 1945 1991 1946 1992 Pi = 4.0*atan(1.0); 1947 lolim = 0 ;1993 lolim = 0.0; 1948 1994 uplim = Pi/2.0; 1949 1995 halfheight = dumLen/2.0; 1950 1996 nord=20; 1951 1997 summ = 0.0; 1952 1998 1953 1999 //do the cylinder orientational average 1954 2000 for(i=0;i<nord;i++) { … … 1964 2010 Vcyl = Pi*radius*radius*dumLen; 1965 2011 answer *= Vcyl*Vcyl; 1966 2012 1967 2013 dr = SchulzPoint_cpr(dumLen,len_avg,zz); 1968 2014 return(dr*answer); … … 1972 2018 double 1973 2019 Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N) 1974 { 2020 { 1975 2021 // qq is the q-value for the calculation (1/A) 1976 2022 // rcore is the core radius of the cylinder (A) … … 1983 2029 double Pi; 1984 2030 int kk; 1985 2031 1986 2032 Pi = 4.0*atan(1.0); 1987 2033 1988 2034 dr1 = rhoc-rhosolv; 1989 2035 dr2 = rhol-rhosolv; 1990 2036 area = Pi*rcore*rcore; 1991 2037 totald=2.0*(thick+length); 1992 2038 1993 2039 besarg1 = qq*rcore*sin(dum); 1994 2040 besarg2 = qq*rcore*sin(dum); 1995 2041 1996 2042 sinarg1 = qq*length*cos(dum); 1997 2043 sinarg2 = qq*(length+thick)*cos(dum); … … 2012 2058 si1 = sin(sinarg1)/sinarg1; 2013 2059 } 2014 if ( besarg2 == 0.0){2060 if (sinarg2 == 0.0){ 2015 2061 si2 = 1.0; 2016 2062 }else{ … … 2022 2068 2023 2069 retval =((t1+t2)*(t1+t2))*sin(dum); 2024 2070 2025 2071 // loop for the structure facture S(q) 2026 2072 sqq=0.0; … … 2028 2074 dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0; 2029 2075 sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt); 2030 } 2031 2076 } 2077 2032 2078 // end of loop for S(q) 2033 2079 sqq=1.0+2.0*sqq/N; 2034 2035 2080 retval *= sqq; 2036 2081 2037 2082 return(retval); 2038 2083 } … … 2042 2087 Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad) 2043 2088 { 2044 2089 2045 2090 double halfheight,uplim,lolim,zi,summ,yyy,Pi; 2046 2091 double answer,dr,Vcyl; 2047 2092 int i,nord; 2048 2093 2049 2094 Pi = 4.0*atan(1.0); 2050 lolim = 0 ;2095 lolim = 0.0; 2051 2096 uplim = Pi/2.0; 2052 2097 halfheight = length/2.0; … … 2054 2099 nord=76; 2055 2100 summ = 0.0; 2056 2101 2057 2102 //do the cylinder orientational average 2058 2103 // for(i=0;i<nord;i++) { … … 2073 2118 Vcyl = Pi*dumRad*dumRad*length; 2074 2119 answer *= Vcyl*Vcyl; 2075 2120 2076 2121 dr = SchulzPoint_cpr(dumRad,radius,zz); 2077 2122 return(dr*answer); … … 2082 2127 { 2083 2128 double dr; 2084 2129 2085 2130 dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0)); 2086 2131 return(exp(dr)); … … 2095 2140 0.1208650973866179e-2,-0.5395239384953e-5}; 2096 2141 int j; 2097 2142 2098 2143 y=x=xx; 2099 2144 tmp=x+5.5; … … 2109 2154 { 2110 2155 double arg,nu,retval; //local variables 2111 2156 2112 2157 nu = nua/a; 2113 arg = qq*a*sqrt(1 +dum*dum*(nu*nu-1));2158 arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0)); 2114 2159 if (arg == 0.0){ 2115 2160 retval =1.0/3.0; … … 2127 2172 { 2128 2173 double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables 2129 2174 2130 2175 gamma = rcore/rshell; 2131 arg1 = qq*rshell*sqrt(1 -dum*dum); //1=shell (outer radius)2132 arg2 = qq*rcore*sqrt(1 -dum*dum); //2=core (inner radius)2176 arg1 = qq*rshell*sqrt(1.0-dum*dum); //1=shell (outer radius) 2177 arg2 = qq*rcore*sqrt(1.0-dum*dum); //2=core (inner radius) 2133 2178 if (arg1 == 0.0){ 2134 2179 lam1 = 1.0; … … 2151 2196 2152 2197 retval = psi*psi*t2*t2; 2153 2198 2154 2199 return(retval); 2155 2200 }//Function HolCylKernel() … … 2160 2205 // mu passed in is really mu*sqrt(1-sig^2) 2161 2206 double arg1,arg2,Pi,tmp1,tmp2; //local variables 2162 2207 2163 2208 Pi = 4.0*atan(1.0); 2164 2209 2165 2210 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 2166 2211 arg1 = (mu/2.0)*cos(Pi*uu/2.0); … … 2177 2222 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 2178 2223 } 2179 2224 2180 2225 return (tmp1*tmp2); 2181 2226 2182 2227 }//Function PPKernel() 2183 2228 … … 2186 2231 TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy) 2187 2232 { 2188 2233 2189 2234 double arg,val,pi; //local variables 2190 2235 2191 2236 pi = 4.0*atan(1.0); 2192 2237 2193 2238 arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2); 2194 2239 arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy); … … 2201 2246 } 2202 2247 return (val); 2203 2248 2204 2249 }//Function TriaxialKernel() 2205 2250 … … 2208 2253 CylKernel(double qq, double rr,double h, double theta) 2209 2254 { 2210 2255 2211 2256 // qq is the q-value for the calculation (1/A) 2212 2257 // rr is the radius of the cylinder (A) … … 2219 2264 siarg = qq * h * cos(theta); 2220 2265 bj =NR_BessJ1(besarg); 2221 2266 2222 2267 //* Computing 2nd power */ 2223 2268 d1 = sin(siarg); … … 2245 2290 2246 2291 return (retval); 2247 2292 2248 2293 }//Function CylKernel() 2249 2294 … … 2251 2296 EllipCylKernel(double qq, double ra,double nu, double theta) 2252 2297 { 2253 //this is the function LAMBDA1^2 in Feigin's notation 2298 //this is the function LAMBDA1^2 in Feigin's notation 2254 2299 // qq is the q-value for the calculation (1/A) 2255 2300 // ra is the transformed radius"a" in Feigin's notation 2256 2301 // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- 2257 2302 // theta is the dummy variable of the integration 2258 2259 double retval,arg; //Local variables 2260 2261 arg = qq*ra*sqrt((1 +nu*nu)/2+(1-nu*nu)*cos(theta)/2);2303 2304 double retval,arg; //Local variables 2305 2306 arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2); 2262 2307 if (arg == 0.0){ 2263 2308 retval = 1.0; … … 2268 2313 //square it 2269 2314 retval *= retval; 2270 2315 2271 2316 return(retval); 2272 2317 2273 2318 }//Function EllipCylKernel() 2274 2319 … … 2277 2322 double ax,z; 2278 2323 double xx,y,ans,ans1,ans2; 2279 2324 2280 2325 if ((ax=fabs(x)) < 8.0) { 2281 2326 y=x*x; … … 2297 2342 if (x < 0.0) ans = -ans; 2298 2343 } 2299 2344 2300 2345 return(ans); 2301 2346 } 2347 2348 /* Lamellar_ParaCrystal - Pedersen's model 2349 2350 */ 2351 double 2352 Lamellar_ParaCrystal(double w[], double q) 2353 { 2354 // Input (fitting) variables are: 2355 //[0] scale factor 2356 //[1] thickness 2357 //[2] number of layers 2358 //[3] spacing between layers 2359 //[4] polydispersity of spacing 2360 //[5] SLD lamellar 2361 //[6] SLD solvent 2362 //[7] incoherent background 2363 // give them nice names 2364 double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn; 2365 double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi; 2366 long n1,n2; 2367 2368 pi = 4.0*atan(1.0); 2369 scale = w[0]; 2370 th = w[1]; 2371 nl = w[2]; 2372 davg = w[3]; 2373 pd = w[4]; 2374 sldLayer = w[5]; 2375 sldSolvent = w[6]; 2376 bkg = w[7]; 2377 2378 contr = w[5] - w[6]; 2379 qval = q; 2380 2381 //get the fractional part of nl, to determine the "mixing" of N's 2382 2383 n1 = trunc(nl); //rounds towards zero 2384 n2 = n1 + 1; 2385 xn = (double)n2 - nl; //fractional contribution of n1 2386 2387 ww = exp(-qval*qval*pd*pd*davg*davg/2.0); 2388 2389 //calculate the n1 contribution 2390 an = paraCryst_an(ww,qval,davg,n1); 2391 Snq = paraCryst_sn(ww,qval,davg,n1,an); 2392 2393 Znq = xn*Snq; 2394 2395 //calculate the n2 contribution 2396 an = paraCryst_an(ww,qval,davg,n2); 2397 Snq = paraCryst_sn(ww,qval,davg,n2,an); 2398 2399 Znq += (1.0-xn)*Snq; 2400 2401 //and the independent contribution 2402 Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg)); 2403 2404 //the limit when NL approaches infinity 2405 // Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 2406 2407 xi = th/2.0; //use 1/2 the bilayer thickness 2408 Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 2409 2410 inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval); 2411 inten *= 1.0e8; 2412 2413 return(scale*inten+bkg); 2414 } 2415 2416 // functions for the lamellar paracrystal model 2417 double 2418 paraCryst_sn(double ww, double qval, double davg, long nl, double an) { 2419 2420 double Snq; 2421 2422 Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) ); 2423 2424 return(Snq); 2425 } 2426 2427 2428 double 2429 paraCryst_an(double ww, double qval, double davg, long nl) { 2430 2431 double an; 2432 2433 an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg); 2434 an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg); 2435 an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg); 2436 an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg); 2437 2438 return(an); 2439 } 2440 2441 2442 /* Spherocylinder : 2443 2444 Uses 76 pt Gaussian quadrature for both integrals 2445 */ 2446 double 2447 Spherocylinder(double w[], double x) 2448 { 2449 int i,j; 2450 double Pi; 2451 double scale,contr,bkg,sldc,slds; 2452 double len,rad,hDist,endRad; 2453 int nordi=76; //order of integration 2454 int nordj=76; 2455 double va,vb; //upper and lower integration limits 2456 double summ,zi,yyy,answer; //running tally of integration 2457 double summj,vaj,vbj,zij; //for the inner integration 2458 double SphCyl_tmp[7],arg1,arg2,inner,be; 2459 2460 2461 scale = w[0]; 2462 rad = w[1]; 2463 len = w[2]; 2464 sldc = w[3]; 2465 slds = w[4]; 2466 bkg = w[5]; 2467 2468 SphCyl_tmp[0] = w[0]; 2469 SphCyl_tmp[1] = w[1]; 2470 SphCyl_tmp[2] = w[2]; 2471 SphCyl_tmp[3] = w[1]; //end radius is same as cylinder radius 2472 SphCyl_tmp[4] = w[3]; 2473 SphCyl_tmp[5] = w[4]; 2474 SphCyl_tmp[6] = w[5]; 2475 2476 hDist = 0; //by definition for this model 2477 endRad = rad; 2478 2479 contr = sldc-slds; 2480 2481 Pi = 4.0*atan(1.0); 2482 va = 0.0; 2483 vb = Pi/2.0; //orintational average, outer integral 2484 vaj = -1.0*hDist/endRad; 2485 vbj = 1.0; //endpoints of inner integral 2486 2487 summ = 0.0; //initialize intergral 2488 2489 for(i=0;i<nordi;i++) { 2490 //setup inner integral over the ellipsoidal cross-section 2491 summj=0.0; 2492 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy 2493 2494 for(j=0;j<nordj;j++) { 2495 //20 gauss points for the inner integral 2496 zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy 2497 yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi); 2498 summj += yyy; 2499 } 2500 //now calculate the value of the inner integral 2501 inner = (vbj-vaj)/2.0*summj; 2502 inner *= 4.0*Pi*endRad*endRad*endRad; 2503 2504 //now calculate outer integrand 2505 arg1 = x*len/2.0*cos(zi); 2506 arg2 = x*rad*sin(zi); 2507 yyy = inner; 2508 2509 if(arg2 == 0) { 2510 be = 0.5; 2511 } else { 2512 be = NR_BessJ1(arg2)/arg2; 2513 } 2514 2515 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2516 yyy += Pi*rad*rad*len*2.0*be; 2517 } else { 2518 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2519 } 2520 yyy *= yyy; 2521 yyy *= sin(zi); // = |A(q)|^2*sin(theta) 2522 yyy *= Gauss76Wt[i]; 2523 summ += yyy; 2524 } //final scaling is done at the end of the function, after the NT_FP64 case 2525 2526 answer = (vb-va)/2.0*summ; 2527 2528 answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0; //divide by volume 2529 answer *= 1.0e8; //convert to cm^-1 2530 answer *= contr*contr; 2531 answer *= scale; 2532 answer += bkg; 2533 2534 return answer; 2535 } 2536 2537 2538 // inner integral of the sphereocylinder model, special case of hDist = 0 2539 // 2540 double 2541 SphCyl_kernel(double w[], double x, double tt, double theta) { 2542 2543 double val,arg1,arg2; 2544 double scale,bkg,sldc,slds; 2545 double len,rad,hDist,endRad,be; 2546 scale = w[0]; 2547 rad = w[1]; 2548 len = w[2]; 2549 endRad = w[3]; 2550 sldc = w[4]; 2551 slds = w[5]; 2552 bkg = w[6]; 2553 2554 hDist = 0.0; 2555 2556 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 2557 arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 2558 2559 if(arg2 == 0) { 2560 be = 0.5; 2561 } else { 2562 be = NR_BessJ1(arg2)/arg2; 2563 } 2564 val = cos(arg1)*(1.0-tt*tt)*be; 2565 2566 return(val); 2567 } 2568 2569 2570 /* Convex Lens : special case where L ~ 0 and hDist < 0 2571 2572 Uses 76 pt Gaussian quadrature for both integrals 2573 */ 2574 double 2575 ConvexLens(double w[], double x) 2576 { 2577 int i,j; 2578 double Pi; 2579 double scale,contr,bkg,sldc,slds; 2580 double len,rad,hDist,endRad; 2581 int nordi=76; //order of integration 2582 int nordj=76; 2583 double va,vb; //upper and lower integration limits 2584 double summ,zi,yyy,answer; //running tally of integration 2585 double summj,vaj,vbj,zij; //for the inner integration 2586 double CLens_tmp[7],arg1,arg2,inner,hh,be; 2587 2588 2589 scale = w[0]; 2590 rad = w[1]; 2591 // len = w[2] 2592 endRad = w[2]; 2593 sldc = w[3]; 2594 slds = w[4]; 2595 bkg = w[5]; 2596 2597 len = 0.01; 2598 2599 CLens_tmp[0] = w[0]; 2600 CLens_tmp[1] = w[1]; 2601 CLens_tmp[2] = len; //length is some small number, essentially zero 2602 CLens_tmp[3] = w[2]; 2603 CLens_tmp[4] = w[3]; 2604 CLens_tmp[5] = w[4]; 2605 CLens_tmp[6] = w[5]; 2606 2607 hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model 2608 2609 contr = sldc-slds; 2610 2611 Pi = 4.0*atan(1.0); 2612 va = 0.0; 2613 vb = Pi/2.0; //orintational average, outer integral 2614 vaj = -1.0*hDist/endRad; 2615 vbj = 1.0; //endpoints of inner integral 2616 2617 summ = 0.0; //initialize intergral 2618 2619 for(i=0;i<nordi;i++) { 2620 //setup inner integral over the ellipsoidal cross-section 2621 summj=0.0; 2622 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy 2623 2624 for(j=0;j<nordj;j++) { 2625 //20 gauss points for the inner integral 2626 zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy 2627 yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi); 2628 summj += yyy; 2629 } 2630 //now calculate the value of the inner integral 2631 inner = (vbj-vaj)/2.0*summj; 2632 inner *= 4.0*Pi*endRad*endRad*endRad; 2633 2634 //now calculate outer integrand 2635 arg1 = x*len/2.0*cos(zi); 2636 arg2 = x*rad*sin(zi); 2637 yyy = inner; 2638 2639 if(arg2 == 0) { 2640 be = 0.5; 2641 } else { 2642 be = NR_BessJ1(arg2)/arg2; 2643 } 2644 2645 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2646 yyy += Pi*rad*rad*len*2.0*be; 2647 } else { 2648 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2649 } 2650 yyy *= yyy; 2651 yyy *= sin(zi); // = |A(q)|^2*sin(theta) 2652 yyy *= Gauss76Wt[i]; 2653 summ += yyy; 2654 } //final scaling is done at the end of the function, after the NT_FP64 case 2655 2656 answer = (vb-va)/2.0*summ; 2657 2658 hh = fabs(hDist); //need positive value for spherical cap volume 2659 answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh)); //divide by volume 2660 answer *= 1.0e8; //convert to cm^-1 2661 answer *= contr*contr; 2662 answer *= scale; 2663 answer += bkg; 2664 2665 return answer; 2666 } 2667 2668 /* Capped Cylinder : special case where L is nonzero and hDist < 0 2669 2670 -- uses the same Kernel as the Convex Lens 2671 2672 Uses 76 pt Gaussian quadrature for both integrals 2673 */ 2674 double 2675 CappedCylinder(double w[], double x) 2676 { 2677 int i,j; 2678 double Pi; 2679 double scale,contr,bkg,sldc,slds; 2680 double len,rad,hDist,endRad; 2681 int nordi=76; //order of integration 2682 int nordj=76; 2683 double va,vb; //upper and lower integration limits 2684 double summ,zi,yyy,answer; //running tally of integration 2685 double summj,vaj,vbj,zij; //for the inner integration 2686 double arg1,arg2,inner,hh,be; 2687 2688 2689 scale = w[0]; 2690 rad = w[1]; 2691 len = w[2]; 2692 endRad = w[3]; 2693 sldc = w[4]; 2694 slds = w[5]; 2695 bkg = w[6]; 2696 2697 hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model 2698 2699 contr = sldc-slds; 2700 2701 Pi = 4.0*atan(1.0); 2702 va = 0.0; 2703 vb = Pi/2.0; //orintational average, outer integral 2704 vaj = -1.0*hDist/endRad; 2705 vbj = 1.0; //endpoints of inner integral 2706 2707 summ = 0.0; //initialize intergral 2708 2709 for(i=0;i<nordi;i++) { 2710 //setup inner integral over the ellipsoidal cross-section 2711 summj=0.0; 2712 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy 2713 2714 for(j=0;j<nordj;j++) { 2715 //20 gauss points for the inner integral 2716 zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy 2717 yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi); //uses the same kernel as ConvexLens, except here L != 0 2718 summj += yyy; 2719 } 2720 //now calculate the value of the inner integral 2721 inner = (vbj-vaj)/2.0*summj; 2722 inner *= 4.0*Pi*endRad*endRad*endRad; 2723 2724 //now calculate outer integrand 2725 arg1 = x*len/2.0*cos(zi); 2726 arg2 = x*rad*sin(zi); 2727 yyy = inner; 2728 2729 if(arg2 == 0) { 2730 be = 0.5; 2731 } else { 2732 be = NR_BessJ1(arg2)/arg2; 2733 } 2734 2735 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2736 yyy += Pi*rad*rad*len*2.0*be; 2737 } else { 2738 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2739 } 2740 2741 2742 2743 yyy *= yyy; 2744 yyy *= sin(zi); // = |A(q)|^2*sin(theta) 2745 yyy *= Gauss76Wt[i]; 2746 summ += yyy; 2747 } //final scaling is done at the end of the function, after the NT_FP64 case 2748 2749 answer = (vb-va)/2.0*summ; 2750 2751 hh = fabs(hDist); //need positive value for spherical cap volume 2752 answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh)); //divide by volume 2753 answer *= 1.0e8; //convert to cm^-1 2754 answer *= contr*contr; 2755 answer *= scale; 2756 answer += bkg; 2757 2758 return answer; 2759 } 2760 2761 2762 2763 // inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0 2764 // 2765 double 2766 ConvLens_kernel(double w[], double x, double tt, double theta) { 2767 2768 double val,arg1,arg2; 2769 double scale,bkg,sldc,slds; 2770 double len,rad,hDist,endRad,be; 2771 scale = w[0]; 2772 rad = w[1]; 2773 len = w[2]; 2774 endRad = w[3]; 2775 sldc = w[4]; 2776 slds = w[5]; 2777 bkg = w[6]; 2778 2779 hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); 2780 2781 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 2782 arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 2783 2784 if(arg2 == 0) { 2785 be = 0.5; 2786 } else { 2787 be = NR_BessJ1(arg2)/arg2; 2788 } 2789 2790 val = cos(arg1)*(1.0-tt*tt)*be; 2791 2792 return(val); 2793 } 2794 2795 2796 /* Dumbbell : special case where L ~ 0 and hDist > 0 2797 2798 Uses 76 pt Gaussian quadrature for both integrals 2799 */ 2800 double 2801 Dumbbell(double w[], double x) 2802 { 2803 int i,j; 2804 double Pi; 2805 double scale,contr,bkg,sldc,slds; 2806 double len,rad,hDist,endRad; 2807 int nordi=76; //order of integration 2808 int nordj=76; 2809 double va,vb; //upper and lower integration limits 2810 double summ,zi,yyy,answer; //running tally of integration 2811 double summj,vaj,vbj,zij; //for the inner integration 2812 double Dumb_tmp[7],arg1,arg2,inner,be; 2813 2814 2815 scale = w[0]; 2816 rad = w[1]; 2817 // len = w[2] 2818 endRad = w[2]; 2819 sldc = w[3]; 2820 slds = w[4]; 2821 bkg = w[5]; 2822 2823 len = 0.01; 2824 2825 Dumb_tmp[0] = w[0]; 2826 Dumb_tmp[1] = w[1]; 2827 Dumb_tmp[2] = len; //length is some small number, essentially zero 2828 Dumb_tmp[3] = w[2]; 2829 Dumb_tmp[4] = w[3]; 2830 Dumb_tmp[5] = w[4]; 2831 Dumb_tmp[6] = w[5]; 2832 2833 hDist = sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model 2834 2835 contr = sldc-slds; 2836 2837 Pi = 4.0*atan(1.0); 2838 va = 0.0; 2839 vb = Pi/2.0; //orintational average, outer integral 2840 vaj = -1.0*hDist/endRad; 2841 vbj = 1.0; //endpoints of inner integral 2842 2843 summ = 0.0; //initialize intergral 2844 2845 for(i=0;i<nordi;i++) { 2846 //setup inner integral over the ellipsoidal cross-section 2847 summj=0.0; 2848 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy 2849 2850 for(j=0;j<nordj;j++) { 2851 //20 gauss points for the inner integral 2852 zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy 2853 yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi); 2854 summj += yyy; 2855 } 2856 //now calculate the value of the inner integral 2857 inner = (vbj-vaj)/2.0*summj; 2858 inner *= 4.0*Pi*endRad*endRad*endRad; 2859 2860 //now calculate outer integrand 2861 arg1 = x*len/2.0*cos(zi); 2862 arg2 = x*rad*sin(zi); 2863 yyy = inner; 2864 2865 if(arg2 == 0) { 2866 be = 0.5; 2867 } else { 2868 be = NR_BessJ1(arg2)/arg2; 2869 } 2870 2871 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2872 yyy += Pi*rad*rad*len*2.0*be; 2873 } else { 2874 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2875 } 2876 yyy *= yyy; 2877 yyy *= sin(zi); // = |A(q)|^2*sin(theta) 2878 yyy *= Gauss76Wt[i]; 2879 summ += yyy; 2880 } //final scaling is done at the end of the function, after the NT_FP64 case 2881 2882 answer = (vb-va)/2.0*summ; 2883 2884 answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0); //divide by volume 2885 answer *= 1.0e8; //convert to cm^-1 2886 answer *= contr*contr; 2887 answer *= scale; 2888 answer += bkg; 2889 2890 return answer; 2891 } 2892 2893 2894 /* Barbell : "normal" case where L is nonzero 0 and hDist > 0 2895 2896 -- uses the same kernel as the Dumbbell case 2897 2898 Uses 76 pt Gaussian quadrature for both integrals 2899 */ 2900 double 2901 Barbell(double w[], double x) 2902 { 2903 int i,j; 2904 double Pi; 2905 double scale,contr,bkg,sldc,slds; 2906 double len,rad,hDist,endRad; 2907 int nordi=76; //order of integration 2908 int nordj=76; 2909 double va,vb; //upper and lower integration limits 2910 double summ,zi,yyy,answer; //running tally of integration 2911 double summj,vaj,vbj,zij; //for the inner integration 2912 double arg1,arg2,inner,be; 2913 2914 2915 scale = w[0]; 2916 rad = w[1]; 2917 len = w[2]; 2918 endRad = w[3]; 2919 sldc = w[4]; 2920 slds = w[5]; 2921 bkg = w[6]; 2922 2923 hDist = sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model 2924 2925 contr = sldc-slds; 2926 2927 Pi = 4.0*atan(1.0); 2928 va = 0.0; 2929 vb = Pi/2.0; //orintational average, outer integral 2930 vaj = -1.0*hDist/endRad; 2931 vbj = 1.0; //endpoints of inner integral 2932 2933 summ = 0.0; //initialize intergral 2934 2935 for(i=0;i<nordi;i++) { 2936 //setup inner integral over the ellipsoidal cross-section 2937 summj=0.0; 2938 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy 2939 2940 for(j=0;j<nordj;j++) { 2941 //20 gauss points for the inner integral 2942 zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy 2943 yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi); //uses the same Kernel as the Dumbbell, here L>0 2944 summj += yyy; 2945 } 2946 //now calculate the value of the inner integral 2947 inner = (vbj-vaj)/2.0*summj; 2948 inner *= 4.0*Pi*endRad*endRad*endRad; 2949 2950 //now calculate outer integrand 2951 arg1 = x*len/2.0*cos(zi); 2952 arg2 = x*rad*sin(zi); 2953 yyy = inner; 2954 2955 if(arg2 == 0) { 2956 be = 0.5; 2957 } else { 2958 be = NR_BessJ1(arg2)/arg2; 2959 } 2960 2961 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2962 yyy += Pi*rad*rad*len*2.0*be; 2963 } else { 2964 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2965 } 2966 yyy *= yyy; 2967 yyy *= sin(zi); // = |A(q)|^2*sin(theta) 2968 yyy *= Gauss76Wt[i]; 2969 summ += yyy; 2970 } //final scaling is done at the end of the function, after the NT_FP64 case 2971 2972 answer = (vb-va)/2.0*summ; 2973 2974 answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0); //divide by volume 2975 answer *= 1.0e8; //convert to cm^-1 2976 answer *= contr*contr; 2977 answer *= scale; 2978 answer += bkg; 2979 2980 return answer; 2981 } 2982 2983 2984 2985 // inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0 2986 // 2987 // inner integral of the Barbell model if L is nonzero 2988 // 2989 double 2990 Dumb_kernel(double w[], double x, double tt, double theta) { 2991 2992 double val,arg1,arg2; 2993 double scale,bkg,sldc,slds; 2994 double len,rad,hDist,endRad,be; 2995 scale = w[0]; 2996 rad = w[1]; 2997 len = w[2]; 2998 endRad = w[3]; 2999 sldc = w[4]; 3000 slds = w[5]; 3001 bkg = w[6]; 3002 3003 hDist = sqrt(fabs(endRad*endRad-rad*rad)); 3004 3005 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 3006 arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 3007 3008 if(arg2 == 0) { 3009 be = 0.5; 3010 } else { 3011 be = NR_BessJ1(arg2)/arg2; 3012 } 3013 val = cos(arg1)*(1.0-tt*tt)*be; 3014 3015 return(val); 3016 } 3017 3018 double PolyCoreBicelle(double dp[], double q) 3019 { 3020 int i; 3021 int nord = 20; 3022 double scale, length, sigma, bkg, radius, radthick, facthick; 3023 double rhoc, rhoh, rhor, rhosolv; 3024 double answer, Vpoly; 3025 double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy; 3026 3027 scale = dp[0]; 3028 radius = dp[1]; 3029 sigma = dp[2]; //sigma is the standard mean deviation 3030 length = dp[3]; 3031 radthick = dp[4]; 3032 facthick= dp[5]; 3033 rhoc = dp[6]; 3034 rhoh = dp[7]; 3035 rhor=dp[8]; 3036 rhosolv = dp[9]; 3037 bkg = dp[10]; 3038 3039 Pi = 4.0*atan(1.0); 3040 3041 lolim = exp(log(radius)-(4.*sigma)); 3042 if (lolim<0.0) { 3043 lolim=0.0; //to avoid numerical error when va<0 (-ve r value) 3044 } 3045 uplim = exp(log(radius)+(4.*sigma)); 3046 3047 summ = 0.0; 3048 Rsqrsumm = 0.0; 3049 3050 for(i=0;i<nord;i++) { 3051 rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 3052 AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); 3053 yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length); 3054 Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions 3055 summ += yyy; 3056 Rsqrsumm += Rsqryyy; 3057 } 3058 3059 answer = (uplim-lolim)/2.0*summ; 3060 Rsqr = (uplim-lolim)/2.0*Rsqrsumm; 3061 //normalize by average cylinder volume 3062 Vpoly = Pi*Rsqr*(length+2*facthick); 3063 answer /= Vpoly; 3064 //convert to [cm-1] 3065 answer *= 1.0e8; 3066 //Scale 3067 answer *= scale; 3068 // add in the background 3069 answer += bkg; 3070 3071 return answer; 3072 3073 } 3074 3075 double 3076 BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){ 3077 3078 double answer,halfheight,Pi; 3079 double lolim,uplim,summ,yyy,zi; 3080 int nord,i; 3081 3082 // set up the integration end points 3083 Pi = 4.0*atan(1.0); 3084 nord = 76; 3085 lolim = 0.0; 3086 uplim = Pi/2; 3087 halfheight = length/2.0; 3088 3089 summ = 0.0; // initialize integral 3090 i=0; 3091 for(i=0;i<nord;i++) { 3092 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 3093 yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi); 3094 summ += yyy; 3095 } 3096 3097 // calculate value of integral to return 3098 answer = (uplim-lolim)/2.0*summ; 3099 return(answer); 3100 } 3101 3102 double 3103 BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum) 3104 { 3105 double dr1,dr2,dr3; 3106 double besarg1,besarg2; 3107 double vol1,vol2,vol3; 3108 double sinarg1,sinarg2; 3109 double t1,t2,t3; 3110 double retval,si1,si2,be1,be2; 3111 3112 double Pi = 4.0*atan(1.0); 3113 3114 dr1 = rhoc-rhoh; 3115 dr2 = rhor-rhosolv; 3116 dr3= rhoh-rhor; 3117 vol1 = Pi*rad*rad*(2.0*length); 3118 vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 3119 vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick); 3120 besarg1 = qq*rad*sin(dum); 3121 besarg2 = qq*(rad+radthick)*sin(dum); 3122 sinarg1 = qq*length*cos(dum); 3123 sinarg2 = qq*(length+facthick)*cos(dum); 3124 3125 if(besarg1 == 0) { 3126 be1 = 0.5; 3127 } else { 3128 be1 = NR_BessJ1(besarg1)/besarg1; 3129 } 3130 if(besarg2 == 0) { 3131 be2 = 0.5; 3132 } else { 3133 be2 = NR_BessJ1(besarg2)/besarg2; 3134 } 3135 if(sinarg1 == 0) { 3136 si1 = 1.0; 3137 } else { 3138 si1 = sin(sinarg1)/sinarg1; 3139 } 3140 if(sinarg2 == 0) { 3141 si2 = 1.0; 3142 } else { 3143 si2 = sin(sinarg2)/sinarg2; 3144 } 3145 t1 = 2.0*vol1*dr1*si1*be1; 3146 t2 = 2.0*vol2*dr2*si2*be2; 3147 t3 = 2.0*vol3*dr3*si2*be1; 3148 3149 retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 3150 return(retval); 3151 3152 } 3153 3154 3155 double 3156 CSPPKernel(double dp[], double mu, double uu) 3157 { 3158 double aa,bb,cc, ta,tb,tc; 3159 double Vin,Vot,V1,V2; 3160 double rhoA,rhoB,rhoC, rhoP, rhosolv; 3161 double dr0, drA,drB, drC; 3162 double arg1,arg2,arg3,arg4,t1,t2, t3, t4; 3163 double Pi,retVal; 3164 3165 aa = dp[1]; 3166 bb = dp[2]; 3167 cc = dp[3]; 3168 ta = dp[4]; 3169 tb = dp[5]; 3170 tc = dp[6]; 3171 rhoA=dp[7]; 3172 rhoB=dp[8]; 3173 rhoC=dp[9]; 3174 rhoP=dp[10]; 3175 rhosolv=dp[11]; 3176 dr0=rhoP-rhosolv; 3177 drA=rhoA-rhosolv; 3178 drB=rhoB-rhosolv; 3179 drC=rhoC-rhosolv; 3180 Vin=(aa*bb*cc); 3181 Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc); 3182 V1=(2.0*ta*bb*cc); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 3183 V2=(2.0*aa*tb*cc); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 3184 aa = aa/bb; 3185 ta=(aa+2.0*ta)/bb; 3186 tb=(aa+2.0*tb)/bb; 3187 3188 Pi = 4.0*atan(1.0); 3189 3190 arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0); 3191 arg2 = (mu/2.0)*cos(Pi*uu/2.0); 3192 arg3= (mu*ta/2.0)*sin(Pi*uu/2.0); 3193 arg4= (mu*tb/2.0)*cos(Pi*uu/2.0); 3194 3195 if(arg1==0.0){ 3196 t1 = 1.0; 3197 } else { 3198 t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1) 3199 } 3200 if(arg2==0.0){ 3201 t2 = 1.0; 3202 } else { 3203 t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2) 3204 } 3205 if(arg3==0.0){ 3206 t3 = 1.0; 3207 } else { 3208 t3 = sin(arg3)/arg3; 3209 } 3210 if(arg4==0.0){ 3211 t4 = 1.0; 3212 } else { 3213 t4 = sin(arg4)/arg4; 3214 } 3215 retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 ); // correct FF : square of sum of phase factors 3216 return(retVal); 3217 3218 } 3219 3220 /* CSParallelepiped : calculates the form factor of a Parallelepiped with a core-shell structure 3221 -- different SLDs can be used for the face and rim 3222 3223 Uses 76 pt Gaussian quadrature for both integrals 3224 */ 3225 double 3226 CSParallelepiped(double dp[], double q) 3227 { 3228 int i,j; 3229 double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg; //local variables of coefficient wave 3230 int nordi=76; //order of integration 3231 int nordj=76; 3232 double va,vb; //upper and lower integration limits 3233 double summ,yyy,answer; //running tally of integration 3234 double summj,vaj,vbj; //for the inner integration 3235 double mu,mudum,arg,sigma,uu,vol; 3236 3237 3238 // Pi = 4.0*atan(1.0); 3239 va = 0.0; 3240 vb = 1.0; //orintational average, outer integral 3241 vaj = 0.0; 3242 vbj = 1.0; //endpoints of inner integral 3243 3244 summ = 0.0; //initialize intergral 3245 3246 scale = dp[0]; 3247 aa = dp[1]; 3248 bb = dp[2]; 3249 cc = dp[3]; 3250 ta = dp[4]; 3251 tb = dp[5]; 3252 tc = dp[6]; // is 0 at the moment 3253 rhoA=dp[7]; //rim A SLD 3254 rhoB=dp[8]; //rim B SLD 3255 rhoC=dp[9]; //rim C SLD 3256 rhoP = dp[10]; //Parallelpiped core SLD 3257 rhosolv=dp[11]; // Solvent SLD 3258 bkg = dp[12]; 3259 3260 mu = q*bb; 3261 vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc; //calculate volume before rescaling 3262 3263 // do the rescaling here, not in the kernel 3264 // normalize all WRT bb 3265 aa = aa/bb; 3266 cc = cc/bb; 3267 3268 for(i=0;i<nordi;i++) { 3269 //setup inner integral over the ellipsoidal cross-section 3270 summj=0.0; 3271 sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy 3272 3273 for(j=0;j<nordj;j++) { 3274 //76 gauss points for the inner integral 3275 uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy 3276 mudum = mu*sqrt(1.0-sigma*sigma); 3277 yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu); 3278 summj += yyy; 3279 } 3280 //now calculate the value of the inner integral 3281 answer = (vbj-vaj)/2.0*summj; 3282 3283 //finish the outer integral cc already scaled 3284 arg = mu*cc*sigma/2.0; 3285 if ( arg == 0.0 ) { 3286 answer *= 1.0; 3287 } else { 3288 answer *= sin(arg)*sin(arg)/arg/arg; 3289 } 3290 3291 //now sum up the outer integral 3292 yyy = Gauss76Wt[i] * answer; 3293 summ += yyy; 3294 } //final scaling is done at the end of the function, after the NT_FP64 case 3295 3296 answer = (vb-va)/2.0*summ; 3297 3298 //normalize by volume 3299 answer /= vol; 3300 //convert to [cm-1] 3301 answer *= 1.0e8; 3302 //Scale 3303 answer *= scale; 3304 // add in the background 3305 answer += bkg; 3306 3307 return answer; 3308 } 3309
Note: See TracChangeset
for help on using the changeset viewer.