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 0 : TabularSpectrum::TabularSpectrum()
52 : :SpectralModel(),
53 0 : tabFreqVal_p(0),
54 0 : 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 0 : freqRef_p=MFrequency::Ref(MFrequency::LSRK);
57 :
58 0 : DebugAssert(ok(), AipsError);
59 0 : }
60 :
61 0 : TabularSpectrum::TabularSpectrum(const MFrequency& refFreq,
62 : const Vector<MFrequency::MVType>& freq,
63 : const Vector<Flux<Double> >& flux,
64 0 : const MFrequency::Ref& refFrame)
65 0 : :SpectralModel(refFreq)
66 : {
67 :
68 0 : Bool stupidTransform = (refFrame.getType() == MFrequency::REST) || (refFrame.getType() == MFrequency::N_Types) || (refFreq.getRef().getType() == MFrequency::REST) || (refFreq.getRef().getType() == MFrequency::N_Types);
69 :
70 0 : if (refFrame.getType() != refFreq.getRef().getType() && !stupidTransform) {
71 0 : referenceFreq_p = MFrequency::Convert(refFreq, refFrame)().getValue().get("Hz").getValue();
72 : } else {
73 0 : referenceFreq_p = refFreq.getValue().get("Hz").getValue();
74 : }
75 0 : setValues(freq, flux, refFrame);
76 0 : DebugAssert(ok(), AipsError);
77 0 : }
78 :
79 0 : TabularSpectrum::TabularSpectrum(const TabularSpectrum& other)
80 0 : :SpectralModel(other)
81 : {
82 0 : operator=(other);
83 0 : DebugAssert(ok(), AipsError);
84 0 : }
85 :
86 0 : TabularSpectrum::~TabularSpectrum() {
87 0 : DebugAssert(ok(), AipsError);
88 0 : }
89 :
90 0 : TabularSpectrum& TabularSpectrum::operator=(const TabularSpectrum& other) {
91 0 : if (this != &other) {
92 0 : SpectralModel::operator=(other);
93 0 : freqRef_p=other.freqRef_p;
94 0 : tabFreqVal_p.resize();
95 0 : tabFreqVal_p=other.tabFreqVal_p;
96 0 : flux_p.resize();
97 0 : flux_p=other.flux_p;
98 0 : referenceFreq_p=other.referenceFreq_p;
99 0 : maxFreq_p=other.maxFreq_p;
100 0 : minFreq_p=other.minFreq_p;
101 0 : ival_p=other.ival_p;
102 0 : qval_p=other.qval_p;
103 0 : uval_p=other.uval_p;
104 0 : vval_p=other.vval_p;
105 : }
106 0 : DebugAssert(ok(), AipsError);
107 0 : return *this;
108 : }
109 :
110 0 : ComponentType::SpectralShape TabularSpectrum::type() const {
111 0 : 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 0 : void TabularSpectrum::setValues(const Vector<MFrequency::MVType>& frequencies, const Vector<Flux<Double> >& flux, const MFrequency::Ref& refFrame) {
124 0 : if(flux.nelements() != frequencies.nelements()){
125 0 : throw(AipsError("frequencies length is not equal to flux length in TabularSpectrum::setValues"));
126 : }
127 :
128 0 : referenceFreq_p=refFreqInFrame(refFrame);
129 :
130 0 : freqRef_p=refFrame;
131 0 : tabFreqVal_p.resize(frequencies.nelements());
132 0 : flux_p.resize();
133 0 : flux_p=flux;
134 0 : ival_p.resize(frequencies.nelements());
135 0 : qval_p.resize(frequencies.nelements());
136 0 : uval_p.resize(frequencies.nelements());
137 0 : vval_p.resize(frequencies.nelements());
138 :
139 0 : for (uInt k=0; k < frequencies.nelements(); ++k){
140 0 : tabFreqVal_p(k)=frequencies(k).get("Hz").getValue();
141 : //IQUV
142 0 : flux_p(k).convertPol(ComponentType::STOKES);
143 0 : ival_p(k)=flux_p(k).value(Stokes::I).getValue();
144 0 : qval_p(k)=flux_p(k).value(Stokes::Q).getValue();
145 0 : uval_p(k)=flux_p(k).value(Stokes::U).getValue();
146 0 : vval_p(k)=flux_p(k).value(Stokes::V).getValue();
147 : }
148 0 : maxFreq_p=max(tabFreqVal_p);
149 0 : minFreq_p=min(tabFreqVal_p);
150 : //Just make sure the refVal_p is calculated
151 0 : this->setRefFrequency(refFrequency());
152 :
153 0 : }
154 0 : void TabularSpectrum::setRefFrequency(const MFrequency& newRefFreq) {
155 0 : SpectralModel::setRefFrequency(newRefFreq);
156 0 : referenceFreq_p=refFreqInFrame(freqRef_p);
157 0 : Vector<Double> xout(1, referenceFreq_p);
158 0 : Vector<Double> scale(1,0.0);
159 0 : refVal_p.resize(4);
160 0 : refVal_p=0.0;
161 0 : Vector<Vector<Double> > iquv(4);
162 0 : iquv[0].reference(ival_p);
163 0 : iquv[1].reference(qval_p);
164 0 : iquv[2].reference(uval_p);
165 0 : iquv[3].reference(vval_p);
166 0 : 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 0 : for (uInt k=0; k < 4; ++k){
169 0 : InterpolateArray1D<Double, Double>::interpolate(scale, xout, tabFreqVal_p, iquv[k], InterpolateArray1D<Double, Double>::linear);
170 0 : refVal_p[k]=scale[0] != 0.0 ? scale[0] : max(iquv[k]);
171 :
172 : }
173 0 : }
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 0 : void TabularSpectrum::sampleStokes(
275 : Matrix<Double>& retvals, const Vector<MFrequency::MVType>& frequencies,
276 : const MFrequency::Ref& refFrame
277 : ) const {
278 0 : ThrowIf(
279 : retvals.shape() != IPosition(2, frequencies.size(), 4),
280 : "Incorrect Matrix shape"
281 : );
282 0 : const auto nSamples = frequencies.nelements();
283 0 : retvals.set(0.0);
284 0 : MFrequency::Convert toThisFrame(refFrame, freqRef_p);
285 0 : 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 0 : refFrame.getType() == MFrequency::REST)
290 0 : || (refFrame.getType() == MFrequency::N_Types)
291 0 : || (freqRef_p.getType() == MFrequency::REST)
292 0 : || (freqRef_p.getType() == MFrequency::N_Types
293 0 : );
294 0 : 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 0 : for(uInt k=0; k<nSamples; ++k) {
301 0 : 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 0 : Vector<Double> scaleone(nSamples);
310 0 : Vector<Vector<Double> > iquv(4);
311 0 : iquv[0].reference(ival_p);
312 0 : iquv[1].reference(qval_p);
313 0 : iquv[2].reference(uval_p);
314 0 : iquv[3].reference(vval_p);
315 0 : for (uInt k=0; k<4; ++k){
316 0 : InterpolateArray1D<Double, Double>::interpolate(
317 0 : scaleone, nu, tabFreqVal_p, iquv[k],
318 : InterpolateArray1D<Double, Double>::linear
319 : );
320 0 : for (uInt i=0; i<nSamples; ++i) {
321 0 : retvals(i, k) = scaleone[i];
322 : }
323 : }
324 0 : }
325 :
326 0 : SpectralModel* TabularSpectrum::clone() const {
327 0 : DebugAssert(ok(), AipsError);
328 0 : SpectralModel* tmpPtr = new TabularSpectrum(*this);
329 0 : AlwaysAssert(tmpPtr != 0, AipsError);
330 0 : return tmpPtr;
331 : }
332 :
333 0 : uInt TabularSpectrum::nParameters() const {
334 0 : return 0;
335 : }
336 :
337 0 : void TabularSpectrum::setParameters(const Vector<Double>& newSpectralParms) {
338 0 : AlwaysAssert(newSpectralParms.nelements() == nParameters(), AipsError);
339 0 : }
340 :
341 0 : Vector<Double> TabularSpectrum::parameters() const {
342 0 : return Vector<Double>(0);
343 : }
344 :
345 0 : void TabularSpectrum::setErrors(const Vector<Double>& newSpectralErrs) {
346 0 : AlwaysAssert(newSpectralErrs.nelements() == nParameters(), AipsError);
347 0 : }
348 :
349 0 : Vector<Double> TabularSpectrum::errors() const {
350 0 : return Vector<Double>(0);
351 : }
352 :
353 0 : Bool TabularSpectrum::fromRecord(String& errorMessage,
354 : const RecordInterface& record) {
355 0 : if (!SpectralModel::fromRecord(errorMessage, record)) return false;
356 : //freqRef
357 0 : 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 0 : Record theTmpMF(record.asRecord("freqRef"));
363 0 : MeasureHolder mh;
364 0 : if(!mh.fromRecord(errorMessage, theTmpMF))
365 0 : return false;
366 0 : if(mh.isMFrequency())
367 0 : freqRef_p=(mh.asMFrequency().getRef());
368 : else
369 0 : return false;
370 0 : }
371 :
372 :
373 : //tabFreqVal
374 0 : 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 0 : tabFreqVal_p.resize();
380 0 : tabFreqVal_p=Vector<Double> (record.asArrayDouble("tabFreqVal"));
381 : }
382 : ////ival
383 0 : 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 0 : ival_p.resize();
389 0 : ival_p=Vector<Double> (record.asArrayDouble("ival"));
390 :
391 0 : qval_p=record.isDefined(String("qval")) ? Vector<Double> (record.asArrayDouble("qval")) : Vector<Double>(ival_p.nelements(), 0.0);
392 0 : uval_p=record.isDefined(String("uval")) ? Vector<Double> (record.asArrayDouble("uval")) : Vector<Double>(ival_p.nelements(), 0.0);
393 0 : vval_p=record.isDefined(String("vval")) ? Vector<Double> (record.asArrayDouble("vval")) : Vector<Double>(ival_p.nelements(), 0.0);
394 : }
395 :
396 : //referenceFreq
397 0 : 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 0 : referenceFreq_p=record.asDouble("referenceFreq");
403 : }
404 0 : setRefFrequency(MFrequency(Quantity(referenceFreq_p, "Hz"), freqRef_p));
405 : //maxFreq and minFreq
406 0 : 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 0 : maxFreq_p=record.asDouble("maxFreq");
412 0 : minFreq_p=record.asDouble("minFreq");
413 : }
414 :
415 0 : return true;
416 : }
417 :
418 0 : Bool TabularSpectrum::toRecord(String& errorMessage,
419 : RecordInterface& record) const {
420 0 : DebugAssert(ok(), AipsError);
421 0 : if (!SpectralModel::toRecord(errorMessage, record)) return false;
422 : //save frame in a temporary MFrequency object
423 0 : MFrequency tmpMF(Quantity(0, "Hz"), freqRef_p);
424 0 : MeasureHolder mh(tmpMF);
425 0 : Record outRec;
426 0 : if(!mh.toRecord(errorMessage, outRec))
427 0 : return false;
428 0 : record.defineRecord("freqRef",outRec);
429 0 : record.define("tabFreqVal", tabFreqVal_p);
430 0 : record.define("ival", ival_p);
431 0 : record.define("qval", qval_p);
432 0 : record.define("uval", uval_p);
433 0 : record.define("vval", vval_p);
434 0 : record.define("referenceFreq", referenceFreq_p);
435 0 : record.define("maxFreq", maxFreq_p);
436 0 : record.define("minFreq", minFreq_p);
437 0 : return true;
438 0 : }
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 0 : Bool TabularSpectrum::ok() const {
450 0 : if (!SpectralModel::ok()) return false;
451 0 : 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 0 : return true;
458 : }
459 :
460 : // Local Variables:
461 : // compile-command: "gmake SpectralIndex"
462 : // End:
463 :
464 : } //# NAMESPACE CASA - END
465 :
|
