Line data Source code
1 : //# TabularSpectrum.cc:
2 : //# Copyright (C) 2010
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: TabularSpectrum.cc 21292 2012-11-28 14:58:19Z gervandiepen $
27 :
28 : #include <components/ComponentModels/TabularSpectrum.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/measures/Measures/MeasureHolder.h>
40 : #include <casacore/casa/Containers/Record.h>
41 : #include <casacore/casa/Quanta/MVFrequency.h>
42 : #include <casacore/casa/Quanta/Quantum.h>
43 : #include <casacore/casa/Utilities/Assert.h>
44 : #include <casacore/casa/Utilities/DataType.h>
45 : #include <casacore/casa/BasicSL/String.h>
46 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
47 :
48 : using namespace casacore;
49 : namespace casa { //# NAMESPACE CASA - BEGIN
50 :
51 33802 : TabularSpectrum::TabularSpectrum()
52 : :SpectralModel(),
53 33802 : tabFreqVal_p(0),
54 67604 : flux_p(0), ival_p(0),qval_p(0), uval_p(0), vval_p(0), referenceFreq_p(0.0), maxFreq_p(0.0), minFreq_p(0.0)
55 : {
56 33802 : freqRef_p=MFrequency::Ref(MFrequency::LSRK);
57 :
58 33802 : DebugAssert(ok(), AipsError);
59 33802 : }
60 :
61 70 : TabularSpectrum::TabularSpectrum(const MFrequency& refFreq,
62 : const Vector<MFrequency::MVType>& freq,
63 : const Vector<Flux<Double> >& flux,
64 70 : const MFrequency::Ref& refFrame)
65 70 : :SpectralModel(refFreq)
66 : {
67 :
68 70 : Bool stupidTransform = (refFrame.getType() == MFrequency::REST) || (refFrame.getType() == MFrequency::N_Types) || (refFreq.getRef().getType() == MFrequency::REST) || (refFreq.getRef().getType() == MFrequency::N_Types);
69 :
70 70 : if (refFrame.getType() != refFreq.getRef().getType() && !stupidTransform) {
71 0 : referenceFreq_p = MFrequency::Convert(refFreq, refFrame)().getValue().get("Hz").getValue();
72 : } else {
73 70 : referenceFreq_p = refFreq.getValue().get("Hz").getValue();
74 : }
75 70 : setValues(freq, flux, refFrame);
76 70 : DebugAssert(ok(), AipsError);
77 70 : }
78 :
79 37710 : TabularSpectrum::TabularSpectrum(const TabularSpectrum& other)
80 37710 : :SpectralModel(other)
81 : {
82 37710 : operator=(other);
83 37710 : DebugAssert(ok(), AipsError);
84 37710 : }
85 :
86 143111 : TabularSpectrum::~TabularSpectrum() {
87 71582 : DebugAssert(ok(), AipsError);
88 143111 : }
89 :
90 37710 : TabularSpectrum& TabularSpectrum::operator=(const TabularSpectrum& other) {
91 37710 : if (this != &other) {
92 37710 : SpectralModel::operator=(other);
93 37710 : freqRef_p=other.freqRef_p;
94 37710 : tabFreqVal_p.resize();
95 37710 : tabFreqVal_p=other.tabFreqVal_p;
96 37710 : flux_p.resize();
97 37710 : flux_p=other.flux_p;
98 37710 : referenceFreq_p=other.referenceFreq_p;
99 37710 : maxFreq_p=other.maxFreq_p;
100 37710 : minFreq_p=other.minFreq_p;
101 37710 : ival_p=other.ival_p;
102 37710 : qval_p=other.qval_p;
103 37710 : uval_p=other.uval_p;
104 37710 : vval_p=other.vval_p;
105 : }
106 37710 : DebugAssert(ok(), AipsError);
107 37710 : return *this;
108 : }
109 :
110 4422 : ComponentType::SpectralShape TabularSpectrum::type() const {
111 4422 : return ComponentType::TABULAR_SPECTRUM;
112 : }
113 :
114 0 : void TabularSpectrum::values(Vector<MFrequency::MVType>& freq, Vector<Flux<Double> >& flux) const {
115 0 : freq.resize(tabFreqVal_p.nelements());
116 0 : flux.resize(flux_p.nelements());
117 0 : flux=flux_p;
118 0 : for (uInt k=0; k < tabFreqVal_p.nelements(); ++k){
119 0 : freq(k)=MVFrequency(Quantity(tabFreqVal_p(k), "Hz"));
120 : }
121 0 : }
122 :
123 312 : void TabularSpectrum::setValues(const Vector<MFrequency::MVType>& frequencies, const Vector<Flux<Double> >& flux, const MFrequency::Ref& refFrame) {
124 312 : if(flux.nelements() != frequencies.nelements()){
125 0 : throw(AipsError("frequencies length is not equal to flux length in TabularSpectrum::setValues"));
126 : }
127 :
128 312 : referenceFreq_p=refFreqInFrame(refFrame);
129 :
130 312 : freqRef_p=refFrame;
131 312 : tabFreqVal_p.resize(frequencies.nelements());
132 312 : flux_p.resize();
133 312 : flux_p=flux;
134 312 : ival_p.resize(frequencies.nelements());
135 312 : qval_p.resize(frequencies.nelements());
136 312 : uval_p.resize(frequencies.nelements());
137 312 : vval_p.resize(frequencies.nelements());
138 :
139 43921 : for (uInt k=0; k < frequencies.nelements(); ++k){
140 43609 : tabFreqVal_p(k)=frequencies(k).get("Hz").getValue();
141 : //IQUV
142 43609 : flux_p(k).convertPol(ComponentType::STOKES);
143 43609 : ival_p(k)=flux_p(k).value(Stokes::I).getValue();
144 43609 : qval_p(k)=flux_p(k).value(Stokes::Q).getValue();
145 43609 : uval_p(k)=flux_p(k).value(Stokes::U).getValue();
146 43609 : vval_p(k)=flux_p(k).value(Stokes::V).getValue();
147 : }
148 312 : maxFreq_p=max(tabFreqVal_p);
149 312 : minFreq_p=min(tabFreqVal_p);
150 : //Just make sure the refVal_p is calculated
151 312 : this->setRefFrequency(refFrequency());
152 :
153 312 : }
154 34266 : void TabularSpectrum::setRefFrequency(const MFrequency& newRefFreq) {
155 34266 : SpectralModel::setRefFrequency(newRefFreq);
156 34266 : referenceFreq_p=refFreqInFrame(freqRef_p);
157 34266 : Vector<Double> xout(1, referenceFreq_p);
158 34266 : Vector<Double> scale(1,0.0);
159 34266 : refVal_p.resize(4);
160 34266 : refVal_p=0.0;
161 34266 : Vector<Vector<Double> > iquv(4);
162 34266 : iquv[0].reference(ival_p);
163 34266 : iquv[1].reference(qval_p);
164 34266 : iquv[2].reference(uval_p);
165 34266 : iquv[3].reference(vval_p);
166 34266 : if(ival_p.nelements() < 1 || tabFreqVal_p.nelements() != ival_p.nelements())
167 0 : throw(AipsError("Values have to be set before referenceFrequency in TabularSpectrum"));
168 171330 : for (uInt k=0; k < 4; ++k){
169 137064 : InterpolateArray1D<Double, Double>::interpolate(scale, xout, tabFreqVal_p, iquv[k], InterpolateArray1D<Double, Double>::linear);
170 137064 : refVal_p[k]=scale[0] != 0.0 ? scale[0] : max(iquv[k]);
171 :
172 : }
173 34266 : }
174 :
175 0 : Double TabularSpectrum::sample(const MFrequency& centerFreq) const {
176 0 : const MFrequency& refFreq(refFrequency());
177 0 : const MFrequency::Ref& centerFreqFrame(centerFreq.getRef());
178 : Double nu;
179 :
180 0 : Bool stupidTransform = (centerFreqFrame.getType() == MFrequency::REST) || (centerFreqFrame.getType() == MFrequency::N_Types) || (freqRef_p.getType() == MFrequency::REST) || (freqRef_p.getType() == MFrequency::N_Types);
181 0 : if (centerFreqFrame.getType() != freqRef_p.getType() && !stupidTransform) {
182 0 : nu = MFrequency::Convert(centerFreq, freqRef_p)().getValue().get("Hz").getValue();
183 : } else {
184 0 : nu = refFreq.getValue().get("Hz").getValue();
185 : }
186 0 : if (nu < minFreq_p || nu > maxFreq_p) {
187 0 : throw(AipsError("TabularSpectrun::sample(...) - "
188 0 : "the frequency requested out of range"));
189 : }
190 :
191 0 : Vector<Double> xout(1, referenceFreq_p);
192 0 : Vector<Double> scale(1,0.0);
193 0 : Double refy=refVal_p[0];
194 0 : xout[0]=nu;
195 0 : InterpolateArray1D<Double, Double>::interpolate(scale, xout, tabFreqVal_p, ival_p, InterpolateArray1D<Double, Double>::linear);
196 :
197 :
198 0 : if(refy !=0.0){
199 0 : return scale[0]/refy;
200 : }
201 :
202 0 : return 0.0 ;
203 0 : }
204 :
205 0 : void TabularSpectrum::sampleStokes(const MFrequency& centerFreq, Vector<Double>& retval) const {
206 0 : const MFrequency& refFreq(refFrequency());
207 0 : const MFrequency::Ref& centerFreqFrame(centerFreq.getRef());
208 : Double nu;
209 0 : retval.resize(4);
210 0 : retval.set(0.0);
211 0 : Bool stupidTransform = (centerFreqFrame.getType() == MFrequency::REST) || (centerFreqFrame.getType() == MFrequency::N_Types) || (freqRef_p.getType() == MFrequency::REST) || (freqRef_p.getType() == MFrequency::N_Types);
212 0 : if (centerFreqFrame.getType() != freqRef_p.getType() && !stupidTransform) {
213 0 : nu = MFrequency::Convert(centerFreq, freqRef_p)().getValue().get("Hz").getValue();
214 : } else {
215 0 : nu = refFreq.getValue().get("Hz").getValue();
216 : }
217 0 : if (nu < minFreq_p || nu > maxFreq_p) {
218 0 : throw(AipsError("TabularSpectrun::sample(...) - "
219 0 : "the frequency requested out of range"));
220 : }
221 :
222 0 : Vector<Double> xout(1, referenceFreq_p);
223 0 : Vector<Double> scale(1,0.0);
224 0 : xout[0]=nu;
225 0 : Vector<Vector<Double> > iquv(4);
226 0 : iquv[0].reference(ival_p);
227 0 : iquv[1].reference(qval_p);
228 0 : iquv[2].reference(uval_p);
229 0 : iquv[3].reference(vval_p);
230 0 : for (uInt k=0; k < 4; ++k){
231 0 : InterpolateArray1D<Double, Double>::interpolate(scale, xout, tabFreqVal_p, iquv[k], InterpolateArray1D<Double, Double>::linear);
232 0 : retval(k)=scale(0);
233 : }
234 :
235 0 : }
236 :
237 0 : void TabularSpectrum::sample(Vector<Double>& scale,
238 : const Vector<MFrequency::MVType>& frequencies,
239 : const MFrequency::Ref& refFrame) const {
240 0 : const uInt nSamples = frequencies.nelements();
241 :
242 0 : MFrequency::Convert toThisFrame(refFrame, freqRef_p);
243 0 : Vector<Double> nu(frequencies.nelements());
244 : //try frame conversion only if it is not something stupid...
245 : //if it is then assume the frequencies are fine as is.
246 0 : Bool stupidTransform = (refFrame.getType() == MFrequency::REST) || (refFrame.getType() == MFrequency::N_Types) || (freqRef_p.getType() == MFrequency::REST) || (freqRef_p.getType() == MFrequency::N_Types);
247 0 : if ((refFrame.getType() != freqRef_p.getType()) && !stupidTransform) {
248 0 : for(uInt k=0; k < nSamples; ++k){
249 0 : nu(k) = toThisFrame(frequencies(k).getValue()).getValue().getValue();
250 : }
251 : } else {
252 0 : for(uInt k=0; k< nSamples; ++k){
253 0 : nu(k) = frequencies(k).getValue();
254 : }
255 : }
256 : /* Vector<Double> xout(1, referenceFreq_p);
257 : Vector<Double> refVal(1,0.0);
258 : InterpolateArray1D<Double, Double>::interpolate(refVal, xout, tabFreqVal_p, ival_p, InterpolateArray1D<Double, Double>::linear);
259 : scale.resize(nSamples);
260 : */
261 0 : InterpolateArray1D<Double, Double>::interpolate(scale, nu, tabFreqVal_p, ival_p, InterpolateArray1D<Double, Double>::linear);
262 0 : if(refVal_p(0) !=0.0){
263 0 : for (uInt i = 0; i < nSamples; i++) {
264 0 : scale(i) = scale(i)/refVal_p(0);
265 : }
266 : }
267 : else{
268 0 : if(max(scale) != 0.0)
269 0 : scale /= max(scale);
270 : }
271 :
272 0 : }
273 :
274 33446 : void TabularSpectrum::sampleStokes(
275 : Matrix<Double>& retvals, const Vector<MFrequency::MVType>& frequencies,
276 : const MFrequency::Ref& refFrame
277 : ) const {
278 33446 : ThrowIf(
279 : retvals.shape() != IPosition(2, frequencies.size(), 4),
280 : "Incorrect Matrix shape"
281 : );
282 33446 : const auto nSamples = frequencies.nelements();
283 33446 : retvals.set(0.0);
284 33446 : MFrequency::Convert toThisFrame(refFrame, freqRef_p);
285 33446 : Vector<Double> nu(frequencies.nelements());
286 : //try frame conversion only if it is not something stupid...
287 : //if it is then assume the frequencies are fine as is.
288 : Bool stupidTransform = (
289 33446 : refFrame.getType() == MFrequency::REST)
290 33446 : || (refFrame.getType() == MFrequency::N_Types)
291 33446 : || (freqRef_p.getType() == MFrequency::REST)
292 66892 : || (freqRef_p.getType() == MFrequency::N_Types
293 33446 : );
294 33446 : if ((refFrame.getType() != freqRef_p.getType()) && !stupidTransform) {
295 0 : for(uInt k=0; k < nSamples; ++k){
296 0 : nu(k) = toThisFrame(frequencies(k).getValue()).getValue().getValue();
297 : }
298 : }
299 : else {
300 13232678 : for(uInt k=0; k<nSamples; ++k) {
301 13199232 : nu(k) = frequencies(k).getValue();
302 : }
303 : }
304 : /* Vector<Double> xout(1, referenceFreq_p);
305 : Vector<Double> refVal(1,0.0);
306 : InterpolateArray1D<Double, Double>::interpolate(refVal, xout, tabFreqVal_p, ival_p, InterpolateArray1D<Double, Double>::linear);
307 : scale.resize(nSamples);
308 : */
309 33446 : Vector<Double> scaleone(nSamples);
310 33446 : Vector<Vector<Double> > iquv(4);
311 33446 : iquv[0].reference(ival_p);
312 33446 : iquv[1].reference(qval_p);
313 33446 : iquv[2].reference(uval_p);
314 33446 : iquv[3].reference(vval_p);
315 167230 : for (uInt k=0; k<4; ++k){
316 133784 : InterpolateArray1D<Double, Double>::interpolate(
317 133784 : scaleone, nu, tabFreqVal_p, iquv[k],
318 : InterpolateArray1D<Double, Double>::linear
319 : );
320 52930712 : for (uInt i=0; i<nSamples; ++i) {
321 52796928 : retvals(i, k) = scaleone[i];
322 : }
323 : }
324 33446 : }
325 :
326 37710 : SpectralModel* TabularSpectrum::clone() const {
327 37710 : DebugAssert(ok(), AipsError);
328 37710 : SpectralModel* tmpPtr = new TabularSpectrum(*this);
329 37710 : AlwaysAssert(tmpPtr != 0, AipsError);
330 37710 : return tmpPtr;
331 : }
332 :
333 164 : uInt TabularSpectrum::nParameters() const {
334 164 : return 0;
335 : }
336 :
337 82 : void TabularSpectrum::setParameters(const Vector<Double>& newSpectralParms) {
338 82 : AlwaysAssert(newSpectralParms.nelements() == nParameters(), AipsError);
339 82 : }
340 :
341 103 : Vector<Double> TabularSpectrum::parameters() const {
342 103 : return Vector<Double>(0);
343 : }
344 :
345 82 : void TabularSpectrum::setErrors(const Vector<Double>& newSpectralErrs) {
346 82 : AlwaysAssert(newSpectralErrs.nelements() == nParameters(), AipsError);
347 82 : }
348 :
349 103 : Vector<Double> TabularSpectrum::errors() const {
350 103 : return Vector<Double>(0);
351 : }
352 :
353 33560 : Bool TabularSpectrum::fromRecord(String& errorMessage,
354 : const RecordInterface& record) {
355 33560 : if (!SpectralModel::fromRecord(errorMessage, record)) return false;
356 : //freqRef
357 33560 : if (!record.isDefined(String("freqRef"))) {
358 0 : errorMessage += "The 'TabularSpectrum' record must have an 'freqRef' field\n";
359 0 : return false;
360 : }
361 : else{
362 33560 : Record theTmpMF(record.asRecord("freqRef"));
363 33560 : MeasureHolder mh;
364 33560 : if(!mh.fromRecord(errorMessage, theTmpMF))
365 0 : return false;
366 33560 : if(mh.isMFrequency())
367 33560 : freqRef_p=(mh.asMFrequency().getRef());
368 : else
369 0 : return false;
370 33560 : }
371 :
372 :
373 : //tabFreqVal
374 33560 : if (!record.isDefined(String("tabFreqVal"))) {
375 0 : errorMessage += "The 'TabularSpectrum' record must have an 'tabFreqVal' field\n";
376 0 : return false;
377 : }
378 : else{
379 33560 : tabFreqVal_p.resize();
380 33560 : tabFreqVal_p=Vector<Double> (record.asArrayDouble("tabFreqVal"));
381 : }
382 : ////ival
383 33560 : if (!record.isDefined(String("ival"))) {
384 0 : errorMessage += "The 'TabularSpectrum' record must have an 'ival' field\n";
385 0 : return false;
386 : }
387 : else{
388 33560 : ival_p.resize();
389 33560 : ival_p=Vector<Double> (record.asArrayDouble("ival"));
390 :
391 33560 : qval_p=record.isDefined(String("qval")) ? Vector<Double> (record.asArrayDouble("qval")) : Vector<Double>(ival_p.nelements(), 0.0);
392 33560 : uval_p=record.isDefined(String("uval")) ? Vector<Double> (record.asArrayDouble("uval")) : Vector<Double>(ival_p.nelements(), 0.0);
393 33560 : vval_p=record.isDefined(String("vval")) ? Vector<Double> (record.asArrayDouble("vval")) : Vector<Double>(ival_p.nelements(), 0.0);
394 : }
395 :
396 : //referenceFreq
397 33560 : if (!record.isDefined(String("referenceFreq"))) {
398 0 : errorMessage += "The 'TabularSpectrum' record must have an 'referenceFreq' field\n";
399 0 : return false;
400 : }
401 : else{
402 33560 : referenceFreq_p=record.asDouble("referenceFreq");
403 : }
404 33560 : setRefFrequency(MFrequency(Quantity(referenceFreq_p, "Hz"), freqRef_p));
405 : //maxFreq and minFreq
406 33560 : if (!record.isDefined(String("maxFreq")) || !record.isDefined(String("minFreq")) ) {
407 0 : errorMessage += "The 'TabularSpectrum' record must have a 'maxFreq' and a 'minFreq' field\n";
408 0 : return false;
409 : }
410 : else{
411 33560 : maxFreq_p=record.asDouble("maxFreq");
412 33560 : minFreq_p=record.asDouble("minFreq");
413 : }
414 :
415 33560 : return true;
416 : }
417 :
418 4319 : Bool TabularSpectrum::toRecord(String& errorMessage,
419 : RecordInterface& record) const {
420 4319 : DebugAssert(ok(), AipsError);
421 4319 : if (!SpectralModel::toRecord(errorMessage, record)) return false;
422 : //save frame in a temporary MFrequency object
423 8638 : MFrequency tmpMF(Quantity(0, "Hz"), freqRef_p);
424 4319 : MeasureHolder mh(tmpMF);
425 4319 : Record outRec;
426 4319 : if(!mh.toRecord(errorMessage, outRec))
427 0 : return false;
428 4319 : record.defineRecord("freqRef",outRec);
429 4319 : record.define("tabFreqVal", tabFreqVal_p);
430 4319 : record.define("ival", ival_p);
431 4319 : record.define("qval", qval_p);
432 4319 : record.define("uval", uval_p);
433 4319 : record.define("vval", vval_p);
434 4319 : record.define("referenceFreq", referenceFreq_p);
435 4319 : record.define("maxFreq", maxFreq_p);
436 4319 : record.define("minFreq", minFreq_p);
437 4319 : return true;
438 4319 : }
439 :
440 0 : Bool TabularSpectrum::convertUnit(String& errorMessage,
441 : const RecordInterface& record) {
442 0 : if (!record.isDefined("freqRef")){
443 0 : errorMessage+="Not a tabularSpectrum object";
444 0 : return false;
445 : }
446 0 : return true;
447 : }
448 :
449 226958 : Bool TabularSpectrum::ok() const {
450 226958 : if (!SpectralModel::ok()) return false;
451 226958 : if (refFrequency().getValue().getValue() <= 0.0) {
452 0 : LogIO logErr(LogOrigin("TabularSpectrum", "ok()"));
453 : logErr << LogIO::SEVERE << "The reference frequency is zero or negative!"
454 0 : << LogIO::POST;
455 0 : return false;
456 0 : }
457 226958 : return true;
458 : }
459 :
460 : // Local Variables:
461 : // compile-command: "gmake SpectralIndex"
462 : // End:
463 :
464 : } //# NAMESPACE CASA - END
465 :
|
