30 #include "Eigen/Sparse"
49 _astrometryModel(astrometryModel),
50 _refractionCoefficient(0),
53 _nParRefrac(_associations->getNFilters()),
64 _referenceColor += i->color;
69 _referenceColor /= double(
count);
72 LOGLS_INFO(
_log,
"Reference Color: " << _referenceColor <<
" sig " << _sigCol);
80 Point AstrometryFit::transformFittedStar(FittedStar
const &fittedStar, AstrometryTransform
const &sky2TP,
81 Point
const &refractionVector,
double refractionCoeff,
83 Point fittedStarInTP = sky2TP.apply(fittedStar);
84 if (fittedStar.mightMove) {
85 fittedStarInTP.x += fittedStar.pmx * mjd;
86 fittedStarInTP.y += fittedStar.pmy * mjd;
91 double color = fittedStar.color - _referenceColor;
92 fittedStarInTP.x += refractionVector.x * color * refractionCoeff;
93 fittedStarInTP.y += refractionVector.y * color * refractionCoeff;
94 return fittedStarInTP;
100 static void tweakAstromMeasurementErrors(FatPoint &P, MeasuredStar
const &Ms,
double error) {
101 static bool called =
false;
102 static double increment = 0;
113 void AstrometryFit::leastSquareDerivativesMeasurement(CcdImage
const &ccdImage, TripletList &tripletList,
114 Eigen::VectorXd &fullGrad,
115 MeasuredStarList
const *msList)
const {
125 if (msList) assert(&(msList->front()->getCcdImage()) == &ccdImage);
128 const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
130 std::size_t npar_mapping = (_fittingDistortions) ? mapping->getNpar() : 0;
132 std::size_t npar_refrac = (_fittingRefrac) ? 1 : 0;
134 std::size_t npar_tot = npar_mapping + npar_pos + npar_refrac + npar_pm;
137 if (npar_tot == 0)
return;
139 if (_fittingDistortions) mapping->getMappingIndices(indices);
142 double mjd = ccdImage.getMjd() - _JDRef;
144 Point refractionVector = ccdImage.getRefractionVector();
146 auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
148 AstrometryTransformLinear dypdy;
152 Eigen::Matrix2d transW(2, 2);
153 Eigen::Matrix2d alpha(2, 2);
154 Eigen::VectorXd grad(npar_tot);
156 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
157 const MeasuredStarList &catalog = (msList) ? *msList : ccdImage.getCatalogForFit();
159 for (
auto &i : catalog) {
160 const MeasuredStar &ms = *i;
161 if (!ms.isValid())
continue;
164 tweakAstromMeasurementErrors(inPos, ms, _posError);
168 if (_fittingDistortions)
169 mapping->computeTransformAndDerivatives(inPos, outPos, H);
171 mapping->transformPosAndErrors(inPos, outPos);
174 double det = outPos.vx * outPos.vy -
std::pow(outPos.vxy, 2);
175 if (
det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
176 LOGLS_WARN(
_log,
"Inconsistent measurement errors: drop measurement at "
177 << Point(ms) <<
" in image " << ccdImage.getName());
180 transW(0, 0) = outPos.vy /
det;
181 transW(1, 1) = outPos.vx /
det;
182 transW(0, 1) = transW(1, 0) = -outPos.vxy /
det;
185 alpha(0, 0) =
sqrt(transW(0, 0));
187 alpha(1, 0) = transW(0, 1) / alpha(0, 0);
190 alpha(1, 1) = 1. /
sqrt(
det * transW(0, 0));
195 Point fittedStarInTP =
196 transformFittedStar(*
fs, *sky2TP, refractionVector, _refractionCoefficient, mjd);
201 sky2TP->computeDerivative(*
fs, dypdy, 1e-3);
204 H(npar_mapping, 0) = -dypdy.A11();
205 H(npar_mapping + 1, 0) = -dypdy.A12();
206 H(npar_mapping, 1) = -dypdy.A21();
207 H(npar_mapping + 1, 1) = -dypdy.A22();
208 indices[npar_mapping] =
fs->getIndexInMatrix();
209 indices.at(npar_mapping + 1) =
fs->getIndexInMatrix() + 1;
214 if (_fittingPM &&
fs->mightMove) {
216 H(ipar + 1, 1) = -mjd;
217 indices[ipar] =
fs->getIndexInMatrix() + 2;
218 indices[ipar + 1] =
fs->getIndexInMatrix() + 3;
221 if (_fittingRefrac) {
224 double color =
fs->color - _referenceColor;
226 H(ipar, 0) = -refractionVector.x * color;
227 H(ipar, 1) = -refractionVector.y * color;
228 indices[ipar] = _refracPosInMatrix;
233 Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
241 for (
std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
243 double val = halpha(ipar, ic);
244 if (
val == 0)
continue;
245 tripletList.addTriplet(indices[ipar], kTriplets + ic,
val);
247 fullGrad(indices[ipar]) += grad(ipar);
251 tripletList.setNextFreeIndex(kTriplets);
254 void AstrometryFit::leastSquareDerivativesReference(FittedStarList
const &fittedStarList,
255 TripletList &tripletList,
256 Eigen::VectorXd &fullGrad)
const {
265 if (!_fittingPos)
return;
269 Eigen::Matrix2d W(2, 2);
270 Eigen::Matrix2d alpha(2, 2);
271 Eigen::Matrix2d H(2, 2), halpha(2, 2), HW(2, 2);
272 AstrometryTransformLinear der;
273 Eigen::Vector2d res, grad;
274 Eigen::Index indices[2 +
NPAR_PM];
275 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
282 TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
283 for (
auto const &i : fittedStarList) {
284 const FittedStar &
fs = *i;
285 const RefStar *rs =
fs.getRefStar();
286 if (rs ==
nullptr)
continue;
287 proj.setTangentPoint(
fs);
290 proj.transformPosAndErrors(*rs, rsProj);
292 proj.computeDerivative(
fs, der, 1e-4);
294 H(0, 0) = -der.A11();
295 H(1, 0) = -der.A12();
296 H(0, 1) = -der.A21();
297 H(1, 1) = -der.A22();
299 double det = rsProj.vx * rsProj.vy -
std::pow(rsProj.vxy, 2);
300 if (rsProj.vx <= 0 || rsProj.vy <= 0 ||
det <= 0) {
301 LOGLS_WARN(
_log,
"RefStar error matrix not positive definite for: " << *rs);
304 W(0, 0) = rsProj.vy /
det;
305 W(0, 1) = W(1, 0) = -rsProj.vxy /
det;
306 W(1, 1) = rsProj.vx /
det;
309 alpha(0, 0) =
sqrt(W(0, 0));
311 alpha(1, 0) = W(0, 1) / alpha(0, 0);
312 alpha(1, 1) = 1. /
sqrt(
det * W(0, 0));
314 indices[0] =
fs.getIndexInMatrix();
315 indices[1] =
fs.getIndexInMatrix() + 1;
316 unsigned npar_tot = 2;
332 for (
std::size_t ipar = 0; ipar < npar_tot; ++ipar) {
333 for (
unsigned ic = 0; ic < 2; ++ic) {
334 double val = halpha(ipar, ic);
335 if (
val == 0)
continue;
336 tripletList.addTriplet(indices[ipar], kTriplets + ic,
val);
338 fullGrad(indices[ipar]) += grad(ipar);
342 tripletList.setNextFreeIndex(kTriplets);
345 void AstrometryFit::accumulateStatImage(CcdImage
const &ccdImage, Chi2Accumulator &accum)
const {
352 const AstrometryMapping *mapping = _astrometryModel->getMapping(ccdImage);
354 double mjd = ccdImage.getMjd() - _JDRef;
356 Point refractionVector = ccdImage.getRefractionVector();
358 auto sky2TP = _astrometryModel->getSkyToTangentPlane(ccdImage);
360 Eigen::Matrix2Xd transW(2, 2);
362 auto &catalog = ccdImage.getCatalogForFit();
363 for (
auto const &ms : catalog) {
364 if (!ms->isValid())
continue;
366 FatPoint inPos = *ms;
367 tweakAstromMeasurementErrors(inPos, *ms, _posError);
371 mapping->transformPosAndErrors(inPos, outPos);
372 double det = outPos.vx * outPos.vy -
std::pow(outPos.vxy, 2);
373 if (
det <= 0 || outPos.vx <= 0 || outPos.vy <= 0) {
374 LOGLS_WARN(
_log,
" Inconsistent measurement errors :drop measurement at "
375 << Point(*ms) <<
" in image " << ccdImage.getName());
378 transW(0, 0) = outPos.vy /
det;
379 transW(1, 1) = outPos.vx /
det;
380 transW(0, 1) = transW(1, 0) = -outPos.vxy /
det;
383 Point fittedStarInTP =
384 transformFittedStar(*
fs, *sky2TP, refractionVector, _refractionCoefficient, mjd);
386 Eigen::Vector2d res(fittedStarInTP.x - outPos.x, fittedStarInTP.y - outPos.y);
387 double chi2Val = res.transpose() * transW * res;
389 accum.addEntry(chi2Val, 2, ms);
393 void AstrometryFit::accumulateStatImageList(
CcdImageList const &ccdImageList, Chi2Accumulator &accum)
const {
394 for (
auto const &ccdImage : ccdImageList) {
395 accumulateStatImage(*ccdImage, accum);
399 void AstrometryFit::accumulateStatRefStars(Chi2Accumulator &accum)
const {
407 FittedStarList &fittedStarList =
_associations->fittedStarList;
408 TanRaDecToPixel proj(AstrometryTransformLinear(), Point(0., 0.));
409 for (
auto const &
fs : fittedStarList) {
410 const RefStar *rs =
fs->getRefStar();
411 if (rs ==
nullptr)
continue;
412 proj.setTangentPoint(*
fs);
415 proj.transformPosAndErrors(*rs, rsProj);
417 double rx = rsProj.x;
418 double ry = rsProj.y;
419 double det = rsProj.vx * rsProj.vy -
std::pow(rsProj.vxy, 2);
420 double wxx = rsProj.vy /
det;
421 double wyy = rsProj.vx /
det;
422 double wxy = -rsProj.vxy /
det;
423 double chi2 = wxx *
std::pow(rx, 2) + 2 * wxy * rx * ry + wyy *
std::pow(ry, 2);
424 accum.addEntry(chi2, 2,
fs);
431 void AstrometryFit::getIndicesOfMeasuredStar(MeasuredStar
const &measuredStar,
433 if (_fittingDistortions) {
434 const AstrometryMapping *mapping = _astrometryModel->getMapping(measuredStar.getCcdImage());
435 mapping->getMappingIndices(indices);
438 Eigen::Index fsIndex =
fs->getIndexInMatrix();
454 _fittingDistortions = (
_whatToFit.
find(
"Distortions") != std::string::npos);
455 _fittingPos = (
_whatToFit.
find(
"Positions") != std::string::npos);
456 _fittingRefrac = (
_whatToFit.
find(
"Refrac") != std::string::npos);
457 if (_sigCol == 0 && _fittingRefrac) {
459 "Cannot fit refraction coefficients without a color lever arm. Ignoring refraction.");
460 _fittingRefrac =
false;
465 _nParDistortions = 0;
466 if (_fittingDistortions) _nParDistortions = _astrometryModel->assignIndices(
_whatToFit, 0);
471 for (
auto &fittedStar : fittedStarList) {
476 fittedStar->setIndexInMatrix(ipar);
478 if ((_fittingPM)&fittedStar->mightMove) ipar +=
NPAR_PM;
481 _nParPositions = ipar - _nParDistortions;
482 if (_fittingRefrac) {
483 _refracPosInMatrix = ipar;
492 "AstrometryFit::offsetParams : the provided vector length is not compatible with "
493 "the current whatToFit setting");
494 if (_fittingDistortions) _astrometryModel->offsetParams(delta);
498 for (
auto const &i : fittedStarList) {
503 Eigen::Index index =
fs.getIndexInMatrix();
504 fs.x += delta(index);
505 fs.y += delta(index + 1);
506 if ((_fittingPM)&
fs.mightMove) {
507 fs.pmx += delta(index + 2);
508 fs.pmy += delta(index + 3);
512 if (_fittingRefrac) {
513 _refractionCoefficient += delta(_refracPosInMatrix);
519 const char *what2fit[] = {
"Positions",
522 "Positions Distortions",
524 "Distortions Refrac",
525 "Positions Distortions Refrac"};
527 const char *what2fit[] = {
"Positions",
"Distortions",
"Positions Distortions"};
529 for (
unsigned k = 0; k <
sizeof(what2fit) /
sizeof(what2fit[0]); ++k) {
536 jacobian.setFromTriplets(tripletList.
begin(), tripletList.
end());
546 ofile <<
"#id" << separator <<
"xccd" << separator <<
"yccd " << separator;
547 ofile <<
"rx" << separator <<
"ry" << separator;
548 ofile <<
"xtp" << separator <<
"ytp" << separator;
549 ofile <<
"mag" << separator <<
"mjd" << separator;
550 ofile <<
"xErr" << separator <<
"yErr" << separator <<
"xyCov" << separator;
551 ofile <<
"xtpi" << separator <<
"ytpi" << separator;
552 ofile <<
"rxi" << separator <<
"ryi" << separator;
553 ofile <<
"color" << separator <<
"fsindex" << separator;
554 ofile <<
"ra" << separator <<
"dec" << separator;
555 ofile <<
"chi2" << separator <<
"nm" << separator;
556 ofile <<
"chip" << separator <<
"visit" <<
std::endl;
558 ofile <<
"#id in source catalog" << separator <<
"coordinates in CCD (pixels)" << separator << separator;
559 ofile <<
"residual on TP (degrees)" << separator << separator;
560 ofile <<
"transformed coordinate in TP (degrees)" << separator << separator;
561 ofile <<
"rough magnitude" << separator <<
"Modified Julian Date of the measurement" << separator;
562 ofile <<
"transformed measurement uncertainty (degrees)" << separator << separator << separator;
563 ofile <<
"as-read position on TP (degrees)" << separator << separator;
564 ofile <<
"as-read residual on TP (degrees)" << separator << separator;
565 ofile <<
"currently unused" << separator <<
"unique index of the fittedStar" << separator;
566 ofile <<
"on sky position of fittedStar" << separator << separator;
567 ofile <<
"contribution to Chi2 (2D dofs)" << separator <<
"number of measurements of this fittedStar"
569 ofile <<
"chip id" << separator <<
"visit id" <<
std::endl;
572 for (
auto const &ccdImage : ccdImageList) {
575 const auto readTransform = ccdImage->getReadWcs();
576 const Point &refractionVector = ccdImage->getRefractionVector();
577 double mjd = ccdImage->getMjd() - _JDRef;
578 for (
auto const &ms : cat) {
579 if (!ms->isValid())
continue;
582 tweakAstromMeasurementErrors(inPos, *ms, _posError);
584 auto sky2TP = _astrometryModel->getSkyToTangentPlane(*ccdImage);
586 compose(*sky2TP, *readTransform);
587 FatPoint inputTpPos = readPixToTangentPlane->apply(inPos);
590 Point fittedStarInTP =
591 transformFittedStar(*
fs, *sky2TP, refractionVector, _refractionCoefficient, mjd);
592 Point res = tpPos - fittedStarInTP;
593 Point inputRes = inputTpPos - fittedStarInTP;
595 double wxx = tpPos.
vy /
det;
596 double wyy = tpPos.
vx /
det;
597 double wxy = -tpPos.
vxy /
det;
598 double chi2 = wxx * res.
x * res.
x + wyy * res.
y * res.
y + 2 * wxy * res.
x * res.
y;
601 ofile << ms->getId() << separator << ms->x << separator << ms->y << separator;
602 ofile << res.
x << separator << res.
y << separator;
603 ofile << tpPos.
x << separator << tpPos.
y << separator;
604 ofile <<
fs->getMag() << separator << mjd << separator;
605 ofile << tpPos.
vx << separator << tpPos.
vy << separator << tpPos.
vxy << separator;
606 ofile << inputTpPos.
x << separator << inputTpPos.
y << separator;
607 ofile << inputRes.
x << separator << inputRes.
y << separator;
608 ofile <<
fs->color << separator <<
fs->getIndexInMatrix() << separator;
609 ofile <<
fs->x << separator <<
fs->y << separator;
610 ofile << chi2 << separator <<
fs->getMeasurementCount() << separator;
611 ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() <<
std::endl;
620 ofile <<
"#ra" << separator <<
"dec " << separator;
621 ofile <<
"rx" << separator <<
"ry" << separator;
622 ofile <<
"mag" << separator;
623 ofile <<
"xErr" << separator <<
"yErr" << separator <<
"xyCov" << separator;
624 ofile <<
"color" << separator <<
"fsindex" << separator;
625 ofile <<
"chi2" << separator <<
"nm" <<
std::endl;
627 ofile <<
"#coordinates of fittedStar (degrees)" << separator << separator;
628 ofile <<
"residual on TP (degrees)" << separator << separator;
629 ofile <<
"magnitude" << separator;
630 ofile <<
"refStar transformed measurement uncertainty (degrees)" << separator << separator << separator;
631 ofile <<
"currently unused" << separator <<
"unique index of the fittedStar" << separator;
632 ofile <<
"refStar contribution to Chi2 (2D dofs)" << separator
633 <<
"number of measurements of this FittedStar" <<
std::endl;
638 for (
auto const &i : fittedStarList) {
641 if (rs ==
nullptr)
continue;
646 double rx = rsProj.
x;
647 double ry = rsProj.
y;
649 double wxx = rsProj.
vy /
det;
650 double wyy = rsProj.
vx /
det;
651 double wxy = -rsProj.
vxy /
det;
652 double chi2 = wxx *
std::pow(rx, 2) + 2 * wxy * rx * ry + wyy *
std::pow(ry, 2);
655 ofile <<
fs.x << separator <<
fs.y << separator;
656 ofile << rx << separator << ry << separator;
657 ofile <<
fs.getMag() << separator;
658 ofile << rsProj.
vx << separator << rsProj.
vy << separator << rsProj.
vxy << separator;
659 ofile <<
fs.color << separator <<
fs.getIndexInMatrix() << separator;
660 ofile << chi2 << separator <<
fs.getMeasurementCount() <<
std::endl;