Changeset 634 for sans/XOP_Dev
 Timestamp:
 Mar 3, 2010 9:08:42 AM (13 years ago)
 Location:
 sans/XOP_Dev/SANSAnalysis/lib
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
r632 r634 1919 1919 si1 = sin(sinarg1)/sinarg1; 1920 1920 } 1921 if ( besarg2 == 0.0){1921 if (sinarg2 == 0.0){ 1922 1922 si2 = 1.0; 1923 1923 }else{ … … 1968 1968 si1 = sin(sinarg1)/sinarg1; 1969 1969 } 1970 if ( besarg2 == 0.0){1970 if (sinarg2 == 0.0){ 1971 1971 si2 = 1.0; 1972 1972 }else{ … … 2058 2058 si1 = sin(sinarg1)/sinarg1; 2059 2059 } 2060 if ( besarg2 == 0.0){2060 if (sinarg2 == 0.0){ 2061 2061 si2 = 1.0; 2062 2062 }else{ … … 2456 2456 double summ,zi,yyy,answer; //running tally of integration 2457 2457 double summj,vaj,vbj,zij; //for the inner integration 2458 double SphCyl_tmp[7],arg1,arg2,inner ;2458 double SphCyl_tmp[7],arg1,arg2,inner,be; 2459 2459 2460 2460 … … 2506 2506 arg2 = x*rad*sin(zi); 2507 2507 yyy = inner; 2508 2509 if(arg2 == 0) { 2510 be = 0.5; 2511 } else { 2512 be = NR_BessJ1(arg2)/arg2; 2513 } 2508 2514 2509 2515 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2510 yyy += Pi*rad*rad*len*2.0* NR_BessJ1(arg2)/arg2;2516 yyy += Pi*rad*rad*len*2.0*be; 2511 2517 } else { 2512 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0* NR_BessJ1(arg2)/arg2;2518 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2513 2519 } 2514 2520 yyy *= yyy; … … 2537 2543 double val,arg1,arg2; 2538 2544 double scale,bkg,sldc,slds; 2539 double len,rad,hDist,endRad ;2545 double len,rad,hDist,endRad,be; 2540 2546 scale = w[0]; 2541 2547 rad = w[1]; … … 2551 2557 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2552 2558 2553 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 2559 if(arg2 == 0) { 2560 be = 0.5; 2561 } else { 2562 be = NR_BessJ1(arg2)/arg2; 2563 } 2564 val = cos(arg1)*(1.0tt*tt)*be; 2554 2565 2555 2566 return(val); … … 2573 2584 double summ,zi,yyy,answer; //running tally of integration 2574 2585 double summj,vaj,vbj,zij; //for the inner integration 2575 double CLens_tmp[7],arg1,arg2,inner,hh ;2586 double CLens_tmp[7],arg1,arg2,inner,hh,be; 2576 2587 2577 2588 … … 2626 2637 yyy = inner; 2627 2638 2639 if(arg2 == 0) { 2640 be = 0.5; 2641 } else { 2642 be = NR_BessJ1(arg2)/arg2; 2643 } 2644 2628 2645 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2629 yyy += Pi*rad*rad*len*2.0* NR_BessJ1(arg2)/arg2;2646 yyy += Pi*rad*rad*len*2.0*be; 2630 2647 } else { 2631 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0* NR_BessJ1(arg2)/arg2;2648 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2632 2649 } 2633 2650 yyy *= yyy; … … 2667 2684 double summ,zi,yyy,answer; //running tally of integration 2668 2685 double summj,vaj,vbj,zij; //for the inner integration 2669 double arg1,arg2,inner,hh ;2686 double arg1,arg2,inner,hh,be; 2670 2687 2671 2688 … … 2710 2727 yyy = inner; 2711 2728 2729 if(arg2 == 0) { 2730 be = 0.5; 2731 } else { 2732 be = NR_BessJ1(arg2)/arg2; 2733 } 2734 2712 2735 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2713 yyy += Pi*rad*rad*len*2.0* NR_BessJ1(arg2)/arg2;2736 yyy += Pi*rad*rad*len*2.0*be; 2714 2737 } else { 2715 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0* NR_BessJ1(arg2)/arg2;2738 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2716 2739 } 2740 2741 2742 2717 2743 yyy *= yyy; 2718 2744 yyy *= sin(zi); // = A(q)^2*sin(theta) … … 2742 2768 double val,arg1,arg2; 2743 2769 double scale,bkg,sldc,slds; 2744 double len,rad,hDist,endRad ;2770 double len,rad,hDist,endRad,be; 2745 2771 scale = w[0]; 2746 2772 rad = w[1]; … … 2756 2782 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2757 2783 2758 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 2784 if(arg2 == 0) { 2785 be = 0.5; 2786 } else { 2787 be = NR_BessJ1(arg2)/arg2; 2788 } 2789 2790 val = cos(arg1)*(1.0tt*tt)*be; 2759 2791 2760 2792 return(val); … … 2778 2810 double summ,zi,yyy,answer; //running tally of integration 2779 2811 double summj,vaj,vbj,zij; //for the inner integration 2780 double Dumb_tmp[7],arg1,arg2,inner ;2812 double Dumb_tmp[7],arg1,arg2,inner,be; 2781 2813 2782 2814 … … 2831 2863 yyy = inner; 2832 2864 2865 if(arg2 == 0) { 2866 be = 0.5; 2867 } else { 2868 be = NR_BessJ1(arg2)/arg2; 2869 } 2870 2833 2871 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2834 yyy += Pi*rad*rad*len*2.0* NR_BessJ1(arg2)/arg2;2872 yyy += Pi*rad*rad*len*2.0*be; 2835 2873 } else { 2836 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0* NR_BessJ1(arg2)/arg2;2874 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2837 2875 } 2838 2876 yyy *= yyy; … … 2872 2910 double summ,zi,yyy,answer; //running tally of integration 2873 2911 double summj,vaj,vbj,zij; //for the inner integration 2874 double arg1,arg2,inner ;2912 double arg1,arg2,inner,be; 2875 2913 2876 2914 … … 2915 2953 yyy = inner; 2916 2954 2955 if(arg2 == 0) { 2956 be = 0.5; 2957 } else { 2958 be = NR_BessJ1(arg2)/arg2; 2959 } 2960 2917 2961 if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2918 yyy += Pi*rad*rad*len*2.0* NR_BessJ1(arg2)/arg2;2962 yyy += Pi*rad*rad*len*2.0*be; 2919 2963 } else { 2920 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0* NR_BessJ1(arg2)/arg2;2964 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 2921 2965 } 2922 2966 yyy *= yyy; … … 2948 2992 double val,arg1,arg2; 2949 2993 double scale,bkg,sldc,slds; 2950 double len,rad,hDist,endRad ;2994 double len,rad,hDist,endRad,be; 2951 2995 scale = w[0]; 2952 2996 rad = w[1]; … … 2962 3006 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2963 3007 2964 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 3008 if(arg2 == 0) { 3009 be = 0.5; 3010 } else { 3011 be = NR_BessJ1(arg2)/arg2; 3012 } 3013 val = cos(arg1)*(1.0tt*tt)*be; 2965 3014 2966 3015 return(val); … … 3059 3108 double sinarg1,sinarg2; 3060 3109 double t1,t2,t3; 3061 double retval ;3110 double retval,si1,si2,be1,be2; 3062 3111 3063 3112 double Pi = 4.0*atan(1.0); … … 3074 3123 sinarg2 = qq*(length+facthick)*cos(dum); 3075 3124 3076 t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 3077 t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 3078 t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; 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; 3079 3148 3080 3149 retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 
sans/XOP_Dev/SANSAnalysis/lib/libSphere.c
r632 r634 1287 1287 qr=x*rcore; 1288 1288 contr = rhocorerhoshel; 1289 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1289 if(qr == 0){ 1290 bes = 1.0; 1291 }else{ 1292 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1293 } 1290 1294 vol = 4.0*pi/3.0*rcore*rcore*rcore; 1291 1295 f = vol*bes*contr; … … 1293 1297 qr=x*(rcore+thick); 1294 1298 contr = rhoshelrhosolv; 1295 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1299 if(qr == 0){ 1300 bes = 1.0; 1301 }else{ 1302 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1303 } 1296 1304 vol = 4.0*pi/3.0*pow((rcore+thick),3); 1297 1305 f += vol*bes*contr; … … 1342 1350 qr=x*rcore; 1343 1351 contr = rhocorerhoshel1; 1344 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1352 if(qr == 0){ 1353 bes = 1.0; 1354 }else{ 1355 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1356 } 1345 1357 vol = 4.0*pi/3.0*rcore*rcore*rcore; 1346 1358 f = vol*bes*contr; … … 1348 1360 qr=x*(rcore+thick1); 1349 1361 contr = rhoshel1rhoshel2; 1350 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1362 if(qr == 0){ 1363 bes = 1.0; 1364 }else{ 1365 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1366 } 1351 1367 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1352 1368 f += vol*bes*contr; … … 1354 1370 qr=x*(rcore+thick1+thick2); 1355 1371 contr = rhoshel2rhosolv; 1356 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1372 if(qr == 0){ 1373 bes = 1.0; 1374 }else{ 1375 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1376 } 1357 1377 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1358 1378 f += vol*bes*contr; … … 1409 1429 qr=x*rcore; 1410 1430 contr = rhocorerhoshel1; 1411 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1431 if(qr == 0){ 1432 bes = 1.0; 1433 }else{ 1434 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1435 } 1412 1436 vol = 4.0*pi/3.0*rcore*rcore*rcore; 1413 1437 f = vol*bes*contr; … … 1415 1439 qr=x*(rcore+thick1); 1416 1440 contr = rhoshel1rhoshel2; 1417 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1441 if(qr == 0){ 1442 bes = 1.0; 1443 }else{ 1444 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1445 } 1418 1446 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1419 1447 f += vol*bes*contr; … … 1421 1449 qr=x*(rcore+thick1+thick2); 1422 1450 contr = rhoshel2rhoshel3; 1423 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1451 if(qr == 0){ 1452 bes = 1.0; 1453 }else{ 1454 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1455 } 1424 1456 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1425 1457 f += vol*bes*contr; … … 1427 1459 qr=x*(rcore+thick1+thick2+thick3); 1428 1460 contr = rhoshel3rhosolv; 1429 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1461 if(qr == 0){ 1462 bes = 1.0; 1463 }else{ 1464 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1465 } 1430 1466 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 1431 1467 f += vol*bes*contr; … … 1485 1521 qr=x*rcore; 1486 1522 contr = rhocorerhoshel1; 1487 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1523 if(qr == 0){ 1524 bes = 1.0; 1525 }else{ 1526 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1527 } 1488 1528 vol = 4.0*pi/3.0*rcore*rcore*rcore; 1489 1529 f = vol*bes*contr; … … 1491 1531 qr=x*(rcore+thick1); 1492 1532 contr = rhoshel1rhoshel2; 1493 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1533 if(qr == 0){ 1534 bes = 1.0; 1535 }else{ 1536 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1537 } 1494 1538 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1495 1539 f += vol*bes*contr; … … 1497 1541 qr=x*(rcore+thick1+thick2); 1498 1542 contr = rhoshel2rhoshel3; 1499 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1543 if(qr == 0){ 1544 bes = 1.0; 1545 }else{ 1546 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1547 } 1500 1548 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1501 1549 f += vol*bes*contr; … … 1503 1551 qr=x*(rcore+thick1+thick2+thick3); 1504 1552 contr = rhoshel3rhoshel4; 1505 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1553 if(qr == 0){ 1554 bes = 1.0; 1555 }else{ 1556 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1557 } 1506 1558 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 1507 1559 f += vol*bes*contr; … … 1509 1561 qr=x*(rcore+thick1+thick2+thick3+thick4); 1510 1562 contr = rhoshel4rhosolv; 1511 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1563 if(qr == 0){ 1564 bes = 1.0; 1565 }else{ 1566 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1567 } 1512 1568 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 1513 1569 f += vol*bes*contr; 
sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c
r632 r634 154 154 ans = G1*exp(x*x*Rg1*Rg1/3.0); 155 155 ans += B1*pow((erf1*erf1*erf1/x),Pow1); 156 157 if(x == 0) { 158 ans = G1; 159 } 156 160 157 161 ans *= scale; … … 190 194 ans += G2*exp(x*x*Rg2*Rg2/3.0); 191 195 ans += B2*pow((erf2*erf2*erf2/x),Pow2); 196 197 if(x == 0) { 198 ans = G1+G2; 199 } 192 200 193 201 ans *= scale; … … 232 240 ans += G2*exp(x*x*Rg2*Rg2/3.0) + B2*exp(x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2); 233 241 ans += G3*exp(x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3); 242 243 if(x == 0) { 244 ans = G1+G2+G3; 245 } 234 246 235 247 ans *= scale; … … 280 292 ans += G3*exp(x*x*Rg3*Rg3/3.0) + B3*exp(x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3); 281 293 ans += G4*exp(x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4); 294 295 if(x == 0) { 296 ans = G1+G2+G3+G4; 297 } 282 298 283 299 ans *= scale;
Note: See TracChangeset
for help on using the changeset viewer.