Commits
1294 1294 | |
1295 1295 | |
1296 1296 | |
1297 1297 | void deleteEphemeris(map<AtmPhaseCorrection, Table*>& apc2EphemTable_m) { |
1298 1298 | for ( map<AtmPhaseCorrection, Table*>::iterator iter = apc2EphemTable_m.begin(); iter!=apc2EphemTable_m.end(); ++iter) |
1299 1299 | if (apc2EphemTable_m[iter->first] != NULL) delete apc2EphemTable_m[iter->first]; |
1300 1300 | } |
1301 1301 | |
1302 1302 | std::map<int, double>ephemStartTime_m; |
1303 1303 | |
1304 - | void fillEphemeris(ASDM* ds_p, uint64_t timeStepInNanoSecond, bool interpolate_ephemeris, bool tabulate_ephemeris_polynomials, string telescopeName) { |
1304 + | void fillEphemeris(ASDM* ds_p, uint64_t timeStepInNanoSecond, bool interpolate_ephemeris, string telescopeName) { |
1305 1305 | LOGENTER("fillEphemeris"); |
1306 1306 | |
1307 1307 | try { |
1308 1308 | // Retrieve the Ephemeris table's content. |
1309 1309 | EphemerisTable & eT = ds_p->getEphemeris(); |
1310 1310 | const vector<EphemerisRow *>& eR_v = eT.get(); |
1311 1311 | |
1312 1312 | infostream.str(""); |
1313 1313 | infostream << "The dataset has " << eT.size() << " ephemeris row(s)..."; |
1314 1314 | info(infostream.str()); |
1359 1359 | * 1) MJD0 |
1360 1360 | * 2) DMJD0 |
1361 1361 | * 3) NAME |
1362 1362 | * 4) GeoLong in deg. |
1363 1363 | * 5) GeoLat in deg. |
1364 1364 | * 6) GeoDist in km. |
1365 1365 | * 7) origin (a textual summary of the origin of this ephemeris) |
1366 1366 | * |
1367 1367 | * We derive these values from the informations found in the first element of i2e_m[ephemerisId] |
1368 1368 | */ |
1369 - | vector<EphemerisRow *>& v = i2e_m[ephemerisId]; |
1370 - | vector<double> observerLocation = v[0]->getObserverLocation(); |
1369 + | vector<EphemerisRow *>& ephRow_v = i2e_m[ephemerisId]; |
1370 + | vector<double> observerLocation = ephRow_v[0]->getObserverLocation(); |
1371 1371 | double geoLong = radian2degree(observerLocation[0]); // in order to get degrees. |
1372 1372 | double geoLat = radian2degree(observerLocation[1]); // in order to get degrees. |
1373 1373 | double geoDist = observerLocation[2] / 1000.0; // in order to get km (supposedly above the reference ellipsoid) |
1374 1374 | |
1375 - | int64_t t0ASDM = v[0]->getTimeInterval().getStart().get(); // The first time recorded for this ephemerisId. |
1375 + | int64_t t0ASDM = ephRow_v[0]->getTimeInterval().getStart().get(); // The first time recorded for this ephemerisId. |
1376 1376 | int64_t q = t0ASDM / timeStepInNanoSecond; |
1377 1377 | int64_t r = t0ASDM % timeStepInNanoSecond; |
1378 1378 | int64_t t0MS = t0ASDM; |
1379 1379 | if ( r != 0 ) { |
1380 1380 | q = q + 1; |
1381 1381 | t0MS = q * timeStepInNanoSecond; |
1382 1382 | } |
1383 1383 | |
1384 1384 | double mjd0 = ArrayTime(t0MS).getMJD(); |
1385 1385 | |
1386 - | double dmjd = interpolate_ephemeris ? 0.001 : v[0]->getTimeInterval().getDuration().get() / 1000000000LL / 86400.0; // Grid time step == 0.001 if ephemeris interpolation requested |
1386 + | double dmjd = interpolate_ephemeris ? 0.001 : ephRow_v[0]->getTimeInterval().getDuration().get() / 1000000000LL / 86400.0; // Grid time step == 0.001 if ephemeris interpolation requested |
1387 1387 | // otherwise == the interval of time of the first element of ephemeris converted in days. |
1388 1388 | // *SUPPOSEDLY* constant over all the ephemeris. |
1389 1389 | |
1390 1390 | // determine the position reference system |
1391 - | double equator = v[0]->getEquinoxEquator(); |
1391 + | double equator = ephRow_v[0]->getEquinoxEquator(); |
1392 1392 | string posref = "unknown"; |
1393 1393 | if (equator == 2000.) { // the Ephemeris table presently only stores the equator |
1394 1394 | posref = "ICRF/ICRS"; |
1395 1395 | } |
1396 1396 | |
1397 1397 | // Prepare the table keywords with the values computed above. |
1398 1398 | TableDesc tableDesc; |
1399 1399 | |
1400 - | tableDesc.comment() = v[0]->getOrigin(); |
1400 + | tableDesc.comment() = ephRow_v[0]->getOrigin(); |
1401 1401 | time_t now = time(0); |
1402 1402 | struct tm tstruct; |
1403 1403 | char buf[80]; |
1404 1404 | |
1405 1405 | tstruct = *gmtime(&now); |
1406 1406 | |
1407 1407 | strftime(buf, sizeof(buf), "%Y/%m/%d/%H:%M", &tstruct); |
1408 1408 | string creationDate(buf); |
1409 1409 | |
1410 1410 | tableDesc.rwKeywordSet().define("VS_CREATE", creationDate); |
1492 1492 | // Now it's time to fill the EPHEM table(s). |
1493 1493 | // |
1494 1494 | // Below the vectors which will be used to populate slices of Ephemeris tables in the MSes. |
1495 1495 | // |
1496 1496 | vector<double> mjdMS_v; |
1497 1497 | vector<double> raMS_v; |
1498 1498 | vector<double> decMS_v; |
1499 1499 | vector<double> distanceMS_v; |
1500 1500 | vector<double> radVelMS_v; |
1501 1501 | |
1502 - | bool numPolyDirIsOne = v[0]->getNumPolyDir() == 1; |
1503 - | bool numPolyDistIsOne = v[0]->getNumPolyDist() == 1; |
1504 - | bool radVelExists = v[0]->isRadVelExists() && v[0]->isNumPolyRadVelExists(); |
1505 - | bool numPolyRadVelIsOne = radVelExists ? v[0]->getNumPolyRadVel() == 1 : false; |
1502 + | bool numPolyDirIsOne = ephRow_v[0]->getNumPolyDir() == 1; |
1503 + | bool numPolyDistIsOne = ephRow_v[0]->getNumPolyDist() == 1; |
1504 + | bool radVelExists = ephRow_v[0]->isRadVelExists() && ephRow_v[0]->isNumPolyRadVelExists(); |
1505 + | bool numPolyRadVelIsOne = radVelExists ? ephRow_v[0]->getNumPolyRadVel() == 1 : false; |
1506 1506 | |
1507 + | LOG ("numPolyDirIsOne = " + TO_STRING(numPolyDirIsOne)); |
1508 + | LOG ("numPolyDistIsOne = " + TO_STRING(numPolyDistIsOne)); |
1509 + | LOG ("radVelExists = " + TO_STRING(radVelExists)); |
1510 + | LOG ("numPolyRadVelIsOne = " + TO_STRING(numPolyRadVelIsOne)); |
1507 1511 | |
1508 - | |
1509 - | if (interpolate_ephemeris) { |
1512 + | // the single row case |
1513 + | if (ephRow_v.size()==1) { |
1514 + | |
1515 + | if (timeStepInNanoSecond <= 0.0) { |
1516 + | timeStepInNanoSecond = 0.001; |
1517 + | } |
1518 + | |
1519 + | infostream.str(""); |
1520 + | infostream << "The MS Ephemeris table for ephemerisId = '" << ephemerisId |
1521 + | << "' will be produced by tabulating the polynomials found in the single row for 'dir', 'distance' and optionally 'radVel' with a timestep of '" |
1522 + | << timeStepInNanoSecond / 1.e9 / 86400. << "' days"; |
1523 + | info(infostream.str()); |
1524 + | |
1525 + | dmjd = timeStepInNanoSecond / 1.e9 / 86400.; // to be written in the keyword DMJD |
1526 + | |
1510 1527 | // |
1511 - | // Check that for each polynomial column the degree is always null or never null and that the optional fields are always present or always absent. |
1512 - | // And also verify that there is no "hole" in the time range covered by the sequence of ArrayTime intervals when the degree is == 0. |
1528 + | // Calculate the grid of times where the polynomials will be tabulated. |
1529 + | // This grid contains all the times which : |
1530 + | // - are multiple of the tabulation time steps, |
1531 + | // - contained in the arraytime interval of validity of the current ASDM Ephemeris row. |
1532 + | // |
1533 + | // At this point times are expressed in nanoseconds. |
1513 1534 | // |
1514 - | vector<double> duration_v; // In seconds |
1515 - | vector<double> time_v; // " " |
1516 - | vector<double> timeOffset_v; // (midpoint - origin) in seconds |
1517 1535 | |
1518 - | time_v.push_back(1.0e-09*v[0]->getTimeInterval().getStart().get()); |
1519 - | timeOffset_v.push_back(1.0e-09*(v[0]->getTimeInterval().getMidPoint().get()-v[0]->getTimeOrigin().get())); |
1536 + | // this case has just one element in ephRow_v |
1537 + | |
1538 + | int64_t tstartASDM = ephRow_v[0]->getTimeInterval().getStart().get(); |
1539 + | int64_t tendASDM = tstartASDM + ephRow_v[0]->getTimeInterval().getDuration().get(); |
1540 + | int64_t q = tstartASDM / timeStepInNanoSecond; |
1541 + | int64_t r = tstartASDM % timeStepInNanoSecond; |
1542 + | int64_t tstartMS = tstartASDM; |
1543 + | |
1544 + | if ( r!= 0 ) { |
1545 + | q = q + 1; |
1546 + | tstartMS = q * timeStepInNanoSecond; |
1547 + | } |
1548 + | |
1549 + | vector<int64_t> tabulation_time_v; |
1550 + | int64_t t = tstartMS - timeStepInNanoSecond; // One extra timestep before the beginning of the interval of validity. |
1551 + | do { |
1552 + | tabulation_time_v.push_back(t); |
1553 + | t += timeStepInNanoSecond; |
1554 + | } |
1555 + | while (t <= tendASDM); |
1556 + | tabulation_time_v.push_back(t); // One extra timestep after the end of the interval of validity. |
1557 + | |
1558 + | // |
1559 + | // Tabulate the MS Ephemeris columns for each tabulation time, in s. |
1560 + | // |
1561 + | double timeOrigin = ephRow_v[0]->getTimeOrigin().get() * 1.0e-09 / 86400. ; // to "days" |
1562 + | |
1563 + | // Convert from `radians to degrees. |
1564 + | const vector<vector<double> >& dir_v = ephRow_v[0]->getDir(); |
1565 + | vector<double> ra_coeff_v; |
1566 + | vector<double> dec_coeff_v; |
1567 + | for (unsigned int idir = 0; idir < dir_v.size(); idir++) { |
1568 + | ra_coeff_v.push_back(dir_v[idir][0]); |
1569 + | dec_coeff_v.push_back(dir_v[idir][1]); |
1570 + | } |
1571 + | |
1572 + | LOG(" There are " + TO_STRING(ra_coeff_v.size()) + " RA coeffs and " + TO_STRING(dec_coeff_v.size()) + " DEC coeffs"); |
1573 + | |
1574 + | // Convert radian to degree. |
1575 + | std::transform(ra_coeff_v.begin(), ra_coeff_v.end(), ra_coeff_v.begin(), radian2degree); |
1576 + | std::transform(dec_coeff_v.begin(), dec_coeff_v.end(), dec_coeff_v.begin(), radian2degree); |
1577 + | |
1578 + | // Convert from m to AU. |
1579 + | vector<double> distance_coeff_v = ephRow_v[0]->getDistance(); |
1580 + | std::transform(distance_coeff_v.begin(), distance_coeff_v.end(), distance_coeff_v.begin(), m2au); |
1581 + | |
1582 + | vector<double> radvel_coeff_v; |
1583 + | if (ephRow_v[0]->isRadVelExists()) { |
1584 + | radvel_coeff_v = ephRow_v[0]->getRadVel(); |
1585 + | // Convert from m per s to AU per day. |
1586 + | std::transform(radvel_coeff_v.begin(), radvel_coeff_v.end(), radvel_coeff_v.begin(), mpers2auperd); |
1587 + | } |
1588 + | |
1589 + | // And proceed... |
1590 + | LOG ("There will be " + TO_STRING(tabulation_time_v.size()) + " time steps used to tabulate the polynomials."); |
1591 + | for (unsigned int itab = 0; itab < tabulation_time_v.size(); itab++) { |
1592 + | double tabulation_time = tabulation_time_v[itab] * 1.0e-09 / 86400.0 ; // It appeared that times should be expressed in "day" !! |
1593 + | |
1594 + | // MJD |
1595 + | mjdMS_v.push_back(ArrayTime(tabulation_time_v[itab]).getMJD()); |
1596 + | |
1597 + | // RA / DEC |
1598 + | raMS_v.push_back(evalPoly(ra_coeff_v.size(), |
1599 + | ra_coeff_v, |
1600 + | timeOrigin, |
1601 + | tabulation_time)); |
1520 1602 | |
1603 + | decMS_v.push_back(evalPoly(dec_coeff_v.size(), |
1604 + | dec_coeff_v, |
1605 + | timeOrigin, |
1606 + | tabulation_time)); |
1607 + | |
1608 + | // DISTANCE |
1609 + | distanceMS_v.push_back(evalPoly(distance_coeff_v.size(), |
1610 + | distance_coeff_v, |
1611 + | timeOrigin, |
1612 + | tabulation_time)); |
1613 + | |
1614 + | // RADVEL |
1615 + | if (radVelExists) { |
1616 + | radVelMS_v.push_back(evalPoly(radvel_coeff_v.size(), |
1617 + | radvel_coeff_v, |
1618 + | timeOrigin, |
1619 + | tabulation_time)); |
1620 + | } |
1621 + | } |
1622 + | } else { |
1623 + | // both of these cases require that the numPoly* is either all 1 or always greater than 1. Check here. |
1521 1624 | errstream.str(""); |
1522 - | for (unsigned int i = 1; i < v.size(); i++) { |
1523 - | if (numPolyDirIsOne != (v[i]->getNumPolyDir() == 1)) { |
1625 + | for (unsigned int i = 1; i < ephRow_v.size(); i++) { |
1626 + | if (numPolyDirIsOne != (ephRow_v[i]->getNumPolyDir() == 1)) { |
1524 1627 | errstream << "In the table Ephemeris the value of the field 'numPolyDir' is expected to be always equal to 1 or always greater than 1. This rule is violated at line #" << i <<"."; |
1525 1628 | error(errstream.str()); |
1526 1629 | } |
1527 - | |
1528 - | if (numPolyDistIsOne != (v[i]->getNumPolyDist() == 1)) { |
1630 + | |
1631 + | if (numPolyDistIsOne != (ephRow_v[i]->getNumPolyDist() == 1)) { |
1529 1632 | errstream << "In the table Ephemeris the value of the field 'numPolyDist' is expected to be always equal to 1 or always greater than 1. This rule is violated at line #" << i <<"."; |
1530 1633 | error(errstream.str()); |
1531 1634 | } |
1532 - | |
1533 - | if (radVelExists != (v[i]->isRadVelExists() && v[i]->isNumPolyRadVelExists())) { |
1635 + | |
1636 + | if (radVelExists != (ephRow_v[i]->isRadVelExists() && ephRow_v[i]->isNumPolyRadVelExists())) { |
1534 1637 | errstream << "In the table Ephemeris the fields 'radVel' and 'numPolyRadVel' are expected to be always absent or always present. This rule is violated at line #" << i <<"."; |
1535 1638 | error(errstream.str()); |
1536 1639 | } |
1537 - | |
1640 + | |
1538 1641 | if (radVelExists) { |
1539 - | if (numPolyRadVelIsOne != (v[i]->getNumPolyRadVel() == 1)) { |
1642 + | if (numPolyRadVelIsOne != (ephRow_v[i]->getNumPolyRadVel() == 1)) { |
1540 1643 | errstream << "In the table Ephemeris the value of the field 'numPolyRadVel' is expected to be always equal to 1 or always greater than 1. This rule is violated at line #" << i <<"."; |
1541 1644 | error(errstream.str()); |
1542 1645 | } |
1543 1646 | } |
1544 - | |
1647 + | } |
1648 + | if (!interpolate_ephemeris && (numPolyDirIsOne && numPolyDistIsOne && (!radVelExists || numPolyRadVelIsOne))) { |
1649 + | // interpolation is NOT requested and all possible polynomial columns are simple scalars, numPoly==1 |
1650 + | // Just copy ephemeris without any interpolation. Just adapt the units. |
1651 + | infostream.str(""); |
1652 + | infostream << "The MS Ephemeris table for ephemerisId = '" << ephemerisId |
1653 + | << "' will be produced by copying the values found in the ASDM with no interpolation"; |
1654 + | info(infostream.str()); |
1655 + | for(const EphemerisRow *eR_p: ephRow_v) { |
1656 + | mjdMS_v.push_back(eR_p->getTimeInterval().getMidPoint().getMJD()); // MJD |
1657 + | vector<vector<double> > dir = eR_p->getDir(); |
1658 + | raMS_v.push_back(radian2degree(dir[0][0])); // deg |
1659 + | decMS_v.push_back(radian2degree(dir[0][1])); // deg |
1660 + | distanceMS_v.push_back(m2au(eR_p->getDistance()[0])); // AU |
1661 + | if (radVelExists) { |
1662 + | radVelMS_v.push_back(mpers2auperd(eR_p->getRadVel()[0])); // AU/d |
1663 + | } |
1664 + | } |
1665 + | } else { |
1666 + | // the general case, polynomials are evaluated and scalars are interpolated to the same timesteps. |
1667 + | infostream.str(""); |
1668 + | infostream << "The MS Ephemeris table for ephemerisId = '" << ephemerisId |
1669 + | << "' will be produced by tabulating the polynomials found in the 'dir', 'distance' and optionally 'radVel' columns with a timestep of '" |
1670 + | << timeStepInNanoSecond / 1.e9 / 86400. << "' days"; |
1671 + | info(infostream.str()); |
1545 1672 | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1546 - | int64_t start_i = v[i]->getTimeInterval().getStart().get() ; |
1547 - | int64_t start_i_1 = v[i-1]->getTimeInterval().getStart().get(); |
1548 - | int64_t duration_i_1 = v[i-1]->getTimeInterval().getDuration().get(); |
1549 - | if (start_i != (start_i_1 + duration_i_1)) { |
1550 - | infostream.str(""); |
1551 - | infostream << "The value of 'timeInterval' at row #" << i-1 << " does not cover the time range up to the start time of the next row. The polynomial will be evaluated despite the presence of this 'hole'"; |
1552 - | info(infostream.str()); |
1673 + | infostream.str(""); |
1674 + | infostream << "the polynomials of order 0 will be interpolated between mid-interval values using the same timestep."; |
1675 + | info(infostream.str()); |
1676 + | } |
1677 + | |
1678 + | // Check that for each polynomial column the degree is always null or never null and that the optional fields are always present or always absent. |
1679 + | // And also verify that there is no "hole" in the time range covered by the sequence of ArrayTime intervals when the degree is == 0. |
1680 + | |
1681 + | vector<double> duration_v; // In seconds |
1682 + | vector<double> time_v; // " " |
1683 + | vector<double> timeOffset_v; // (midpoint - origin) in seconds |
1684 + | |
1685 + | time_v.push_back(1.0e-09*ephRow_v[0]->getTimeInterval().getStart().get()); |
1686 + | timeOffset_v.push_back(1.0e-09*(ephRow_v[0]->getTimeInterval().getMidPoint().get()-ephRow_v[0]->getTimeOrigin().get())); |
1687 + | |
1688 + | // look for holes in the intervals for the degree == 0 case |
1689 + | for (unsigned int i = 1; i < ephRow_v.size(); i++) { |
1690 + | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1691 + | int64_t start_i = ephRow_v[i]->getTimeInterval().getStart().get() ; |
1692 + | int64_t start_i_1 = ephRow_v[i-1]->getTimeInterval().getStart().get(); |
1693 + | int64_t duration_i_1 = ephRow_v[i-1]->getTimeInterval().getDuration().get(); |
1694 + | if (start_i != (start_i_1 + duration_i_1)) { |
1695 + | infostream.str(""); |
1696 + | infostream << "The value of 'timeInterval' at row #" << i-1 << " does not cover the time range up to the start time of the next row. The polynomial will be evaluated despite the presence of this 'hole'"; |
1697 + | info(infostream.str()); |
1698 + | } |
1699 + | duration_v.push_back(1.0e-09*(start_i - start_i_1)); |
1700 + | time_v.push_back(1.0e-09*start_i); |
1553 1701 | } |
1554 - | duration_v.push_back(1.0e-09*(start_i - start_i_1)); |
1555 - | time_v.push_back(1.0e-09*start_i); |
1556 1702 | } |
1557 - | } |
1558 - | |
1559 - | LOG ("numPolyDirIsOne = " + TO_STRING(numPolyDirIsOne)); |
1560 - | LOG ("numPolyDistIsOne = " + TO_STRING(numPolyDistIsOne)); |
1561 - | LOG ("radVelExists = " + TO_STRING(radVelExists)); |
1562 - | LOG ("numPolyRadVelIsOne = " + TO_STRING(numPolyRadVelIsOne)); |
1563 1703 | |
1704 + | // The number of tabulated values (i.e. the number of rows in the MS Ephemerides) table depend on the |
1705 + | // degrees of each polynomial column. If there is at least one such degree which is equal to 1 then |
1706 + | // we exclude the last element of the vector ephRow_v, i.e. the last ArrayTimeInterval, since on this time interval |
1707 + | // we would miss one end value for the interpolation (so far extrapolation is excluded). |
1564 1708 | |
1565 - | // |
1566 - | // The number of tabulated values (i.e. the number of rows in the MS Ephemerides) table depend on the |
1567 - | // degrees of each polynomial column. If there is at least one such degree which is equal to 1 then |
1568 - | // we exclude the last element of the vector v, i.e. the last ArrayTimeInterval, since on this time interval |
1569 - | // we would miss one end value for the interpolation (so far extrapolation is excluded). |
1570 - | // |
1571 - | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1572 - | // Then just "forget" the last element. |
1573 - | LOG("Erasing the last element of v (size before = '" + TO_STRING(v.size()) + "')"); |
1574 - | v.erase(v.begin() + v.size() - 1); |
1575 - | LOG("Erasing the last element of v (size after = '" + TO_STRING(v.size()) + "')"); |
1576 - | |
1577 - | LOG("Erasing the last element of duration_v (size before = '" + TO_STRING(duration_v.size()) + "')"); |
1578 - | duration_v.erase(duration_v.begin() + duration_v.size() - 1); |
1579 - | LOG("Erasing the last element of duration_v (size after = '" + TO_STRING(duration_v.size()) + "')"); |
1580 - | } |
1709 + | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1710 + | // Then just "forget" the last element. |
1711 + | LOG("Erasing the last element of ephRow_v (size before = '" + TO_STRING(ephRow_v.size()) + "')"); |
1712 + | ephRow_v.erase(ephRow_v.begin() + ephRow_v.size() - 1); |
1581 1713 | |
1582 - | // |
1583 - | // Determine the timely ordered sequence of indexes in v which will be used to tabulate the ephemeris data to be put into the MS table. |
1584 - | // |
1585 - | LOG("Prepare the time ordered sequence of indexes used to tabulate the ephemeris data to be written in the MS table."); |
1586 - | typedef pair<uint32_t, int64_t> atiIdxMStime_pair; |
1587 - | vector<atiIdxMStime_pair> atiIdxMStime_v; |
1714 + | LOG("Erasing the last element of duration_v (size before = '" + TO_STRING(duration_v.size()) + "')"); |
1715 + | duration_v.erase(duration_v.begin() + duration_v.size() - 1); |
1716 + | } |
1588 1717 | |
1589 - | uint32_t index = 0; |
1590 - | int64_t tMS = t0MS; |
1718 + | // Determine the timely ordered sequence of indexes in ephRow_v which will be used to tabulate the ephemeris data to be put into the MS table. |
1719 + | LOG("Prepare the time ordered sequence of indexes used to tabulate the ephemeris data to be written in the MS table."); |
1720 + | typedef pair<uint32_t, int64_t> atiIdxMStime_pair; |
1721 + | vector<atiIdxMStime_pair> atiIdxMStime_v; |
1591 1722 | |
1592 - | atiIdxMStime_v.push_back(atiIdxMStime_pair(index, tMS)); |
1593 - | LOG ("size of atiIdxMStime_v="+TO_STRING(atiIdxMStime_v.size())+", index = "+TO_STRING(index)+", tMS = "+TO_STRING(tMS)); |
1594 - | tMS += timeStepInNanoSecond; |
1723 + | uint32_t index = 0; |
1724 + | int64_t tMS = t0MS; |
1595 1725 | |
1596 - | int64_t start = v[index]->getTimeInterval().getStart().get(); |
1597 - | int64_t end = start + v[index]->getTimeInterval().getDuration().get(); |
1598 - | do { |
1599 - | if (tMS < end) { |
1600 - | atiIdxMStime_v.push_back(atiIdxMStime_pair(index, tMS)); |
1601 - | tMS += timeStepInNanoSecond; |
1602 - | LOG ("size of atiIdxMStime_v="+TO_STRING(atiIdxMStime_v.size())+", index = "+TO_STRING(index)+", tMS = "+TO_STRING(tMS)); |
1603 - | } |
1604 - | else { |
1605 - | index++; |
1606 - | end = v[index]->getTimeInterval().getStart().get() + v[index]->getTimeInterval().getDuration().get(); |
1607 - | } |
1726 + | atiIdxMStime_v.push_back(atiIdxMStime_pair(index, tMS)); |
1727 + | LOG ("size of atiIdxMStime_v="+TO_STRING(atiIdxMStime_v.size())+", index = "+TO_STRING(index)+", tMS = "+TO_STRING(tMS)); |
1728 + | tMS += timeStepInNanoSecond; |
1608 1729 | |
1609 - | } while (index < v.size()-1); |
1730 + | int64_t start = ephRow_v[index]->getTimeInterval().getStart().get(); |
1731 + | int64_t end = start + ephRow_v[index]->getTimeInterval().getDuration().get(); |
1732 + | do { |
1733 + | if (tMS < end) { |
1734 + | atiIdxMStime_v.push_back(atiIdxMStime_pair(index, tMS)); |
1735 + | tMS += timeStepInNanoSecond; |
1736 + | LOG ("size of atiIdxMStime_v="+TO_STRING(atiIdxMStime_v.size())+", index = "+TO_STRING(index)+", tMS = "+TO_STRING(tMS)); |
1737 + | } |
1738 + | else { |
1739 + | index++; |
1740 + | end = ephRow_v[index]->getTimeInterval().getStart().get() + ephRow_v[index]->getTimeInterval().getDuration().get(); |
1741 + | } |
1742 + | |
1743 + | } while (index < ephRow_v.size()-1); |
1610 1744 | |
1611 - | LOG("atiIdxMStime_v has " + TO_STRING(atiIdxMStime_v.size()) + " elements."); |
1745 + | LOG("atiIdxMStime_v has " + TO_STRING(atiIdxMStime_v.size()) + " elements."); |
1612 1746 | |
1613 - | // |
1614 - | // Prepare the coefficients which will be used for the tabulation. |
1747 + | // Prepare the coefficients which will be used for the tabulation. |
1615 1748 | |
1616 - | LOG("Prepare the coefficients which will be used for the tabulations."); |
1617 - | vector<vector<double> > raASDM_vv; |
1618 - | vector<double> raASDM_v; |
1619 - | vector<vector<double> > decASDM_vv; |
1620 - | vector<double> decASDM_v; |
1621 - | vector<vector<double> > distanceASDM_vv; |
1622 - | vector<double> distanceASDM_v; |
1623 - | vector<vector<double> > radVelASDM_vv; |
1624 - | vector<double> radVelASDM_v; |
1625 - | vector<double> empty_v; |
1626 - | vector<double> temp_v; |
1627 - | |
1628 - | cout.precision(10); |
1629 - | for (unsigned int i = 0; i < v.size(); i++) { |
1630 - | LOG_EPHEM("original " + TO_STRING (ArrayTime(v[i]->getTimeInterval().getStart().get()).getMJD())); |
1631 - | vector<vector<double> > temp_vv = v[i]->getDir(); |
1632 - | if (numPolyDirIsOne) { |
1633 - | raASDM_v.push_back(radian2degree(temp_vv[0][0])); |
1634 - | decASDM_v.push_back(radian2degree(temp_vv[0][1])); |
1635 - | LOG_EPHEM (" " + TO_STRING(raASDM_v.back()) + " " + TO_STRING(decASDM_v.back())); |
1636 - | } |
1637 - | else { |
1638 - | raASDM_vv.push_back(empty_v); |
1639 - | decASDM_vv.push_back(empty_v); |
1640 - | for (int j = 0; j < v[i]->getNumPolyDir(); j++) { |
1641 - | raASDM_vv.back().push_back(radian2degree(temp_vv[j][0])); |
1642 - | decASDM_vv.back().push_back(radian2degree(temp_vv[j][1])); ; |
1749 + | LOG("Prepare the coefficients which will be used for the tabulations."); |
1750 + | vector<vector<double> > raASDM_vv; |
1751 + | vector<double> raASDM_v; |
1752 + | vector<vector<double> > decASDM_vv; |
1753 + | vector<double> decASDM_v; |
1754 + | vector<vector<double> > distanceASDM_vv; |
1755 + | vector<double> distanceASDM_v; |
1756 + | vector<vector<double> > radVelASDM_vv; |
1757 + | vector<double> radVelASDM_v; |
1758 + | vector<double> empty_v; |
1759 + | vector<double> temp_v; |
1760 + | |
1761 + | cout.precision(10); |
1762 + | for (unsigned int i = 0; i < ephRow_v.size(); i++) { |
1763 + | LOG_EPHEM("original " + TO_STRING (ArrayTime(ephRow_v[i]->getTimeInterval().getStart().get()).getMJD())); |
1764 + | vector<vector<double> > temp_vv = ephRow_v[i]->getDir(); |
1765 + | if (numPolyDirIsOne) { |
1766 + | raASDM_v.push_back(radian2degree(temp_vv[0][0])); |
1767 + | decASDM_v.push_back(radian2degree(temp_vv[0][1])); |
1768 + | LOG_EPHEM (" " + TO_STRING(raASDM_v.back()) + " " + TO_STRING(decASDM_v.back())); |
1643 1769 | } |
1644 - | } |
1645 - | |
1646 - | temp_v = v[i]->getDistance(); |
1647 - | if (numPolyDistIsOne) { |
1648 - | distanceASDM_v.push_back(m2au(temp_v[0])); // AU |
1649 - | LOG_EPHEM (" " + TO_STRING(distanceASDM_v.back())); |
1650 - | } |
1651 - | else { |
1652 - | distanceASDM_vv.push_back(empty_v); |
1653 - | for (int j = 0; j < v[i]->getNumPolyDist(); j++) |
1654 - | distanceASDM_vv.back().push_back(m2au(temp_v[j])); // AU |
1655 - | } |
1656 - | |
1657 - | if (radVelExists) { |
1658 - | temp_v = v[i]->getRadVel(); |
1659 - | if (numPolyRadVelIsOne) { |
1660 - | radVelASDM_v.push_back(mpers2auperd(temp_v[0])); // AU/d |
1661 - | LOG_EPHEM(" " + TO_STRING(radVelASDM_v.back())); |
1770 + | else { |
1771 + | raASDM_vv.push_back(empty_v); |
1772 + | decASDM_vv.push_back(empty_v); |
1773 + | for (int j = 0; j < ephRow_v[i]->getNumPolyDir(); j++) { |
1774 + | raASDM_vv.back().push_back(radian2degree(temp_vv[j][0])); |
1775 + | decASDM_vv.back().push_back(radian2degree(temp_vv[j][1])); ; |
1776 + | } |
1777 + | } |
1778 + | |
1779 + | temp_v = ephRow_v[i]->getDistance(); |
1780 + | if (numPolyDistIsOne) { |
1781 + | distanceASDM_v.push_back(m2au(temp_v[0])); // AU |
1782 + | LOG_EPHEM (" " + TO_STRING(distanceASDM_v.back())); |
1662 1783 | } |
1663 1784 | else { |
1664 - | radVelASDM_vv.push_back(empty_v); |
1665 - | for (int j = 0; j < v[i]->getNumPolyRadVel(); j++) |
1666 - | radVelASDM_vv.back().push_back(mpers2auperd(temp_v[j])); // AU/d |
1667 - | } |
1668 - | } |
1669 - | LOG_EPHEM("\n"); |
1670 - | } |
1671 - | |
1672 - | |
1673 - | // Preparing the coefficients of piecewise polynomial of degree 1. |
1674 - | // |
1675 - | // The calculations below are done only once and will be useful for all the columns |
1676 - | // requiring the cubic spline interpolation. |
1677 - | // |
1678 - | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1679 - | if (numPolyDirIsOne) { |
1680 - | LOG("Compute the linear interpolation coefficients for RAD"); |
1681 - | linearInterpCoeff(v.size(), |
1682 - | time_v, |
1683 - | timeOffset_v, |
1684 - | raASDM_v, |
1685 - | raASDM_vv); |
1686 - | |
1687 - | LOG("Compute the linear interpolation coefficients for DEC"); |
1688 - | linearInterpCoeff(v.size(), |
1689 - | time_v, |
1690 - | timeOffset_v, |
1691 - | decASDM_v, |
1692 - | decASDM_vv); |
1693 - | } |
1694 - | |
1695 - | if (numPolyDistIsOne) { |
1696 - | LOG("Compute the linear interolation coefficients for Dist"); |
1697 - | linearInterpCoeff(v.size(), |
1698 - | time_v, |
1699 - | timeOffset_v, |
1700 - | distanceASDM_v, |
1701 - | distanceASDM_vv); |
1702 - | } |
1703 - | |
1704 - | if (radVelExists && numPolyRadVelIsOne) { |
1705 - | LOG("Compute the linear interpolation coefficients for RadVel"); |
1706 - | linearInterpCoeff(v.size(), |
1707 - | time_v, |
1708 - | timeOffset_v, |
1709 - | radVelASDM_v, |
1710 - | radVelASDM_vv); |
1711 - | } |
1712 - | } |
1713 - | // End of interpolating with piecewise polynomial of degree 1. |
1714 - | |
1715 - | |
1716 - | for (atiIdxMStime_pair atiIdxMStime: atiIdxMStime_v) { |
1717 - | // |
1718 - | // MJD |
1719 - | mjdMS_v.push_back(ArrayTime(atiIdxMStime.second).getMJD()); |
1720 - | LOG_EPHEM( "resampled " + TO_STRING(mjdMS_v.back())); |
1721 - | LOG("mjdMS_v -> "+TO_STRING(mjdMS_v.back())); |
1722 - | |
1723 - | double timeOrigin = 1.0e-09 * v[atiIdxMStime.first]->getTimeOrigin().get(); |
1724 - | double time = 1.0e-09 * atiIdxMStime.second; |
1725 - | |
1726 - | LOG("timeOrigin="+TO_STRING(timeOrigin)+", time="+TO_STRING(time)); |
1727 - | |
1728 - | // |
1729 - | // RA / DEC |
1730 - | LOG("Eval poly for RA"); |
1731 - | LOG("atiIdxMStime.first = " + TO_STRING(atiIdxMStime.first)); |
1732 - | raMS_v.push_back(evalPoly(raASDM_vv[atiIdxMStime.first].size(), |
1733 - | raASDM_vv[atiIdxMStime.first], |
1734 - | timeOrigin, |
1735 - | time)); |
1736 - | LOG_EPHEM(" " + TO_STRING(raMS_v.back())); |
1737 - | LOG("raMS_v -> "+TO_STRING(raMS_v.back())); |
1738 - | |
1739 - | LOG("Eval poly for DEC"); |
1740 - | decMS_v.push_back(evalPoly(decASDM_vv[atiIdxMStime.first].size(), |
1741 - | decASDM_vv[atiIdxMStime.first], |
1742 - | timeOrigin, |
1743 - | time)); |
1744 - | LOG_EPHEM(" " + TO_STRING(decMS_v.back())); |
1745 - | |
1746 - | // |
1747 - | // Distance |
1748 - | LOG("Eval poly for distance"); |
1749 - | distanceMS_v.push_back(evalPoly(distanceASDM_vv[atiIdxMStime.first].size(), |
1750 - | distanceASDM_vv[atiIdxMStime.first], |
1751 - | timeOrigin, |
1752 - | time)); |
1753 - | LOG_EPHEM(" " + TO_STRING(distanceMS_v.back())); |
1754 - | // |
1755 - | // Radvel |
1756 - | if (radVelExists) { |
1757 - | LOG("Eval poly for radvel"); |
1758 - | radVelMS_v.push_back(evalPoly(radVelASDM_vv[atiIdxMStime.first].size(), |
1759 - | radVelASDM_vv[atiIdxMStime.first], |
1760 - | timeOrigin, |
1761 - | time)); |
1762 - | LOG_EPHEM(" " + TO_STRING(radVelMS_v.back())); |
1763 - | } |
1764 - | LOG_EPHEM("\n"); |
1765 - | } |
1766 - | } // end if interpolate_ephemeris |
1767 - | else if (tabulate_ephemeris_polynomials && |
1768 - | v.size()==1 && |
1769 - | timeStepInNanoSecond > 0 |
1770 - | ) { |
1771 - | |
1772 - | |
1773 - | infostream << "The MS Ephemeris table for ephemerisId = '" << ephemerisId |
1774 - | << "' will be produced by tabulating the polynomials found in the columns 'dir', 'distance' and optionally 'radVel' with a timestep of '" |
1775 - | << timeStepInNanoSecond / 1.e9 / 86400. << "' days"; |
1776 - | |
1777 - | |
1778 - | dmjd = timeStepInNanoSecond / 1.e9 / 86400.; // to be written in the keyword DMJD |
1779 - | |
1780 - | for (unsigned int i = 0; i < v.size(); i++) { |
1781 - | // |
1782 - | // Calculate the grid of times where the polynomials will be tabulated. |
1783 - | // This grid contains all the times which : |
1784 - | // - are multiple of the tabulation time steps, |
1785 - | // - contained in the arraytime interval of validity of the current ASDM Ephemeris row. |
1786 - | // |
1787 - | // At this point times are expressed in nanoseconds. |
1788 - | // |
1789 - | |
1790 - | int64_t tstartASDM = v[i]->getTimeInterval().getStart().get(); |
1791 - | int64_t tendASDM = tstartASDM + v[i]->getTimeInterval().getDuration().get(); |
1792 - | int64_t q = tstartASDM / timeStepInNanoSecond; |
1793 - | int64_t r = tstartASDM % timeStepInNanoSecond; |
1794 - | int64_t tstartMS = tstartASDM; |
1795 - | |
1796 - | if ( r!= 0 ) { |
1797 - | q = q + 1; |
1798 - | tstartMS = q * timeStepInNanoSecond; |
1785 + | distanceASDM_vv.push_back(empty_v); |
1786 + | for (int j = 0; j < ephRow_v[i]->getNumPolyDist(); j++) |
1787 + | distanceASDM_vv.back().push_back(m2au(temp_v[j])); // AU |
1788 + | } |
1789 + | |
1790 + | if (radVelExists) { |
1791 + | temp_v = ephRow_v[i]->getRadVel(); |
1792 + | if (numPolyRadVelIsOne) { |
1793 + | radVelASDM_v.push_back(mpers2auperd(temp_v[0])); // AU/d |
1794 + | LOG_EPHEM(" " + TO_STRING(radVelASDM_v.back())); |
1795 + | } |
1796 + | else { |
1797 + | radVelASDM_vv.push_back(empty_v); |
1798 + | for (int j = 0; j < ephRow_v[i]->getNumPolyRadVel(); j++) |
1799 + | radVelASDM_vv.back().push_back(mpers2auperd(temp_v[j])); // AU/d |
1800 + | } |
1801 + | } |
1802 + | LOG_EPHEM("\n"); |
1799 1803 | } |
1800 1804 | |
1801 - | vector<int64_t> tabulation_time_v; |
1802 - | int64_t t = tstartMS - timeStepInNanoSecond; // One extra timestep before the beginning of the interval of validity. |
1803 - | do { |
1804 - | tabulation_time_v.push_back(t); |
1805 - | t += timeStepInNanoSecond; |
1806 - | } |
1807 - | while (t <= tendASDM); |
1808 - | tabulation_time_v.push_back(t); // One extra timestep after the end of the interval of validity. |
1809 - | |
1810 1805 | |
1806 + | // Preparing the coefficients of piecewise polynomial of degree 1. |
1811 1807 | // |
1812 - | // Let's tabulate the MS Ephemeris column for each tabulation time, in s. |
1808 + | // The calculations below are done only once and will be useful for all the columns |
1809 + | // requiring the cubic spline interpolation. |
1813 1810 | // |
1814 - | double timeOrigin = v[i]->getTimeOrigin().get() * 1.0e-09 / 86400. ; // It appeared that time must be expressed in "days" !! |
1815 - | |
1816 - | // Convert from radians to degrees. |
1817 - | const vector<vector<double> >& dir_v = v[i]->getDir(); |
1818 - | vector<double> ra_coeff_v; |
1819 - | vector<double> dec_coeff_v; |
1820 - | for (unsigned int idir = 0; idir < dir_v.size(); idir++) { |
1821 - | ra_coeff_v.push_back(dir_v[idir][0]); |
1822 - | dec_coeff_v.push_back(dir_v[idir][1]); |
1811 + | if (numPolyDirIsOne || numPolyDistIsOne || (radVelExists && numPolyRadVelIsOne)) { |
1812 + | if (numPolyDirIsOne) { |
1813 + | LOG("Compute the linear interpolation coefficients for RAD"); |
1814 + | linearInterpCoeff(ephRow_v.size(), |
1815 + | time_v, |
1816 + | timeOffset_v, |
1817 + | raASDM_v, |
1818 + | raASDM_vv); |
1819 + | |
1820 + | LOG("Compute the linear interpolation coefficients for DEC"); |
1821 + | linearInterpCoeff(ephRow_v.size(), |
1822 + | time_v, |
1823 + | timeOffset_v, |
1824 + | decASDM_v, |
1825 + | decASDM_vv); |
1826 + | } |
1827 + | |
1828 + | if (numPolyDistIsOne) { |
1829 + | LOG("Compute the linear interolation coefficients for Dist"); |
1830 + | linearInterpCoeff(ephRow_v.size(), |
1831 + | time_v, |
1832 + | timeOffset_v, |
1833 + | distanceASDM_v, |
1834 + | distanceASDM_vv); |
1835 + | } |
1836 + | |
1837 + | if (radVelExists && numPolyRadVelIsOne) { |
1838 + | LOG("Compute the linear interpolation coefficients for RadVel"); |
1839 + | linearInterpCoeff(ephRow_v.size(), |
1840 + | time_v, |
1841 + | timeOffset_v, |
1842 + | radVelASDM_v, |
1843 + | radVelASDM_vv); |
1844 + | } |
1823 1845 | } |
1824 - | |
1825 - | LOG(" There are " + TO_STRING(ra_coeff_v.size()) + " RA coeffs and " + TO_STRING(dec_coeff_v.size()) + " DEC coeffs"); |
1826 - | |
1827 - | // Convert radian to degree. |
1828 - | std::transform(ra_coeff_v.begin(), ra_coeff_v.end(), ra_coeff_v.begin(), radian2degree); |
1829 - | std::transform(dec_coeff_v.begin(), dec_coeff_v.end(), dec_coeff_v.begin(), radian2degree); |
1830 - | |
1831 - | // Convert from m to AU. |
1832 - | vector<double> distance_coeff_v = v[i]->getDistance(); |
1833 - | std::transform(distance_coeff_v.begin(), distance_coeff_v.end(), distance_coeff_v.begin(), m2au); |
1846 + | // End of interpolating with piecewise polynomial of degree 1. |
1834 1847 | |
1835 - | |
1836 - | vector<double> radvel_coeff_v; |
1837 - | if (v[i]->isRadVelExists()) { |
1838 - | radvel_coeff_v = v[i]->getRadVel(); |
1839 - | // Convert from m per s to AU per day. |
1840 - | std::transform(radvel_coeff_v.begin(), radvel_coeff_v.end(), radvel_coeff_v.begin(), mpers2auperd); |
1841 - | } else { |
1842 - | // Let's derivate the polynomial in distance. |
1843 - | for (unsigned int iDistance=1; iDistance < distance_coeff_v.size(); iDistance++) |
1844 - | radvel_coeff_v.push_back(distance_coeff_v[iDistance]*iDistance); |
1845 - | } |
1846 1848 | |
1847 - | // And proceed... |
1848 - | LOG ("There will be " + TO_STRING(tabulation_time_v.size()) + " time steps used to tabulate the polynomials."); |
1849 - | for (unsigned int itab = 0; itab < tabulation_time_v.size(); itab++) { |
1850 - | double tabulation_time = tabulation_time_v[itab] * 1.0e-09 / 86400.0 ; // It appeared that times should be expressed in "day" !! |
1851 - | |
1852 - | |
1849 + | for (atiIdxMStime_pair atiIdxMStime: atiIdxMStime_v) { |
1853 1850 | // |
1854 1851 | // MJD |
1855 - | // |
1856 - | mjdMS_v.push_back(ArrayTime(tabulation_time_v[itab]).getMJD()); |
1852 + | mjdMS_v.push_back(ArrayTime(atiIdxMStime.second).getMJD()); |
1853 + | LOG_EPHEM( "resampled " + TO_STRING(mjdMS_v.back())); |
1854 + | LOG("mjdMS_v -> "+TO_STRING(mjdMS_v.back())); |
1855 + | |
1856 + | double timeOrigin = 1.0e-09 * ephRow_v[atiIdxMStime.first]->getTimeOrigin().get(); |
1857 + | double time = 1.0e-09 * atiIdxMStime.second; |
1858 + | |
1859 + | LOG("timeOrigin="+TO_STRING(timeOrigin)+", time="+TO_STRING(time)); |
1857 1860 | |
1858 - | // |
1859 1861 | // RA / DEC |
1860 - | // |
1861 - | raMS_v.push_back(evalPoly(ra_coeff_v.size(), |
1862 - | ra_coeff_v, |
1862 + | LOG("Eval poly for RA"); |
1863 + | LOG("atiIdxMStime.first = " + TO_STRING(atiIdxMStime.first)); |
1864 + | raMS_v.push_back(evalPoly(raASDM_vv[atiIdxMStime.first].size(), |
1865 + | raASDM_vv[atiIdxMStime.first], |
1863 1866 | timeOrigin, |
1864 - | tabulation_time)); |
1867 + | time)); |
1868 + | LOG_EPHEM(" " + TO_STRING(raMS_v.back())); |
1869 + | LOG("raMS_v -> "+TO_STRING(raMS_v.back())); |
1865 1870 | |
1866 - | decMS_v.push_back(evalPoly(dec_coeff_v.size(), |
1867 - | dec_coeff_v, |
1871 + | LOG("Eval poly for DEC"); |
1872 + | decMS_v.push_back(evalPoly(decASDM_vv[atiIdxMStime.first].size(), |
1873 + | decASDM_vv[atiIdxMStime.first], |
1868 1874 | timeOrigin, |
1869 - | tabulation_time)); |
1875 + | time)); |
1876 + | LOG_EPHEM(" " + TO_STRING(decMS_v.back())); |
1870 1877 | |
1871 - | |
1872 - | // |
1873 - | // DISTANCE |
1874 - | // |
1875 - | distanceMS_v.push_back(evalPoly(distance_coeff_v.size(), |
1876 - | distance_coeff_v, |
1878 + | // Distance |
1879 + | LOG("Eval poly for distance"); |
1880 + | distanceMS_v.push_back(evalPoly(distanceASDM_vv[atiIdxMStime.first].size(), |
1881 + | distanceASDM_vv[atiIdxMStime.first], |
1877 1882 | timeOrigin, |
1878 - | tabulation_time)); |
1879 - | |
1880 - | // |
1881 - | // RADVEL |
1882 - | // |
1883 - | if (radVelExists) { |
1884 - | radVelMS_v.push_back(evalPoly(radvel_coeff_v.size(), |
1885 - | radvel_coeff_v, |
1883 + | time)); |
1884 + | LOG_EPHEM(" " + TO_STRING(distanceMS_v.back())); |
1885 + | |
1886 + | // Radvel |
1887 + | if (radVelExists) { |
1888 + | LOG("Eval poly for radvel"); |
1889 + | radVelMS_v.push_back(evalPoly(radVelASDM_vv[atiIdxMStime.first].size(), |
1890 + | radVelASDM_vv[atiIdxMStime.first], |
1886 1891 | timeOrigin, |
1887 - | tabulation_time)); |
1892 + | time)); |
1893 + | LOG_EPHEM(" " + TO_STRING(radVelMS_v.back())); |
1888 1894 | } |
1895 + | LOG_EPHEM("\n"); |
1889 1896 | } |
1890 - | } |
1891 - | } |
1892 - | else { |
1893 - | // Just copy ephemeris without any interpolation. Just adapt the units. |
1894 - | for(const EphemerisRow *eR_p: v) { |
1895 - | mjdMS_v.push_back(eR_p->getTimeInterval().getMidPoint().getMJD()); // MJD |
1896 - | vector<vector<double> > dir = eR_p->getDir(); |
1897 - | raMS_v.push_back(radian2degree(dir[0][0])); // deg |
1898 - | decMS_v.push_back(radian2degree(dir[0][1])); // deg |
1899 - | distanceMS_v.push_back(m2au(eR_p->getDistance()[0])); // AU |
1900 - | if (radVelExists) |
1901 - | radVelMS_v.push_back(mpers2auperd(eR_p->getRadVel()[0])); // AU/d |
1902 - | } |
1903 - | } // end not if interpolate_ephemeris |
1897 + | } // end of the general interpolation case |
1898 + | } // end of the non-single row case, everything is ready to push to the MS now |
1904 1899 | |
1905 1900 | // |
1906 1901 | // Record the starting time + one time step so that Field can use it if needed. |
1907 1902 | // |
1908 1903 | ephemStartTime_m[ephemerisId] = ArrayTime(mjdMS_v[1]).get()*1.e-09; |
1909 1904 | |
1910 1905 | // Now the data are ready to be written to the MS Ephemeris table. |
1911 1906 | // Let's proceed, using Slicers. |
1912 - | // |
1913 1907 | |
1914 1908 | unsigned int numRows = raMS_v.size(); |
1915 1909 | Slicer slicer(IPosition(1, 0), |
1916 1910 | IPosition(1, numRows-1), |
1917 1911 | Slicer::endIsLast); |
1918 1912 | |
1919 1913 | for (map<AtmPhaseCorrection, ASDM2MSFiller*>::iterator iter = msFillers.begin(); |
1920 1914 | iter != msFillers.end(); |
1921 1915 | ++iter) { |
1922 1916 | Table * table_p = apc2EphemTable_m[iter->first]; |
4311 4305 | |
4312 4306 | uint64_t bdfSliceSizeInMb = 0; // The default size of the BDF slice hold in memory. |
4313 4307 | |
4314 4308 | bool mute = false; |
4315 4309 | |
4316 4310 | bool doparallel = false; |
4317 4311 | |
4318 4312 | bool ac_xc_per_timestamp = false; // for the time being the option is 'preserve the old order' |
4319 4313 | |
4320 4314 | bool interpolate_ephemeris = false ; |
4321 - | bool tabulate_ephemeris_polynomials = false; |
4322 4315 | double polyephem_tabtimestep = 0.001; |
4323 4316 | |
4324 4317 | // Process command line options and parameters. |
4325 4318 | po::variables_map vm; |
4326 4319 | |
4327 4320 | try { |
4328 4321 | |
4329 4322 | |
4330 4323 | // Declare the supported options. |
4331 4324 | po::options_description generic("Converts an ASDM dataset into a CASA measurement set.\n" |
4355 4348 | ("no-caldevice", "The CalDevice table will be ignored.") |
4356 4349 | ("no-ephemeris", "The ephemeris table will be ignored.") |
4357 4350 | ("no-syspower", "the SysPower table will be ignored.") |
4358 4351 | ("no-pointing", "The Pointing table will be ignored.") |
4359 4352 | ("check-row-uniqueness", "The row uniqueness constraint will be checked in the tables where it's defined") |
4360 4353 | ("bdf-slice-size", po::value<uint64_t>(&bdfSliceSizeInMb)->default_value(500), "The maximum amount of memory expressed as an integer in units of megabytes (1024*1024) allocated for BDF data. The default is 500 (megabytes)") |
4361 4354 | //("parallel", "run with multithreading mode.") |
4362 4355 | ("lazy", "defers the production of the observational data in the MS Main table (DATA column) - Purely experimental, don't use in production !") |
4363 4356 | ("with-pointing-correction", "add (ASDM::Pointing::encoder - ASDM::Pointing::pointingDirection) to the value to be written in MS::Pointing::direction - (related with JIRA tickets CSV-2878 and ICT-1532))") |
4364 4357 | ("ac-xc-per-timestamp", po::value<string>()->default_value("no"), "if set to yes, then the filler writes in that order autocorrelations and cross correlations rows for one given data description and timestamp. Otherwise auto correlations data are grouped for a sequence of time stamps and then come the cross correlations data for the same sequence of timestamps.") |
4365 - | ("polyephem-tabtimestep", po::value<double>(&polyephem_tabtimestep), "Defines the time step used to tabulate the polynomials found in the columns 'dir', 'distance' and optionally 'radVel' of the ASDM Ephemeris table. The unit to express the time step is the day and the default value is 0.001. If 'radvel' is not present then the radial velocity will be obtained by tabulating the derivative the polynomial found in 'distance'.") |
4358 + | ("polyephem-tabtimestep", po::value<double>(&polyephem_tabtimestep)->default_value(0.001), "Defines the time step used to tabulate the polynomials found in the columns 'dir', 'distance' and optionally 'radVel' of the ASDM Ephemeris table. The unit to express the time step is the day and the default value is 0.001. If 'radvel' is not present then the radial velocity will be obtained by tabulating the derivative the polynomial found in 'distance'.") |
4366 4359 | ("interpolate-ephemeris", po::value<string>()->default_value("no"), "if set to 'yes' then the filler will resample the sequence of times found in the ASDM Ephemeris table into an evenly spaced sequence of times on which the ephemeris paarameters will obtained by an interpolation of degree 1. Otherwise (!= 'yes') the ephemeris parameters will be copies of what's in the ASDM Ephemeris table on the same sequence of times"); |
4367 4360 | |
4368 4361 | // Hidden options, will be allowed both on command line and |
4369 4362 | // in config file, but will not be shown to the user. |
4370 4363 | po::options_description hidden("Hidden options"); |
4371 4364 | hidden.add_options() |
4372 4365 | ("asdm-directory", po::value< string >(), "asdm directory") |
4373 4366 | ; |
4374 4367 | hidden.add_options() |
4375 4368 | ("ms-directory-prefix", po::value< string >(), "ms directory prefix") |
4610 4603 | infostream.str(""); |
4611 4604 | infostream << "The lazy filler can not be used with any input data selection (icm, isrt, its)." << endl; |
4612 4605 | infostream << "The non-lazy version of the filler will be used." << endl; |
4613 4606 | // an option is being ignored, warn the user even if not verbose |
4614 4607 | warning(infostream.str()); |
4615 4608 | } |
4616 4609 | |
4617 4610 | // Do we consider another order than ac_xc_per_timestamp ? |
4618 4611 | ac_xc_per_timestamp = boost::algorithm::to_lower_copy(vm["ac-xc-per-timestamp"].as<string>()) == "yes"; |
4619 4612 | |
4620 - | // Do we want to tabulate polynomial present in the ephemeris table ? |
4621 - | tabulate_ephemeris_polynomials = (vm.count("polyephem-tabtimestep") != 0); |
4622 - | |
4623 - | if (tabulate_ephemeris_polynomials) { |
4624 - | // If we tabluate then ignore all the other options about ephemeris |
4625 - | // infostream.str(); |
4626 - | // infostream << "The MS Ephemeris table(s) will be produced by tabulating the polynomials found in the columns 'dir', 'distance' and optionally 'radVel' with a timestep of '" |
4627 - | // << polyephem_tabtimestep << "' day, i.e. '"<< ((uint64_t) (polyephem_tabtimestep * 86400 * 1.e09)) <<"' nanoseconds."; |
4628 - | // info(infostream.str()); |
4629 - | } else { |
4630 - | // Do we want interpolate the values found in the ASDM Ephemeris table or not ? |
4631 - | interpolate_ephemeris = boost::algorithm::to_lower_copy(vm["interpolate-ephemeris"].as<string>()) == "yes"; |
4632 - | infostream.str(""); |
4633 - | if (interpolate_ephemeris) { |
4634 - | infostream << "the MS Ephemeris table(s) will be produced by interpolation of the values present in the ASDM Ephemeris table on a resampled time grid."; |
4635 - | } |
4636 - | else { |
4637 - | infostream << "the MS Ephemeris tables(s) will be produced by simple copies of the values found in the ASDM Ephemeris table with just units conversion."; |
4638 - | } |
4639 - | info(infostream.str()); |
4640 - | } |
4613 + | // Do we want interpolate the values found in the ASDM Ephemeris table or not ? |
4614 + | interpolate_ephemeris = boost::algorithm::to_lower_copy(vm["interpolate-ephemeris"].as<string>()) == "yes"; |
4641 4615 | } |
4642 4616 | catch (std::exception& e) { |
4643 4617 | errstream.str(""); |
4644 4618 | errstream << e.what(); |
4645 4619 | error(errstream.str()); |
4646 4620 | } |
4647 4621 | // this just a dummy number for now (innthread may |
4648 4622 | // come from user input in future..) |
4649 4623 | // Also setting environment variable, OMP_NUM_THREADS=1 |
4650 4624 | // one can excute multiwrite part in a single thread. |
5523 5497 | errstream.str(""); |
5524 5498 | errstream << e.what(); |
5525 5499 | error(errstream.str()); |
5526 5500 | } |
5527 5501 | |
5528 5502 | // Process the Ephemeris table. |
5529 5503 | // |
5530 5504 | // Create and fill the MS ephemeris table(s) with a time interpolation time step set to 86400000000 nanoseconds ( 1/1000 day). |
5531 5505 | if (processEphemeris) { |
5532 5506 | uint64_t timeStepInNanoSeconds = polyephem_tabtimestep * 86400 * 1.e09; |
5533 - | fillEphemeris(ds, timeStepInNanoSeconds, interpolate_ephemeris, tabulate_ephemeris_polynomials, telescopeName); |
5507 + | fillEphemeris(ds, timeStepInNanoSeconds, interpolate_ephemeris, telescopeName); |
5534 5508 | } |
5535 5509 | |
5536 5510 | // Process the Field table. |
5537 5511 | // Now it respects the degree of the polynomials but it ignores the ephemerisId. |
5538 5512 | // The ephemerisId will be processed during the call to fillEphemeris. |
5539 5513 | // |
5540 5514 | fillField(ds, processEphemeris); |
5541 5515 | |
5542 5516 | // Process the FlagCmd table. |
5543 5517 | // |