Line data    Source code

       1             : //# SkyComponent.cc:  this defines SkyComponent
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2003
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: SkyComponent.cc 21292 2012-11-28 14:58:19Z gervandiepen $
      27             : 
      28             : #include <casacore/casa/Quanta/QMath.h>
      29             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      30             : #include <components/ComponentModels/SkyComponent.h>
      31             : #include <components/ComponentModels/ComponentShape.h>
      32             : #include <components/ComponentModels/Flux.h>
      33             : #include <components/ComponentModels/SkyCompRep.h>
      34             : #include <components/ComponentModels/SpectralModel.h>
      35             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      36             : #include <casacore/casa/Arrays/Vector.h>
      37             : #include <casacore/casa/Containers/RecordInterface.h>
      38             : #include <casacore/casa/Exceptions/Error.h>
      39             : #include <iomanip>
      40             : #include <casacore/casa/Logging/LogIO.h>
      41             : #include <casacore/casa/Logging/LogOrigin.h>
      42             : #include <casacore/casa/Quanta/MVTime.h>
      43             : #include <casacore/casa/Utilities/Precision.h>
      44             : 
      45             : #include <casacore/measures/Measures/MDirection.h>
      46             : #include <casacore/measures/Measures/MFrequency.h>
      47             : 
      48             : #include <casacore/casa/Quanta/MVAngle.h>
      49             : #include <casacore/casa/Utilities/Assert.h>
      50             : #include <casacore/casa/BasicSL/String.h>
      51             : 
      52             : #include <iomanip>
      53             : 
      54             : using namespace casacore;
      55             : namespace casa { //# NAMESPACE CASA - BEGIN
      56             : 
      57        9976 : SkyComponent::SkyComponent()
      58        9976 :   :itsCompPtr(new SkyCompRep) 
      59        9976 : {}
      60             : 
      61          24 : SkyComponent::SkyComponent(const ComponentType::Shape& shape)   
      62          24 :   :itsCompPtr(new SkyCompRep(shape))
      63          24 : {}
      64             : 
      65         308 : SkyComponent::SkyComponent(const ComponentType::Shape& shape,
      66         308 :                            const ComponentType::SpectralShape& spectralModel) 
      67         308 :   :itsCompPtr(new SkyCompRep(shape, spectralModel))
      68         308 : {}
      69             : 
      70       11204 : SkyComponent::SkyComponent(const Flux<Double>& flux,
      71             :                            const ComponentShape& shape, 
      72       11204 :                            const SpectralModel& spectrum)
      73       11204 :   :itsCompPtr(new SkyCompRep(flux, shape, spectrum))
      74       11204 : {}
      75             : 
      76        6028 : SkyComponent::SkyComponent(const SkyComponent& other) 
      77             :   :SkyCompBase(other),
      78        6028 :    itsCompPtr (other.itsCompPtr)
      79        6028 : {}
      80             : 
      81       27540 : SkyComponent::~SkyComponent() {}
      82             : 
      83        8720 : SkyComponent& SkyComponent::operator=(const SkyComponent& other) {
      84        8720 :   if (this != &other)
      85        8720 :     itsCompPtr = other.itsCompPtr;
      86        8720 :   return *this;
      87             : }
      88             : 
      89       43388 : Flux<Double>& SkyComponent::flux() {
      90       43388 :   return itsCompPtr->flux();
      91             : }
      92             : 
      93       11083 : const Flux<Double>& SkyComponent::flux() const {
      94       11083 :   return itsCompPtr->flux();
      95             : }
      96             : 
      97       14211 : const ComponentShape& SkyComponent::shape() const {
      98       14211 :   return itsCompPtr->shape();
      99             : }
     100             : 
     101      154842 : ComponentShape& SkyComponent::shape() {
     102      154842 :   return itsCompPtr->shape();
     103             : }
     104             : 
     105        3076 : void SkyComponent::setShape(const ComponentShape& newShape) {
     106        3076 :   itsCompPtr->setShape(newShape);
     107        3076 : }
     108             : 
     109       11476 : const SpectralModel& SkyComponent::spectrum() const {
     110       11476 :   return itsCompPtr->spectrum();
     111             : }
     112             : 
     113        1049 : SpectralModel& SkyComponent::spectrum() {
     114        1049 :   return itsCompPtr->spectrum();
     115             : }
     116             : 
     117         451 : void SkyComponent::setSpectrum(const SpectralModel& newSpectrum) {
     118         451 :   itsCompPtr->setSpectrum(newSpectrum);
     119         451 : }
     120             : 
     121       11813 : String& SkyComponent::label() {
     122       11813 :   return itsCompPtr->label();
     123             : }
     124             : 
     125       11083 : const String& SkyComponent::label() const {
     126       11083 :   return itsCompPtr->label();
     127             : }
     128             : 
     129       11083 : Vector<Double>& SkyComponent::optionalParameters() {
     130       11083 :    return itsCompPtr->optionalParameters();
     131             : }
     132             : 
     133       11083 : const Vector<Double>& SkyComponent::optionalParameters() const {
     134       11083 :    return itsCompPtr->optionalParameters();
     135             : }
     136             : 
     137           0 : Bool SkyComponent::isPhysical() const {
     138           0 :   return itsCompPtr->isPhysical();
     139             : }
     140             : 
     141        2491 : Flux<Double> SkyComponent::sample(const MDirection& direction, 
     142             :                               const MVAngle& pixelLatSize, 
     143             :                               const MVAngle& pixelLongSize, 
     144             :                               const MFrequency& centerFrequency) const {
     145        2491 :   return itsCompPtr->sample(direction, pixelLatSize, pixelLongSize, 
     146        2491 :                             centerFrequency);
     147             : }
     148             : 
     149        4439 : void SkyComponent::sample(Cube<Double>& samples, const Unit& reqUnit,
     150             :                           const Vector<MVDirection>& directions, 
     151             :                           const MeasRef<MDirection>& dirRef, 
     152             :                           const MVAngle& pixelLatSize, 
     153             :                           const MVAngle& pixelLongSize, 
     154             :                           const Vector<MVFrequency>& frequencies,
     155             :                           const MeasRef<MFrequency>& freqRef) const {
     156        4439 :   itsCompPtr->sample(samples, reqUnit,
     157             :                      directions, dirRef, pixelLatSize, pixelLongSize,
     158        4439 :                      frequencies, freqRef);
     159        4439 : }
     160             : 
     161     2128905 : Flux<Double> SkyComponent::visibility(const Vector<Double>& uvw,
     162             :                                       const Double& frequency) const {
     163     2128905 :   return itsCompPtr->visibility(uvw, frequency);
     164             : }
     165             : 
     166        5435 : void SkyComponent::visibility(Cube<DComplex>& visibilities,
     167             :                               const Matrix<Double>& uvws,
     168             :                               const Vector<Double>& frequencies) const {
     169        5435 :   itsCompPtr->visibility(visibilities, uvws, frequencies);
     170        5435 : }
     171             : 
     172        5840 : Bool SkyComponent::fromRecord(String& errorMessage, 
     173             :                               const RecordInterface& record) {
     174        5840 :   return itsCompPtr->fromRecord(errorMessage, record);
     175             : }
     176             : 
     177        6537 : Bool SkyComponent::toRecord(String& errorMessage,
     178             :                             RecordInterface& record) const {
     179        6537 :   return itsCompPtr->toRecord(errorMessage, record);
     180             : }
     181             : 
     182       11083 : SkyComponent SkyComponent::copy() const {
     183       11083 :   SkyComponent newComp(flux().copy(), shape(), spectrum());
     184       11083 :   newComp.label() = label();
     185       11083 :   newComp.optionalParameters() = optionalParameters();
     186       11083 :   return newComp;
     187             : }
     188             : 
     189         393 : void SkyComponent::fromPixel (Double& fluxRatio, const Vector<Double>& parameters,
     190             :                               const Unit& brightnessUnitIn,
     191             :                               const GaussianBeam& restoringBeam,
     192             :                               const CoordinateSystem& cSys,
     193             :                               ComponentType::Shape componentShape,
     194             :                               Stokes::StokesTypes stokes)
     195             : {
     196         393 :    itsCompPtr->fromPixel(fluxRatio, parameters, brightnessUnitIn, restoringBeam,
     197             :                          cSys, componentShape, stokes);
     198         393 : }
     199             : 
     200          58 : Vector<Double> SkyComponent::toPixel (const Unit& brightnessUnitIn,
     201             :                                       const GaussianBeam& restoringBeam,
     202             :                                       const CoordinateSystem& cSys,
     203             :                                       Stokes::StokesTypes stokes) const
     204             : {
     205             :    return itsCompPtr->toPixel(brightnessUnitIn, restoringBeam,
     206          58 :                               cSys, stokes);
     207             : }
     208             : 
     209             : 
     210             : 
     211             : 
     212           0 : Bool SkyComponent::ok() const {
     213           0 :   if (itsCompPtr.null()) {
     214           0 :     LogIO logErr(LogOrigin("SkyComponent", "ok()"));
     215             :     logErr << LogIO::SEVERE << "Internal pointer is not pointing to anything"
     216           0 :            << LogIO::POST;
     217           0 :     return false;
     218             :   }
     219           0 :   if (!itsCompPtr->ok()) {
     220           0 :     LogIO logErr(LogOrigin("SkyComponent", "ok()"));
     221             :     logErr << LogIO::SEVERE << "Component representation is not ok"
     222           0 :            << LogIO::POST;
     223           0 :     return false;
     224             :   }
     225           0 :   return true;
     226             : }
     227             : 
     228           0 : String SkyComponent::summarize(
     229             :         const DirectionCoordinate *const &dc, Bool longErrOnGreatCircle
     230             : ) const {
     231           0 :         ostringstream ldpar;
     232           0 :         if (shape().type()==ComponentType::LDISK) {
     233           0 :                 ldpar << " (limb-darkening exponent: "<<optionalParameters()(0) <<" )"<<endl;
     234             :         }
     235           0 :         ostringstream summary;
     236           0 :         summary << "SUMMARY OF COMPONENT " << label() << endl;
     237           0 :         summary << "Shape: " << shape().ident() << ldpar.str() <<endl;
     238           0 :         const Flux<Double> myFlux = flux();
     239           0 :         Quantum<Vector<std::complex<double> > > fluxValue;
     240           0 :         myFlux.value(fluxValue);
     241           0 :         summary << "Flux density: " << fluxValue << " +/- " << myFlux.errors() << endl;
     242           0 :         summary << "Spectral model: " << spectrum().ident() << endl;
     243           0 :         std::shared_ptr<Vector<Double>> pixCoords;
     244           0 :         summary << "Position: " <<  positionToString(pixCoords, dc, longErrOnGreatCircle) << endl;
     245           0 :         summary << "Size: " << endl << shape().sizeToString() << endl;
     246           0 :         return summary.str();
     247             : }
     248             : 
     249         393 : String SkyComponent::positionToString(
     250             :     std::shared_ptr<casacore::Vector<casacore::Double>>& pixelCoords,
     251             :         const DirectionCoordinate *const &dc ,Bool longErrOnGreatCircle
     252             : ) const {
     253         393 :     pixelCoords.reset();
     254             :         // FIXME essentially cut and paste of Gareth's python code. Needs work.
     255         786 :         ostringstream position;
     256         786 :         MDirection mdir = shape().refDirection();
     257             : 
     258         786 :         Quantity lat = mdir.getValue().getLat("rad");
     259         786 :         String latString = MVAngle(lat).string(MVAngle::ANGLE_CLEAN, 8);
     260             : 
     261         786 :         Quantity longitude = mdir.getValue().getLong("rad");
     262         786 :         String longString = MVTime(longitude).string(MVTime::TIME, 9);
     263             : 
     264         786 :         Quantity dLat = shape().refDirectionErrorLat();
     265         393 :         dLat.convert("rad");
     266             : 
     267         786 :         Quantity dLong = shape().refDirectionErrorLong();
     268         393 :         dLong.convert("rad");
     269             : 
     270             :         // Add error estimates to ra/dec strings if an error is given (either >0)
     271             : 
     272         393 :         uInt precision = 1;
     273         393 :         if ( dLong.getValue() != 0 || dLat.getValue() != 0 ) {
     274         393 :                 dLong.convert("s");
     275         393 :                 dLat.convert("arcsec");
     276         393 :                 Double drasec  = roundDouble(dLong.getValue());
     277         393 :                 Double ddecarcsec = roundDouble(dLat.getValue());
     278         786 :                 Vector<Double> dravec(2), ddecvec(2);
     279         393 :                 dravec.set(drasec);
     280         393 :                 ddecvec.set(ddecarcsec);
     281         393 :                 precision = precisionForValueErrorPairs(dravec,ddecvec);
     282         393 :                 longString = MVTime(longitude).string(MVTime::TIME, 6+precision);
     283         393 :                 latString =  MVAngle(lat).string(MVAngle::ANGLE, 6+precision);
     284             :         }
     285         393 :         std::pair<String, String> labels = _axisLabels(dc);
     286         393 :         position << "Position ---" << endl;
     287         393 :         position << "       --- " << labels.first << "   " << longString;
     288         393 :         if (dLong.getValue() == 0) {
     289           0 :                 position << " (fixed)" << endl;
     290             :         }
     291             :         else {
     292         786 :                 Quantity timeError = longErrOnGreatCircle ? dLong/cos(lat) : dLong;
     293         393 :                 position << " +/- " << std::fixed
     294         393 :                         << std::setprecision(precision) << timeError << " ("
     295         786 :                         << dLong.getValue("arcsec") << " arcsec"
     296             :                         << (longErrOnGreatCircle ? " along great circle" : "")
     297         393 :                         << ")" << endl;
     298             :         }
     299         393 :         position << "       --- " << labels.second << " " << latString;
     300         393 :         if (dLat.getValue() == 0) {
     301           0 :                 position << " (fixed)" << endl;
     302             :         }
     303             :         else {
     304         393 :                 position << " +/- " << dLat << endl;
     305             :         }
     306         393 :         if (dc) {
     307         786 :                 const Vector<String> units = dc->worldAxisUnits();
     308         786 :                 Vector<Double> world(dc->nWorldAxes(), 0);
     309         393 :                 world[0] = longitude.getValue(units[0]);
     310         393 :                 world[1] = lat.getValue(units[1]);
     311             :                 // TODO do the pixel computations in another method
     312         393 :                 pixelCoords.reset(new Vector<Double>());
     313         393 :                 if (dc->toPixel(*pixelCoords, world)) {
     314         786 :                         Vector<Double> increment = dc->increment();
     315         393 :                         Double longPixErr = dLong.getValue() == 0
     316         393 :                                 ? 0 : abs(dLong.getValue(units[0])/increment[0]);
     317         393 :                         Double latPixErr = dLat.getValue() == 0
     318         393 :                                 ? 0 : abs(dLat.getValue(units[1])/increment[1]);
     319         786 :                         Vector<Double> longPix(2), latPix(2);
     320         393 :                         longPix.set(roundDouble(longPixErr));
     321         393 :                         latPix.set(roundDouble(latPixErr));
     322         393 :                         precision = precisionForValueErrorPairs(longPix, latPix);
     323         393 :                         position << std::fixed <<  std::setprecision(precision);
     324         393 :                         position << "       --- " << labels.first << " " << (*pixelCoords)[0];
     325         393 :                         if (dLong.getValue() == 0) {
     326           0 :                                 position << " (fixed)" << endl;
     327             :                         }
     328             :                         else {
     329         393 :                                 position << " +/- " << longPixErr << " pixels" << endl;
     330             :                         }
     331         393 :                         position << "       --- " << labels.second << " " << (*pixelCoords)[1];
     332         393 :                         if (dLat.getValue() == 0) {
     333           0 :                                 position << " (fixed)" << endl;
     334             :                         }
     335             :                         else {
     336         393 :                                 position << " +/- " << latPixErr << " pixels" << endl;
     337             :                         }
     338             :                 }
     339             :                 else {
     340           0 :                         position << "unable to determine position in pixels:" << dc->errorMessage() << endl;
     341             :                 }
     342             :         }
     343         786 :         return position.str();
     344             : }
     345             : 
     346         393 : std::pair<String, String> SkyComponent::_axisLabels(
     347             :         const DirectionCoordinate *const &dc
     348             : ) {
     349         393 :         std::pair<String, String> labels;
     350         393 :         if (dc) {
     351         786 :                 Vector<String> names = dc->worldAxisNames();
     352        1179 :                 for (uInt i=0; i<2; i++) {
     353        1572 :                         String label;
     354         786 :                         names[i].downcase();
     355         786 :                         if (names[i] == "right ascension") {
     356         338 :                                 label = "ra:";
     357             :                         }
     358         448 :                         else if (names[i] == "declination") {
     359         338 :                                 label = "dec:";
     360             :                         }
     361         110 :                         else if (names[i] == "longitude") {
     362          55 :                                 label = "long:";
     363             :                         }
     364          55 :                         else if (names[i] == "latitude") {
     365          55 :                                 label = "lat:";
     366             :                         }
     367             :                         else {
     368           0 :                                 label = names[i] + ":";
     369             :                         }
     370         786 :                         if (i == 0) {
     371         393 :                                 labels.first = label;
     372             :                         }
     373             :                         else {
     374         393 :                                 labels.second = label;
     375             :                         }
     376             :                 }
     377             :         }
     378             :         else {
     379           0 :                 labels.first  = "long:";
     380           0 :                 labels.second = "lat: ";
     381             :         }
     382             :         typedef std::string::size_type size_type;
     383         393 :         size_type f = labels.first.size();
     384         393 :         size_type s = labels.second.size();
     385         393 :         if (f > s) {
     386          55 :                 size_type d = f - s;
     387         110 :                 for (size_type i=0; i<d ;i++) {
     388          55 :                         labels.second += " ";
     389             :                 }
     390             :         }
     391         338 :         else if (f < s) {
     392         338 :                 size_type d = s - f;
     393         676 :                 for (size_type i=0; i<d ;i++) {
     394         338 :                         labels.first += " ";
     395             :                 }
     396             :         }
     397         393 :         return labels;
     398             : }
     399             : 
     400             : }
     401             :