- Timestamp:
- 10/19/12 17:25:35 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STGrid.cpp
r2671 r2673 1598 1598 #define grdjinc1 grdjinc1_ 1599 1599 #endif 1600 #if defined(USE_CASAPY) 1600 1601 extern "C" { 1601 1602 void grdsf(Double*, Double*); … … 1603 1604 void grdjinc1(Double*, Double*, Int*, Double*); 1604 1605 } 1606 #else 1607 extern "C" { 1608 void grdsf(Double*, Double*); 1609 } 1610 void grdgauss(Double *hwhm, Double *val, Double *out) 1611 { 1612 *out = exp(-log(2.0) * (*val / *hwhm) * (*val / *hwhm)); 1613 } 1614 void grdjinc1(Double *c, Double *val, Int *normalize, Double *out) 1615 { 1616 LogIO os(LogOrigin("","grdjinc1",WHERE)); 1617 os << LogIO::SEVERE 1618 << "Jinc gridding is not supported" << LogIO::POST; 1619 1620 // Calculate J_1(x) using approximate formula 1621 Double x = C::pi * *val / *c; 1622 Double ax = fabs(x); 1623 Double ans; 1624 if ( ax < 8.0 ) { 1625 Double y = x * x; 1626 Double ans1 = x * (72362614232.0 + y * (-7895059235.0 1627 + y * (242396853.1 + y * (-2972611.439 1628 + y * (15704.48260 + y * (-30.16036606)))))); 1629 Double ans2 = 144725228442.0 + y * (2300535178.0 1630 + y * (18583304.74 + y * (99447.43394 1631 + y * (376.9991397 + y * 1.0)))); 1632 ans = ans1 / ans2; 1633 } 1634 else { 1635 Double z = 8.0 / ax; 1636 Double y = z * z; 1637 Double xx = ax - 2.356194491; 1638 Double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4 1639 + y * (0.2457520174e-5 + y * (-0.240337019e-6)))); 1640 Double ans2 = 0.04687499995 + y * (-0.2002690873e-3 1641 + y * (0.8449199096e-5 + y * (-0.88228987e-6 1642 + y * (0.105787412e-6)))); 1643 ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 1644 - z * sin(xx) * ans2); 1645 if (x < 0.0) 1646 ans = -ans; 1647 } 1648 1649 // Then, calculate Jinc 1650 if (x == 0.0) { 1651 *out = 0.5; 1652 } 1653 else { 1654 *out = ans / x; 1655 } 1656 1657 if (*normalize == 1) 1658 *out = *out / 0.5; 1659 } 1660 #endif 1605 1661 void STGrid::spheroidalFunc( Vector<Float> &convFunc ) 1606 1662 { … … 1767 1823 << "convType_ = " << convType_ << endl 1768 1824 << "convSupport_ = " << convSupport_ << endl 1769 << "truncate_ = " << truncate << endl1825 << "truncate_ = " << truncate << "pixel" << endl 1770 1826 << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST; 1771 1827 } … … 1815 1871 << "convType_ = " << convType_ << endl 1816 1872 << "convSupport_ = " << convSupport_ << endl 1817 << "truncate_ = " << truncate << endl1873 << "truncate_ = " << truncate << "pixel" << endl 1818 1874 << "gwidth_ = " << pixelGW << "pixel" << endl 1819 1875 << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
Note:
See TracChangeset
for help on using the changeset viewer.