source: branches/casa-release-4_3-test/external-alma/components/SpectralComponents/SpectralElementFactory.cc@ 3137

Last change on this file since 3137 was 2980, checked in by Malte Marquarding, 10 years ago

Add a copy of casacore/components/SpectralComponents to external-alma directory to prepare for its removal from casacore-trunk. DOn;t activate in SConscript yet.
File size: 7.2 KB
Line 
1//# Spectral2Element.cc: Record conversion for SpectralElement
2//# Copyright (C) 2001,2004
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: SpectralElementFactory.cc 21465 2014-06-19 05:56:56Z gervandiepen $
27
28//# Includes
29
30#include <casa/Containers/Record.h>
31#include <casa/Utilities/PtrHolder.h>
32#include <components/SpectralComponents/CompiledSpectralElement.h>
33#include <components/SpectralComponents/GaussianSpectralElement.h>
34#include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
35#include <components/SpectralComponents/LogTransformedPolynomialSpectralElement.h>
36#include <components/SpectralComponents/LorentzianSpectralElement.h>
37#include <components/SpectralComponents/PolynomialSpectralElement.h>
38#include <components/SpectralComponents/PowerLogPolynomialSpectralElement.h>
39#include <components/SpectralComponents/SpectralElementFactory.h>
40
41namespace casa { //# NAMESPACE CASA - BEGIN
42
43SpectralElement* SpectralElementFactory::fromRecord(
44 const RecordInterface &in
45) {
46 PtrHolder<SpectralElement> specEl;
47 String origin = "SpectralElementFactory::fromRecord: ";
48 if (
49 ! in.isDefined("type")
50 || in.type(in.idToNumber(RecordFieldId("type"))) != TpString
51 ) {
52 throw AipsError("Record does not represent a SpectralElement");
53 }
54 String stp;
55 SpectralElement::Types tp;
56 in.get(RecordFieldId("type"), stp);
57 if (!SpectralElement::toType(tp, stp)) {
58 throw AipsError("Unknown spectral type in SpectralElement::fromRecord\n");
59 }
60
61 Vector<Double> errs;
62
63 // Get the errors if defined in record
64
65 if (in.isDefined("errors")) {
66 DataType dataType = in.dataType("errors");
67 if (dataType == TpArrayDouble) {
68 in.get("errors", errs);
69 }
70 else if (dataType == TpArrayFloat) {
71 Vector<Float> v;
72 in.get("errors", v);
73 errs.resize(v.nelements());
74 convertArray(errs, v);
75 }
76 else if (dataType == TpArrayInt) {
77 Vector<Int> v;
78 in.get("errors", v);
79 errs.resize(v.nelements());
80 convertArray(errs, v);
81 }
82 else {
83 throw AipsError(
84 "SpectralElement::fromRecord: errors field "
85 "must be double, float or int\n"
86 );
87 }
88 }
89
90 Vector<Double> param;
91 if (in.isDefined(("parameters"))) {
92 DataType dataType = in.dataType("parameters");
93 if (dataType == TpArrayDouble) {
94 in.get("parameters", param);
95 }
96 else if (dataType == TpArrayFloat) {
97 Vector<Float> v;
98 in.get("parameters", v);
99 param.resize(v.nelements());
100 convertArray(param, v);
101 }
102 else if (dataType == TpArrayInt) {
103 Vector<Int> v;
104 in.get("parameters", v);
105 param.resize(v.nelements());
106 convertArray(param, v);
107 }
108 else {
109 throw AipsError(
110 origin +
111 "SpectralElement::fromRecord: parameters field "
112 "must be double, float or int\n"
113 );
114 }
115 }
116
117 // Make sizes of errors and params equal
118 if (errs.nelements() == 0) {
119 errs.resize(param.nelements());
120 errs = 0.0;
121 }
122 if (errs.nelements() != param.nelements()) {
123 throw AipsError(
124 origin +
125 "SpectralElement::fromRecord must have equal lengths "
126 "for parameters and errors fields"
127 );
128 }
129 switch (tp) {
130 case SpectralElement::GAUSSIAN:
131 if (param.nelements() != 3) {
132 throw AipsError(
133 origin + "Illegal number of parameters for Gaussian element"
134 );
135 }
136 if (param(2) <= 0.0) {
137 throw AipsError(
138 origin +
139 "The width of a Gaussian element must be positive"
140 );
141 }
142 param(2) = GaussianSpectralElement::sigmaFromFWHM (param(2));
143 errs(2) = GaussianSpectralElement::sigmaFromFWHM (errs(2));
144 specEl.set(new GaussianSpectralElement(param(0), param(1), param(2)));
145 specEl->setError(errs);
146 break;
147 case SpectralElement::LORENTZIAN:
148 if (param.nelements() != 3) {
149 throw AipsError(
150 "Illegal number of parameters for Lorentzian element"
151 );
152 }
153 if (param(2) <= 0.0) {
154 throw AipsError(
155 "The width of a Lorentzian element must be positive"
156 );
157 }
158 specEl.set(new LorentzianSpectralElement(param(0), param(1), param(2)));
159 specEl->setError(errs);
160 break;
161 case SpectralElement::POLYNOMIAL:
162 if (param.nelements() == 0) {
163 throw AipsError(
164 "Polynomial spectral element must have order "
165 "of at least zero"
166 );
167 }
168 specEl.set(new PolynomialSpectralElement(param.nelements() - 1));
169 specEl->set(param);
170 specEl->setError(errs);
171 break;
172 case SpectralElement::COMPILED:
173 if (
174 in.isDefined("compiled")
175 && in.type(in.idToNumber(RecordFieldId("compiled"))) == TpString
176 ) {
177 String function;
178 in.get(RecordFieldId("compiled"), function);
179 specEl.set(new CompiledSpectralElement(function, param));
180 specEl->setError(errs);
181 }
182 else {
183 throw AipsError(
184 "No compiled string in SpectralElement::fromRecord\n"
185 );
186 }
187 break;
188
189 case SpectralElement::GMULTIPLET:
190 {
191 if (! in.isDefined("gaussians")) {
192 throw AipsError("gaussians not defined in record");
193 }
194 if (! in.isDefined("fixedMatrix")) {
195 throw AipsError("fixed matrix not defined in record");
196 }
197 Record gaussians = in.asRecord("gaussians");
198 uInt i = 0;
199 vector<GaussianSpectralElement> comps(0);
200 while(True) {
201 String id = "*" + String::toString(i);
202 if (gaussians.isDefined(id)) {
203 PtrHolder<SpectralElement> gauss(fromRecord(gaussians.asRecord(id)));
204 comps.push_back(
205 *dynamic_cast<GaussianSpectralElement*>(
206 gauss.ptr()
207 )
208 );
209 i++;
210 }
211 else {
212 break;
213 }
214 }
215 Matrix<Double> fixedMatrix = in.asArrayDouble("fixedMatrix");
216 fixedMatrix.reform(IPosition(2, comps.size()-1, 3));
217 specEl.set(new GaussianMultipletSpectralElement(comps, fixedMatrix));
218 }
219 break;
220
221 case SpectralElement::POWERLOGPOLY: {
222 specEl.set(new PowerLogPolynomialSpectralElement(param));
223 specEl->set(param);
224 specEl->setError(errs);
225 }
226 break;
227
228 case SpectralElement::LOGTRANSPOLY: {
229 specEl.set(new LogTransformedPolynomialSpectralElement(param));
230 specEl->set(param);
231 specEl->setError(errs);
232 }
233 break;
234
235 default:
236 throw AipsError(
237 "Unhandled or illegal spectral element record in "
238 "SpectralElementFactory::fromRecord\n"
239 );
240 }
241
242 if (in.isDefined("fixed")) {
243 specEl->fix(in.asArrayBool("fixed"));
244 }
245 // ready to return, fish out the pointer and return it without deleting it
246 SpectralElement *sp = specEl.ptr();
247 specEl.clear(False);
248 return sp;
249}
250
251
252} //# NAMESPACE CASA - END
253
Note: See TracBrowser for help on using the repository browser.