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 16249 : SpectralModel::SpectralModel()
51 32498 : :itsRefFreq(Quantum<Double>(1, "GHz"), MFrequency::DEFAULT),
52 : itsFreqUnit("GHz"),
53 48747 : itsFreqErr(0, itsFreqUnit)
54 : {
55 16249 : DebugAssert(SpectralModel::ok(), AipsError);
56 16249 : }
57 :
58 92 : SpectralModel::SpectralModel(const MFrequency& refFreq, const Unit& freqUnit)
59 : :itsRefFreq(refFreq),
60 : itsFreqUnit(freqUnit),
61 92 : itsFreqErr(0, itsFreqUnit)
62 : {
63 92 : DebugAssert(SpectralModel::ok(), AipsError);
64 92 : }
65 :
66 17131 : SpectralModel::SpectralModel(const SpectralModel& other)
67 : :RecordTransformable(),
68 17131 : itsRefFreq(other.itsRefFreq),
69 17131 : itsFreqUnit(other.itsFreqUnit),
70 17131 : itsFreqErr(other.itsFreqErr)
71 : {
72 17131 : DebugAssert(SpectralModel::ok(), AipsError);
73 17131 : }
74 :
75 8430 : SpectralModel& SpectralModel::operator=(const SpectralModel& other) {
76 8430 : if (this != &other) {
77 8430 : itsRefFreq = other.itsRefFreq;
78 8430 : itsFreqUnit = other.itsFreqUnit;
79 8430 : itsFreqErr = other.itsFreqErr;
80 : }
81 8430 : DebugAssert(SpectralModel::ok(), AipsError);
82 8430 : return *this;
83 : }
84 :
85 33472 : SpectralModel::~SpectralModel() {
86 33472 : }
87 :
88 385 : const String& SpectralModel::ident() const {
89 385 : DebugAssert(SpectralModel::ok(), AipsError);
90 385 : static String typeString;
91 385 : typeString = ComponentType::name(type());
92 385 : return typeString;
93 : }
94 :
95 91197 : const MFrequency& SpectralModel::refFrequency() const {
96 91197 : DebugAssert(SpectralModel::ok(), AipsError);
97 91197 : return itsRefFreq;
98 : }
99 :
100 6002 : Double SpectralModel::refFreqInFrame(const MFrequency::Ref& elframe) const {
101 6002 : const MFrequency& refFreq(refFrequency());
102 : Double nu0;
103 6002 : Bool stupidTransform = (elframe.getType() == MFrequency::REST) || (elframe.getType() == MFrequency::N_Types) || (refFreq.getRef().getType() == MFrequency::REST) || (refFreq.getRef().getType() == MFrequency::N_Types);
104 6002 : if (elframe != refFreq.getRef() && !stupidTransform) {
105 : try{
106 859 : 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 5143 : nu0 = refFreq.getValue().getValue();
113 : }
114 :
115 6002 : return nu0;
116 : }
117 :
118 11898 : void SpectralModel::setRefFrequency(const MFrequency& newRefFreq) {
119 11898 : itsRefFreq = newRefFreq;
120 11898 : DebugAssert(SpectralModel::ok(), AipsError);
121 11898 : }
122 :
123 6937 : const Unit& SpectralModel::frequencyUnit() const {
124 6937 : DebugAssert(SpectralModel::ok(), AipsError);
125 6937 : return itsFreqUnit;
126 : }
127 :
128 37 : void SpectralModel::convertFrequencyUnit(const Unit& freqUnit) {
129 37 : itsFreqUnit = freqUnit;
130 37 : DebugAssert(SpectralModel::ok(), AipsError);
131 37 : }
132 :
133 308 : void SpectralModel::setRefFrequencyError(const Quantum<Double>& newRefFreqErr){
134 308 : 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 308 : itsFreqErr = newRefFreqErr;
140 308 : DebugAssert(SpectralModel::ok(), AipsError);
141 308 : }
142 :
143 385 : const Quantum<Double>& SpectralModel::refFrequencyError() const {
144 385 : 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 6148 : Bool SpectralModel::fromRecord(String& errorMessage,
160 : const RecordInterface& record) {
161 12296 : const String freqString("frequency");
162 6148 : if (!record.isDefined(freqString)) {
163 0 : errorMessage += "The 'frequency' field does not exist\n";
164 0 : return false;
165 : }
166 12296 : const RecordFieldId frequency(freqString);
167 6148 : 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 12296 : const Record& freqRecord = record.asRecord(frequency);
179 12296 : MeasureHolder mh;
180 6148 : if (!mh.fromRecord(errorMessage, freqRecord)) {
181 0 : errorMessage += "Could not parse the reference frequency\n";
182 0 : return false;
183 : }
184 6148 : if (!mh.isMFrequency()) {
185 0 : errorMessage += "The reference frequency is not a frequency measure\n";
186 0 : return false;
187 : }
188 6148 : SpectralModel::setRefFrequency(mh.asMFrequency());
189 6148 : return true;
190 : }
191 :
192 6937 : Bool SpectralModel::toRecord(String& errorMessage,
193 : RecordInterface& record) const {
194 6937 : DebugAssert(SpectralModel::ok(), AipsError);
195 6937 : record.define(RecordFieldId("type"), ComponentType::name(type()));
196 13874 : Record freqRecord;
197 13874 : const MeasureHolder mh(refFrequency());
198 6937 : if (!mh.toRecord(errorMessage, freqRecord)) {
199 0 : errorMessage += "Could not convert the reference frequency to a record\n";
200 0 : return false;
201 : }
202 6937 : const String m0String("m0");
203 6937 : if (freqRecord.isDefined(m0String)) {
204 13874 : const RecordFieldId m0(m0String);
205 6937 : if (freqRecord.dataType(m0) == TpRecord) {
206 13874 : Record m0Rec = freqRecord.asRecord(m0);
207 13874 : QuantumHolder qh;
208 13874 : String error;
209 6937 : if (qh.fromRecord(error, m0Rec) && qh.isQuantumDouble()) {
210 13874 : Quantum<Double> q = qh.asQuantumDouble();
211 6937 : q.convert(frequencyUnit());
212 6937 : qh = QuantumHolder(q);
213 6937 : if (qh.toRecord(error, m0Rec)) {
214 6937 : freqRecord.defineRecord(m0, m0Rec);
215 : }
216 : }
217 : }
218 : }
219 6937 : record.defineRecord(RecordFieldId("frequency"), freqRecord);
220 6937 : return true;
221 : }
222 :
223 5840 : ComponentType::SpectralShape SpectralModel::
224 : getType(String& errorMessage, const RecordInterface& record) {
225 11680 : const String typeString("type");
226 5840 : 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 11680 : const RecordFieldId type(typeString);
232 5840 : 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 5840 : 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 5840 : const String& typeVal = record.asString(type);
243 5840 : return ComponentType::spectralShape(typeVal);
244 : }
245 :
246 2410747 : Bool SpectralModel::ok() const {
247 2410747 : 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 2410747 : return true;
255 : }
256 :
257 308 : Bool SpectralModel::badError(const Quantum<Double>& quantum) {
258 308 : const UnitVal hz(1, "Hz");
259 616 : 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 :
|
