Line data Source code
1 : //# SpectralModel.cc:
2 : //# Copyright (C) 1998,1999,2000,2002,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: SpectralModel.cc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
27 :
28 : #include <components/ComponentModels/SpectralModel.h>
29 : #include <casacore/casa/Containers/Record.h>
30 : #include <casacore/casa/Containers/RecordFieldId.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/measures/Measures/MeasureHolder.h>
35 : #include <casacore/measures/Measures/MFrequency.h>
36 : #include <casacore/measures/Measures/MCFrequency.h>
37 : #include <casacore/measures/Measures/MeasConvert.h>
38 : #include <casacore/casa/Quanta/QuantumHolder.h>
39 : #include <casacore/casa/Quanta/Quantum.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Utilities/DataType.h>
42 : #include <casacore/casa/BasicSL/String.h>
43 : #include <casacore/casa/Logging/LogIO.h>
44 : #include <casacore/casa/Logging/LogOrigin.h>
45 : #include <iostream>
46 :
47 : using namespace casacore;
48 : namespace casa { //# NAMESPACE CASA - BEGIN
49 :
50 404 : SpectralModel::SpectralModel()
51 808 : :itsRefFreq(Quantum<Double>(1, "GHz"), MFrequency::DEFAULT),
52 : itsFreqUnit("GHz"),
53 1212 : itsFreqErr(0, itsFreqUnit)
54 : {
55 404 : DebugAssert(SpectralModel::ok(), AipsError);
56 404 : }
57 :
58 0 : SpectralModel::SpectralModel(const MFrequency& refFreq, const Unit& freqUnit)
59 : :itsRefFreq(refFreq),
60 : itsFreqUnit(freqUnit),
61 0 : itsFreqErr(0, itsFreqUnit)
62 : {
63 0 : DebugAssert(SpectralModel::ok(), AipsError);
64 0 : }
65 :
66 1 : SpectralModel::SpectralModel(const SpectralModel& other)
67 : :RecordTransformable(),
68 1 : itsRefFreq(other.itsRefFreq),
69 1 : itsFreqUnit(other.itsFreqUnit),
70 1 : itsFreqErr(other.itsFreqErr)
71 : {
72 1 : DebugAssert(SpectralModel::ok(), AipsError);
73 1 : }
74 :
75 0 : SpectralModel& SpectralModel::operator=(const SpectralModel& other) {
76 0 : if (this != &other) {
77 0 : itsRefFreq = other.itsRefFreq;
78 0 : itsFreqUnit = other.itsFreqUnit;
79 0 : itsFreqErr = other.itsFreqErr;
80 : }
81 0 : DebugAssert(SpectralModel::ok(), AipsError);
82 0 : return *this;
83 : }
84 :
85 404 : SpectralModel::~SpectralModel() {
86 404 : }
87 :
88 0 : const String& SpectralModel::ident() const {
89 0 : DebugAssert(SpectralModel::ok(), AipsError);
90 0 : static String typeString;
91 0 : typeString = ComponentType::name(type());
92 0 : return typeString;
93 : }
94 :
95 400 : const MFrequency& SpectralModel::refFrequency() const {
96 400 : DebugAssert(SpectralModel::ok(), AipsError);
97 400 : return itsRefFreq;
98 : }
99 :
100 0 : Double SpectralModel::refFreqInFrame(const MFrequency::Ref& elframe) const {
101 0 : const MFrequency& refFreq(refFrequency());
102 : Double nu0;
103 0 : Bool stupidTransform = (elframe.getType() == MFrequency::REST) || (elframe.getType() == MFrequency::N_Types) || (refFreq.getRef().getType() == MFrequency::REST) || (refFreq.getRef().getType() == MFrequency::N_Types);
104 0 : if (elframe != refFreq.getRef() && !stupidTransform) {
105 : try{
106 0 : nu0 = (MFrequency::Convert(refFreq, elframe))().getValue().getValue();
107 : }
108 0 : catch(...){
109 0 : throw(AipsError(String("Cannot transform frequency from ") + MFrequency::showType(elframe.getType()) + String(" to ") + MFrequency::showType(refFreq.getRef().getType())));
110 : }
111 : } else {
112 0 : nu0 = refFreq.getValue().getValue();
113 : }
114 :
115 0 : return nu0;
116 : }
117 :
118 200 : void SpectralModel::setRefFrequency(const MFrequency& newRefFreq) {
119 200 : itsRefFreq = newRefFreq;
120 200 : DebugAssert(SpectralModel::ok(), AipsError);
121 200 : }
122 :
123 200 : const Unit& SpectralModel::frequencyUnit() const {
124 200 : DebugAssert(SpectralModel::ok(), AipsError);
125 200 : return itsFreqUnit;
126 : }
127 :
128 0 : void SpectralModel::convertFrequencyUnit(const Unit& freqUnit) {
129 0 : itsFreqUnit = freqUnit;
130 0 : DebugAssert(SpectralModel::ok(), AipsError);
131 0 : }
132 :
133 0 : void SpectralModel::setRefFrequencyError(const Quantum<Double>& newRefFreqErr){
134 0 : if (badError(newRefFreqErr)) {
135 0 : LogIO logErr(LogOrigin("SpectralModel", "setRefFrequencyError"));
136 : logErr << "The errors must be positive values with units like Hz"
137 0 : << LogIO::EXCEPTION;
138 : }
139 0 : itsFreqErr = newRefFreqErr;
140 0 : DebugAssert(SpectralModel::ok(), AipsError);
141 0 : }
142 :
143 0 : const Quantum<Double>& SpectralModel::refFrequencyError() const {
144 0 : return itsFreqErr;
145 : }
146 :
147 0 : void SpectralModel::sample(Vector<Double>& scale,
148 : const Vector<MFrequency::MVType>& frequencies,
149 : const MFrequency::Ref& refFrame) const {
150 0 : DebugAssert(SpectralModel::ok(), AipsError);
151 0 : const uInt nSamples = frequencies.nelements();
152 0 : DebugAssert(scale.nelements() == nSamples, AipsError);
153 :
154 0 : for (uInt i = 0; i < nSamples; i++) {
155 0 : scale(i) = sample(MFrequency(frequencies(i), refFrame));
156 : }
157 0 : }
158 :
159 200 : Bool SpectralModel::fromRecord(String& errorMessage,
160 : const RecordInterface& record) {
161 400 : const String freqString("frequency");
162 200 : if (!record.isDefined(freqString)) {
163 0 : errorMessage += "The 'frequency' field does not exist\n";
164 0 : return false;
165 : }
166 400 : const RecordFieldId frequency(freqString);
167 200 : if (record.dataType(frequency) != TpRecord) {
168 0 : if ((record.dataType(frequency) == TpString) &&
169 0 : (record.shape(frequency) == IPosition(1,1)) &&
170 0 : (record.asString(frequency) == String("current"))) {
171 0 : return true;
172 : } else {
173 0 : errorMessage += "The 'frequency' field must be a record\n";
174 0 : errorMessage += "or the string 'current' (to use the current value)\n";
175 0 : return false;
176 : }
177 : }
178 400 : const Record& freqRecord = record.asRecord(frequency);
179 400 : MeasureHolder mh;
180 200 : if (!mh.fromRecord(errorMessage, freqRecord)) {
181 0 : errorMessage += "Could not parse the reference frequency\n";
182 0 : return false;
183 : }
184 200 : if (!mh.isMFrequency()) {
185 0 : errorMessage += "The reference frequency is not a frequency measure\n";
186 0 : return false;
187 : }
188 200 : SpectralModel::setRefFrequency(mh.asMFrequency());
189 200 : return true;
190 : }
191 :
192 200 : Bool SpectralModel::toRecord(String& errorMessage,
193 : RecordInterface& record) const {
194 200 : DebugAssert(SpectralModel::ok(), AipsError);
195 200 : record.define(RecordFieldId("type"), ComponentType::name(type()));
196 400 : Record freqRecord;
197 400 : const MeasureHolder mh(refFrequency());
198 200 : if (!mh.toRecord(errorMessage, freqRecord)) {
199 0 : errorMessage += "Could not convert the reference frequency to a record\n";
200 0 : return false;
201 : }
202 200 : const String m0String("m0");
203 200 : if (freqRecord.isDefined(m0String)) {
204 400 : const RecordFieldId m0(m0String);
205 200 : if (freqRecord.dataType(m0) == TpRecord) {
206 400 : Record m0Rec = freqRecord.asRecord(m0);
207 400 : QuantumHolder qh;
208 400 : String error;
209 200 : if (qh.fromRecord(error, m0Rec) && qh.isQuantumDouble()) {
210 400 : Quantum<Double> q = qh.asQuantumDouble();
211 200 : q.convert(frequencyUnit());
212 200 : qh = QuantumHolder(q);
213 200 : if (qh.toRecord(error, m0Rec)) {
214 200 : freqRecord.defineRecord(m0, m0Rec);
215 : }
216 : }
217 : }
218 : }
219 200 : record.defineRecord(RecordFieldId("frequency"), freqRecord);
220 200 : return true;
221 : }
222 :
223 200 : ComponentType::SpectralShape SpectralModel::
224 : getType(String& errorMessage, const RecordInterface& record) {
225 400 : const String typeString("type");
226 200 : if (!record.isDefined(typeString)) {
227 : errorMessage +=
228 0 : String("\nThe record does not have a 'type' field.");
229 0 : return ComponentType::UNKNOWN_SPECTRAL_SHAPE;
230 : }
231 400 : const RecordFieldId type(typeString);
232 200 : if (record.dataType(type) != TpString) {
233 0 : errorMessage += String("\nThe 'type' field in the spectrum record,")
234 0 : + String("must be a String");
235 0 : return ComponentType::UNKNOWN_SPECTRAL_SHAPE;
236 : }
237 200 : if (record.shape(type) != IPosition(1,1)) {
238 0 : errorMessage += String("The 'type' field, in the spectrum record,") +
239 0 : String(" must have only 1 element\n");
240 0 : return ComponentType::UNKNOWN_SPECTRAL_SHAPE;
241 : }
242 200 : const String& typeVal = record.asString(type);
243 200 : return ComponentType::spectralShape(typeVal);
244 : }
245 :
246 3419 : Bool SpectralModel::ok() const {
247 3419 : if (itsFreqUnit != Unit("GHz")) {
248 0 : LogIO logErr(LogOrigin("SpectralModel", "ok()"));
249 : logErr << LogIO::SEVERE << "The reference frequency has units with "
250 : << endl << " different dimensions than the Hz."
251 0 : << LogIO::POST;
252 0 : return false;
253 : }
254 3419 : return true;
255 : }
256 :
257 0 : Bool SpectralModel::badError(const Quantum<Double>& quantum) {
258 0 : const UnitVal hz(1, "Hz");
259 0 : return !(quantum.check(hz)) || (quantum.getValue() < 0.0);
260 : }
261 :
262 : // Local Variables:
263 : // compile-command: "gmake SpectralModel"
264 : // End:
265 :
266 : } //# NAMESPACE CASA - END
267 :
|
