00001 //# SpectralIndex.h: Models the spectral variation with a spectral index 00002 //# Copyright (C) 1998-2014 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# $Id: SpectralIndex.h 21229 2012-04-02 12:00:20Z gervandiepen $ 00027 00028 #ifndef COMPONENTS_SPECTRALINDEX_H 00029 #define COMPONENTS_SPECTRALINDEX_H 00030 00031 #include <casa/aips.h> 00032 #include <components/ComponentModels/ComponentType.h> 00033 #include <components/ComponentModels/SpectralModel.h> 00034 00035 namespace casa { //# NAMESPACE CASA - BEGIN 00036 00037 class MFrequency; 00038 class RecordInterface; 00039 class String; 00040 template <class T> class Vector; 00041 00042 // <summary>Models the spectral variation with a spectral index</summary> 00043 00044 // <use visibility=export> 00045 00046 // <reviewed reviewer="" date="yyyy/mm/dd" tests="tSpectralIndex" demos="dSpectralModel"> 00047 // </reviewed> 00048 00049 // <prerequisite> 00050 // <li> <linkto class="SpectralModel">SpectralModel</linkto> 00051 // </prerequisite> 00052 // 00053 // <synopsis> 00054 00055 // This class models the spectral variation of a component with a spectral 00056 // index. 00057 00058 // This class like the other spectral models becomes more useful when used 00059 // through the <linkto class=SkyComponent>SkyComponent</linkto> class, which 00060 // incorporates the flux and spatial variation of the emission, or through the 00061 // <linkto class=ComponentList>ComponentList</linkto> class, which handles 00062 // groups of SkyComponent objects. 00063 00064 // A spectral index is the exponent in a power law model for the variation flux 00065 // with frequency. It is mathematically is defined as: 00066 // <srcblock> 00067 // (nu / nu_0)^alpha 00068 // </srcblock> 00069 // Where: 00070 // <dl compact> 00071 // <dt><src>nu_0</src><dd> is the reference frequency 00072 // <dt><src>alpha</src><dd> is the spectral index 00073 // <dt><src>nu</src><dd> is the user specified frequency 00074 // </dl> 00075 00076 // As with all classes derived from SpectralModel the basic operation of this 00077 // class is to model the flux as a function of frequency. This class does not 00078 // know what the flux is at the reference frequency. Instead the sample 00079 // functions return factors that are used to scale the flux and calculate the 00080 // amount of flux at a specified frequency. 00081 00082 // Besides the reference frequency this class has one parameter; the spectral 00083 // index. This parameter can be set & queried using the general purpose 00084 // <src>parameters</src> functions or the class specific <src>index</src> 00085 // functions. 00086 00087 //<Dec-2013> Have added a full stokes spectral variation ... 00088 // This is done via setting setStokesIndex. 00089 // If setIndex is used then only the first element of 4 parameters is set 00090 // The 4 elements are 00091 // 0) alpha for stokes I ..such that 00092 // 1) alpha for linear pol fraction i.e make sure that Q and U is such that sqrt(Q_0^2+U_0^2)/I_0 * (nu/nu_0)^alpha(1)) is obeyed. 00093 // 2) Rot measure (RM=alpha(2)) value to rotate linear pol by angle RM*(lambda^2- lambda_0^2) 00094 // 3) alpha for circular pol fraction i.e to make V such that V/I=V_0/I_0 *(nu/nu_0)^alpha(3) 00095 //</Dec-2013> 00096 // This class also contains functions (<src>toRecord</src> & 00097 // <src>fromRecord</src>) which perform the conversion between Records and 00098 // SpectralIndex objects. These functions define how a SpectralIndex 00099 // object is represented in glish. The format of the record that is generated 00100 // and accepted by these functions is: 00101 // <srcblock> 00102 // c := [type = 'spectral index', 00103 // frequency = [type = 'frequency', 00104 // refer = 'lsr', 00105 // m0 = [value = 1, unit = 'GHz'] 00106 // ], 00107 // index = 0.7 00108 // ] 00109 // </srcblock> 00110 // The frequency field contains a record representation of a frequency measure 00111 // and its format is defined in the Measures module. Its refer field defines 00112 // the reference frame for the direction and the m0 field defines the value of 00113 // the reference frequency. The parsing of the type field is case 00114 // insensitive. The index field contains the spectral index. 00115 // </synopsis> 00116 00117 // 00118 // <example> 00119 // These examples are coded in the tSpectralModel.h file. 00120 // <h4>Example 1:</h4> 00121 // In this example a SpectralIndex object is created and used to calculate the 00122 // flux at a number of frequencies. 00123 // <srcblock> 00124 //<Dec-2013> The example below was NEVER implemented... 00125 // Could well have written the following and call it documentation: 00126 // Sous un arbre, vos laitues naissent-elles ? 00127 // Si vos laitues naissent, vos navets aussi naissent ! 00128 // Leaving it as is for archeological purposes 00129 // </Dec-2013> 00130 // SpectralIndex siModel; 00131 // siModel.setRefFrequency(MFrequency(Quantity(1.0, "GHz"))); 00132 // siModel.setIndex(1.0, Stokes::I); 00133 // siModel.setIndex(0.5, Stokes::Q); 00134 // siModel.setIndex(0.5, Stokes::U); 00135 // siModel.setIndex(-1.0, Stokes::V); 00136 // const Flux<Double> LBandFlux(1.0, 1.0, 1.0, 1.0); 00137 // const MVFrequency step(Quantity(100.0, "MHz")); 00138 // MVFrequency sampleFreq = siModel.refFrequency().getValue(); 00139 // Flux<Double> sampleFlux; 00140 // cout << "Frequency\t I-Flux\t Q-Flux\t U-Flux\t V-Flux\n"; 00141 // for (uInt i = 0; i < 11; i++) { 00142 // sampleFlux = LBandFlux.copy(); 00143 // sampleFlux.convertPol(ComponentType::LINEAR); 00144 // sampleFlux.convertUnit(Unit("WU")); 00145 // siModel.sample(sampleFlux, 00146 // MFrequency(sampleFreq, siModel.refFrequency().getRef())); 00147 // cout << setprecision(3) << sampleFreq.get("GHz") 00148 // << "\t\t " << sampleFlux.value(0u).re 00149 // << "\t " << sampleFlux.value(1u).re 00150 // << "\t " << sampleFlux.value(2u).re 00151 // << "\t " << sampleFlux.value(3u).re 00152 // << " " << sampleFlux.unit().getName() << endl; 00153 // sampleFreq += step; 00154 // } 00155 //<Dec-2013> 00156 // Now for an example 00158 // const MFrequency f1(Quantity(1.0, "GHz"), MFrequency::LSRK); 00159 // const MFrequency f2(Quantity(2.0, "GHz"), MFrequency::LSRK); 00160 // SpectralIndex siModel; 00161 // siModel.setIndex(1.0); 00162 // cout << "scale value at 1 GHz for setIndex 1.0 " << siModel.sample(f1) << endl; 00163 // Vector<Double> indices(4); 00164 // indices(0)=1.0; indices(1)=0.2; indices(2)=0.0005; indices(3)=0.1; 00165 // siModel.setStokesIndex(indices); 00166 // Vector<Double> iquv(4); 00167 // iquv(0)=10.0; iquv(1)=0.2; iquv(2)=0.4; iquv(3)=0.1; 00168 // cerr << "iquv in " << iquv << " indices " << indices << endl; 00169 // siModel.sampleStokes(f1, iquv); 00170 // cerr << "scale value of I at 1.0 GHz " << siModel.sample(f1) << " iquv out " << iquv << endl; 00171 // siModel.sampleStokes(f2, iquv); 00172 // cerr << "scale value of I at 2.0 GHz " << siModel.sample(f2) << " iquv out " << iquv << endl; 00173 //</Dec-2013> 00174 // </srcblock> 00175 // </example> 00176 // 00177 // <motivation> A Spectral Index frequency variation is the most widely used 00178 // model in radio astronomy. In particular the NFRA package 'newstar' uses it 00179 // extensively. 00180 // </motivation> 00181 // 00182 // <todo asof="1999/11/23"> 00183 // <li> Nothing I hope 00184 // </todo> 00185 00186 // <linkfrom anchor="SpectralIndex" classes="SpectralModel ConstantSpectrum"> 00187 // <here>SpectralIndex</here> - Uses a spectral index to model the spectrum 00188 // </linkfrom> 00189 00190 class SpectralIndex: public SpectralModel 00191 { 00192 public: 00193 // The default SpectralIndex has a reference frequency of 1 GHz in the LSR 00194 // frame and a spectral index of zero. As such it is no different from the 00195 // ConstantSpectrum class (except slower). 00196 SpectralIndex(); 00197 00198 // Construct a SpectralIndex with specified reference frequency and 00199 // exponent. 00200 SpectralIndex(const MFrequency& refFreq, Double exponent = 0.0); 00201 00202 // The copy constructor uses copy semantics 00203 SpectralIndex(const SpectralIndex& other); 00204 00205 // The destructor does nothing special. 00206 virtual ~SpectralIndex(); 00207 00208 // The assignment operator uses copy semantics. 00209 SpectralIndex& operator=(const SpectralIndex& other); 00210 00211 // return the actual spectral type ie., ComponentType::SPECTRAL_INDEX 00212 virtual ComponentType::SpectralShape type() const; 00213 00214 // set/get the spectral index. 00215 // <group> 00216 const Double& index() const; 00217 void setIndex(const Double& newIndex); 00218 const Vector<Double>& stokesIndex() const; 00219 void setStokesIndex(const Vector<Double>& newIndex); 00220 // </group> 00221 00222 // Return the scaling factor that indicates what proportion of the flux is at 00223 // the specified frequency. ie. if the centreFrequency argument is the 00224 // reference frequency then this function will always return one. At other 00225 // frequencies it will return a non-negative number. 00226 virtual Double sample(const MFrequency& centerFrequency) const; 00227 00228 virtual void sampleStokes(const MFrequency& centerFrequency, Vector<Double>& iquv) const; 00229 // Same as the previous function except that many frequencies can be sampled 00230 // at once. The reference frame must be the same for all the specified 00231 // frequencies. Uses a customised implementation for improved speed. 00232 virtual void sample(Vector<Double>& scale, 00233 const Vector<MFrequency::MVType>& frequencies, 00234 const MFrequency::Ref& refFrame) const; 00235 00236 virtual void sampleStokes(Vector<Vector<Double> >& scale, 00237 const Vector<MFrequency::MVType>& frequencies, 00238 const MFrequency::Ref& refFrame) const; 00239 00240 // Return a pointer to a copy of this object upcast to a SpectralModel 00241 // object. The class that uses this function is responsible for deleting the 00242 // pointer. This is used to implement a virtual copy constructor. 00243 virtual SpectralModel* clone() const; 00244 00245 // return the number of parameters. There is one parameter or 4 for this spectral 00246 // model, namely the spectral index for I or I,Q,U,V. So you supply a unit length vector 00247 // or one 4 element long when 00248 // using these functions. Otherwise an exception (AipsError) may be thrown. 00249 // <group> 00250 virtual uInt nParameters() const; 00251 virtual void setParameters(const Vector<Double>& newSpectralParms); 00252 virtual Vector<Double> parameters() const; 00253 virtual void setErrors(const Vector<Double>& newSpectralErrs); 00254 virtual Vector<Double> errors() const; 00255 // </group> 00256 00257 // These functions convert between a Record and a SpectralIndex. These 00258 // functions define how a SpectralIndex object is represented in glish and 00259 // this is detailed in the synopsis above. These functions return False if 00260 // the record is malformed and append an error message to the supplied string 00261 // giving the reason. 00262 // <group> 00263 virtual Bool fromRecord(String& errorMessage, const RecordInterface& record); 00264 virtual Bool toRecord(String& errorMessage, RecordInterface& record) const; 00265 // </group> 00266 00267 // Convert the parameters of the spectral index object to the specified 00268 // units. Only one field of the supplied record is used, namely 'index'. This 00269 // field is optional as the spectral index is a unitless quantity. If the 00270 // index field is specified it must have the empty string as its value. This 00271 // function always returns True unless the index field is specified and does 00272 // not contain an empty string. 00273 virtual Bool convertUnit(String& errorMessage, 00274 const RecordInterface& record); 00275 00276 // Function which checks the internal data of this class for consistant 00277 // values. Returns True if everything is fine otherwise returns False. 00278 virtual Bool ok() const; 00279 00280 private: 00281 Double itsIndex; 00282 Vector<Double> itsStokesIndex; 00283 Double itsError; 00284 }; 00285 00286 } //# NAMESPACE CASA - END 00287 00288 #endif