Changeset 453 for sans/XOP_Dev/SANSAnalysis/lib/libSphere.c
 Timestamp:
 Nov 18, 2008 3:26:32 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/SANSAnalysis/lib/libSphere.c
r235 r453 1224 1224 } 1225 1225 1226 double 1227 OneShell(double dp[], double q) 1228 { 1229 // variables are: 1230 //[0] scale factor 1231 //[1] radius of core [] 1232 //[2] SLD of the core [2] 1233 //[3] thickness of the shell [] 1234 //[4] SLD of the shell 1235 //[5] SLD of the solvent 1236 //[6] background [cm1] 1237 1238 double x,pi; 1239 double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg; //my local names 1240 double bes,f,vol,qr,contr,f2; 1241 1242 pi = 4.0*atan(1.0); 1243 x=q; 1244 1245 scale = dp[0]; 1246 rcore = dp[1]; 1247 rhocore = dp[2]; 1248 thick = dp[3]; 1249 rhoshel = dp[4]; 1250 rhosolv = dp[5]; 1251 bkg = dp[6]; 1252 1253 // core first, then add in shell 1254 qr=x*rcore; 1255 contr = rhocorerhoshel; 1256 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1257 vol = 4.0*pi/3.0*rcore*rcore*rcore; 1258 f = vol*bes*contr; 1259 //now the shell 1260 qr=x*(rcore+thick); 1261 contr = rhoshelrhosolv; 1262 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1263 vol = 4.0*pi/3.0*pow((rcore+thick),3); 1264 f += vol*bes*contr; 1265 1266 // normalize to particle volume and rescale from [1] to [cm1] 1267 f2 = f*f/vol*1.0e8; 1268 1269 //scale if desired 1270 f2 *= scale; 1271 // then add in the background 1272 f2 += bkg; 1273 1274 return(f2); 1275 } 1276 1277 double 1278 TwoShell(double dp[], double q) 1279 { 1280 // variables are: 1281 //[0] scale factor 1282 //[1] radius of core [] 1283 //[2] SLD of the core [2] 1284 //[3] thickness of shell 1 [] 1285 //[4] SLD of shell 1 1286 //[5] thickness of shell 2 [] 1287 //[6] SLD of shell 2 1288 //[7] SLD of the solvent 1289 //[8] background [cm1] 1290 1291 double x,pi; 1292 double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names 1293 double bes,f,vol,qr,contr,f2; 1294 double rhoshel2,thick2; 1295 1296 pi = 4.0*atan(1.0); 1297 x=q; 1298 1299 scale = dp[0]; 1300 rcore = dp[1]; 1301 rhocore = dp[2]; 1302 thick1 = dp[3]; 1303 rhoshel1 = dp[4]; 1304 thick2 = dp[5]; 1305 rhoshel2 = dp[6]; 1306 rhosolv = dp[7]; 1307 bkg = dp[8]; 1308 // core first, then add in shells 1309 qr=x*rcore; 1310 contr = rhocorerhoshel1; 1311 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1312 vol = 4.0*pi/3*rcore*rcore*rcore; 1313 f = vol*bes*contr; 1314 //now the shell (1) 1315 qr=x*(rcore+thick1); 1316 contr = rhoshel1rhoshel2; 1317 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1318 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1319 f += vol*bes*contr; 1320 //now the shell (2) 1321 qr=x*(rcore+thick1+thick2); 1322 contr = rhoshel2rhosolv; 1323 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1324 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1325 f += vol*bes*contr; 1326 1327 1328 // normalize to particle volume and rescale from [1] to [cm1] 1329 f2 = f*f/vol*1.0e8; 1330 1331 //scale if desired 1332 f2 *= scale; 1333 // then add in the background 1334 f2 += bkg; 1335 1336 return(f2); 1337 } 1338 1339 double 1340 ThreeShell(double dp[], double q) 1341 { 1342 // variables are: 1343 //[0] scale factor 1344 //[1] radius of core [] 1345 //[2] SLD of the core [2] 1346 //[3] thickness of shell 1 [] 1347 //[4] SLD of shell 1 1348 //[5] thickness of shell 2 [] 1349 //[6] SLD of shell 2 1350 //[7] thickness of shell 3 1351 //[8] SLD of shell 3 1352 //[9] SLD of solvent 1353 //[10] background [cm1] 1354 1355 double x,pi; 1356 double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names 1357 double bes,f,vol,qr,contr,f2; 1358 double rhoshel2,thick2,rhoshel3,thick3; 1359 1360 pi = 4.0*atan(1.0); 1361 x=q; 1362 1363 scale = dp[0]; 1364 rcore = dp[1]; 1365 rhocore = dp[2]; 1366 thick1 = dp[3]; 1367 rhoshel1 = dp[4]; 1368 thick2 = dp[5]; 1369 rhoshel2 = dp[6]; 1370 thick3 = dp[7]; 1371 rhoshel3 = dp[8]; 1372 rhosolv = dp[9]; 1373 bkg = dp[10]; 1374 1375 // core first, then add in shells 1376 qr=x*rcore; 1377 contr = rhocorerhoshel1; 1378 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1379 vol = 4.0*pi/3*rcore*rcore*rcore; 1380 f = vol*bes*contr; 1381 //now the shell (1) 1382 qr=x*(rcore+thick1); 1383 contr = rhoshel1rhoshel2; 1384 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1385 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1386 f += vol*bes*contr; 1387 //now the shell (2) 1388 qr=x*(rcore+thick1+thick2); 1389 contr = rhoshel2rhoshel3; 1390 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1391 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1392 f += vol*bes*contr; 1393 //now the shell (3) 1394 qr=x*(rcore+thick1+thick2+thick3); 1395 contr = rhoshel3rhosolv; 1396 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1397 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 1398 f += vol*bes*contr; 1399 1400 // normalize to particle volume and rescale from [1] to [cm1] 1401 f2 = f*f/vol*1.0e8; 1402 1403 //scale if desired 1404 f2 *= scale; 1405 // then add in the background 1406 f2 += bkg; 1407 1408 return(f2); 1409 } 1410 1411 double 1412 FourShell(double dp[], double q) 1413 { 1414 // variables are: 1415 //[0] scale factor 1416 //[1] radius of core [] 1417 //[2] SLD of the core [2] 1418 //[3] thickness of shell 1 [] 1419 //[4] SLD of shell 1 1420 //[5] thickness of shell 2 [] 1421 //[6] SLD of shell 2 1422 //[7] thickness of shell 3 1423 //[8] SLD of shell 3 1424 //[9] thickness of shell 3 1425 //[10] SLD of shell 3 1426 //[11] SLD of solvent 1427 //[12] background [cm1] 1428 1429 double x,pi; 1430 double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names 1431 double bes,f,vol,qr,contr,f2; 1432 double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; 1433 1434 pi = 4.0*atan(1.0); 1435 x=q; 1436 1437 scale = dp[0]; 1438 rcore = dp[1]; 1439 rhocore = dp[2]; 1440 thick1 = dp[3]; 1441 rhoshel1 = dp[4]; 1442 thick2 = dp[5]; 1443 rhoshel2 = dp[6]; 1444 thick3 = dp[7]; 1445 rhoshel3 = dp[8]; 1446 thick4 = dp[9]; 1447 rhoshel4 = dp[10]; 1448 rhosolv = dp[11]; 1449 bkg = dp[12]; 1450 1451 // core first, then add in shells 1452 qr=x*rcore; 1453 contr = rhocorerhoshel1; 1454 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1455 vol = 4.0*pi/3*rcore*rcore*rcore; 1456 f = vol*bes*contr; 1457 //now the shell (1) 1458 qr=x*(rcore+thick1); 1459 contr = rhoshel1rhoshel2; 1460 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1461 vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 1462 f += vol*bes*contr; 1463 //now the shell (2) 1464 qr=x*(rcore+thick1+thick2); 1465 contr = rhoshel2rhoshel3; 1466 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1467 vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 1468 f += vol*bes*contr; 1469 //now the shell (3) 1470 qr=x*(rcore+thick1+thick2+thick3); 1471 contr = rhoshel3rhoshel4; 1472 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1473 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 1474 f += vol*bes*contr; 1475 //now the shell (4) 1476 qr=x*(rcore+thick1+thick2+thick3+thick4); 1477 contr = rhoshel4rhosolv; 1478 bes = 3.0*(sin(qr)qr*cos(qr))/(qr*qr*qr); 1479 vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 1480 f += vol*bes*contr; 1481 1482 1483 // normalize to particle volume and rescale from [1] to [cm1] 1484 f2 = f*f/vol*1.0e8; 1485 1486 //scale if desired 1487 f2 *= scale; 1488 // then add in the background 1489 f2 += bkg; 1490 1491 return(f2); 1492 } 1493 1494 double 1495 PolyOneShell(double dp[], double x) 1496 { 1497 double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz; //my local names 1498 double va,vb,summ,yyy,zi; 1499 double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi; 1500 int nord=76,ii; 1501 1502 pi = 4.0*atan(1.0); 1503 1504 scale = dp[0]; 1505 rcore = dp[1]; 1506 pd = dp[2]; 1507 rhocore = dp[3]; 1508 thick = dp[4]; 1509 rhoshel = dp[5]; 1510 rhosolv = dp[6]; 1511 bkg = dp[7]; 1512 1513 zz = (1.0/pd)*(1.0/pd)1.0; //polydispersity of the core only 1514 1515 range = 8.0; //std deviations for the integration 1516 va = rcore*(1.0range*pd); 1517 if (va<0) { 1518 va=0; //otherwise numerical error when pd >= 0.3, making a<0 1519 } 1520 if (pd>0.3) { 1521 range = range + (pd0.3)*18.0; //stretch upper range to account for skewed tail 1522 } 1523 vb = rcore*(1.0+range*pd); // is this far enough past avg radius? 1524 1525 //temp set scale=1 and bkg=0 for quadrature calc 1526 temp_1sf[0] = 1.0; 1527 temp_1sf[1] = dp[1]; //the core radius will be changed in the loop 1528 temp_1sf[2] = dp[3]; 1529 temp_1sf[3] = dp[4]; 1530 temp_1sf[4] = dp[5]; 1531 temp_1sf[5] = dp[6]; 1532 temp_1sf[6] = 0.0; 1533 1534 summ = 0.0; // initialize integral 1535 for(ii=0;ii<nord;ii+=1) { 1536 // calculate Gauss points on integration interval (rvalue for evaluation) 1537 zi = ( Gauss76Z[ii]*(vbva) + vb + va )/2.0; 1538 temp_1sf[1] = zi; 1539 yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x); 1540 //unnormalize by volume 1541 yyy *= 4.0*pi/3.0*pow((zi+thick),3); 1542 summ += yyy; //add to the running total of the quadrature 1543 } 1544 // calculate value of integral to return 1545 answer = (vbva)/2.0*summ; 1546 1547 //renormalize by the average volume 1548 zp1 = zz + 1.0; 1549 zp2 = zz + 2.0; 1550 zp3 = zz + 3.0; 1551 vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3); 1552 answer /= vpoly; 1553 //scale 1554 answer *= scale; 1555 // add in the background 1556 answer += bkg; 1557 1558 return(answer); 1559 } 1560 1561 double 1562 PolyTwoShell(double dp[], double x) 1563 { 1564 double scale,rcore,rhocore,rhosolv,bkg,pd,zz; //my local names 1565 double va,vb,summ,yyy,zi; 1566 double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi; 1567 int nord=76,ii; 1568 double thick1,thick2; 1569 double rhoshel1,rhoshel2; 1570 1571 scale = dp[0]; 1572 rcore = dp[1]; 1573 pd = dp[2]; 1574 rhocore = dp[3]; 1575 thick1 = dp[4]; 1576 rhoshel1 = dp[5]; 1577 thick2 = dp[6]; 1578 rhoshel2 = dp[7]; 1579 rhosolv = dp[8]; 1580 bkg = dp[9]; 1581 1582 pi = 4.0*atan(1.0); 1583 1584 zz = (1.0/pd)*(1.0/pd)1.0; //polydispersity of the core only 1585 1586 range = 8.0; //std deviations for the integration 1587 va = rcore*(1.0range*pd); 1588 if (va<0) { 1589 va=0; //otherwise numerical error when pd >= 0.3, making a<0 1590 } 1591 if (pd>0.3) { 1592 range = range + (pd0.3)*18.0; //stretch upper range to account for skewed tail 1593 } 1594 vb = rcore*(1.0+range*pd); // is this far enough past avg radius? 1595 1596 //temp set scale=1 and bkg=0 for quadrature calc 1597 temp_2sf[0] = 1.0; 1598 temp_2sf[1] = dp[1]; //the core radius will be changed in the loop 1599 temp_2sf[2] = dp[3]; 1600 temp_2sf[3] = dp[4]; 1601 temp_2sf[4] = dp[5]; 1602 temp_2sf[5] = dp[6]; 1603 temp_2sf[6] = dp[7]; 1604 temp_2sf[7] = dp[8]; 1605 temp_2sf[8] = 0.0; 1606 1607 summ = 0.0; // initialize integral 1608 for(ii=0;ii<nord;ii+=1) { 1609 // calculate Gauss points on integration interval (rvalue for evaluation) 1610 zi = ( Gauss76Z[ii]*(vbva) + vb + va )/2.0; 1611 temp_2sf[1] = zi; 1612 yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x); 1613 //unnormalize by volume 1614 yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3); 1615 summ += yyy; //add to the running total of the quadrature 1616 } 1617 // calculate value of integral to return 1618 answer = (vbva)/2.0*summ; 1619 1620 //renormalize by the average volume 1621 zp1 = zz + 1.0; 1622 zp2 = zz + 2.0; 1623 zp3 = zz + 3.0; 1624 vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3); 1625 answer /= vpoly; 1626 //scale 1627 answer *= scale; 1628 // add in the background 1629 answer += bkg; 1630 1631 return(answer); 1632 } 1633 1634 double 1635 PolyThreeShell(double dp[], double x) 1636 { 1637 double scale,rcore,rhocore,rhosolv,bkg,pd,zz; //my local names 1638 double va,vb,summ,yyy,zi; 1639 double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi; 1640 int nord=76,ii; 1641 double thick1,thick2,thick3; 1642 double rhoshel1,rhoshel2,rhoshel3; 1643 1644 scale = dp[0]; 1645 rcore = dp[1]; 1646 pd = dp[2]; 1647 rhocore = dp[3]; 1648 thick1 = dp[4]; 1649 rhoshel1 = dp[5]; 1650 thick2 = dp[6]; 1651 rhoshel2 = dp[7]; 1652 thick3 = dp[8]; 1653 rhoshel3 = dp[9]; 1654 rhosolv = dp[10]; 1655 bkg = dp[11]; 1656 1657 pi = 4.0*atan(1.0); 1658 1659 zz = (1.0/pd)*(1.0/pd)1.0; //polydispersity of the core only 1660 1661 range = 8.0; //std deviations for the integration 1662 va = rcore*(1.0range*pd); 1663 if (va<0) { 1664 va=0; //otherwise numerical error when pd >= 0.3, making a<0 1665 } 1666 if (pd>0.3) { 1667 range = range + (pd0.3)*18.0; //stretch upper range to account for skewed tail 1668 } 1669 vb = rcore*(1.0+range*pd); // is this far enough past avg radius? 1670 1671 //temp set scale=1 and bkg=0 for quadrature calc 1672 temp_3sf[0] = 1.0; 1673 temp_3sf[1] = dp[1]; //the core radius will be changed in the loop 1674 temp_3sf[2] = dp[3]; 1675 temp_3sf[3] = dp[4]; 1676 temp_3sf[4] = dp[5]; 1677 temp_3sf[5] = dp[6]; 1678 temp_3sf[6] = dp[7]; 1679 temp_3sf[7] = dp[8]; 1680 temp_3sf[8] = dp[9]; 1681 temp_3sf[9] = dp[10]; 1682 temp_3sf[10] = 0.0; 1683 1684 summ = 0.0; // initialize integral 1685 for(ii=0;ii<nord;ii+=1) { 1686 // calculate Gauss points on integration interval (rvalue for evaluation) 1687 zi = ( Gauss76Z[ii]*(vbva) + vb + va )/2.0; 1688 temp_3sf[1] = zi; 1689 yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x); 1690 //unnormalize by volume 1691 yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3); 1692 summ += yyy; //add to the running total of the quadrature 1693 } 1694 // calculate value of integral to return 1695 answer = (vbva)/2.0*summ; 1696 1697 //renormalize by the average volume 1698 zp1 = zz + 1.0; 1699 zp2 = zz + 2.0; 1700 zp3 = zz + 3.0; 1701 vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3); 1702 answer /= vpoly; 1703 //scale 1704 answer *= scale; 1705 // add in the background 1706 answer += bkg; 1707 1708 return(answer); 1709 } 1710 1711 double 1712 PolyFourShell(double dp[], double x) 1713 { 1714 double scale,rcore,rhocore,rhosolv,bkg,pd,zz; //my local names 1715 double va,vb,summ,yyy,zi; 1716 double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi; 1717 int nord=76,ii; 1718 double thick1,thick2,thick3,thick4; 1719 double rhoshel1,rhoshel2,rhoshel3,rhoshel4; 1720 1721 scale = dp[0]; 1722 rcore = dp[1]; 1723 pd = dp[2]; 1724 rhocore = dp[3]; 1725 thick1 = dp[4]; 1726 rhoshel1 = dp[5]; 1727 thick2 = dp[6]; 1728 rhoshel2 = dp[7]; 1729 thick3 = dp[8]; 1730 rhoshel3 = dp[9]; 1731 thick4 = dp[10]; 1732 rhoshel4 = dp[11]; 1733 rhosolv = dp[12]; 1734 bkg = dp[13]; 1735 1736 pi = 4.0*atan(1.0); 1737 1738 zz = (1.0/pd)*(1.0/pd)1.0; //polydispersity of the core only 1739 1740 range = 8.0; //std deviations for the integration 1741 va = rcore*(1.0range*pd); 1742 if (va<0) { 1743 va=0; //otherwise numerical error when pd >= 0.3, making a<0 1744 } 1745 if (pd>0.3) { 1746 range = range + (pd0.3)*18.0; //stretch upper range to account for skewed tail 1747 } 1748 vb = rcore*(1.0+range*pd); // is this far enough past avg radius? 1749 1750 //temp set scale=1 and bkg=0 for quadrature calc 1751 temp_4sf[0] = 1.0; 1752 temp_4sf[1] = dp[1]; //the core radius will be changed in the loop 1753 temp_4sf[2] = dp[3]; 1754 temp_4sf[3] = dp[4]; 1755 temp_4sf[4] = dp[5]; 1756 temp_4sf[5] = dp[6]; 1757 temp_4sf[6] = dp[7]; 1758 temp_4sf[7] = dp[8]; 1759 temp_4sf[8] = dp[9]; 1760 temp_4sf[9] = dp[10]; 1761 temp_4sf[10] = dp[11]; 1762 temp_4sf[11] = dp[12]; 1763 temp_4sf[12] = 0.0; 1764 1765 summ = 0.0; // initialize integral 1766 for(ii=0;ii<nord;ii+=1) { 1767 // calculate Gauss points on integration interval (rvalue for evaluation) 1768 zi = ( Gauss76Z[ii]*(vbva) + vb + va )/2.0; 1769 temp_4sf[1] = zi; 1770 yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x); 1771 //unnormalize by volume 1772 yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3); 1773 summ += yyy; //add to the running total of the quadrature 1774 } 1775 // calculate value of integral to return 1776 answer = (vbva)/2.0*summ; 1777 1778 //renormalize by the average volume 1779 zp1 = zz + 1.0; 1780 zp2 = zz + 2.0; 1781 zp3 = zz + 3.0; 1782 vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3); 1783 answer /= vpoly; 1784 //scale 1785 answer *= scale; 1786 // add in the background 1787 answer += bkg; 1788 1789 return(answer); 1790 } 1791 1792 1793 /* BCC_ParaCrystal : calculates the form factor of a Triaxial Ellipsoid at the given xvalue p>x 1794 1795 Uses 150 pt Gaussian quadrature for both integrals 1796 1797 */ 1798 double 1799 BCC_ParaCrystal(double w[], double x) 1800 { 1801 int i,j; 1802 double Pi; 1803 double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv; //local variables of coefficient wave 1804 int nordi=150; //order of integration 1805 int nordj=150; 1806 double va,vb; //upper and lower integration limits 1807 double summ,zi,yyy,answer; //running tally of integration 1808 double summj,vaj,vbj,zij; //for the inner integration 1809 1810 Pi = 4.0*atan(1.0); 1811 va = 0.0; 1812 vb = 2.0*Pi; //orintational average, outer integral 1813 vaj = 0.0; 1814 vbj = Pi; //endpoints of inner integral 1815 1816 summ = 0.0; //initialize intergral 1817 1818 scale = w[0]; 1819 Dnn = w[1]; //Nearest neighbor distance A 1820 gg = w[2]; //Paracrystal distortion factor 1821 Rad = w[3]; //Sphere radius 1822 sld = w[4]; 1823 sldSolv = w[5]; 1824 background = w[6]; 1825 1826 contrast = sld  sldSolv; 1827 1828 //Volume fraction calculated from lattice symmetry and sphere radius 1829 latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3); 1830 1831 for(i=0;i<nordi;i++) { 1832 //setup inner integral over the ellipsoidal crosssection 1833 summj=0.0; 1834 zi = ( Gauss150Z[i]*(vbva) + va + vb )/2.0; //the outer dummy is phi 1835 for(j=0;j<nordj;j++) { 1836 //20 gauss points for the inner integral 1837 zij = ( Gauss150Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the inner dummy is theta 1838 yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij); 1839 summj += yyy; 1840 } 1841 //now calculate the value of the inner integral 1842 answer = (vbjvaj)/2.0*summj; 1843 1844 //now calculate outer integral 1845 yyy = Gauss150Wt[i] * answer; 1846 summ += yyy; 1847 } //final scaling is done at the end of the function, after the NT_FP64 case 1848 1849 answer = (vbva)/2.0*summ; 1850 // Multiply by contrast^2 1851 answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 1852 // add in the background 1853 answer += background; 1854 1855 return answer; 1856 } 1857 1858 // xx is phi (outer) 1859 // yy is theta (inner) 1860 double 1861 BCC_Integrand(double w[], double qq, double xx, double yy) { 1862 1863 double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 1864 1865 Dnn = w[1]; //Nearest neighbor distance A 1866 gg = w[2]; //Paracrystal distortion factor 1867 aa = Dnn; 1868 Da = gg*aa; 1869 1870 Pi = 4.0*atan(1.0); 1871 temp1 = qq*qq*Da*Da; 1872 temp3 = qq*aa; 1873 1874 retVal = BCCeval(yy,xx,temp1,temp3); 1875 retVal /=4.0*Pi; 1876 1877 return(retVal); 1878 } 1879 1880 double 1881 BCCeval(double Theta, double Phi, double temp1, double temp3) { 1882 1883 double temp6,temp7,temp8,temp9,temp10; 1884 double result; 1885 1886 temp6 = sin(Theta); 1887 temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta); 1888 temp8 = 1.0*sin(Theta)*cos(Phi)sin(Theta)*sin(Phi)+cos(Theta); 1889 temp9 = 1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)cos(Theta); 1890 temp10 = exp((1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 1891 result = pow(1.0(temp10*temp10),3)*temp6/((1.02.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.02.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.02.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 1892 1893 return (result); 1894 } 1895 1896 double 1897 SphereForm_Paracrystal(double radius, double delrho, double x) { 1898 1899 double bes,f,vol,f2,pi; 1900 pi = 4.0*atan(1.0); 1901 // 1902 //handle q==0 separately 1903 if(x==0) { 1904 f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8; 1905 return(f); 1906 } 1907 1908 bes = 3.0*(sin(x*radius)x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius); 1909 vol = 4.0*pi/3.0*radius*radius*radius; 1910 f = vol*bes*delrho ; // [=] 1911 // normalize to single particle volume, convert to 1/cm 1912 f2 = f * f / vol * 1.0e8; // [=] 1/cm 1913 1914 return (f2); 1915 } 1916 1917 /* FCC_ParaCrystal : calculates the form factor of a Triaxial Ellipsoid at the given xvalue p>x 1918 1919 Uses 150 pt Gaussian quadrature for both integrals 1920 1921 */ 1922 double 1923 FCC_ParaCrystal(double w[], double x) 1924 { 1925 int i,j; 1926 double Pi; 1927 double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv; //local variables of coefficient wave 1928 int nordi=150; //order of integration 1929 int nordj=150; 1930 double va,vb; //upper and lower integration limits 1931 double summ,zi,yyy,answer; //running tally of integration 1932 double summj,vaj,vbj,zij; //for the inner integration 1933 1934 Pi = 4.0*atan(1.0); 1935 va = 0.0; 1936 vb = 2.0*Pi; //orintational average, outer integral 1937 vaj = 0.0; 1938 vbj = Pi; //endpoints of inner integral 1939 1940 summ = 0.0; //initialize intergral 1941 1942 scale = w[0]; 1943 Dnn = w[1]; //Nearest neighbor distance A 1944 gg = w[2]; //Paracrystal distortion factor 1945 Rad = w[3]; //Sphere radius 1946 sld = w[4]; 1947 sldSolv = w[5]; 1948 background = w[6]; 1949 1950 contrast = sld  sldSolv; 1951 //Volume fraction calculated from lattice symmetry and sphere radius 1952 latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3); 1953 1954 for(i=0;i<nordi;i++) { 1955 //setup inner integral over the ellipsoidal crosssection 1956 summj=0.0; 1957 zi = ( Gauss150Z[i]*(vbva) + va + vb )/2.0; //the outer dummy is phi 1958 for(j=0;j<nordj;j++) { 1959 //20 gauss points for the inner integral 1960 zij = ( Gauss150Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the inner dummy is theta 1961 yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij); 1962 summj += yyy; 1963 } 1964 //now calculate the value of the inner integral 1965 answer = (vbjvaj)/2.0*summj; 1966 1967 //now calculate outer integral 1968 yyy = Gauss150Wt[i] * answer; 1969 summ += yyy; 1970 } //final scaling is done at the end of the function, after the NT_FP64 case 1971 1972 answer = (vbva)/2.0*summ; 1973 // Multiply by contrast^2 1974 answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 1975 // add in the background 1976 answer += background; 1977 1978 return answer; 1979 } 1980 1981 1982 // xx is phi (outer) 1983 // yy is theta (inner) 1984 double 1985 FCC_Integrand(double w[], double qq, double xx, double yy) { 1986 1987 double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 1988 1989 Pi = 4.0*atan(1.0); 1990 Dnn = w[1]; //Nearest neighbor distance A 1991 gg = w[2]; //Paracrystal distortion factor 1992 aa = Dnn; 1993 Da = gg*aa; 1994 1995 temp1 = qq*qq*Da*Da; 1996 temp3 = qq*aa; 1997 1998 retVal = FCCeval(yy,xx,temp1,temp3); 1999 retVal /=4*Pi; 2000 2001 return(retVal); 2002 } 2003 2004 double 2005 FCCeval(double Theta, double Phi, double temp1, double temp3) { 2006 2007 double temp6,temp7,temp8,temp9,temp10; 2008 double result; 2009 2010 temp6 = sin(Theta); 2011 temp7 = sin(Theta)*sin(Phi)+cos(Theta); 2012 temp8 = 1.0*sin(Theta)*cos(Phi)+cos(Theta); 2013 temp9 = 1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi); 2014 temp10 = exp((1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 2015 result = pow((1.0(temp10*temp10)),3)*temp6/((1.02.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.02.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.02.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 2016 2017 return (result); 2018 } 2019 2020 2021 /* SC_ParaCrystal : calculates the form factor of a Triaxial Ellipsoid at the given xvalue p>x 2022 2023 Uses 150 pt Gaussian quadrature for both integrals 2024 2025 */ 2026 double 2027 SC_ParaCrystal(double w[], double x) 2028 { 2029 int i,j; 2030 double Pi; 2031 double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv; //local variables of coefficient wave 2032 int nordi=150; //order of integration 2033 int nordj=150; 2034 double va,vb; //upper and lower integration limits 2035 double summ,zi,yyy,answer; //running tally of integration 2036 double summj,vaj,vbj,zij; //for the inner integration 2037 2038 Pi = 4.0*atan(1.0); 2039 va = 0.0; 2040 vb = 2.0*Pi; //orintational average, outer integral 2041 vaj = 0.0; 2042 vbj = Pi; //endpoints of inner integral 2043 2044 summ = 0.0; //initialize intergral 2045 2046 scale = w[0]; 2047 Dnn = w[1]; //Nearest neighbor distance A 2048 gg = w[2]; //Paracrystal distortion factor 2049 Rad = w[3]; //Sphere radius 2050 sld = w[4]; 2051 sldSolv = w[5]; 2052 background = w[6]; 2053 2054 contrast = sld  sldSolv; 2055 //Volume fraction calculated from lattice symmetry and sphere radius 2056 latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3); 2057 2058 for(i=0;i<nordi;i++) { 2059 //setup inner integral over the ellipsoidal crosssection 2060 summj=0.0; 2061 zi = ( Gauss150Z[i]*(vbva) + va + vb )/2.0; //the outer dummy is phi 2062 for(j=0;j<nordj;j++) { 2063 //20 gauss points for the inner integral 2064 zij = ( Gauss150Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the inner dummy is theta 2065 yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij); 2066 summj += yyy; 2067 } 2068 //now calculate the value of the inner integral 2069 answer = (vbjvaj)/2.0*summj; 2070 2071 //now calculate outer integral 2072 yyy = Gauss150Wt[i] * answer; 2073 summ += yyy; 2074 } //final scaling is done at the end of the function, after the NT_FP64 case 2075 2076 answer = (vbva)/2.0*summ; 2077 // Multiply by contrast^2 2078 answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 2079 // add in the background 2080 answer += background; 2081 2082 return answer; 2083 } 2084 2085 // xx is phi (outer) 2086 // yy is theta (inner) 2087 double 2088 SC_Integrand(double w[], double qq, double xx, double yy) { 2089 2090 double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi; 2091 2092 Pi = 4.0*atan(1.0); 2093 Dnn = w[1]; //Nearest neighbor distance A 2094 gg = w[2]; //Paracrystal distortion factor 2095 aa = Dnn; 2096 Da = gg*aa; 2097 2098 temp1 = qq*qq*Da*Da; 2099 temp2 = pow( 1.0exp(1.0*temp1) ,3); 2100 temp3 = qq*aa; 2101 temp4 = 2.0*exp(0.5*temp1); 2102 temp5 = exp(1.0*temp1); 2103 2104 2105 retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5); 2106 retVal /= 4*Pi; 2107 2108 return(retVal); 2109 } 2110 2111 double 2112 SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure 2113 2114 double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation 2115 double result; 2116 2117 temp6 = sin(Theta); 2118 temp7 = 1.0*temp3*sin(Theta)*cos(Phi); 2119 temp8 = temp3*sin(Theta)*sin(Phi); 2120 temp9 = temp3*cos(Theta); 2121 result = temp6/((1.0temp4*cos((temp7))+temp5)*(1.0temp4*cos((temp8))+temp5)*(1.0temp4*cos((temp9))+temp5)); 2122 2123 return (result); 2124 } 2125
Note: See TracChangeset
for help on using the changeset viewer.