Line data    Source code

       1             : //# SpectralIndex.cc:
       2             : //# Copyright (C) 1998,1999,2000,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: SpectralIndex.cc 21292 2012-11-28 14:58:19Z gervandiepen $
      27             : 
      28             : #include <components/ComponentModels/SpectralIndex.h>
      29             : #include <casacore/casa/Arrays/Vector.h>
      30             : #include <casacore/casa/Containers/RecordInterface.h>
      31             : #include <casacore/casa/Exceptions/Error.h>
      32             : #include <casacore/casa/Arrays/IPosition.h>
      33             : #include <casacore/casa/Logging/LogIO.h>
      34             : #include <casacore/casa/Logging/LogOrigin.h>
      35             : #include <casacore/casa/BasicMath/Math.h>
      36             : #include <casacore/measures/Measures/MFrequency.h>
      37             : #include <casacore/measures/Measures/MCFrequency.h>
      38             : #include <casacore/measures/Measures/MeasConvert.h>
      39             : #include <casacore/casa/Quanta/MVFrequency.h>
      40             : #include <casacore/casa/Quanta/Quantum.h>
      41             : #include <casacore/casa/Utilities/Assert.h>
      42             : #include <casacore/casa/Utilities/DataType.h>
      43             : #include <casacore/casa/BasicSL/String.h>
      44             : 
      45             : using namespace casacore;
      46             : namespace casa { //# NAMESPACE CASA - BEGIN
      47             : 
      48       12347 : SpectralIndex::SpectralIndex()
      49             :   :SpectralModel(),
      50       12347 :    itsIndex(0.0),
      51           0 :    itsStokesIndex(Vector<Double>(1, 0.0)),
      52       12347 :    itsError(0.0)
      53             : {
      54       12347 :   DebugAssert(ok(), AipsError);
      55       12347 : }
      56             : 
      57          21 : SpectralIndex::SpectralIndex(const MFrequency& refFreq, Double exponent)
      58             :   :SpectralModel(refFreq),
      59          21 :    itsIndex(exponent),
      60          21 :    itsStokesIndex(Vector<Double>(1, exponent)),
      61          21 :    itsError(0.0)
      62             : {
      63          21 :   DebugAssert(ok(), AipsError);
      64          21 : }
      65             : 
      66       11902 : SpectralIndex::SpectralIndex(const SpectralIndex& other) 
      67             :   :SpectralModel(other),
      68       11902 :    itsIndex(other.itsIndex),
      69       11902 :    itsStokesIndex(other.itsStokesIndex),
      70       11902 :    itsError(other.itsError)
      71             : {
      72       11902 :   DebugAssert(ok(), AipsError);
      73       11902 : }
      74             : 
      75       48523 : SpectralIndex::~SpectralIndex() {
      76       24270 :   DebugAssert(ok(), AipsError);
      77       48523 : }
      78             : 
      79           0 : SpectralIndex& SpectralIndex::operator=(const SpectralIndex& other) {
      80           0 :   if (this != &other) {
      81           0 :     SpectralModel::operator=(other);
      82           0 :     itsIndex = other.itsIndex;
      83           0 :     itsStokesIndex.resize(other.itsStokesIndex.nelements());
      84           0 :     itsStokesIndex=other.itsStokesIndex;
      85           0 :     itsError = other.itsError;
      86             :   }
      87           0 :   DebugAssert(ok(), AipsError);
      88           0 :   return *this;
      89             : }
      90             : 
      91        5361 : ComponentType::SpectralShape SpectralIndex::type() const {
      92        5361 :   DebugAssert(ok(), AipsError);
      93        5361 :   return ComponentType::SPECTRAL_INDEX;
      94             : }
      95             : 
      96        3034 : const Double& SpectralIndex::index() const {
      97        3034 :    DebugAssert(ok(), AipsError);
      98        3034 :    return itsIndex;
      99             : }
     100           0 : const Vector<Double>& SpectralIndex::stokesIndex() const {
     101           0 :   DebugAssert(ok(), AipsError);
     102           0 :   return itsStokesIndex;
     103             : 
     104             : }
     105       12330 :   void SpectralIndex::setIndex(const Double& newIndex) { 
     106       12330 :   itsIndex = newIndex;
     107       12330 :   itsStokesIndex.resize(1);
     108       12330 :   itsStokesIndex[0]=newIndex;
     109       12330 :   DebugAssert(ok(), AipsError);
     110       12330 : }
     111             : 
     112             : 
     113          14 :   void SpectralIndex::setStokesIndex(const Vector<Double>& newIndex) { 
     114          14 :   itsIndex = newIndex[0];
     115          14 :   itsStokesIndex.resize();
     116          14 :   itsStokesIndex=newIndex;
     117          14 :   DebugAssert(ok(), AipsError);
     118          14 : }
     119             : 
     120         514 : Double SpectralIndex::sample(const MFrequency& centerFreq) const {
     121         514 :   DebugAssert(ok(), AipsError);
     122         514 :   const MFrequency& refFreq(refFrequency());
     123         514 :   const MFrequency::Ref& centerFreqFrame(centerFreq.getRef());
     124             :   Double nu0;
     125         514 :   if (centerFreqFrame != refFreq.getRef()) {
     126         514 :     nu0 = MFrequency::Convert(refFreq,centerFreqFrame)().getValue().getValue();
     127             :   } else {
     128           0 :     nu0 = refFreq.getValue().getValue();
     129             :   }
     130         514 :   if (nu0 <= 0.0) {
     131           0 :     throw(AipsError("SpectralIndex::sample(...) - "
     132           0 :                     "the reference frequency is zero or negative"));
     133             :   }
     134         514 :   const Double nu = centerFreq.getValue().getValue();
     135        1028 :   return pow(nu/nu0, itsIndex);
     136         514 : }
     137             : 
     138           0 : void SpectralIndex::sampleStokes(
     139             :     const MFrequency& centerFreq, Vector<Double>& iquv
     140             : ) const {
     141           0 :     Vector<Double> scale(4);
     142           0 :     if(itsStokesIndex.nelements()==1){
     143           0 :       scale.resize(1);
     144           0 :       scale(0)=sample(centerFreq);
     145           0 :       iquv(0) *= scale(0);
     146           0 :       return;
     147             :     }
     148           0 :     if(iquv.nelements() != 4)
     149           0 :       throw(AipsError("SpectralIndex::sampleStokes(...) 4 stokes parameters expected"));
     150           0 :     Double nu0=refFreqInFrame(centerFreq.getRef());
     151             :     
     152           0 :     if (nu0 <= 0.0) {
     153           0 :       throw(AipsError("SpectralIndex::sampleStokes(...) - "
     154           0 :                     "the reference frequency is zero or negative"));
     155             :     }
     156           0 :     const Double nu = centerFreq.getValue().getValue();
     157           0 :     iquv[0] *= pow(nu/nu0, itsStokesIndex[0]);
     158             :     //u and v get scaled so that fractional linear pol 
     159             :     // is scaled by pow(nu/nu0, alpha[1]);
     160             :     // as I is scaled by alpha[0]
     161             :     //easily shown then that u and q are scaled by pow(nu/nu0, alpha[0]+alpha[1])
     162             :     //similarly for scaling of fractional  circular pol
     163             :     
     164             : 
     165           0 :     for (uInt k=1; k < 3; ++k){
     166           0 :       iquv[k]*=pow(nu/nu0, itsStokesIndex[0]+itsStokesIndex[1]);
     167             :     }
     168             :     ///RM type rotation of linear pol. 
     169           0 :     if(itsStokesIndex[2] != 0.0){
     170             :       //angle= RM x lambda^2 : itsStokesIndex[2]=RM
     171             :       //if
     172             :       // Q'= cos(2*angle)*Q - sin(2*angle)*U
     173             :       // U'= sin(2*angle) *Q + cos(2*angle)*U
     174             :       // Q' and U' preserves fractional linear pol..and rotates it by angle
     175           0 :       Double twoalpha=2*itsStokesIndex[2]*C::c*C::c*(nu0*nu0-nu*nu)/(nu*nu*nu0*nu0);
     176           0 :       Double cos2alph=cos(twoalpha); 
     177           0 :       Double sin2alph=sin(twoalpha);
     178             :       //Q'
     179           0 :       Double tempQ=cos2alph*iquv[1]-sin2alph*iquv[2];
     180             :       //U'
     181           0 :       iquv[2]=sin2alph*iquv[1]+cos2alph*iquv[2];
     182           0 :       iquv[1]=tempQ;
     183             :     }
     184           0 :     iquv[3]*=pow(nu/nu0, itsStokesIndex[0]+itsStokesIndex[3]);
     185             :     
     186           0 : }
     187           0 : void SpectralIndex::sample(Vector<Double>& scale, 
     188             :                            const Vector<MFrequency::MVType>& frequencies, 
     189             :                            const MFrequency::Ref& refFrame) const {
     190           0 :   DebugAssert(ok(), AipsError);
     191           0 :   const uInt nSamples = frequencies.nelements();
     192           0 :   DebugAssert(scale.nelements() == nSamples, AipsError);
     193             : 
     194           0 :   const MFrequency& refFreq(refFrequency());
     195             :   Double nu0;
     196           0 :   if (refFrame != refFreq.getRef()) {
     197           0 :     nu0 = MFrequency::Convert(refFreq, refFrame)().getValue().getValue();
     198             :   } else {
     199           0 :     nu0 = refFreq.getValue().getValue();
     200             :   }
     201           0 :   if (nu0 <= 0.0) {
     202           0 :     throw(AipsError("SpectralIndex::sample(...) - "
     203           0 :                     "the reference frequency is zero or negative"));
     204             :   }
     205             : 
     206             :   Double nu;
     207           0 :   for (uInt i = 0; i < nSamples; i++) {
     208           0 :     nu = frequencies(i).getValue();
     209           0 :     scale(i) = pow(nu/nu0, itsIndex);
     210             :   }
     211           0 : }
     212             : 
     213       10010 : void SpectralIndex::sampleStokes(
     214             :     Matrix<Double>& iquv, const Vector<MFrequency::MVType>& frequencies, 
     215             :     const MFrequency::Ref& refFrame
     216             : ) const {
     217       10010 :     ThrowIf(
     218             :         iquv.shape() != IPosition(2, frequencies.size(), 4),
     219             :         "Incorrect Matrix shape"
     220             :     );
     221       10010 :     const auto nSamples = frequencies.nelements();
     222       10010 :     const auto nu0 = refFreqInFrame(refFrame); 
     223       10010 :     ThrowIf(nu0 <= 0.0, "the reference frequency is zero or negative");
     224             : 
     225             :     Double nu;
     226             :     //fractional linear and circular pol variation
     227             :     //u and v get scaled so that fractional linear pol 
     228             :     // is scaled by pow(nu/nu0, alpha[1]);
     229             :     // as I is scaled by alpha[0]
     230             :     //easily shown then that u and q are scaled by pow(nu/nu0, alpha[0]+alpha[1])
     231             :     //similarly for scaling of fractional  circular pol
     232      938904 :     for (uInt i=0; i<nSamples; ++i) {
     233      928894 :         nu = frequencies(i).getValue();
     234      928894 :         iquv(i, 0) *= pow(nu/nu0, itsStokesIndex(0));
     235      928894 :         iquv(i, 3) *= pow(nu/nu0, itsStokesIndex(0) + itsStokesIndex(3));
     236      928894 :         iquv(i, 1) *= pow(nu/nu0, itsStokesIndex[0] + itsStokesIndex[1]);
     237      928894 :         iquv(i, 2) *= pow(nu/nu0, itsStokesIndex[0] + itsStokesIndex[1]);
     238             :     }
     239             :     ///RM type rotation of linear pol. 
     240       10010 :     if(itsStokesIndex[2] != 0.0) {
     241             :         //angle= RM x lambda^2 : itsStokesIndex[2]=RM
     242             :         //if
     243             :         // Q'= cos(2*angle)*Q - sin(2*angle)*U
     244             :         // U'= sin(2*angle) *Q + cos(2*angle)*U
     245             :         // Q' and U' preserves fractional linear pol..and rotates it by angle
     246      472188 :         for (uInt i = 0; i < nSamples; i++) {
     247      468587 :             nu = frequencies(i).getValue();
     248      468587 :             Double twoalpha = 2*itsStokesIndex[2]*C::c*C::c*(nu0*nu0-nu*nu)/(nu*nu*nu0*nu0);
     249      468587 :             Double cos2alph = cos(twoalpha); 
     250      468587 :             Double sin2alph = sin(twoalpha);
     251             :             //Q'
     252      468587 :             Double tempQ = cos2alph*iquv(i, 1) - sin2alph*iquv(i, 2);
     253             :             //U'
     254      468587 :             iquv(i, 2) = sin2alph*iquv(i, 1) + cos2alph*iquv(i, 2);
     255      468587 :             iquv(i, 1) = tempQ;
     256             :         }
     257             :     }
     258       10010 : }
     259             :    
     260       11902 : SpectralModel* SpectralIndex::clone() const {
     261       11902 :   DebugAssert(ok(), AipsError);
     262       11902 :   SpectralModel* tmpPtr = new SpectralIndex(*this);
     263       11902 :   AlwaysAssert(tmpPtr != 0, AipsError);
     264       11902 :   return tmpPtr;
     265             : }
     266             : 
     267       14660 : uInt SpectralIndex::nParameters() const {
     268       14660 :   DebugAssert(ok(), AipsError);
     269       14660 :   return 1;
     270             : }
     271             : 
     272        1165 : void SpectralIndex::setParameters(const Vector<Double>& newSpectralParms) {
     273        1165 :   DebugAssert(newSpectralParms.nelements() == nParameters(), AipsError);
     274        1165 :   itsIndex = newSpectralParms(0);
     275        1165 :   DebugAssert(ok(), AipsError);
     276        1165 : }
     277             : 
     278        2327 : Vector<Double> SpectralIndex::parameters() const {
     279        2327 :   DebugAssert(ok(), AipsError);
     280        2327 :   return Vector<Double>(1, itsIndex);
     281             : }
     282             : 
     283       13495 : void SpectralIndex::setErrors(const Vector<Double>& newSpectralErrs) {
     284       13495 :   DebugAssert(newSpectralErrs.nelements() == nParameters(), AipsError);
     285       13495 :   if (newSpectralErrs(0) < 0.0) {
     286           0 :     LogIO logErr(LogOrigin("SpectralIndex", "setErrors(...)"));
     287             :     logErr << "The errors must be non-negative."
     288           0 :            << LogIO::EXCEPTION;
     289           0 :   }
     290       13495 :   itsError = newSpectralErrs(0);
     291       13495 :   DebugAssert(ok(), AipsError);
     292       13495 : }
     293             : 
     294        5361 : Vector<Double> SpectralIndex::errors() const {
     295        5361 :   DebugAssert(ok(), AipsError);
     296        5361 :   return Vector<Double>(1, itsError);
     297             : }
     298             : 
     299       12330 : Bool SpectralIndex::fromRecord(String& errorMessage, 
     300             :                                const RecordInterface& record) {
     301       12330 :   if (!SpectralModel::fromRecord(errorMessage, record)) return false;
     302       12330 :   if (!record.isDefined(String("index"))) {
     303           0 :     errorMessage += "The 'spectrum' record must have an 'index' field\n";
     304           0 :     return false;
     305             :   }
     306             : //
     307             :   {
     308       12330 :      const RecordFieldId index("index");
     309       12330 :      const IPosition shape(1,1);
     310       12330 :      if (record.shape(index) != shape) {
     311           0 :        errorMessage += "The 'index' field must be a scalar\n";
     312           0 :        return false;
     313             :      }
     314             :      Double indexVal;
     315       12330 :      switch (record.dataType(index)) {
     316       12330 :      case TpDouble:
     317             :      case TpFloat:
     318             :      case TpInt:
     319       12330 :        indexVal = record.asDouble(index);
     320       12330 :        break;
     321           0 :      default:
     322           0 :        errorMessage += "The 'index' field must be a real number\n";
     323           0 :        return false;
     324             :      }
     325       12330 :      setIndex(indexVal);
     326       12330 :   }
     327       12330 :   if(record.isDefined("stokesindex")){
     328       12330 :     Vector<Double> tempstokes=record.asArrayDouble("stokesindex");
     329       12330 :     if((tempstokes.nelements() != 4) && (tempstokes.nelements() != 1) ){
     330           0 :       errorMessage += "Stokes indices is not of length 1 or 4\n";
     331           0 :       return false;
     332             :     }
     333       12330 :     itsStokesIndex.resize();
     334       12330 :     itsStokesIndex=tempstokes;
     335       12330 :   }
     336             : //
     337             :   {
     338       12330 :       Vector<Double> errorVals(1, 0.0);
     339       12330 :       if (record.isDefined("error")) {
     340       12330 :         const RecordFieldId error("error");
     341       12330 :         const IPosition shape(1,1);
     342       12330 :         if (record.shape(error) != shape) {
     343           0 :           errorMessage += "The 'error' field must be a scalar\n";
     344           0 :           return false;
     345             :         }
     346       12330 :         switch (record.dataType(error)) {
     347       12330 :         case TpDouble:
     348             :         case TpFloat:
     349             :         case TpInt:
     350       12330 :             errorVals[0] = record.asDouble(error);
     351       12330 :           break;
     352           0 :         default:
     353           0 :           errorMessage += "The 'error' field must be a real number\n";
     354           0 :           return false;
     355             :         }
     356       12330 :      }
     357             : //
     358       12330 :      setErrors(errorVals);
     359       12330 :   }
     360             : //
     361       12330 :   DebugAssert(ok(), AipsError);
     362       12330 :   return true;
     363             : }
     364             : 
     365        3034 : Bool SpectralIndex::toRecord(String& errorMessage,
     366             :                              RecordInterface& record) const {
     367        3034 :   DebugAssert(ok(), AipsError);
     368        3034 :   if (!SpectralModel::toRecord(errorMessage, record)) return false;
     369        3034 :   record.define("index", index());
     370        3034 :   if(itsStokesIndex.nelements() != 0)
     371        3034 :     record.define("stokesindex", itsStokesIndex);
     372        3034 :   record.define("error", errors()(0));
     373        3034 :   return true;
     374             : }
     375             : 
     376           0 : Bool SpectralIndex::convertUnit(String& errorMessage,
     377             :                                 const RecordInterface& record) {
     378           0 :   const String fieldString("index");
     379           0 :   if (!record.isDefined(fieldString)) {
     380           0 :     return true;
     381             :   }
     382           0 :   const RecordFieldId field(fieldString);
     383           0 :   if (!((record.dataType(field) == TpString) && 
     384           0 :         (record.shape(field) == IPosition(1,1)) &&
     385           0 :         (record.asString(field) == ""))) {
     386           0 :     errorMessage += "The 'index' field must be an empty string\n";
     387           0 :     errorMessage += "(and not a vector of strings)\n";
     388           0 :     return false;
     389             :   }
     390           0 :   return true;
     391           0 : }
     392             : 
     393      135948 : Bool SpectralIndex::ok() const {
     394      135948 :     if (!SpectralModel::ok()) return false;
     395      135948 :     ThrowIf(
     396             :         refFrequency().getValue().getValue() <= 0.0,
     397             :         "The reference frequency is zero or negative!"
     398             :     ); 
     399      135948 :     ThrowIf(abs(
     400             :         itsIndex) > 100,
     401             :         "The absolute value of the spectral index is greater than 100!"
     402             :     ); 
     403      135948 :     return true;
     404             : }
     405             : 
     406             : }
     407             :