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: casa-feedback@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       51202 : SkyComponent::SkyComponent()
      58       51202 :   :itsCompPtr(new SkyCompRep) 
      59       51202 : {}
      60             : 
      61           2 : SkyComponent::SkyComponent(const ComponentType::Shape& shape)   
      62           2 :   :itsCompPtr(new SkyCompRep(shape))
      63           2 : {}
      64             : 
      65           1 : SkyComponent::SkyComponent(const ComponentType::Shape& shape,
      66           1 :                            const ComponentType::SpectralShape& spectralModel) 
      67           1 :   :itsCompPtr(new SkyCompRep(shape, spectralModel))
      68           1 : {}
      69             : 
      70           1 : SkyComponent::SkyComponent(const Flux<Double>& flux,
      71             :                            const ComponentShape& shape, 
      72           1 :                            const SpectralModel& spectrum)
      73           1 :   :itsCompPtr(new SkyCompRep(flux, shape, spectrum))
      74           1 : {}
      75             : 
      76       51202 : SkyComponent::SkyComponent(const SkyComponent& other) 
      77             :   :SkyCompBase(other),
      78       51202 :    itsCompPtr (other.itsCompPtr)
      79       51202 : {}
      80             : 
      81      102406 : SkyComponent::~SkyComponent() {}
      82             : 
      83       25602 : SkyComponent& SkyComponent::operator=(const SkyComponent& other) {
      84       25602 :   if (this != &other)
      85       25602 :     itsCompPtr = other.itsCompPtr;
      86       25602 :   return *this;
      87             : }
      88             : 
      89       51203 : Flux<Double>& SkyComponent::flux() {
      90       51203 :   return itsCompPtr->flux();
      91             : }
      92             : 
      93           0 : const Flux<Double>& SkyComponent::flux() const {
      94           0 :   return itsCompPtr->flux();
      95             : }
      96             : 
      97         800 : const ComponentShape& SkyComponent::shape() const {
      98         800 :   return itsCompPtr->shape();
      99             : }
     100             : 
     101       51207 : ComponentShape& SkyComponent::shape() {
     102       51207 :   return itsCompPtr->shape();
     103             : }
     104             : 
     105           0 : void SkyComponent::setShape(const ComponentShape& newShape) {
     106           0 :   itsCompPtr->setShape(newShape);
     107           0 : }
     108             : 
     109           0 : const SpectralModel& SkyComponent::spectrum() const {
     110           0 :   return itsCompPtr->spectrum();
     111             : }
     112             : 
     113           0 : SpectralModel& SkyComponent::spectrum() {
     114           0 :   return itsCompPtr->spectrum();
     115             : }
     116             : 
     117           0 : void SkyComponent::setSpectrum(const SpectralModel& newSpectrum) {
     118           0 :   itsCompPtr->setSpectrum(newSpectrum);
     119           0 : }
     120             : 
     121           0 : String& SkyComponent::label() {
     122           0 :   return itsCompPtr->label();
     123             : }
     124             : 
     125           0 : const String& SkyComponent::label() const {
     126           0 :   return itsCompPtr->label();
     127             : }
     128             : 
     129           0 : Vector<Double>& SkyComponent::optionalParameters() {
     130           0 :    return itsCompPtr->optionalParameters();
     131             : }
     132             : 
     133           0 : const Vector<Double>& SkyComponent::optionalParameters() const {
     134           0 :    return itsCompPtr->optionalParameters();
     135             : }
     136             : 
     137           0 : Bool SkyComponent::isPhysical() const {
     138           0 :   return itsCompPtr->isPhysical();
     139             : }
     140             : 
     141           0 : Flux<Double> SkyComponent::sample(const MDirection& direction, 
     142             :                               const MVAngle& pixelLatSize, 
     143             :                               const MVAngle& pixelLongSize, 
     144             :                               const MFrequency& centerFrequency) const {
     145           0 :   return itsCompPtr->sample(direction, pixelLatSize, pixelLongSize, 
     146           0 :                             centerFrequency);
     147             : }
     148             : 
     149           0 : 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           0 :   itsCompPtr->sample(samples, reqUnit,
     157             :                      directions, dirRef, pixelLatSize, pixelLongSize,
     158             :                      frequencies, freqRef);
     159           0 : }
     160             : 
     161           0 : Flux<Double> SkyComponent::visibility(const Vector<Double>& uvw,
     162             :                                       const Double& frequency) const {
     163           0 :   return itsCompPtr->visibility(uvw, frequency);
     164             : }
     165             : 
     166       25600 : void SkyComponent::visibility(Cube<DComplex>& visibilities,
     167             :                               const Matrix<Double>& uvws,
     168             :                               const Vector<Double>& frequencies) const {
     169       25600 :   itsCompPtr->visibility(visibilities, uvws, frequencies);
     170       25600 : }
     171             : 
     172       25600 : Bool SkyComponent::fromRecord(String& errorMessage, 
     173             :                               const RecordInterface& record) {
     174       25600 :   return itsCompPtr->fromRecord(errorMessage, record);
     175             : }
     176             : 
     177         800 : Bool SkyComponent::toRecord(String& errorMessage,
     178             :                             RecordInterface& record) const {
     179         800 :   return itsCompPtr->toRecord(errorMessage, record);
     180             : }
     181             : 
     182           0 : SkyComponent SkyComponent::copy() const {
     183           0 :   SkyComponent newComp(flux().copy(), shape(), spectrum());
     184           0 :   newComp.label() = label();
     185           0 :   newComp.optionalParameters() = optionalParameters();
     186           0 :   return newComp;
     187           0 : }
     188             : 
     189           0 : 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           0 :    itsCompPtr->fromPixel(fluxRatio, parameters, brightnessUnitIn, restoringBeam,
     197             :                          cSys, componentShape, stokes);
     198           0 : }
     199             : 
     200           0 : 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           0 :                               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           0 :   }
     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           0 :   }
     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           0 : }
     248             : 
     249           0 : String SkyComponent::positionToString(
     250             :     std::shared_ptr<casacore::Vector<casacore::Double>>& pixelCoords,
     251             :         const DirectionCoordinate *const &dc ,Bool longErrOnGreatCircle
     252             : ) const {
     253           0 :     pixelCoords.reset();
     254             :         // FIXME essentially cut and paste of Gareth's python code. Needs work.
     255           0 :         ostringstream position;
     256           0 :         MDirection mdir = shape().refDirection();
     257             : 
     258           0 :         Quantity lat = mdir.getValue().getLat("rad");
     259           0 :         String latString = MVAngle(lat).string(MVAngle::ANGLE_CLEAN, 8);
     260             : 
     261           0 :         Quantity longitude = mdir.getValue().getLong("rad");
     262           0 :         String longString = MVTime(longitude).string(MVTime::TIME, 9);
     263             : 
     264           0 :         Quantity dLat = shape().refDirectionErrorLat();
     265           0 :         dLat.convert("rad");
     266             : 
     267           0 :         Quantity dLong = shape().refDirectionErrorLong();
     268           0 :         dLong.convert("rad");
     269             : 
     270             :         // Add error estimates to ra/dec strings if an error is given (either >0)
     271             : 
     272           0 :         uInt precision = 1;
     273           0 :         if ( dLong.getValue() != 0 || dLat.getValue() != 0 ) {
     274           0 :                 dLong.convert("s");
     275           0 :                 dLat.convert("arcsec");
     276           0 :                 Double drasec  = roundDouble(dLong.getValue());
     277           0 :                 Double ddecarcsec = roundDouble(dLat.getValue());
     278           0 :                 Vector<Double> dravec(2), ddecvec(2);
     279           0 :                 dravec.set(drasec);
     280           0 :                 ddecvec.set(ddecarcsec);
     281           0 :                 precision = precisionForValueErrorPairs(dravec,ddecvec);
     282           0 :                 longString = MVTime(longitude).string(MVTime::TIME, 6+precision);
     283           0 :                 latString =  MVAngle(lat).string(MVAngle::ANGLE, 6+precision);
     284           0 :         }
     285           0 :         std::pair<String, String> labels = _axisLabels(dc);
     286           0 :         position << "Position ---" << endl;
     287           0 :         position << "       --- " << labels.first << "   " << longString;
     288           0 :         if (dLong.getValue() == 0) {
     289           0 :                 position << " (fixed)" << endl;
     290             :         }
     291             :         else {
     292           0 :                 Quantity timeError = longErrOnGreatCircle ? dLong/cos(lat) : dLong;
     293           0 :                 position << " +/- " << std::fixed
     294           0 :                         << std::setprecision(precision) << timeError << " ("
     295           0 :                         << dLong.getValue("arcsec") << " arcsec"
     296             :                         << (longErrOnGreatCircle ? " along great circle" : "")
     297           0 :                         << ")" << endl;
     298           0 :         }
     299           0 :         position << "       --- " << labels.second << " " << latString;
     300           0 :         if (dLat.getValue() == 0) {
     301           0 :                 position << " (fixed)" << endl;
     302             :         }
     303             :         else {
     304           0 :                 position << " +/- " << dLat << endl;
     305             :         }
     306           0 :         if (dc) {
     307           0 :                 const Vector<String> units = dc->worldAxisUnits();
     308           0 :                 Vector<Double> world(dc->nWorldAxes(), 0);
     309           0 :                 world[0] = longitude.getValue(units[0]);
     310           0 :                 world[1] = lat.getValue(units[1]);
     311             :                 // TODO do the pixel computations in another method
     312           0 :                 pixelCoords.reset(new Vector<Double>());
     313           0 :                 if (dc->toPixel(*pixelCoords, world)) {
     314           0 :                         Vector<Double> increment = dc->increment();
     315           0 :                         Double longPixErr = dLong.getValue() == 0
     316           0 :                                 ? 0 : abs(dLong.getValue(units[0])/increment[0]);
     317           0 :                         Double latPixErr = dLat.getValue() == 0
     318           0 :                                 ? 0 : abs(dLat.getValue(units[1])/increment[1]);
     319           0 :                         Vector<Double> longPix(2), latPix(2);
     320           0 :                         longPix.set(roundDouble(longPixErr));
     321           0 :                         latPix.set(roundDouble(latPixErr));
     322           0 :                         precision = precisionForValueErrorPairs(longPix, latPix);
     323           0 :                         position << std::fixed <<  std::setprecision(precision);
     324           0 :                         position << "       --- " << labels.first << " " << (*pixelCoords)[0];
     325           0 :                         if (dLong.getValue() == 0) {
     326           0 :                                 position << " (fixed)" << endl;
     327             :                         }
     328             :                         else {
     329           0 :                                 position << " +/- " << longPixErr << " pixels" << endl;
     330             :                         }
     331           0 :                         position << "       --- " << labels.second << " " << (*pixelCoords)[1];
     332           0 :                         if (dLat.getValue() == 0) {
     333           0 :                                 position << " (fixed)" << endl;
     334             :                         }
     335             :                         else {
     336           0 :                                 position << " +/- " << latPixErr << " pixels" << endl;
     337             :                         }
     338           0 :                 }
     339             :                 else {
     340           0 :                         position << "unable to determine position in pixels:" << dc->errorMessage() << endl;
     341             :                 }
     342           0 :         }
     343           0 :         return position.str();
     344           0 : }
     345             : 
     346           0 : std::pair<String, String> SkyComponent::_axisLabels(
     347             :         const DirectionCoordinate *const &dc
     348             : ) {
     349           0 :         std::pair<String, String> labels;
     350           0 :         if (dc) {
     351           0 :                 Vector<String> names = dc->worldAxisNames();
     352           0 :                 for (uInt i=0; i<2; i++) {
     353           0 :                         String label;
     354           0 :                         names[i].downcase();
     355           0 :                         if (names[i] == "right ascension") {
     356           0 :                                 label = "ra:";
     357             :                         }
     358           0 :                         else if (names[i] == "declination") {
     359           0 :                                 label = "dec:";
     360             :                         }
     361           0 :                         else if (names[i] == "longitude") {
     362           0 :                                 label = "long:";
     363             :                         }
     364           0 :                         else if (names[i] == "latitude") {
     365           0 :                                 label = "lat:";
     366             :                         }
     367             :                         else {
     368           0 :                                 label = names[i] + ":";
     369             :                         }
     370           0 :                         if (i == 0) {
     371           0 :                                 labels.first = label;
     372             :                         }
     373             :                         else {
     374           0 :                                 labels.second = label;
     375             :                         }
     376           0 :                 }
     377           0 :         }
     378             :         else {
     379           0 :                 labels.first  = "long:";
     380           0 :                 labels.second = "lat: ";
     381             :         }
     382             :         typedef std::string::size_type size_type;
     383           0 :         size_type f = labels.first.size();
     384           0 :         size_type s = labels.second.size();
     385           0 :         if (f > s) {
     386           0 :                 size_type d = f - s;
     387           0 :                 for (size_type i=0; i<d ;i++) {
     388           0 :                         labels.second += " ";
     389             :                 }
     390             :         }
     391           0 :         else if (f < s) {
     392           0 :                 size_type d = s - f;
     393           0 :                 for (size_type i=0; i<d ;i++) {
     394           0 :                         labels.first += " ";
     395             :                 }
     396             :         }
     397           0 :         return labels;
     398           0 : }
     399             : 
     400             : }
     401             :