Changeset 34c2649 in sasview for sansmodels/src/sans/models/libigor


Ignore:
Timestamp:
Nov 21, 2011 4:52:25 PM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
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:
1810613
Parents:
afea8fe
Message:

Re #4 Fixed warnings

Location:
sansmodels/src/sans/models/libigor
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/libigor/libCylinder.c

    r3a084629 r34c2649  
    13371337         
    13381338        return answer; 
    1339 } 
    1340  
    1341 /////////functions for WRC implementation of flexible cylinders 
    1342 static double 
    1343 Sk_WR(double q, double L, double b) 
    1344 { 
    1345         // 
    1346         double p1,p2,p1short,p2short,q0,qconnect; 
    1347         double C,epsilon,ans,q0short,Sexvmodify,pi; 
    1348          
    1349         pi = 4.0*atan(1.0); 
    1350          
    1351         p1 = 4.12; 
    1352         p2 = 4.42; 
    1353         p1short = 5.36; 
    1354         p2short = 5.62; 
    1355         q0 = 3.1; 
    1356         qconnect = q0/b; 
    1357         //       
    1358         q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
    1359          
    1360         // 
    1361         if(L/b > 10.0) { 
    1362                 C = 3.06/pow((L/b),0.44); 
    1363                 epsilon = 0.176; 
    1364         } else { 
    1365                 C = 1.0; 
    1366                 epsilon = 0.170; 
    1367         } 
    1368         // 
    1369          
    1370         if( L > 4*b ) { // Longer Chains 
    1371                 if (q*b <= 3.1) {               //Modified by Yun on Oct. 15, 
    1372                         Sexvmodify = Sexvnew(q, L, b); 
    1373                         ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
    1374                 } else { //q(i)*b > 3.1 
    1375                         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); 
    1376                 }  
    1377         } else { //L <= 4*b Shorter Chains 
    1378                 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
    1379                         if (q*b<=0.01) { 
    1380                                 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
    1381                         } else { 
    1382                                 ans = Sdebye1(q,L,b); 
    1383                         } 
    1384                 } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
    1385                         ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); 
    1386                 } 
    1387         } 
    1388          
    1389         return(ans); 
    1390         //return(a2long(q, L, b, p1, p2, q0)); 
    1391 } 
    1392  
    1393 //WR named this w (too generic) 
    1394 static double 
    1395 w_WR(double x) 
    1396 { 
    1397     double yy; 
    1398     yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 
    1399          
    1400     return (yy); 
    1401 } 
    1402  
    1403 // 
    1404 static double 
    1405 u1(double q, double L, double b) 
    1406 { 
    1407     double yy; 
    1408          
    1409     yy = Rgsquareshort(q,L,b)*q*q; 
    1410      
    1411     return (yy); 
    1412 } 
    1413  
    1414 // was named u 
    1415 static double 
    1416 u_WR(double q, double L, double b) 
    1417 { 
    1418     double yy; 
    1419     yy = Rgsquare(q,L,b)*q*q; 
    1420     return (yy); 
    1421 } 
    1422  
    1423  
    1424  
    1425 // 
    1426 static double 
    1427 Rgsquarezero(double q, double L, double b) 
    1428 {     
    1429     double yy; 
    1430     yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); 
    1431      
    1432     return (yy); 
    1433 } 
    1434  
    1435 // 
    1436 static double 
    1437 Rgsquareshort(double q, double L, double b) 
    1438 {     
    1439     double yy; 
    1440     yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
    1441      
    1442     return (yy); 
    1443 } 
    1444  
    1445 // 
    1446 static double 
    1447 Rgsquare(double q, double L, double b) 
    1448 { 
    1449     double yy; 
    1450     yy = AlphaSquare(L/b)*L*b/6.0; 
    1451      
    1452     return (yy); 
    1453 } 
    1454  
    1455 // 
    1456 static double 
    1457 AlphaSquare(double x) 
    1458 {     
    1459     double yy; 
    1460     yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
    1461          
    1462     return (yy); 
    1463 } 
    1464  
    1465 // ?? funciton is not used - but should the log actually be log10??? 
    1466 static double 
    1467 miu(double x) 
    1468 {     
    1469     double yy; 
    1470     yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); 
    1471      
    1472     return (yy); 
    1473 } 
    1474  
    1475 // 
    1476 static double 
    1477 Sdebye(double q, double L, double b) 
    1478 {     
    1479     double yy; 
    1480     yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 
    1481          
    1482     return (yy); 
    1483 } 
    1484  
    1485 // 
    1486 static double 
    1487 Sdebye1(double q, double L, double b) 
    1488 {     
    1489     double yy; 
    1490     yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 
    1491      
    1492     return (yy); 
    1493 } 
    1494  
    1495 // 
    1496 static double 
    1497 Sexv(double q, double L, double b) 
    1498 {     
    1499     double yy,C1,C2,C3,miu,Rg2; 
    1500     C1=1.22; 
    1501     C2=0.4288; 
    1502     C3=-1.651; 
    1503     miu = 0.585; 
    1504          
    1505     Rg2 = Rgsquare(q,L,b); 
    1506      
    1507     yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
    1508      
    1509     return (yy); 
    1510 } 
    1511  
    1512 // this must be WR modified version 
    1513 static double 
    1514 Sexvnew(double q, double L, double b) 
    1515 {     
    1516     double yy,C1,C2,C3,miu; 
    1517     double del=1.05,C_star2,Rg2; 
    1518      
    1519     C1=1.22; 
    1520     C2=0.4288; 
    1521     C3=-1.651; 
    1522     miu = 0.585; 
    1523          
    1524     //calculating the derivative to decide on the corection (cutoff) term? 
    1525     // I have modified this from WRs original code 
    1526      
    1527     if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 
    1528         C_star2 = 0.0; 
    1529     } else { 
    1530         C_star2 = 1.0; 
    1531     } 
    1532          
    1533     Rg2 = Rgsquare(q,L,b); 
    1534      
    1535         yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
    1536          
    1537     return (yy); 
    1538 } 
    1539  
    1540 // these are the messy ones 
    1541 static double 
    1542 a2short(double q, double L, double b, double p1short, double p2short, double q0) 
    1543 {     
    1544     double yy,Rg2_sh; 
    1545     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 
    1546     double pi; 
    1547          
    1548     E = 2.718281828459045091; 
    1549         pi = 4.0*atan(1.0); 
    1550     Rg2_sh = Rgsquareshort(q,L,b); 
    1551     Rg2_sh2 = Rg2_sh*Rg2_sh; 
    1552     t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    1553     Et1 = pow(E,t1); 
    1554     Emt1 =pow(E,-t1);  
    1555     q02 = q0*q0; 
    1556     q0p = pow(q0,(-4.0 + p2short) ); 
    1557      
    1558     //E is the number e 
    1559     yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); 
    1560          
    1561     return (yy); 
    1562 } 
    1563  
    1564 // 
    1565 static double 
    1566 a1short(double q, double L, double b, double p1short, double p2short, double q0) 
    1567 {     
    1568     double yy,Rg2_sh; 
    1569     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 
    1570         double pi; 
    1571      
    1572     E = 2.718281828459045091; 
    1573         pi = 4.0*atan(1.0); 
    1574     Rg2_sh = Rgsquareshort(q,L,b); 
    1575     Rg2_sh2 = Rg2_sh*Rg2_sh; 
    1576     b3 = b*b*b; 
    1577     t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    1578     Et1 = pow(E,t1); 
    1579     Emt1 =pow(E,-t1);  
    1580     q02 = q0*q0; 
    1581     q0p = pow(q0,(-4.0 + p1short) ); 
    1582  
    1583     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)))))); 
    1584  
    1585     return(yy); 
    1586 } 
    1587  
    1588 // this one will be lots of trouble 
    1589 static double 
    1590 a2long(double q, double L, double b, double p1, double p2, double q0) 
    1591 { 
    1592     double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 
    1593     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 
    1594     double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 
    1595      
    1596     pi = 4.0*atan(1.0); 
    1597     E = 2.718281828459045091; 
    1598     if( L/b > 10.0) { 
    1599         C = 3.06/pow((L/b),0.44); 
    1600     } else { 
    1601         C = 1.0; 
    1602     } 
    1603          
    1604     C1 = 1.22; 
    1605     C2 = 0.4288; 
    1606     C3 = -1.651; 
    1607     C4 = 1.523; 
    1608     C5 = 0.1477; 
    1609     miu = 0.585; 
    1610          
    1611     Rg2 = Rgsquare(q,L,b); 
    1612     Rg22 = Rg2*Rg2; 
    1613     b2 = b*b; 
    1614     b3 = b*b*b; 
    1615     b4 = b3*b; 
    1616     q02 = q0*q0; 
    1617     q03 = q0*q0*q0; 
    1618     q04 = q03*q0; 
    1619     q05 = q04*q0; 
    1620      
    1621     t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
    1622  
    1623     t2 = (b*C*(((-1.0*((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*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 
    1624  
    1625     t3 = (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.0))/(2.0*C5); 
    1626  
    1627     t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
    1628  
    1629     t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    1630  
    1631     t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
    1632  
    1633     t7 = (((-((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)); 
    1634  
    1635     t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    1636  
    1637     t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 
    1638  
    1639     t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    1640  
    1641  
    1642     yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((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))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
    1643  
    1644     return (yy); 
    1645 } 
    1646  
    1647 //need to define this on my own 
    1648 static double 
    1649 sech_WR(double x) 
    1650 {        
    1651         return(1/cosh(x)); 
    1652 } 
    1653  
    1654 // 
    1655 static double 
    1656 a1long(double q, double L, double b, double p1, double p2, double q0) 
    1657 {     
    1658     double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 
    1659     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 
    1660     double E,pi; 
    1661     double b2,b3,b4,q02,q03,q04,q05,Rg22; 
    1662          
    1663     pi = 4.0*atan(1.0); 
    1664     E = 2.718281828459045091; 
    1665      
    1666     if( L/b > 10.0) { 
    1667         C = 3.06/pow((L/b),0.44); 
    1668     } else { 
    1669         C = 1.0; 
    1670     } 
    1671          
    1672     C1 = 1.22; 
    1673     C2 = 0.4288; 
    1674     C3 = -1.651; 
    1675     C4 = 1.523; 
    1676     C5 = 0.1477; 
    1677     miu = 0.585; 
    1678          
    1679     Rg2 = Rgsquare(q,L,b); 
    1680     Rg22 = Rg2*Rg2; 
    1681     b2 = b*b; 
    1682     b3 = b*b*b; 
    1683     b4 = b3*b; 
    1684     q02 = q0*q0; 
    1685     q03 = q0*q0*q0; 
    1686     q04 = q03*q0; 
    1687     q05 = q04*q0; 
    1688      
    1689     t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
    1690      
    1691     t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    1692      
    1693     t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
    1694      
    1695     t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    1696          
    1697     t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 
    1698      
    1699     t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); 
    1700      
    1701     t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
    1702      
    1703     t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
    1704      
    1705     t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    1706      
    1707     t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    1708      
    1709     t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
    1710      
    1711     t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    1712      
    1713     t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
    1714      
    1715     t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    1716      
    1717     t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
    1718          
    1719      
    1720     yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); 
    1721          
    1722     return (yy); 
    17231339} 
    17241340 
    … …  
    21321748} 
    21331749 
    2134 static double 
    2135 gammaln(double xx) 
    2136 { 
    2137         double x,y,tmp,ser; 
    2138         static double cof[6]={76.18009172947146,-86.50532032941677, 
    2139                 24.01409824083091,-1.231739572450155, 
    2140                 0.1208650973866179e-2,-0.5395239384953e-5}; 
    2141         int j; 
    2142          
    2143         y=x=xx; 
    2144         tmp=x+5.5; 
    2145         tmp -= (x+0.5)*log(tmp); 
    2146         ser=1.000000000190015; 
    2147         for (j=0;j<=5;j++) ser += cof[j]/++y; 
    2148         return -tmp+log(2.5066282746310005*ser/x); 
    2149 } 
    2150  
    21511750 
    21521751double 

  • sansmodels/src/sans/models/libigor/libCylinder.h

    r3a084629 r34c2649  
    6767 
    6868/////////functions for WRC implementation of flexible cylinders 
    69 static double Sk_WR(double q, double L, double b); 
    70 static double w_WR(double x); 
    71 static double u1(double q, double L, double b); 
    72 static double u_WR(double q, double L, double b); 
    73 static double Rgsquarezero(double q, double L, double b); 
    74 static double Rgsquareshort(double q, double L, double b); 
    75 static double Rgsquare(double q, double L, double b); 
    76 static double AlphaSquare(double x); 
    77 static double miu(double x); 
    78 static double Sdebye(double q, double L, double b); 
    79 static double Sdebye1(double q, double L, double b); 
    80 static double Sexv(double q, double L, double b); 
    81 static double Sexvnew(double q, double L, double b); 
    82 static double a2short(double q, double L, double b, double p1short, double p2short, double q0); 
    83 static double a1short(double q, double L, double b, double p1short, double p2short, double q0); 
    84 static double a2long(double q, double L, double b, double p1, double p2, double q0); 
    85 static double sech_WR(double x); 
    86 static double a1long(double q, double L, double b, double p1, double p2, double q0); 
    87 static double gammaln(double xx); 
     69static double 
     70gammaln(double xx) 
     71{ 
     72        double x,y,tmp,ser; 
     73        static double cof[6]={76.18009172947146,-86.50532032941677, 
     74                24.01409824083091,-1.231739572450155, 
     75                0.1208650973866179e-2,-0.5395239384953e-5}; 
     76        int j; 
     77 
     78        y=x=xx; 
     79        tmp=x+5.5; 
     80        tmp -= (x+0.5)*log(tmp); 
     81        ser=1.000000000190015; 
     82        for (j=0;j<=5;j++) ser += cof[j]/++y; 
     83        return -tmp+log(2.5066282746310005*ser/x); 
     84} 
     85 
     86// 
     87static double 
     88AlphaSquare(double x) 
     89{ 
     90    double yy; 
     91    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
     92 
     93    return (yy); 
     94} 
     95 
     96// 
     97static double 
     98Rgsquarezero(double q, double L, double b) 
     99{ 
     100    double yy; 
     101    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)))); 
     102 
     103    return (yy); 
     104} 
     105 
     106// 
     107static double 
     108Rgsquareshort(double q, double L, double b) 
     109{ 
     110    double yy; 
     111    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
     112 
     113    return (yy); 
     114} 
     115 
     116// 
     117static double 
     118Rgsquare(double q, double L, double b) 
     119{ 
     120    double yy; 
     121    yy = AlphaSquare(L/b)*L*b/6.0; 
     122 
     123    return (yy); 
     124} 
     125 
     126// ?? funciton is not used - but should the log actually be log10??? 
     127static double 
     128miu(double x) 
     129{ 
     130    double yy; 
     131    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))); 
     132 
     133    return (yy); 
     134} 
     135 
     136//WR named this w (too generic) 
     137static double 
     138w_WR(double x) 
     139{ 
     140    double yy; 
     141    yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 
     142 
     143    return (yy); 
     144} 
     145 
     146// 
     147static double 
     148u1(double q, double L, double b) 
     149{ 
     150    double yy; 
     151 
     152    yy = Rgsquareshort(q,L,b)*q*q; 
     153 
     154    return (yy); 
     155} 
     156 
     157// was named u 
     158static double 
     159u_WR(double q, double L, double b) 
     160{ 
     161    double yy; 
     162    yy = Rgsquare(q,L,b)*q*q; 
     163    return (yy); 
     164} 
     165 
     166 
     167// 
     168static double 
     169Sdebye1(double q, double L, double b) 
     170{ 
     171    double yy; 
     172    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 
     173 
     174    return (yy); 
     175} 
     176 
     177 
     178// 
     179static double 
     180Sdebye(double q, double L, double b) 
     181{ 
     182    double yy; 
     183    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 
     184 
     185    return (yy); 
     186} 
     187 
     188// 
     189static double 
     190Sexv(double q, double L, double b) 
     191{ 
     192    double yy,C1,C2,C3,miu,Rg2; 
     193    C1=1.22; 
     194    C2=0.4288; 
     195    C3=-1.651; 
     196    miu = 0.585; 
     197 
     198    Rg2 = Rgsquare(q,L,b); 
     199 
     200    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))); 
     201 
     202    return (yy); 
     203} 
     204 
     205// this must be WR modified version 
     206static double 
     207Sexvnew(double q, double L, double b) 
     208{ 
     209    double yy,C1,C2,C3,miu; 
     210    double del=1.05,C_star2,Rg2; 
     211 
     212    C1=1.22; 
     213    C2=0.4288; 
     214    C3=-1.651; 
     215    miu = 0.585; 
     216 
     217    //calculating the derivative to decide on the corection (cutoff) term? 
     218    // I have modified this from WRs original code 
     219 
     220    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 
     221        C_star2 = 0.0; 
     222    } else { 
     223        C_star2 = 1.0; 
     224    } 
     225 
     226    Rg2 = Rgsquare(q,L,b); 
     227 
     228        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))); 
     229 
     230    return (yy); 
     231} 
     232 
     233// these are the messy ones 
     234static double 
     235a2short(double q, double L, double b, double p1short, double p2short, double q0) 
     236{ 
     237    double yy,Rg2_sh; 
     238    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 
     239    double pi; 
     240 
     241    E = 2.718281828459045091; 
     242        pi = 4.0*atan(1.0); 
     243    Rg2_sh = Rgsquareshort(q,L,b); 
     244    Rg2_sh2 = Rg2_sh*Rg2_sh; 
     245    t1 = ((q0*q0*Rg2_sh)/(b*b)); 
     246    Et1 = pow(E,t1); 
     247    Emt1 =pow(E,-t1); 
     248    q02 = q0*q0; 
     249    q0p = pow(q0,(-4.0 + p2short) ); 
     250 
     251    //E is the number e 
     252    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))))))); 
     253 
     254    return (yy); 
     255} 
     256 
     257// 
     258static double 
     259a1short(double q, double L, double b, double p1short, double p2short, double q0) 
     260{ 
     261    double yy,Rg2_sh; 
     262    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 
     263        double pi; 
     264 
     265    E = 2.718281828459045091; 
     266        pi = 4.0*atan(1.0); 
     267    Rg2_sh = Rgsquareshort(q,L,b); 
     268    Rg2_sh2 = Rg2_sh*Rg2_sh; 
     269    b3 = b*b*b; 
     270    t1 = ((q0*q0*Rg2_sh)/(b*b)); 
     271    Et1 = pow(E,t1); 
     272    Emt1 =pow(E,-t1); 
     273    q02 = q0*q0; 
     274    q0p = pow(q0,(-4.0 + p1short) ); 
     275 
     276    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)))))); 
     277 
     278    return(yy); 
     279} 
     280 
     281 
     282//need to define this on my own 
     283static double 
     284sech_WR(double x) 
     285{ 
     286        return(1/cosh(x)); 
     287} 
     288 
     289// this one will be lots of trouble 
     290static double 
     291a2long(double q, double L, double b, double p1, double p2, double q0) 
     292{ 
     293    double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 
     294    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 
     295    double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 
     296 
     297    pi = 4.0*atan(1.0); 
     298    E = 2.718281828459045091; 
     299    if( L/b > 10.0) { 
     300        C = 3.06/pow((L/b),0.44); 
     301    } else { 
     302        C = 1.0; 
     303    } 
     304 
     305    C1 = 1.22; 
     306    C2 = 0.4288; 
     307    C3 = -1.651; 
     308    C4 = 1.523; 
     309    C5 = 0.1477; 
     310    miu = 0.585; 
     311 
     312    Rg2 = Rgsquare(q,L,b); 
     313    Rg22 = Rg2*Rg2; 
     314    b2 = b*b; 
     315    b3 = b*b*b; 
     316    b4 = b3*b; 
     317    q02 = q0*q0; 
     318    q03 = q0*q0*q0; 
     319    q04 = q03*q0; 
     320    q05 = q04*q0; 
     321 
     322    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
     323 
     324    t2 = (b*C*(((-1.0*((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*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 
     325 
     326    t3 = (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.0))/(2.0*C5); 
     327 
     328    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
     329 
     330    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     331 
     332    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
     333 
     334    t7 = (((-((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)); 
     335 
     336    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     337 
     338    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 
     339 
     340    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     341 
     342 
     343    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((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))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
     344 
     345    return (yy); 
     346} 
     347// 
     348static double 
     349a1long(double q, double L, double b, double p1, double p2, double q0) 
     350{ 
     351    double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 
     352    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 
     353    double E,pi; 
     354    double b2,b3,b4,q02,q03,q04,q05,Rg22; 
     355 
     356    pi = 4.0*atan(1.0); 
     357    E = 2.718281828459045091; 
     358 
     359    if( L/b > 10.0) { 
     360        C = 3.06/pow((L/b),0.44); 
     361    } else { 
     362        C = 1.0; 
     363    } 
     364 
     365    C1 = 1.22; 
     366    C2 = 0.4288; 
     367    C3 = -1.651; 
     368    C4 = 1.523; 
     369    C5 = 0.1477; 
     370    miu = 0.585; 
     371 
     372    Rg2 = Rgsquare(q,L,b); 
     373    Rg22 = Rg2*Rg2; 
     374    b2 = b*b; 
     375    b3 = b*b*b; 
     376    b4 = b3*b; 
     377    q02 = q0*q0; 
     378    q03 = q0*q0*q0; 
     379    q04 = q03*q0; 
     380    q05 = q04*q0; 
     381 
     382    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)))); 
     383 
     384    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)))))); 
     385 
     386    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)))); 
     387 
     388    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     389 
     390    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 
     391 
     392    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))); 
     393 
     394    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)); 
     395 
     396    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
     397 
     398    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)))))); 
     399 
     400    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)))))); 
     401 
     402    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)); 
     403 
     404    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     405 
     406    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)))); 
     407 
     408    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)))))); 
     409 
     410    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)))); 
     411 
     412 
     413    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))))))))))); 
     414 
     415    return (yy); 
     416} 
     417 
     418 
     419static double 
     420Sk_WR(double q, double L, double b) 
     421{ 
     422        // 
     423        double p1,p2,p1short,p2short,q0; 
     424        double C,ans,q0short,Sexvmodify,pi; 
     425 
     426        pi = 4.0*atan(1.0); 
     427 
     428        p1 = 4.12; 
     429        p2 = 4.42; 
     430        p1short = 5.36; 
     431        p2short = 5.62; 
     432        q0 = 3.1; 
     433        // 
     434        q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
     435 
     436        // 
     437        if(L/b > 10.0) { 
     438                C = 3.06/pow((L/b),0.44); 
     439        } else { 
     440                C = 1.0; 
     441        } 
     442        // 
     443 
     444        if( L > 4*b ) { // Longer Chains 
     445                if (q*b <= 3.1) {               //Modified by Yun on Oct. 15, 
     446                        Sexvmodify = Sexvnew(q, L, b); 
     447                        ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
     448                } else { //q(i)*b > 3.1 
     449                        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); 
     450                } 
     451        } else { //L <= 4*b Shorter Chains 
     452                if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
     453                        if (q*b<=0.01) { 
     454                                ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
     455                        } else { 
     456                                ans = Sdebye1(q,L,b); 
     457                        } 
     458                } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
     459                        ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); 
     460                } 
     461        } 
     462 
     463        return(ans); 
     464        //return(a2long(q, L, b, p1, p2, q0)); 
     465} 
     466 

  • sansmodels/src/sans/models/libigor/libSphere.c

    r3a084629 r34c2649  
    653653} 
    654654 
     655/* 
    655656static double 
    656657LogNormal_distr(double sig, double mu, double pt) 
    … …  
    672673        return(retval); 
    673674} 
     675*/ 
    674676 
    675677// scattering from a core shell sphere with a (Schulz) polydisperse core and constant ratio (shell thickness)/(core radius) 
    … …  
    12201222} 
    12211223 
     1224/* 
    12221225static double 
    12231226SchulzPoint(double x, double avg, double zz) { 
    … …  
    12451248    return -tmp+log(2.5066282746310005*ser/x); 
    12461249} 
     1250*/ 
    12471251 
    12481252double 
    … …  
    12621266        // variables are: 
    12631267        //[0] scale factor 
    1264         //[1] radius of core [] 
    1265         //[2] SLD of the core   [-2] 
    1266         //[3] thickness of the shell    [] 
     1268        //[1] radius of core [ï¿œ] 
     1269        //[2] SLD of the core   [ï¿œ-2] 
     1270        //[3] thickness of the shell    [ï¿œ] 
    12671271        //[4] SLD of the shell 
    12681272        //[5] SLD of the solvent 
    … …  
    13051309        f += vol*bes*contr; 
    13061310         
    1307         // normalize to particle volume and rescale from [-1] to [cm-1] 
     1311        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1] 
    13081312        f2 = f*f/vol*1.0e8; 
    13091313         
    … …  
    13211325        // variables are: 
    13221326        //[0] scale factor 
    1323         //[1] radius of core [] 
    1324         //[2] SLD of the core   [-2] 
    1325         //[3] thickness of shell 1 [] 
     1327        //[1] radius of core [ï¿œ] 
     1328        //[2] SLD of the core   [ï¿œ-2] 
     1329        //[3] thickness of shell 1 [ï¿œ] 
    13261330        //[4] SLD of shell 1 
    1327         //[5] thickness of shell 2 [] 
     1331        //[5] thickness of shell 2 [ï¿œ] 
    13281332        //[6] SLD of shell 2 
    13291333        //[7] SLD of the solvent 
    … …  
    13791383         
    13801384                 
    1381         // normalize to particle volume and rescale from [-1] to [cm-1] 
     1385        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1] 
    13821386        f2 = f*f/vol*1.0e8; 
    13831387         
    … …  
    13951399        // variables are: 
    13961400        //[0] scale factor 
    1397         //[1] radius of core [] 
    1398         //[2] SLD of the core   [-2] 
    1399         //[3] thickness of shell 1 [] 
     1401        //[1] radius of core [ï¿œ] 
     1402        //[2] SLD of the core   [ï¿œ-2] 
     1403        //[3] thickness of shell 1 [ï¿œ] 
    14001404        //[4] SLD of shell 1 
    1401         //[5] thickness of shell 2 [] 
     1405        //[5] thickness of shell 2 [ï¿œ] 
    14021406        //[6] SLD of shell 2 
    14031407        //[7] thickness of shell 3 
    … …  
    14671471        f += vol*bes*contr; 
    14681472                 
    1469         // normalize to particle volume and rescale from [-1] to [cm-1] 
     1473        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1] 
    14701474        f2 = f*f/vol*1.0e8; 
    14711475         
    … …  
    14831487        // variables are: 
    14841488        //[0] scale factor 
    1485         //[1] radius of core [] 
    1486         //[2] SLD of the core   [-2] 
    1487         //[3] thickness of shell 1 [] 
     1489        //[1] radius of core [ï¿œ] 
     1490        //[2] SLD of the core   [ï¿œ-2] 
     1491        //[3] thickness of shell 1 [ï¿œ] 
    14881492        //[4] SLD of shell 1 
    1489         //[5] thickness of shell 2 [] 
     1493        //[5] thickness of shell 2 [ï¿œ] 
    14901494        //[6] SLD of shell 2 
    14911495        //[7] thickness of shell 3 
    … …  
    15701574         
    15711575                 
    1572         // normalize to particle volume and rescale from [-1] to [cm-1] 
     1576        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1] 
    15731577        f2 = f*f/vol*1.0e8; 
    15741578         
    … …  
    19972001        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius); 
    19982002        vol = 4.0*pi/3.0*radius*radius*radius; 
    1999         f = vol*bes*delrho      ;       // [=]  
     2003        f = vol*bes*delrho      ;       // [=] ï¿œ 
    20002004        // normalize to single particle volume, convert to 1/cm 
    20012005        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 

  • sansmodels/src/sans/models/libigor/libSphere.h

    r3a084629 r34c2649  
    4949 
    5050 
    51 static double SchulzPoint(double x, double avg, double zz); 
    52 static double gammln(double xx); 
    53 static double Gauss_distr(double sig, double avg, double pt); 
    54 static double LogNormal_distr(double sig, double mu, double pt); 
     51static double 
     52gammln(double xx) { 
     53    double x,y,tmp,ser; 
     54    static double cof[6]={76.18009172947146,-86.50532032941677, 
     55                24.01409824083091,-1.231739572450155, 
     56                0.1208650973866179e-2,-0.5395239384953e-5}; 
     57    int j; 
    5558 
     59    y=x=xx; 
     60    tmp=x+5.5; 
     61    tmp -= (x+0.5)*log(tmp); 
     62    ser=1.000000000190015; 
     63    for (j=0;j<=5;j++) ser += cof[j]/++y; 
     64    return -tmp+log(2.5066282746310005*ser/x); 
     65} 
     66 
     67static double 
     68LogNormal_distr(double sig, double mu, double pt) 
     69{ 
     70        double retval,pi; 
     71 
     72        pi = 4.0*atan(1.0); 
     73        retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
     74        return(retval); 
     75} 
     76 
     77static double 
     78Gauss_distr(double sig, double avg, double pt) 
     79{ 
     80        double retval,Pi; 
     81 
     82        Pi = 4.0*atan(1.0); 
     83        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0); 
     84        return(retval); 
     85} 
     86 
     87static double SchulzPoint(double x, double avg, double zz) { 
     88    double dr; 
     89    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 
     90    return (exp(dr)); 
     91}; 
     92 

  • sansmodels/src/sans/models/libigor/libTwoPhase.c

    r3a084629 r34c2649  
    421421{ 
    422422        //w[0] = scale 
    423         //w[1] = radius of gyration [] 
     423        //w[1] = radius of gyration [ï¿œ] 
    424424        //w[2] = polydispersity, ratio of Mw/Mn 
    425425        //w[3] = bkg [cm-1] 
    … …  
    475475 
    476476 
    477 static double 
    478 gammln(double xx) { 
    479          
    480     double x,y,tmp,ser; 
    481     static double cof[6]={76.18009172947146,-86.50532032941677, 
    482                 24.01409824083091,-1.231739572450155, 
    483                 0.1208650973866179e-2,-0.5395239384953e-5}; 
    484     int j; 
    485          
    486     y=x=xx; 
    487     tmp=x+5.5; 
    488     tmp -= (x+0.5)*log(tmp); 
    489     ser=1.000000000190015; 
    490     for (j=0;j<=5;j++) ser += cof[j]/++y; 
    491     return -tmp+log(2.5066282746310005*ser/x); 
    492 } 
    493  
    494477double 
    495478GaussianShell(double w[], double x) 
    … …  
    497480        // variables are:                                                        
    498481        //[0] scale 
    499         //[1] radius () 
    500         //[2] thick () (thickness parameter - this is the std. dev. of the Gaussian width of the shell) 
     482        //[1] radius (ï¿œ) 
     483        //[2] thick (ï¿œ) (thickness parameter - this is the std. dev. of the Gaussian width of the shell) 
    501484        //[3] polydispersity of the radius 
    502         //[4] sld shell (-2) 
     485        //[4] sld shell (ï¿œ-2) 
    503486        //[5] sld solvent 
    504487        //[6] background (cm-1) 

  • sansmodels/src/sans/models/libigor/libTwoPhase.h

    r3a084629 r34c2649  
    2424 
    2525/* internal functions */ 
    26 static double gammln(double xx); 
     26static double 
     27gammln(double xx) { 
    2728 
     29    double x,y,tmp,ser; 
     30    static double cof[6]={76.18009172947146,-86.50532032941677, 
     31                24.01409824083091,-1.231739572450155, 
     32                0.1208650973866179e-2,-0.5395239384953e-5}; 
     33    int j; 
     34 
     35    y=x=xx; 
     36    tmp=x+5.5; 
     37    tmp -= (x+0.5)*log(tmp); 
     38    ser=1.000000000190015; 
     39    for (j=0;j<=5;j++) ser += cof[j]/++y; 
     40    return -tmp+log(2.5066282746310005*ser/x); 
     41} 
     42 
     43 
Note: See TracChangeset for help on using the changeset viewer.