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: 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: 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 159279 : SpectralModel::SpectralModel()
51 159279 : :itsRefFreq(Quantum<Double>(1, "GHz"), MFrequency::DEFAULT),
52 159279 : itsFreqUnit("GHz"),
53 318558 : itsFreqErr(0, itsFreqUnit)
54 : {
55 159279 : DebugAssert(SpectralModel::ok(), AipsError);
56 159279 : }
57 :
58 92 : SpectralModel::SpectralModel(const MFrequency& refFreq, const Unit& freqUnit)
59 92 : :itsRefFreq(refFreq),
60 92 : itsFreqUnit(freqUnit),
61 184 : itsFreqErr(0, itsFreqUnit)
62 : {
63 92 : DebugAssert(SpectralModel::ok(), AipsError);
64 92 : }
65 :
66 56873 : SpectralModel::SpectralModel(const SpectralModel& other)
67 : :RecordTransformable(),
68 56873 : itsRefFreq(other.itsRefFreq),
69 56873 : itsFreqUnit(other.itsFreqUnit),
70 113746 : itsFreqErr(other.itsFreqErr)
71 : {
72 56873 : DebugAssert(SpectralModel::ok(), AipsError);
73 56873 : }
74 :
75 37710 : SpectralModel& SpectralModel::operator=(const SpectralModel& other) {
76 37710 : if (this != &other) {
77 37710 : itsRefFreq = other.itsRefFreq;
78 37710 : itsFreqUnit = other.itsFreqUnit;
79 37710 : itsFreqErr = other.itsFreqErr;
80 : }
81 37710 : DebugAssert(SpectralModel::ok(), AipsError);
82 37710 : return *this;
83 : }
84 :
85 216243 : SpectralModel::~SpectralModel() {
86 216243 : }
87 :
88 2695 : const String& SpectralModel::ident() const {
89 2695 : DebugAssert(SpectralModel::ok(), AipsError);
90 2695 : static String typeString;
91 2695 : typeString = ComponentType::name(type());
92 2695 : return typeString;
93 : }
94 :
95 504887 : const MFrequency& SpectralModel::refFrequency() const {
96 504887 : DebugAssert(SpectralModel::ok(), AipsError);
97 504887 : return itsRefFreq;
98 : }
99 :
100 44588 : Double SpectralModel::refFreqInFrame(const MFrequency::Ref& elframe) const {
101 44588 : const MFrequency& refFreq(refFrequency());
102 : Double nu0;
103 44588 : Bool stupidTransform = (elframe.getType() == MFrequency::REST) || (elframe.getType() == MFrequency::N_Types) || (refFreq.getRef().getType() == MFrequency::REST) || (refFreq.getRef().getType() == MFrequency::N_Types);
104 44588 : 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 0 : }
111 : } else {
112 43729 : nu0 = refFreq.getValue().getValue();
113 : }
114 :
115 44588 : return nu0;
116 : }
117 :
118 115623 : void SpectralModel::setRefFrequency(const MFrequency& newRefFreq) {
119 115623 : itsRefFreq = newRefFreq;
120 115623 : DebugAssert(SpectralModel::ok(), AipsError);
121 115623 : }
122 :
123 10047 : const Unit& SpectralModel::frequencyUnit() const {
124 10047 : DebugAssert(SpectralModel::ok(), AipsError);
125 10047 : 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 1463 : void SpectralModel::setRefFrequencyError(const Quantum<Double>& newRefFreqErr){
134 1463 : 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 0 : }
139 1463 : itsFreqErr = newRefFreqErr;
140 1463 : DebugAssert(SpectralModel::ok(), AipsError);
141 1463 : }
142 :
143 2695 : const Quantum<Double>& SpectralModel::refFrequencyError() const {
144 2695 : 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 79438 : Bool SpectralModel::fromRecord(String& errorMessage,
160 : const RecordInterface& record) {
161 79438 : const String freqString("frequency");
162 79438 : if (!record.isDefined(freqString)) {
163 0 : errorMessage += "The 'frequency' field does not exist\n";
164 0 : return false;
165 : }
166 79438 : const RecordFieldId frequency(freqString);
167 79438 : 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 79438 : const Record& freqRecord = record.asRecord(frequency);
179 79438 : MeasureHolder mh;
180 79438 : if (!mh.fromRecord(errorMessage, freqRecord)) {
181 0 : errorMessage += "Could not parse the reference frequency\n";
182 0 : return false;
183 : }
184 79438 : if (!mh.isMFrequency()) {
185 0 : errorMessage += "The reference frequency is not a frequency measure\n";
186 0 : return false;
187 : }
188 79438 : SpectralModel::setRefFrequency(mh.asMFrequency());
189 79438 : return true;
190 79438 : }
191 :
192 10047 : Bool SpectralModel::toRecord(String& errorMessage,
193 : RecordInterface& record) const {
194 10047 : DebugAssert(SpectralModel::ok(), AipsError);
195 10047 : record.define(RecordFieldId("type"), ComponentType::name(type()));
196 10047 : Record freqRecord;
197 10047 : const MeasureHolder mh(refFrequency());
198 10047 : if (!mh.toRecord(errorMessage, freqRecord)) {
199 0 : errorMessage += "Could not convert the reference frequency to a record\n";
200 0 : return false;
201 : }
202 10047 : const String m0String("m0");
203 10047 : if (freqRecord.isDefined(m0String)) {
204 10047 : const RecordFieldId m0(m0String);
205 10047 : if (freqRecord.dataType(m0) == TpRecord) {
206 10047 : Record m0Rec = freqRecord.asRecord(m0);
207 10047 : QuantumHolder qh;
208 10047 : String error;
209 10047 : if (qh.fromRecord(error, m0Rec) && qh.isQuantumDouble()) {
210 10047 : Quantum<Double> q = qh.asQuantumDouble();
211 10047 : q.convert(frequencyUnit());
212 10047 : qh = QuantumHolder(q);
213 10047 : if (qh.toRecord(error, m0Rec)) {
214 10047 : freqRecord.defineRecord(m0, m0Rec);
215 : }
216 10047 : }
217 10047 : }
218 10047 : }
219 10047 : record.defineRecord(RecordFieldId("frequency"), freqRecord);
220 10047 : return true;
221 10047 : }
222 :
223 77975 : ComponentType::SpectralShape SpectralModel::
224 : getType(String& errorMessage, const RecordInterface& record) {
225 77975 : const String typeString("type");
226 77975 : 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 77975 : const RecordFieldId type(typeString);
232 77975 : 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 77975 : 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 77975 : const String& typeVal = record.asString(type);
243 77975 : return ComponentType::spectralShape(typeVal);
244 77975 : }
245 :
246 3861136 : Bool SpectralModel::ok() const {
247 3861136 : 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 0 : }
254 3861136 : return true;
255 : }
256 :
257 1463 : Bool SpectralModel::badError(const Quantum<Double>& quantum) {
258 1463 : const UnitVal hz(1, "Hz");
259 2926 : return !(quantum.check(hz)) || (quantum.getValue() < 0.0);
260 1463 : }
261 :
262 : // Local Variables:
263 : // compile-command: "gmake SpectralModel"
264 : // End:
265 :
266 : } //# NAMESPACE CASA - END
267 :
|
