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/PowerLogPoly.h>
29 : #include <casacore/casa/Arrays/Vector.h>
30 : #include <casacore/casa/Arrays/ArrayLogical.h>
31 : #include <casacore/casa/Containers/RecordInterface.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 : #include <casacore/casa/Arrays/IPosition.h>
34 : #include <casacore/casa/Logging/LogIO.h>
35 : #include <casacore/casa/Logging/LogOrigin.h>
36 : #include <casacore/casa/BasicMath/Math.h>
37 : #include <casacore/measures/Measures/MFrequency.h>
38 : #include <casacore/measures/Measures/MCFrequency.h>
39 : #include <casacore/measures/Measures/MeasConvert.h>
40 : #include <casacore/casa/Quanta/MVFrequency.h>
41 : #include <casacore/casa/Quanta/Quantum.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 : #include <casacore/casa/Utilities/DataType.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 :
46 : using namespace casacore;
47 : namespace casa { //# NAMESPACE CASA - BEGIN
48 :
49 185 : PowerLogPoly::PowerLogPoly() : SpectralModel() {}
50 :
51 1 : PowerLogPoly::PowerLogPoly(
52 : const MFrequency& refFreq, const Vector<Double> coeffs
53 1 : ) : SpectralModel(refFreq), _coeffs(coeffs.copy()),
54 2 : _errors(Vector<Double>(coeffs.size(), 0)) {}
55 :
56 4764 : PowerLogPoly::PowerLogPoly(const PowerLogPoly& other)
57 4764 : : SpectralModel(other), _coeffs(other._coeffs.copy()),
58 9528 : _errors(other._errors.copy()) {}
59 :
60 9900 : PowerLogPoly::~PowerLogPoly() {}
61 :
62 0 : PowerLogPoly& PowerLogPoly::operator=(const PowerLogPoly& other) {
63 0 : if (this != &other) {
64 0 : SpectralModel::operator=(other);
65 0 : _coeffs.resize(other._coeffs.size());
66 0 : _coeffs = other._coeffs.copy();
67 0 : _errors.resize(other._errors.size());
68 0 : _errors = other._errors.copy();
69 : }
70 0 : return *this;
71 : }
72 :
73 188 : ComponentType::SpectralShape PowerLogPoly::type() const {
74 188 : return ComponentType::PLP;
75 : }
76 :
77 2181 : Double PowerLogPoly::_getNu0(const MFrequency::Ref& refFrame) const {
78 2181 : const MFrequency& refFreq(refFrequency());
79 2181 : const auto nu0 = refFrame == refFreq.getRef()
80 2181 : ? refFreq.getValue().getValue()
81 2181 : : MFrequency::Convert(refFreq, refFrame)().getValue().getValue();
82 2181 : ThrowIf(nu0 <= 0.0, "the reference frequency is zero or negative");
83 2181 : return nu0;
84 : }
85 :
86 174936 : Double PowerLogPoly::_getIntensityRatio(Double x) const {
87 174936 : const auto logx = log(x);
88 174936 : Double exponent = 0;
89 697743 : for(uInt i=0; i<_coeffs.size(); ++i) {
90 522807 : if (i == 0) {
91 174936 : exponent = _coeffs[i];
92 : }
93 347871 : else if (i == 1) {
94 174936 : exponent += _coeffs[i] * logx;
95 : }
96 172935 : else if (i == 2) {
97 172935 : exponent += _coeffs[i] * logx * logx;
98 : }
99 : else {
100 0 : exponent += _coeffs[i] * pow(logx, i);
101 : }
102 : }
103 174936 : return pow(x, exponent);
104 : }
105 :
106 2000 : Double PowerLogPoly::sample(const MFrequency& centerFreq) const {
107 2000 : const auto x = centerFreq.getValue()/_getNu0(centerFreq.getRef());
108 2000 : return _getIntensityRatio(x);
109 : }
110 :
111 2000 : void PowerLogPoly::sampleStokes(
112 : const MFrequency& centerFreq, Vector<Double>& iquv
113 : ) const {
114 2000 : ThrowIf(iquv.size() != 4, "Four stokes parameters in iquv are expected");
115 2000 : iquv[0] *= sample(centerFreq);
116 : // TODO add full stokes support.
117 2000 : }
118 :
119 0 : void PowerLogPoly::sample(
120 : Vector<Double>& scale, const Vector<MFrequency::MVType>& frequencies,
121 : const MFrequency::Ref& refFrame
122 : ) const {
123 0 : const auto nSamples = frequencies.size();
124 0 : ThrowIf(
125 : scale.size() != nSamples,
126 : "A Vector of length " + String::toString(nSamples) + " is required"
127 : );
128 0 : const auto nu0 = _getNu0(refFrame);
129 0 : for (uInt i=0; i<nSamples; ++i) {
130 0 : scale[i] = _getIntensityRatio(frequencies[i].getValue()/nu0);
131 : }
132 0 : }
133 :
134 181 : void PowerLogPoly::sampleStokes(
135 : Matrix<Double>& iquv, const Vector<MFrequency::MVType>& frequencies,
136 : const MFrequency::Ref& refFrame
137 : ) const {
138 181 : ThrowIf(
139 : iquv.shape() != IPosition(2, frequencies.size(), 4),
140 : "Incorrect Matrix shape"
141 : );
142 181 : const auto nSamples = frequencies.size();
143 181 : const auto nu0 = _getNu0(refFrame);
144 173117 : for (uInt i=0; i<nSamples; ++i) {
145 172936 : iquv(i, 0) *= _getIntensityRatio(frequencies[i].getValue()/nu0);
146 : }
147 : // TODO full polarization implementation to come
148 181 : }
149 :
150 4764 : SpectralModel* PowerLogPoly::clone() const {
151 4764 : SpectralModel* tmpPtr = new PowerLogPoly(*this);
152 4764 : AlwaysAssert(tmpPtr, AipsError);
153 4764 : return tmpPtr;
154 : }
155 :
156 370 : uInt PowerLogPoly::nParameters() const {
157 370 : return _coeffs.size();
158 : }
159 :
160 189 : void PowerLogPoly::setParameters(const Vector<Double>& newSpectralParms) {
161 189 : _coeffs.resize(0);
162 189 : _coeffs = newSpectralParms.copy();
163 189 : const auto s = _coeffs.size();
164 189 : if (_errors.size() != s) {
165 185 : _errors.resize(s, True);
166 : }
167 189 : }
168 :
169 188 : Vector<Double> PowerLogPoly::parameters() const {
170 188 : return _coeffs.copy();
171 : }
172 :
173 189 : void PowerLogPoly::setErrors(const Vector<Double>& newSpectralErrs) {
174 189 : ThrowIf(anyLT(newSpectralErrs, 0.0), "The errors must be non-negative.");
175 189 : _errors.resize(newSpectralErrs.size());
176 189 : _errors = newSpectralErrs.copy();
177 189 : const auto s = _errors.size();
178 189 : if (_coeffs.size() != s) {
179 0 : _coeffs.resize(s, True);
180 : }
181 189 : }
182 :
183 188 : Vector<Double> PowerLogPoly::errors() const {
184 188 : return _errors.copy();
185 : }
186 :
187 185 : Bool PowerLogPoly::fromRecord(
188 : String& errorMessage, const RecordInterface& record) {
189 185 : if (! SpectralModel::fromRecord(errorMessage, record)) {
190 0 : return false;
191 : }
192 185 : if (!record.isDefined(String("coeffs"))) {
193 0 : errorMessage += "The record must have an 'coeffs' field";
194 0 : return false;
195 : }
196 : {
197 185 : const RecordFieldId coeffs("coeffs");
198 185 : Vector<Double> coeffVal;
199 185 : switch (record.dataType(coeffs)) {
200 185 : case TpArrayDouble:
201 : case TpArrayFloat:
202 : case TpArrayInt:
203 185 : coeffVal = record.asArrayDouble(coeffs);
204 185 : break;
205 0 : default:
206 0 : errorMessage += "The 'coeff' field must be an array of numbers";
207 0 : return false;
208 : }
209 185 : setParameters(coeffVal);
210 : }
211 : // TODO full stokes implementation to come
212 : {
213 185 : Vector<Double> errorVals(nParameters(), 0.0);
214 185 : if (record.isDefined("error")) {
215 185 : const RecordFieldId error("error");
216 185 : switch (record.dataType(error)) {
217 185 : case TpArrayDouble:
218 : case TpArrayFloat:
219 : case TpArrayInt: {
220 185 : const auto x = record.asArrayDouble(error);
221 185 : if (x.size() != nParameters()) {
222 0 : errorMessage = "Coefficient and error array lengths are not equal";
223 0 : return false;
224 : }
225 : else {
226 185 : setErrors(x);
227 : }
228 : }
229 185 : break;
230 0 : default:
231 0 : errorMessage += "The 'error' field must be an array of numbers";
232 0 : return false;
233 : }
234 : }
235 : }
236 185 : return true;
237 : }
238 :
239 185 : Bool PowerLogPoly::toRecord(
240 : String& errorMessage, RecordInterface& record) const {
241 185 : if (! SpectralModel::toRecord(errorMessage, record)) {
242 0 : return false;
243 : }
244 185 : record.define("coeffs", parameters());
245 : // TODO full stokes support still to come
246 185 : record.define("error", errors());
247 185 : return true;
248 : }
249 :
250 0 : Bool PowerLogPoly::convertUnit(
251 : String&, const RecordInterface&
252 : ) {
253 : // parameters are dimensionless, no conversion can be done
254 0 : return true;
255 : }
256 :
257 4586 : Bool PowerLogPoly::ok() const {
258 4586 : if (! SpectralModel::ok()) {
259 0 : return false;
260 : }
261 4586 : if (refFrequency().getValue().getValue() <= 0.0) {
262 0 : LogIO logErr(LogOrigin("PowerLogPoly", "ok()"));
263 : logErr << LogIO::SEVERE << "The reference frequency is zero or negative!"
264 0 : << LogIO::POST;
265 0 : return false;
266 : }
267 4586 : if (anyGT(abs(_coeffs), 100.0)) {
268 0 : LogIO logErr(LogOrigin("PowerLogPoly", "ok()"));
269 : logErr << LogIO::SEVERE << "At least one coefficient is greater than 100!"
270 0 : << LogIO::POST;
271 0 : return false;
272 : }
273 4586 : return true;
274 : }
275 :
276 : } //# NAMESPACE CASA - END
|
