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: 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: 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         731 : SpectralIndex::SpectralIndex()
      49             :   :SpectralModel(),
      50             :    itsIndex(0.0),
      51             :    itsStokesIndex(Vector<Double>(1, 0.0)),
      52         731 :    itsError(0.0)
      53             : {
      54         731 :   DebugAssert(ok(), AipsError);
      55         731 : }
      56             : 
      57          21 : SpectralIndex::SpectralIndex(const MFrequency& refFreq, Double exponent)
      58             :   :SpectralModel(refFreq),
      59             :    itsIndex(exponent),
      60             :    itsStokesIndex(Vector<Double>(1, exponent)),
      61          21 :    itsError(0.0)
      62             : {
      63          21 :   DebugAssert(ok(), AipsError);
      64          21 : }
      65             : 
      66        1441 : SpectralIndex::SpectralIndex(const SpectralIndex& other) 
      67             :   :SpectralModel(other),
      68        1441 :    itsIndex(other.itsIndex),
      69        1441 :    itsStokesIndex(other.itsStokesIndex),
      70        1441 :    itsError(other.itsError)
      71             : {
      72        1441 :   DebugAssert(ok(), AipsError);
      73        1441 : }
      74             : 
      75        4369 : SpectralIndex::~SpectralIndex() {
      76        2193 :   DebugAssert(ok(), AipsError);
      77        4369 : }
      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         741 : ComponentType::SpectralShape SpectralIndex::type() const {
      92         741 :   DebugAssert(ok(), AipsError);
      93         741 :   return ComponentType::SPECTRAL_INDEX;
      94             : }
      95             : 
      96         724 : const Double& SpectralIndex::index() const {
      97         724 :    DebugAssert(ok(), AipsError);
      98         724 :    return itsIndex;
      99             : }
     100           0 : const Vector<Double>& SpectralIndex::stokesIndex() const {
     101           0 :   DebugAssert(ok(), AipsError);
     102           0 :   return itsStokesIndex;
     103             : 
     104             : }
     105         714 :   void SpectralIndex::setIndex(const Double& newIndex) { 
     106         714 :   itsIndex = newIndex;
     107         714 :   itsStokesIndex.resize(1);
     108         714 :   itsStokesIndex[0]=newIndex;
     109         714 :   DebugAssert(ok(), AipsError);
     110         714 : }
     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        1028 :   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             :     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             : }
     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             :       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             : }
     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             :     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         704 : void SpectralIndex::sampleStokes(
     214             :     Matrix<Double>& iquv, const Vector<MFrequency::MVType>& frequencies, 
     215             :     const MFrequency::Ref& refFrame
     216             : ) const {
     217         704 :     ThrowIf(
     218             :         iquv.shape() != IPosition(2, frequencies.size(), 4),
     219             :         "Incorrect Matrix shape"
     220             :     );
     221         704 :     const auto nSamples = frequencies.nelements();
     222         704 :     const auto nu0 = refFreqInFrame(refFrame); 
     223         704 :     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      367494 :     for (uInt i=0; i<nSamples; ++i) {
     233      366790 :         nu = frequencies(i).getValue();
     234      366790 :         iquv(i, 0) *= pow(nu/nu0, itsStokesIndex(0));
     235      366790 :         iquv(i, 3) *= pow(nu/nu0, itsStokesIndex(0) + itsStokesIndex(3));
     236      366790 :         iquv(i, 1) *= pow(nu/nu0, itsStokesIndex[0] + itsStokesIndex[1]);
     237      366790 :         iquv(i, 2) *= pow(nu/nu0, itsStokesIndex[0] + itsStokesIndex[1]);
     238             :     }
     239             :     ///RM type rotation of linear pol. 
     240         704 :     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      279886 :         for (uInt i = 0; i < nSamples; i++) {
     247      279472 :             nu = frequencies(i).getValue();
     248      279472 :             Double twoalpha = 2*itsStokesIndex[2]*C::c*C::c*(nu0*nu0-nu*nu)/(nu*nu*nu0*nu0);
     249      279472 :             Double cos2alph = cos(twoalpha); 
     250      279472 :             Double sin2alph = sin(twoalpha);
     251             :             //Q'
     252      279472 :             Double tempQ = cos2alph*iquv(i, 1) - sin2alph*iquv(i, 2);
     253             :             //U'
     254      279472 :             iquv(i, 2) = sin2alph*iquv(i, 1) + cos2alph*iquv(i, 2);
     255      279472 :             iquv(i, 1) = tempQ;
     256             :         }
     257             :     }
     258         704 : }
     259             :    
     260        1441 : SpectralModel* SpectralIndex::clone() const {
     261        1441 :   DebugAssert(ok(), AipsError);
     262        1441 :   SpectralModel* tmpPtr = new SpectralIndex(*this);
     263        1441 :   AlwaysAssert(tmpPtr != 0, AipsError);
     264        1441 :   return tmpPtr;
     265             : }
     266             : 
     267         734 : uInt SpectralIndex::nParameters() const {
     268         734 :   DebugAssert(ok(), AipsError);
     269         734 :   return 1;
     270             : }
     271             : 
     272          10 : void SpectralIndex::setParameters(const Vector<Double>& newSpectralParms) {
     273          10 :   DebugAssert(newSpectralParms.nelements() == nParameters(), AipsError);
     274          10 :   itsIndex = newSpectralParms(0);
     275          10 :   DebugAssert(ok(), AipsError);
     276          10 : }
     277             : 
     278          17 : Vector<Double> SpectralIndex::parameters() const {
     279          17 :   DebugAssert(ok(), AipsError);
     280          17 :   return Vector<Double>(1, itsIndex);
     281             : }
     282             : 
     283         724 : void SpectralIndex::setErrors(const Vector<Double>& newSpectralErrs) {
     284         724 :   DebugAssert(newSpectralErrs.nelements() == nParameters(), AipsError);
     285         724 :   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             :   }
     290         724 :   itsError = newSpectralErrs(0);
     291         724 :   DebugAssert(ok(), AipsError);
     292         724 : }
     293             : 
     294         741 : Vector<Double> SpectralIndex::errors() const {
     295         741 :   DebugAssert(ok(), AipsError);
     296         741 :   return Vector<Double>(1, itsError);
     297             : }
     298             : 
     299         714 : Bool SpectralIndex::fromRecord(String& errorMessage, 
     300             :                                const RecordInterface& record) {
     301         714 :   if (!SpectralModel::fromRecord(errorMessage, record)) return false;
     302         714 :   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         714 :      const RecordFieldId index("index");
     309         714 :      const IPosition shape(1,1);
     310         714 :      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         714 :      switch (record.dataType(index)) {
     316         714 :      case TpDouble:
     317             :      case TpFloat:
     318             :      case TpInt:
     319         714 :        indexVal = record.asDouble(index);
     320         714 :        break;
     321           0 :      default:
     322           0 :        errorMessage += "The 'index' field must be a real number\n";
     323           0 :        return false;
     324             :      }
     325         714 :      setIndex(indexVal);
     326             :   }
     327         714 :   if(record.isDefined("stokesindex")){
     328         714 :     Vector<Double> tempstokes=record.asArrayDouble("stokesindex");
     329         714 :     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         714 :     itsStokesIndex.resize();
     334         714 :     itsStokesIndex=tempstokes;
     335             :   }
     336             : //
     337             :   {
     338         714 :       Vector<Double> errorVals(1, 0.0);
     339         714 :       if (record.isDefined("error")) {
     340         714 :         const RecordFieldId error("error");
     341         714 :         const IPosition shape(1,1);
     342         714 :         if (record.shape(error) != shape) {
     343           0 :           errorMessage += "The 'error' field must be a scalar\n";
     344           0 :           return false;
     345             :         }
     346         714 :         switch (record.dataType(error)) {
     347         714 :         case TpDouble:
     348             :         case TpFloat:
     349             :         case TpInt:
     350         714 :             errorVals[0] = record.asDouble(error);
     351         714 :           break;
     352           0 :         default:
     353           0 :           errorMessage += "The 'error' field must be a real number\n";
     354           0 :           return false;
     355             :         }
     356             :      }
     357             : //
     358         714 :      setErrors(errorVals);
     359             :   }
     360             : //
     361         714 :   DebugAssert(ok(), AipsError);
     362         714 :   return true;
     363             : }
     364             : 
     365         724 : Bool SpectralIndex::toRecord(String& errorMessage,
     366             :                              RecordInterface& record) const {
     367         724 :   DebugAssert(ok(), AipsError);
     368         724 :   if (!SpectralModel::toRecord(errorMessage, record)) return false;
     369         724 :   record.define("index", index());
     370         724 :   if(itsStokesIndex.nelements() != 0)
     371         724 :     record.define("stokesindex", itsStokesIndex);
     372         724 :   record.define("error", errors()(0));
     373         724 :   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             : }
     392             : 
     393       12924 : Bool SpectralIndex::ok() const {
     394       12924 :     if (!SpectralModel::ok()) return false;
     395       12924 :     ThrowIf(
     396             :         refFrequency().getValue().getValue() <= 0.0,
     397             :         "The reference frequency is zero or negative!"
     398             :     ); 
     399       12924 :     ThrowIf(abs(
     400             :         itsIndex) > 100,
     401             :         "The absolute value of the spectral index is greater than 100!"
     402             :     ); 
     403       12924 :     return true;
     404             : }
     405             : 
     406             : }
     407             :