Ignore:
Timestamp:
Nov 18, 2008 3:26:32 PM (14 years ago)
Author:
srkline
Message:

Additions to the library of the 2008 model functions. Direct proting of the Igor code, duplicated by these XOPs

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

    r235 r453  
    12241224} 
    12251225 
     1226double 
     1227OneShell(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        [cm-1] 
     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 = rhocore-rhoshel; 
     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 = rhoshel-rhosolv; 
     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 [cm-1] 
     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 
     1277double 
     1278TwoShell(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        [cm-1] 
     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 = rhocore-rhoshel1; 
     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 = rhoshel1-rhoshel2; 
     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 = rhoshel2-rhosolv; 
     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 [cm-1] 
     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 
     1339double 
     1340ThreeShell(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       [cm-1] 
     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 = rhocore-rhoshel1; 
     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 = rhoshel1-rhoshel2; 
     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 = rhoshel2-rhoshel3; 
     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 = rhoshel3-rhosolv; 
     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 [cm-1] 
     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 
     1411double 
     1412FourShell(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       [cm-1] 
     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 = rhocore-rhoshel1; 
     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 = rhoshel1-rhoshel2; 
     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 = rhoshel2-rhoshel3; 
     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 = rhoshel3-rhoshel4; 
     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 = rhoshel4-rhosolv; 
     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 [cm-1] 
     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 
     1494double 
     1495PolyOneShell(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.0-range*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 + (pd-0.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 (r-value for evaluation) 
     1537                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1538                temp_1sf[1] = zi; 
     1539                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x); 
     1540                //un-normalize 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 = (vb-va)/2.0*summ; 
     1546         
     1547        //re-normalize 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 
     1561double 
     1562PolyTwoShell(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.0-range*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 + (pd-0.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 (r-value for evaluation) 
     1610                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1611                temp_2sf[1] = zi; 
     1612                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x); 
     1613                //un-normalize 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 = (vb-va)/2.0*summ; 
     1619         
     1620        //re-normalize 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 
     1634double 
     1635PolyThreeShell(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.0-range*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 + (pd-0.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 (r-value for evaluation) 
     1687                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1688                temp_3sf[1] = zi; 
     1689                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x); 
     1690                //un-normalize 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 = (vb-va)/2.0*summ; 
     1696         
     1697        //re-normalize 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 
     1711double 
     1712PolyFourShell(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.0-range*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 + (pd-0.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 (r-value for evaluation) 
     1768                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1769                temp_4sf[1] = zi; 
     1770                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x); 
     1771                //un-normalize 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 = (vb-va)/2.0*summ; 
     1777         
     1778        //re-normalize 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 x-value p->x 
     1794 
     1795Uses 150 pt Gaussian quadrature for both integrals 
     1796 
     1797*/ 
     1798double 
     1799BCC_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 cross-section 
     1833                summj=0.0; 
     1834                zi = ( Gauss150Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/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) 
     1860double 
     1861BCC_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 
     1880double 
     1881BCCeval(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.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 
     1892         
     1893        return (result); 
     1894} 
     1895 
     1896double 
     1897SphereForm_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 x-value p->x 
     1918 
     1919Uses 150 pt Gaussian quadrature for both integrals 
     1920 
     1921*/ 
     1922double 
     1923FCC_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 cross-section 
     1956                summj=0.0; 
     1957                zi = ( Gauss150Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/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) 
     1984double 
     1985FCC_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 
     2004double 
     2005FCCeval(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.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.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 x-value p->x 
     2022 
     2023Uses 150 pt Gaussian quadrature for both integrals 
     2024 
     2025*/ 
     2026double 
     2027SC_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 cross-section 
     2060                summj=0.0; 
     2061                zi = ( Gauss150Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/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) 
     2087double 
     2088SC_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.0-exp(-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 
     2111double 
     2112SCeval(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.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5)); 
     2122         
     2123        return (result); 
     2124} 
     2125 
Note: See TracChangeset for help on using the changeset viewer.