Changeset 8e91f01 in sasview for sansmodels/src/libigor/libCylinder.c
- Timestamp:
- Aug 4, 2009 1:31:21 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:
- b4679de
- Parents:
- 42f193a
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/libigor/libCylinder.c
r34c3020 r8e91f01 25 25 double uplim,lolim; //upper and lower integration limits 26 26 double summ,zi,yyy,answer,vcyl; //running tally of integration 27 27 28 28 Pi = 4.0*atan(1.0); 29 29 lolim = 0; 30 30 uplim = Pi/2.0; 31 31 32 32 summ = 0.0; //initialize intergral 33 33 34 34 scale = dp[0]; //make local copies in case memory moves 35 35 radius = dp[1]; … … 43 43 summ += yyy; 44 44 } 45 45 46 46 answer = (uplim-lolim)/2.0*summ; 47 47 // Multiply by contrast^2 … … 57 57 // add in the background 58 58 answer += bkg; 59 59 60 60 return answer; 61 61 } … … 81 81 double summ,zi,yyy,answer,vell; //running tally of integration 82 82 double summj,vaj,vbj,zij,arg; //for the inner integration 83 83 84 84 Pi = 4.0*atan(1.0); 85 85 va = 0; … … 87 87 vaj=0; 88 88 vbj=Pi; //endpoints of inner integral 89 89 90 90 summ = 0.0; //initialize intergral 91 91 92 92 scale = dp[0]; //make local copies in case memory moves 93 93 ra = dp[1]; … … 111 111 //divide integral by Pi 112 112 answer /=Pi; 113 113 114 114 //now calculate outer integral 115 115 arg = q*length*zi/2; … … 130 130 // add in the background 131 131 answer += bkg; 132 132 133 133 return answer; 134 134 } … … 156 156 double summ,zi,yyy,answer,vell; //running tally of integration 157 157 double summj,vaj,vbj,zij,arg; //for the inner integration 158 158 159 159 Pi = 4.0*atan(1.0); 160 160 va = 0; … … 162 162 vaj=0; 163 163 vbj=Pi; //endpoints of inner integral 164 164 165 165 summ = 0.0; //initialize intergral 166 166 167 167 scale = dp[0]; //make local copies in case memory moves 168 168 ra = dp[1]; … … 186 186 //divide integral by Pi 187 187 answer /=Pi; 188 188 189 189 //now calculate outer integral 190 190 arg = q*length*zi/2; … … 192 192 summ += yyy; 193 193 } 194 194 195 195 answer = (vb-va)/2.0*summ; 196 196 // Multiply by contrast^2 … … 205 205 answer *= scale; 206 206 // add in the background 207 answer += bkg; 208 207 answer += bkg; 208 209 209 return answer; 210 210 } … … 231 231 double summ,zi,yyy,answer; //running tally of integration 232 232 double summj,vaj,vbj,zij; //for the inner integration 233 233 234 234 Pi = 4.0*atan(1.0); 235 235 va = 0; … … 237 237 vaj = 0; 238 238 vbj = 1; //endpoints of inner integral 239 239 240 240 summ = 0.0; //initialize intergral 241 241 242 242 scale = dp[0]; //make local copies in case memory moves 243 243 aa = dp[1]; … … 258 258 //now calculate the value of the inner integral 259 259 answer = (vbj-vaj)/2.0*summj; 260 260 261 261 //now calculate outer integral 262 262 yyy = Gauss76Wt[i] * answer; 263 263 summ += yyy; 264 264 } //final scaling is done at the end of the function, after the NT_FP64 case 265 265 266 266 answer = (vb-va)/2.0*summ; 267 267 // Multiply by contrast^2 … … 275 275 // add in the background 276 276 answer += bkg; 277 277 278 278 return answer; 279 279 } … … 301 301 double summj,vaj,vbj; //for the inner integration 302 302 double mu,mudum,arg,sigma,uu,vol; 303 304 303 304 305 305 // Pi = 4.0*atan(1.0); 306 306 va = 0; … … 308 308 vaj = 0; 309 309 vbj = 1; //endpoints of inner integral 310 310 311 311 summ = 0.0; //initialize intergral 312 312 313 313 scale = dp[0]; //make local copies in case memory moves 314 314 aa = dp[1]; … … 317 317 delrho = dp[4]; 318 318 bkg = dp[5]; 319 319 320 320 mu = q*bb; 321 321 vol = aa*bb*cc; … … 323 323 aa = aa/bb; 324 324 cc = cc/bb; 325 325 326 326 for(i=0;i<nordi;i++) { 327 327 //setup inner integral over the ellipsoidal cross-section 328 328 summj=0; 329 329 sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy 330 330 331 331 for(j=0;j<nordj;j++) { 332 332 //76 gauss points for the inner integral … … 338 338 //now calculate the value of the inner integral 339 339 answer = (vbj-vaj)/2.0*summj; 340 340 341 341 arg = mu*cc*sigma/2; 342 342 if ( arg == 0 ) { … … 345 345 answer *= sin(arg)*sin(arg)/arg/arg; 346 346 } 347 347 348 348 //now sum up the outer integral 349 349 yyy = Gauss76Wt[i] * answer; 350 350 summ += yyy; 351 351 } //final scaling is done at the end of the function, after the NT_FP64 case 352 352 353 353 answer = (vb-va)/2.0*summ; 354 354 // Multiply by contrast^2 … … 362 362 // add in the background 363 363 answer += bkg; 364 364 365 365 return answer; 366 366 } … … 385 385 double va,vb,zi; //upper and lower integration limits 386 386 double summ,answer,pi; //running tally of integration 387 387 388 388 pi = 4.0*atan(1.0); 389 389 va = 0; 390 390 vb = 1; //limits of numerical integral 391 391 392 392 summ = 0.0; //initialize intergral 393 393 394 394 scale = dp[0]; //make local copies in case memory moves 395 395 rcore = dp[1]; … … 398 398 delrho = dp[4]; 399 399 bkg = dp[5]; 400 400 401 401 for(i=0;i<nord;i++) { 402 402 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 403 403 summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi); 404 404 } 405 405 406 406 answer = (vb-va)/2.0*summ; 407 407 // Multiply by contrast^2 … … 415 415 // add in the background 416 416 answer += bkg; 417 417 418 418 return answer; 419 419 } … … 438 438 double va,vb,zi; //upper and lower integration limits 439 439 double summ,answer,pi; //running tally of integration 440 440 441 441 pi = 4.0*atan(1.0); 442 442 va = 0; 443 443 vb = 1; //limits of numerical integral 444 444 445 445 summ = 0.0; //initialize intergral 446 446 447 447 scale = dp[0]; //make local copies in case memory moves 448 448 nua = dp[1]; … … 450 450 delrho = dp[3]; 451 451 bkg = dp[4]; 452 452 453 453 for(i=0;i<nord;i++) { 454 454 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 455 455 summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi); 456 456 } 457 457 458 458 answer = (vb-va)/2.0*summ; 459 459 // Multiply by contrast^2 … … 467 467 // add in the background 468 468 answer += bkg; 469 469 470 470 return answer; 471 471 } … … 485 485 double summ,zi,yyy,answer,Vpoly; //running tally of integration 486 486 double range,zz,Pi; 487 487 488 488 Pi = 4.0*atan(1.0); 489 489 range = 3.4; 490 490 491 491 summ = 0.0; //initialize intergral 492 492 493 493 scale = dp[0]; //make local copies in case memory moves 494 494 radius = dp[1]; … … 497 497 delrho = dp[4]; 498 498 bkg = dp[5]; 499 499 500 500 zz = (1.0/pd)*(1.0/pd) - 1.0; 501 501 502 502 lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution 503 503 if(lolim<0) { … … 508 508 } 509 509 uplim = radius*(1.0+range*pd); 510 510 511 511 for(i=0;i<nord;i++) { 512 512 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 514 514 summ += yyy; 515 515 } 516 516 517 517 answer = (uplim-lolim)/2.0*summ; 518 518 //normalize by average cylinder volume … … 526 526 // add in the background 527 527 answer += bkg; 528 528 529 529 return answer; 530 530 } … … 543 543 double summ,zi,yyy,answer,Vpoly; //running tally of integration 544 544 double range,zz,Pi; 545 546 545 546 547 547 Pi = 4.0*atan(1.0); 548 548 range = 3.4; 549 549 550 550 summ = 0.0; //initialize intergral 551 551 552 552 scale = dp[0]; //make local copies in case memory moves 553 553 radius = dp[1]; … … 556 556 delrho = dp[4]; 557 557 bkg = dp[5]; 558 558 559 559 zz = (1.0/pd)*(1.0/pd) - 1.0; 560 560 561 561 lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution 562 562 if(lolim<0) { … … 567 567 } 568 568 uplim = length*(1.0+range*pd); 569 569 570 570 for(i=0;i<nord;i++) { 571 571 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 573 573 summ += yyy; 574 574 } 575 575 576 576 answer = (uplim-lolim)/2.0*summ; 577 577 //normalize by average cylinder volume (first moment) … … 585 585 // add in the background 586 586 answer += bkg; 587 587 588 588 return answer; 589 589 } … … 603 603 double summ,zi,yyy,answer,Vcyl; //running tally of integration 604 604 double Pi; 605 606 Pi = 4.0*atan(1.0); 607 605 606 Pi = 4.0*atan(1.0); 607 608 608 lolim = 0.0; 609 609 uplim = Pi/2.0; 610 610 611 611 summ = 0.0; //initialize intergral 612 612 613 613 scale = dp[0]; //make local copies in case memory moves 614 614 rcore = dp[1]; … … 619 619 rhosolv = dp[6]; 620 620 bkg = dp[7]; 621 621 622 622 halfheight = length/2.0; 623 623 624 624 for(i=0;i<nord;i++) { 625 625 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 627 627 summ += yyy; 628 628 } 629 629 630 630 answer = (uplim-lolim)/2.0*summ; 631 // length is the total core length 631 // length is the total core length 632 632 Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick); 633 633 answer /= Vcyl; … … 638 638 // add in the background 639 639 answer += bkg; 640 640 641 641 return answer; 642 642 } … … 657 657 double summ,yyy,answer,Vpoly; //running tally of integration 658 658 double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; 659 660 Pi = 4.0*atan(1.0); 661 659 660 Pi = 4.0*atan(1.0); 661 662 662 summ = 0.0; //initialize intergral 663 663 Rsqrsumm = 0.0; 664 664 665 665 scale = dp[0]; 666 666 radius = dp[1]; … … 673 673 rhosolv = dp[8]; 674 674 bkg = dp[9]; 675 675 676 676 lolim = exp(log(radius)-(4.*sigma)); 677 677 if (lolim<0) { … … 679 679 } 680 680 uplim = exp(log(radius)+(4.*sigma)); 681 681 682 682 for(i=0;i<nord;i++) { 683 683 rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 688 688 Rsqrsumm += Rsqryyy; 689 689 } 690 690 691 691 answer = (uplim-lolim)/2.0*summ; 692 692 Rsqr = (uplim-lolim)/2.0*Rsqrsumm; … … 700 700 // add in the background 701 701 answer += bkg; 702 702 703 703 return answer; 704 704 } … … 717 717 double summ,zi,yyy,answer,oblatevol; //running tally of integration 718 718 double Pi; 719 720 Pi = 4.0*atan(1.0); 721 719 720 Pi = 4.0*atan(1.0); 721 722 722 lolim = 0.0; 723 723 uplim = 1.0; 724 724 725 725 summ = 0.0; //initialize intergral 726 727 726 727 728 728 scale = dp[0]; //make local copies in case memory moves 729 729 crmaj = dp[1]; … … 733 733 delpc = dp[5]; 734 734 delps = dp[6]; 735 bkg = dp[7]; 736 735 bkg = dp[7]; 736 737 737 for(i=0;i<nord;i++) { 738 738 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 740 740 summ += yyy; 741 741 } 742 742 743 743 answer = (uplim-lolim)/2.0*summ; 744 744 // normalize by particle volume 745 745 oblatevol = 4*Pi/3*trmaj*trmaj*trmin; 746 746 answer /= oblatevol; 747 747 748 748 //convert to [cm-1] 749 749 answer *= 1.0e8; … … 752 752 // add in the background 753 753 answer += bkg; 754 754 755 755 return answer; 756 756 } … … 769 769 double summ,zi,yyy,answer,prolatevol; //running tally of integration 770 770 double Pi; 771 772 Pi = 4.0*atan(1.0); 773 771 772 Pi = 4.0*atan(1.0); 773 774 774 lolim = 0.0; 775 775 uplim = 1.0; 776 776 777 777 summ = 0.0; //initialize intergral 778 778 779 779 scale = dp[0]; //make local copies in case memory moves 780 780 crmaj = dp[1]; … … 784 784 delpc = dp[5]; 785 785 delps = dp[6]; 786 bkg = dp[7]; 787 786 bkg = dp[7]; 787 788 788 for(i=0;i<nord;i++) { 789 789 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 791 791 summ += yyy; 792 792 } 793 793 794 794 answer = (uplim-lolim)/2.0*summ; 795 795 // normalize by particle volume 796 796 prolatevol = 4*Pi/3*trmaj*trmin*trmin; 797 797 answer /= prolatevol; 798 798 799 799 //convert to [cm-1] 800 800 answer *= 1.0e8; … … 803 803 // add in the background 804 804 answer += bkg; 805 805 806 806 return answer; 807 807 } … … 820 820 int nord=76; //order of integration 821 821 double Pi; 822 823 824 Pi = 4.0*atan(1.0); 825 822 823 824 Pi = 4.0*atan(1.0); 825 826 826 va = 0.0; 827 827 vb = Pi/2.0; 828 828 829 829 summ = 0.0; //initialize intergral 830 830 831 831 scale = dp[0]; 832 832 rcore = dp[1]; … … 839 839 gsd = dp[8]; 840 840 bkg = dp[9]; 841 841 842 842 d=2.0*thick+length; 843 843 halfheight = length/2.0; 844 844 845 845 for(i=0;i<nord;i++) { 846 846 zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0; … … 848 848 summ += yyy; 849 849 } 850 850 851 851 answer = (vb-va)/2.0*summ; 852 // length is the total core length 852 // length is the total core length 853 853 vcyl=Pi*rcore*rcore*(2.0*thick+length)*N; 854 854 answer /= vcyl; … … 859 859 // add in the background 860 860 answer += bkg; 861 861 862 862 return answer; 863 863 } … … 873 873 double inten, qval,Pq; 874 874 double Pi; 875 876 875 876 877 877 Pi = 4.0*atan(1.0); 878 878 scale = dp[0]; … … 882 882 bkg = dp[4]; 883 883 qval = q; 884 884 885 885 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); 886 886 887 887 inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... 888 888 889 889 inten /= del; //normalize by the thickness (in A) 890 890 891 891 inten *= 1.0e8; // 1/A to 1/cm 892 892 893 893 return(inten+bkg); 894 894 } … … 906 906 double Pi,Euler,dQDefault,fii; 907 907 int ii,NNint; 908 908 909 909 Euler = 0.5772156649; // Euler's constant 910 dQDefault = 0.0 025; //[=] 1/A, q-resolution, default value910 dQDefault = 0.0;//0.0025; //[=] 1/A, q-resolution, default value 911 911 dQ = dQDefault; 912 912 913 913 Pi = 4.0*atan(1.0); 914 914 qval = q; 915 915 916 916 scale = dp[0]; 917 917 dd = dp[1]; … … 922 922 Cp = dp[6]; 923 923 bkg = dp[7]; 924 924 925 925 Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); 926 926 927 927 NNint = (int)NN; //cast to an integer for the loop 928 928 ii=0; 929 929 Sq = 0.0; 930 930 for(ii=1;ii<(NNint-1);ii+=1) { 931 931 932 932 fii = (double)ii; //do I really need to do this? 933 933 934 934 temp = 0.0; 935 935 alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); … … 937 937 t2 = 2.0*qval*qval*dd*dd*alpha; 938 938 t3 = dQ*dQ*dd*dd*ii*ii; 939 939 940 940 temp = 1.0-ii/NN; 941 941 temp *= cos(dd*qval*ii/(1.0+t1)); 942 942 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 943 943 temp /= sqrt(1.0+t1); 944 944 945 945 Sq += temp; 946 946 } 947 947 948 948 Sq *= 2.0; 949 949 Sq += 1.0; 950 950 951 951 inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); 952 952 953 953 inten *= 1.0e8; // 1/A to 1/cm 954 954 955 955 return(inten+bkg); 956 956 } … … 969 969 double Pi,Euler,dQDefault,fii; 970 970 int ii,NNint; 971 972 971 972 973 973 Euler = 0.5772156649; // Euler's constant 974 974 dQDefault = 0.0025; //[=] 1/A, q-resolution, default value 975 975 dQ = dQDefault; 976 976 977 977 Pi = 4.0*atan(1.0); 978 978 qval= q; 979 979 980 980 scale = dp[0]; 981 981 dd = dp[1]; … … 988 988 Cp = dp[8]; 989 989 bkg = dp[9]; 990 991 990 991 992 992 drh = SLD_H - SLD_S; 993 993 drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar 994 994 995 995 Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); 996 996 Pq *= Pq; 997 997 Pq *= 4.0/(qval*qval); 998 998 999 999 NNint = (int)NN; //cast to an integer for the loop 1000 1000 ii=0; 1001 1001 Sq = 0.0; 1002 1002 for(ii=1;ii<(NNint-1);ii+=1) { 1003 1003 1004 1004 fii = (double)ii; //do I really need to do this? 1005 1005 1006 1006 temp = 0.0; 1007 1007 alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); … … 1009 1009 t2 = 2.0*qval*qval*dd*dd*alpha; 1010 1010 t3 = dQ*dQ*dd*dd*ii*ii; 1011 1011 1012 1012 temp = 1.0-ii/NN; 1013 1013 temp *= cos(dd*qval*ii/(1.0+t1)); 1014 1014 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 1015 1015 temp /= sqrt(1.0+t1); 1016 1016 1017 1017 Sq += temp; 1018 1018 } 1019 1019 1020 1020 Sq *= 2.0; 1021 1021 Sq += 1.0; 1022 1022 1023 1023 inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); 1024 1024 1025 1025 inten *= 1.0e8; // 1/A to 1/cm 1026 1026 1027 1027 return(inten+bkg); 1028 1028 } … … 1038 1038 double inten, qval,Pq,drh,drt; 1039 1039 double Pi; 1040 1041 1040 1041 1042 1042 Pi = 4.0*atan(1.0); 1043 1043 qval= q; … … 1049 1049 slds = dp[5]; 1050 1050 bkg = dp[6]; 1051 1052 1051 1052 1053 1053 drh = sldh - slds; 1054 1054 drt = sldt - slds; //correction 13FEB06 by L.Porcar 1055 1055 1056 1056 Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); 1057 1057 Pq *= Pq; 1058 1058 Pq *= 4.0/(qval*qval); 1059 1059 1060 1060 inten = 2.0*Pi*scale*Pq/(qval*qval); //dimensionless... 1061 1061 1062 1062 inten /= 2.0*(delT+delH); //normalize by the bilayer thickness 1063 1063 1064 1064 inten *= 1.0e8; // 1/A to 1/cm 1065 1065 1066 1066 return(inten+bkg); 1067 1067 } … … 1076 1076 double scale,L,B,bkg,rad,qr,cont; 1077 1077 double Pi,flex,crossSect,answer; 1078 1079 1080 Pi = 4.0*atan(1.0); 1081 1078 1079 1080 Pi = 4.0*atan(1.0); 1081 1082 1082 scale = dp[0]; //make local copies in case memory moves 1083 1083 L = dp[1]; … … 1086 1086 cont = dp[4]; 1087 1087 bkg = dp[5]; 1088 1089 1088 1089 1090 1090 qr = q*rad; 1091 1091 1092 1092 flex = Sk_WR(q,L,B); 1093 1093 1094 1094 crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 1095 1095 flex *= crossSect; … … 1098 1098 flex *= 1.0e8; 1099 1099 answer = scale*flex + bkg; 1100 1100 1101 1101 return answer; 1102 1102 } … … 1111 1111 double scale,L,B,bkg,rad,qr,cont,ellRatio; 1112 1112 double Pi,flex,crossSect,answer; 1113 1114 1113 1114 1115 1115 Pi = 4.0*atan(1.0); 1116 1116 scale = dp[0]; //make local copies in case memory moves … … 1121 1121 cont = dp[5]; 1122 1122 bkg = dp[6]; 1123 1123 1124 1124 qr = q*rad; 1125 1125 1126 1126 flex = Sk_WR(q,L,B); 1127 1127 1128 1128 crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio)); 1129 1129 flex *= crossSect; … … 1132 1132 flex *= 1.0e8; 1133 1133 answer = scale*flex + bkg; 1134 1134 1135 1135 return answer; 1136 1136 } … … 1141 1141 double uplim,lolim,Pi,summ,arg,zi,yyy,answer; 1142 1142 int i,nord=76; 1143 1143 1144 1144 Pi = 4.0*atan(1.0); 1145 1145 lolim=0.0; 1146 1146 uplim=Pi/2.0; 1147 1147 summ=0.0; 1148 1148 1149 1149 for(i=0;i<nord;i++) { 1150 1150 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1157 1157 answer *= 2.0/Pi; 1158 1158 return(answer); 1159 1159 1160 1160 } 1161 1161 /* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x … … 1172 1172 double summ,zi,yyy,answer,Vpoly; //running tally of integration 1173 1173 double range,zz,Pi; 1174 1174 1175 1175 Pi = 4.0*atan(1.0); 1176 1176 range = 3.4; 1177 1177 1178 1178 summ = 0.0; //initialize intergral 1179 1179 scale = dp[0]; //make local copies in case memory moves … … 1184 1184 delrho = dp[5]; 1185 1185 bkg = dp[6]; 1186 1186 1187 1187 zz = (1.0/pd)*(1.0/pd) - 1.0; 1188 1188 1189 1189 lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution 1190 1190 if(lolim<0) { … … 1195 1195 } 1196 1196 uplim = length*(1.0+range*pd); 1197 1197 1198 1198 for(i=0;i<nord;i++) { 1199 1199 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1201 1201 summ += yyy; 1202 1202 } 1203 1203 1204 1204 answer = (uplim-lolim)/2.0*summ; 1205 1205 //normalize by average cylinder volume (first moment), using the average length 1206 1206 Vpoly=Pi*radius*radius*length; 1207 1207 answer /= Vpoly; 1208 1208 1209 1209 answer *=delrho*delrho; 1210 1210 1211 1211 //convert to [cm-1] 1212 1212 answer *= 1.0e8; … … 1215 1215 // add in the background 1216 1216 answer += bkg; 1217 1217 1218 1218 return answer; 1219 1219 } … … 1232 1232 double summ,zi,yyy,answer,Vpoly; //running tally of integration 1233 1233 double range,zz,Pi; 1234 1235 1234 1235 1236 1236 Pi = 4.0*atan(1.0); 1237 1237 range = 3.4; 1238 1238 1239 1239 summ = 0.0; //initialize intergral 1240 1240 1241 1241 scale = dp[0]; //make local copies in case memory moves 1242 1242 length = dp[1]; //radius … … 1246 1246 delrho = dp[5]; 1247 1247 bkg = dp[6]; 1248 1248 1249 1249 zz = (1.0/pd)*(1.0/pd) - 1.0; 1250 1250 1251 1251 lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution 1252 1252 if(lolim<0) { … … 1257 1257 } 1258 1258 uplim = radius*(1.0+range*pd); 1259 1259 1260 1260 for(i=0;i<nord;i++) { 1261 1261 //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; … … 1265 1265 summ += yyy; 1266 1266 } 1267 1267 1268 1268 answer = (uplim-lolim)/2.0*summ; 1269 1269 //normalize by average cylinder volume (second moment), using the average radius 1270 1270 Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0); 1271 1271 answer /= Vpoly; 1272 1272 1273 1273 answer *=delrho*delrho; 1274 1274 1275 1275 //convert to [cm-1] 1276 1276 answer *= 1.0e8; … … 1279 1279 // add in the background 1280 1280 answer += bkg; 1281 1281 1282 1282 return answer; 1283 1283 } … … 1290 1290 double p1,p2,p1short,p2short,q0,qconnect; 1291 1291 double C,epsilon,ans,q0short,Sexvmodify,pi; 1292 1292 1293 1293 pi = 4.0*atan(1.0); 1294 1294 1295 1295 p1 = 4.12; 1296 1296 p2 = 4.42; … … 1299 1299 q0 = 3.1; 1300 1300 qconnect = q0/b; 1301 // 1301 // 1302 1302 q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 1303 1303 1304 1304 // 1305 1305 if(L/b > 10.0) { … … 1311 1311 } 1312 1312 // 1313 1313 1314 1314 if( L > 4*b ) { // Longer Chains 1315 1315 if (q*b <= 3.1) { //Modified by Yun on Oct. 15, … … 1318 1318 } else { //q(i)*b > 3.1 1319 1319 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); 1320 } 1320 } 1321 1321 } else { //L <= 4*b Shorter Chains 1322 1322 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { … … 1330 1330 } 1331 1331 } 1332 1332 1333 1333 return(ans); 1334 1334 //return(a2long(q, L, b, p1, p2, q0)); … … 1341 1341 double yy; 1342 1342 yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 1343 1343 1344 1344 return (yy); 1345 1345 } … … 1350 1350 { 1351 1351 double yy; 1352 1352 1353 1353 yy = Rgsquareshort(q,L,b)*q*q; 1354 1354 1355 1355 return (yy); 1356 1356 } … … 1370 1370 static double 1371 1371 Rgsquarezero(double q, double L, double b) 1372 { 1372 { 1373 1373 double yy; 1374 1374 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)))); 1375 1375 1376 1376 return (yy); 1377 1377 } … … 1380 1380 static double 1381 1381 Rgsquareshort(double q, double L, double b) 1382 { 1382 { 1383 1383 double yy; 1384 1384 yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 1385 1385 1386 1386 return (yy); 1387 1387 } … … 1393 1393 double yy; 1394 1394 yy = AlphaSquare(L/b)*L*b/6.0; 1395 1395 1396 1396 return (yy); 1397 1397 } … … 1400 1400 static double 1401 1401 AlphaSquare(double x) 1402 { 1402 { 1403 1403 double yy; 1404 1404 yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 1405 1405 1406 1406 return (yy); 1407 1407 } … … 1410 1410 static double 1411 1411 miu(double x) 1412 { 1412 { 1413 1413 double yy; 1414 1414 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))); 1415 1415 1416 1416 return (yy); 1417 1417 } … … 1420 1420 static double 1421 1421 Sdebye(double q, double L, double b) 1422 { 1422 { 1423 1423 double yy; 1424 1424 yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 1425 1425 1426 1426 return (yy); 1427 1427 } … … 1430 1430 static double 1431 1431 Sdebye1(double q, double L, double b) 1432 { 1432 { 1433 1433 double yy; 1434 1434 yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 1435 1435 1436 1436 return (yy); 1437 1437 } … … 1440 1440 static double 1441 1441 Sexv(double q, double L, double b) 1442 { 1442 { 1443 1443 double yy,C1,C2,C3,miu,Rg2; 1444 1444 C1=1.22; … … 1446 1446 C3=-1.651; 1447 1447 miu = 0.585; 1448 1448 1449 1449 Rg2 = Rgsquare(q,L,b); 1450 1450 1451 1451 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))); 1452 1452 1453 1453 return (yy); 1454 1454 } … … 1457 1457 static double 1458 1458 Sexvnew(double q, double L, double b) 1459 { 1459 { 1460 1460 double yy,C1,C2,C3,miu; 1461 1461 double del=1.05,C_star2,Rg2; 1462 1462 1463 1463 C1=1.22; 1464 1464 C2=0.4288; 1465 1465 C3=-1.651; 1466 1466 miu = 0.585; 1467 1467 1468 1468 //calculating the derivative to decide on the corection (cutoff) term? 1469 1469 // I have modified this from WRs original code 1470 1470 1471 1471 if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 1472 1472 C_star2 = 0.0; … … 1474 1474 C_star2 = 1.0; 1475 1475 } 1476 1476 1477 1477 Rg2 = Rgsquare(q,L,b); 1478 1478 1479 1479 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))); 1480 1480 1481 1481 return (yy); 1482 1482 } … … 1485 1485 static double 1486 1486 a2short(double q, double L, double b, double p1short, double p2short, double q0) 1487 { 1487 { 1488 1488 double yy,Rg2_sh; 1489 1489 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 1490 1490 double pi; 1491 1491 1492 1492 E = 2.718281828459045091; 1493 1493 pi = 4.0*atan(1.0); … … 1496 1496 t1 = ((q0*q0*Rg2_sh)/(b*b)); 1497 1497 Et1 = pow(E,t1); 1498 Emt1 =pow(E,-t1); 1498 Emt1 =pow(E,-t1); 1499 1499 q02 = q0*q0; 1500 1500 q0p = pow(q0,(-4.0 + p2short) ); 1501 1501 1502 1502 //E is the number e 1503 1503 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))))))); 1504 1504 1505 1505 return (yy); 1506 1506 } … … 1509 1509 static double 1510 1510 a1short(double q, double L, double b, double p1short, double p2short, double q0) 1511 { 1511 { 1512 1512 double yy,Rg2_sh; 1513 1513 double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 1514 1514 double pi; 1515 1515 1516 1516 E = 2.718281828459045091; 1517 1517 pi = 4.0*atan(1.0); … … 1521 1521 t1 = ((q0*q0*Rg2_sh)/(b*b)); 1522 1522 Et1 = pow(E,t1); 1523 Emt1 =pow(E,-t1); 1523 Emt1 =pow(E,-t1); 1524 1524 q02 = q0*q0; 1525 1525 q0p = pow(q0,(-4.0 + p1short) ); 1526 1527 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)))))); 1528 1526 1527 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)))))); 1528 1529 1529 return(yy); 1530 1530 } … … 1537 1537 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 1538 1538 double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 1539 1539 1540 1540 pi = 4.0*atan(1.0); 1541 1541 E = 2.718281828459045091; … … 1545 1545 C = 1.0; 1546 1546 } 1547 1547 1548 1548 C1 = 1.22; 1549 1549 C2 = 0.4288; … … 1552 1552 C5 = 0.1477; 1553 1553 miu = 0.585; 1554 1554 1555 1555 Rg2 = Rgsquare(q,L,b); 1556 1556 Rg22 = Rg2*Rg2; … … 1562 1562 q04 = q03*q0; 1563 1563 q05 = q04*q0; 1564 1564 1565 1565 t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 1566 1566 1567 1567 t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L; 1568 1568 1569 1569 t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5); 1570 1570 1571 1571 t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 1572 1572 1573 1573 t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 1574 1574 1575 1575 t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 1576 1576 1577 1577 t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu)); 1578 1578 1579 1579 t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 1580 1580 1581 1581 t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L; 1582 1582 1583 1583 t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 1584 1585 1584 1585 1586 1586 yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 1587 1587 1588 1588 return (yy); 1589 1589 } … … 1592 1592 static double 1593 1593 sech_WR(double x) 1594 { 1594 { 1595 1595 return(1/cosh(x)); 1596 1596 } … … 1599 1599 static double 1600 1600 a1long(double q, double L, double b, double p1, double p2, double q0) 1601 { 1601 { 1602 1602 double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 1603 1603 double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 1604 1604 double E,pi; 1605 1605 double b2,b3,b4,q02,q03,q04,q05,Rg22; 1606 1606 1607 1607 pi = 4.0*atan(1.0); 1608 1608 E = 2.718281828459045091; 1609 1609 1610 1610 if( L/b > 10.0) { 1611 1611 C = 3.06/pow((L/b),0.44); … … 1613 1613 C = 1.0; 1614 1614 } 1615 1615 1616 1616 C1 = 1.22; 1617 1617 C2 = 0.4288; … … 1620 1620 C5 = 0.1477; 1621 1621 miu = 0.585; 1622 1622 1623 1623 Rg2 = Rgsquare(q,L,b); 1624 1624 Rg22 = Rg2*Rg2; … … 1630 1630 q04 = q03*q0; 1631 1631 q05 = q04*q0; 1632 1632 1633 1633 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)))); 1634 1634 1635 1635 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)))))); 1636 1636 1637 1637 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)))); 1638 1638 1639 1639 t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 1640 1640 1641 1641 t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 1642 1642 1643 1643 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))); 1644 1644 1645 1645 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)); 1646 1646 1647 1647 t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 1648 1648 1649 1649 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)))))); 1650 1650 1651 1651 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)))))); 1652 1652 1653 1653 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)); 1654 1654 1655 1655 t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 1656 1656 1657 1657 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)))); 1658 1658 1659 1659 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)))))); 1660 1660 1661 1661 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)))); 1662 1663 1662 1663 1664 1664 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))))))))))); 1665 1665 1666 1666 return (yy); 1667 1667 } … … 1673 1673 // 1674 1674 // FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN 1675 // BY (53) AND (56,57) IN CHEN AND 1675 // BY (53) AND (56,57) IN CHEN AND 1676 1676 // KOTLARCHYK REFERENCE 1677 1677 // … … 1683 1683 // local variables 1684 1684 double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi; 1685 1686 Pi = 4.0*atan(1.0); 1687 1685 1686 Pi = 4.0*atan(1.0); 1687 1688 1688 pi43=4.0/3.0*Pi; 1689 1689 aa = crmaj; … … 1699 1699 gfn = gfnc+gfnt; 1700 1700 gfn2 = gfn*gfn; 1701 1701 1702 1702 return (gfn2); 1703 1703 } … … 1709 1709 // 1710 1710 // <OBLATE ELLIPSOID> 1711 // function gfn4 for oblate ellipsoids 1711 // function gfn4 for oblate ellipsoids 1712 1712 double 1713 1713 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) … … 1715 1715 // local variables 1716 1716 double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 1717 1717 1718 1718 Pi = 4.0*atan(1.0); 1719 1719 pi43=4.0/3.0*Pi; … … 1730 1730 tgfn = gfnc+gfnt; 1731 1731 gfn4 = tgfn*tgfn; 1732 1732 1733 1733 return (gfn4); 1734 1734 } … … 1739 1739 double Pq,vcyl,dl; 1740 1740 double Pi,qr; 1741 1741 1742 1742 Pi = 4.0*atan(1.0); 1743 1743 qr = q*radius; 1744 1744 1745 1745 Pq = Sk_WR(q,zi,lb); //does not have cross section term 1746 1746 Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 1747 1747 1748 1748 vcyl=Pi*radius*radius*zi; 1749 1749 Pq *= vcyl*vcyl; 1750 1750 1751 1751 dl = SchulzPoint_cpr(zi,length,zz); 1752 return (Pq*dl); 1753 1752 return (Pq*dl); 1753 1754 1754 } 1755 1755 … … 1759 1759 double Pq,vcyl,dr; 1760 1760 double Pi,qr; 1761 1761 1762 1762 Pi = 4.0*atan(1.0); 1763 1763 qr = q*zi; 1764 1764 1765 1765 Pq = Sk_WR(q,Lc,Lb); //does not have cross section term 1766 1766 Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 1767 1767 1768 1768 vcyl=Pi*zi*zi*Lc; 1769 1769 Pq *= vcyl*vcyl; 1770 1770 1771 1771 dr = SchulzPoint_cpr(zi,ravg,zz); 1772 return (Pq*dr); 1773 1772 return (Pq*dr); 1773 1774 1774 } 1775 1775 … … 1780 1780 double lolim,uplim,summ,yyy,zi; 1781 1781 int nord,i; 1782 1783 // set up the integration end points 1782 1783 // set up the integration end points 1784 1784 Pi = 4.0*atan(1.0); 1785 1785 nord = 76; … … 1787 1787 uplim = Pi/2; 1788 1788 halfheight = length/2.0; 1789 1789 1790 1790 summ = 0.0; // initialize integral 1791 1791 i=0; … … 1795 1795 summ += yyy; 1796 1796 } 1797 1797 1798 1798 // calculate value of integral to return 1799 1799 answer = (uplim-lolim)/2.0*summ; … … 1803 1803 double 1804 1804 CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum) 1805 { 1805 { 1806 1806 // qq is the q-value for the calculation (1/A) 1807 1807 // radius is the core radius of the cylinder (A) 1808 1808 // radthick and facthick are the radial and face layer thicknesses 1809 1809 // rho(n) are the respective SLD's 1810 // length is the *Half* CORE-LENGTH of the cylinder 1810 // length is the *Half* CORE-LENGTH of the cylinder 1811 1811 // dum is the dummy variable for the integration (theta) 1812 1812 1813 1813 double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 1814 1814 double Pi; 1815 1816 Pi = 4.0*atan(1.0); 1817 1815 1816 Pi = 4.0*atan(1.0); 1817 1818 1818 dr1 = rhoc-rhos; 1819 1819 dr2 = rhos-rhosolv; 1820 1820 vol1 = Pi*rad*rad*(2.0*length); 1821 1821 vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 1822 1822 1823 1823 besarg1 = qq*rad*sin(dum); 1824 1824 besarg2 = qq*(rad+radthick)*sin(dum); 1825 1825 sinarg1 = qq*length*cos(dum); 1826 1826 sinarg2 = qq*(length+facthick)*cos(dum); 1827 1827 1828 1828 t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 1829 1829 t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 1830 1830 1831 1831 retval = ((t1+t2)*(t1+t2))*sin(dum); 1832 1832 return (retval); 1833 1833 1834 1834 } 1835 1835 … … 1838 1838 CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum) 1839 1839 { 1840 1840 1841 1841 double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 1842 1842 double Pi; 1843 1843 1844 1844 Pi = 4.0*atan(1.0); 1845 1845 1846 1846 dr1 = rhoc-rhos; 1847 1847 dr2 = rhos-rhosolv; 1848 1848 vol1 = Pi*rcore*rcore*(2.0*length); 1849 1849 vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick); 1850 1850 1851 1851 besarg1 = qq*rcore*sin(dum); 1852 1852 besarg2 = qq*(rcore+thick)*sin(dum); 1853 1853 sinarg1 = qq*length*cos(dum); 1854 1854 sinarg2 = qq*(length+thick)*cos(dum); 1855 1855 1856 1856 t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 1857 1857 t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 1858 1858 1859 1859 retval = ((t1+t2)*(t1+t2))*sin(dum); 1860 1860 1861 1861 return (retval); 1862 1862 } … … 1865 1865 Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen) 1866 1866 { 1867 1867 1868 1868 double halfheight,uplim,lolim,zi,summ,yyy,Pi; 1869 1869 double answer,dr,Vcyl; 1870 1870 int i,nord; 1871 1871 1872 1872 Pi = 4.0*atan(1.0); 1873 1873 lolim = 0; … … 1876 1876 nord=20; 1877 1877 summ = 0.0; 1878 1878 1879 1879 //do the cylinder orientational average 1880 1880 for(i=0;i<nord;i++) { … … 1890 1890 Vcyl = Pi*radius*radius*dumLen; 1891 1891 answer *= Vcyl*Vcyl; 1892 1892 1893 1893 dr = SchulzPoint_cpr(dumLen,len_avg,zz); 1894 1894 return(dr*answer); … … 1898 1898 double 1899 1899 Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N) 1900 { 1900 { 1901 1901 // qq is the q-value for the calculation (1/A) 1902 1902 // rcore is the core radius of the cylinder (A) … … 1904 1904 // length is the *Half* CORE-LENGTH of the cylinder = L (A) 1905 1905 // dum is the dummy variable for the integration (x in Feigin's notation) 1906 1907 //Local variables 1906 1907 //Local variables 1908 1908 double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; 1909 1909 double Pi; 1910 1910 int kk; 1911 1912 Pi = 4.0*atan(1.0); 1913 1911 1912 Pi = 4.0*atan(1.0); 1913 1914 1914 dr1 = rhoc-rhosolv; 1915 1915 dr2 = rhol-rhosolv; 1916 1916 area = Pi*rcore*rcore; 1917 1917 totald=2.0*(thick+length); 1918 1918 1919 1919 besarg1 = qq*rcore*sin(dum); 1920 1920 besarg2 = qq*rcore*sin(dum); 1921 1921 1922 1922 sinarg1 = qq*length*cos(dum); 1923 1923 sinarg2 = qq*(length+thick)*cos(dum); 1924 1924 1925 1925 t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1); 1926 1926 t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2); 1927 1927 1928 1928 retval =((t1+t2)*(t1+t2))*sin(dum); 1929 1929 1930 1930 // loop for the structure facture S(q) 1931 1931 sqq=0.0; … … 1933 1933 dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0; 1934 1934 sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt); 1935 } 1936 1935 } 1936 1937 1937 // end of loop for S(q) 1938 1938 sqq=1.0+2.0*sqq/N; 1939 1939 retval *= sqq; 1940 1940 1941 1941 return(retval); 1942 1942 } … … 1946 1946 Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad) 1947 1947 { 1948 1948 1949 1949 double halfheight,uplim,lolim,zi,summ,yyy,Pi; 1950 1950 double answer,dr,Vcyl; 1951 1951 int i,nord; 1952 1952 1953 1953 Pi = 4.0*atan(1.0); 1954 1954 lolim = 0; … … 1958 1958 nord=76; 1959 1959 summ = 0.0; 1960 1960 1961 1961 //do the cylinder orientational average 1962 1962 // for(i=0;i<nord;i++) { … … 1977 1977 Vcyl = Pi*dumRad*dumRad*length; 1978 1978 answer *= Vcyl*Vcyl; 1979 1979 1980 1980 dr = SchulzPoint_cpr(dumRad,radius,zz); 1981 1981 return(dr*answer); … … 1986 1986 { 1987 1987 double dr; 1988 1988 1989 1989 dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0)); 1990 1990 return(exp(dr)); … … 1999 1999 0.1208650973866179e-2,-0.5395239384953e-5}; 2000 2000 int j; 2001 2001 2002 2002 y=x=xx; 2003 2003 tmp=x+5.5; … … 2013 2013 { 2014 2014 double arg,nu,retval; //local variables 2015 2015 2016 2016 nu = nua/a; 2017 2017 arg = qq*a*sqrt(1+dum*dum*(nu*nu-1)); 2018 2018 2019 2019 retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); 2020 2020 retval *= retval; 2021 2021 retval *= 9; 2022 2022 2023 2023 return(retval); 2024 2024 }//Function EllipsoidKernel() … … 2028 2028 { 2029 2029 double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables 2030 2030 2031 2031 gamma = rcore/rshell; 2032 2032 arg1 = qq*rshell*sqrt(1-dum*dum); //1=shell (outer radius) … … 2035 2035 lam2 = 2*NR_BessJ1(arg2)/arg2; 2036 2036 psi = 1/(1-gamma*gamma)*(lam1 - gamma*gamma*lam2); //SRK 10/19/00 2037 2037 2038 2038 sinarg = qq*length*dum/2; 2039 2039 t2 = sin(sinarg)/sinarg; 2040 2040 2041 2041 retval = psi*psi*t2*t2; 2042 2042 2043 2043 return(retval); 2044 2044 }//Function HolCylKernel() … … 2049 2049 // mu passed in is really mu*sqrt(1-sig^2) 2050 2050 double arg1,arg2,Pi,tmp1,tmp2; //local variables 2051 2051 2052 2052 Pi = 4.0*atan(1.0); 2053 2053 2054 2054 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 2055 2055 arg1 = (mu/2)*cos(Pi*uu/2); … … 2060 2060 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 2061 2061 } 2062 2062 2063 2063 if (arg2==0) { 2064 2064 tmp2 = 1; … … 2066 2066 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 2067 2067 } 2068 2068 2069 2069 return (tmp1*tmp2); 2070 2070 2071 2071 }//Function PPKernel() 2072 2072 … … 2075 2075 TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy) 2076 2076 { 2077 2077 2078 2078 double arg,val,pi; //local variables 2079 2079 2080 2080 pi = 4.0*atan(1.0); 2081 2081 2082 2082 arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2); 2083 2083 arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy); … … 2085 2085 arg = q*sqrt(arg); 2086 2086 val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); 2087 2087 2088 2088 return (val); 2089 2089 2090 2090 }//Function TriaxialKernel() 2091 2091 … … 2094 2094 CylKernel(double qq, double rr,double h, double theta) 2095 2095 { 2096 2096 2097 2097 // qq is the q-value for the calculation (1/A) 2098 2098 // rr is the radius of the cylinder (A) 2099 2099 // h is the HALF-LENGTH of the cylinder = L/2 (A) 2100 2101 double besarg,bj,retval,d1,t1,b1,t2,b2; //Local variables 2102 2103 2100 2101 double besarg,bj,retval,d1,t1,b1,t2,b2; //Local variables 2102 2103 2104 2104 besarg = qq*rr*sin(theta); 2105 2105 2106 2106 bj =NR_BessJ1(besarg); 2107 2107 2108 2108 //* Computing 2nd power */ 2109 2109 d1 = sin(qq * h * cos(theta)); … … 2119 2119 b2 = d1 * d1; 2120 2120 retval = t1 * t2 / b1 / b2; 2121 2121 2122 2122 return (retval); 2123 2123 2124 2124 }//Function CylKernel() 2125 2125 … … 2127 2127 EllipCylKernel(double qq, double ra,double nu, double theta) 2128 2128 { 2129 //this is the function LAMBDA1^2 in Feigin's notation 2129 //this is the function LAMBDA1^2 in Feigin's notation 2130 2130 // qq is the q-value for the calculation (1/A) 2131 2131 // ra is the transformed radius"a" in Feigin's notation 2132 2132 // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- 2133 2133 // theta is the dummy variable of the integration 2134 2135 double retval,arg; //Local variables 2136 2134 2135 double retval,arg; //Local variables 2136 2137 2137 arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2); 2138 2138 2139 2139 retval = 2*NR_BessJ1(arg)/arg; 2140 2140 2141 2141 //square it 2142 2142 retval *= retval; 2143 2143 2144 2144 return(retval); 2145 2145 2146 2146 }//Function EllipCylKernel() 2147 2147 … … 2150 2150 double ax,z; 2151 2151 double xx,y,ans,ans1,ans2; 2152 2152 2153 2153 if ((ax=fabs(x)) < 8.0) { 2154 2154 y=x*x; … … 2170 2170 if (x < 0.0) ans = -ans; 2171 2171 } 2172 2172 2173 2173 return(ans); 2174 2174 }
Note: See TracChangeset
for help on using the changeset viewer.