Line data Source code
1 : //# MSContinuumSubtractor.cc: Subtract continuum from spectral line data
2 : //# Copyright (C) 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: 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$
27 : //#
28 :
29 : #include <spectrallines/Splatalogue/ListConverter.h>
30 :
31 : #include <casacore/casa/Quanta/Quantum.h>
32 : #include <casacore/casa/Quanta/Unit.h>
33 : #include <casacore/casa/OS/RegularFile.h>
34 : #include <casacore/casa/IO/RegularFileIO.h>
35 : #include <casacore/casa/Utilities/Regex.h>
36 : #include <casacore/tables/Tables/ScaColDesc.h>
37 : #include <casacore/tables/Tables/ScalarColumn.h>
38 : #include <casacore/tables/Tables/SetupNewTab.h>
39 : #include <casacore/tables/DataMan/StandardStMan.h>
40 : #include <casacore/tables/Tables/TableDesc.h>
41 :
42 : using namespace casacore;
43 : namespace casa {
44 :
45 0 : ListConverter::ListConverter(const Vector<String>& filenames, const String& table) :
46 0 : _log(new LogIO()), _listFiles(Vector<RegularFile>(0)), _tableName(table), _freqUnit("") {
47 0 : RegularFile t(_tableName);
48 0 : if (t.exists()) {
49 0 : *_log << "File " << _tableName << " already exists and will not be overwritten" << LogIO::EXCEPTION;
50 : }
51 0 : File::FileWriteStatus tableStatus = t.getWriteStatus();
52 0 : if(tableStatus == File::NOT_CREATABLE ) {
53 0 : *_log << "Table " << _tableName << " cannot be created" << LogIO::EXCEPTION;
54 : }
55 0 : for (Vector<String>::const_iterator iter = filenames.begin(); iter != filenames.end(); iter++) {
56 0 : RegularFile file(*iter);
57 0 : if (! file.exists()) {
58 0 : *_log << "File " << *iter << " does not exist" << LogIO::EXCEPTION;
59 : }
60 0 : if (! file.isReadable()) {
61 0 : *_log << "File " << *iter << " does not have read permission" << LogIO::EXCEPTION;
62 : }
63 0 : _listFiles.resize(_listFiles.size()+1, true);
64 0 : _listFiles[_listFiles.size() - 1] = file;
65 0 : }
66 0 : }
67 :
68 0 : ListConverter::~ListConverter() {
69 0 : delete _log;
70 0 : }
71 :
72 0 : SplatalogueTable* ListConverter::load() {
73 0 : _parseLists();
74 0 : SplatalogueTable *newTable = _defineTable(_species.size());
75 0 : _addData(newTable);
76 0 : newTable->flush();
77 0 : return newTable;
78 : }
79 :
80 0 : void ListConverter::_parseLists() {
81 0 : _log->origin(LogOrigin("ListConverter", __FUNCTION__));
82 :
83 0 : uInt transitionIndex = 0;
84 0 : Regex blankLine("^[ \n\t\r\v\f]+$");
85 0 : Regex lastLine("^Found [0-9]+ rows");
86 0 : Regex debye2("D<sup>2</sup>");
87 0 : Regex k("(K)");
88 0 : Regex cm1("(cm<sup>-1</sup>)");
89 :
90 0 : _frequency.resize(0);
91 0 : for (
92 0 : Vector<RegularFile>::const_iterator iter = _listFiles.begin();
93 0 : iter != _listFiles.end(); iter++
94 : ) {
95 0 : RegularFileIO fileIO(*iter);
96 0 : Int bufSize = 10000;
97 0 : char *buffer = new char[bufSize];
98 : int nRead;
99 0 : String contents;
100 0 : while ((nRead = fileIO.read(bufSize, buffer, false)) == bufSize) {
101 0 : String chunk(buffer, bufSize);
102 0 : contents += chunk;
103 0 : }
104 : // get the last chunk
105 0 : String chunk(buffer, nRead);
106 0 : contents += chunk;
107 :
108 0 : Vector<String> lines = stringToVector(contents, '\n');
109 0 : uInt newSize = _frequency.size() + lines.size();
110 0 : _species.resize(newSize, true);
111 0 : _chemName.resize(newSize, true);
112 0 : _qns.resize(newSize, true);
113 0 : _lineList.resize(newSize, true);
114 0 : _recommended.resize(newSize, true);
115 0 : _frequency.resize(newSize, true);
116 0 : _intensity.resize(newSize, true);
117 0 : _smu2.resize(newSize, true);
118 0 : _logA.resize(newSize, true);
119 0 : _eL.resize(newSize, true);
120 0 : _eU.resize(newSize, true);
121 0 : uInt lineCount = 0;
122 0 : for(Vector<String>::iterator liter=lines.begin(); liter!=lines.end(); liter++) {
123 0 : if (
124 0 : liter->empty() || liter->firstchar() == '#'
125 0 : || liter->matches(blankLine) || liter->matches(lastLine)
126 : ) {
127 : // skip comment
128 0 : lineCount++;
129 0 : continue;
130 : }
131 0 : uInt tabCount = liter->freq('\t');
132 0 : if (tabCount != 10) {
133 0 : *_log << "bad format for line " << (lineCount+1)
134 0 : << " in file " << iter->path().baseName() << LogIO::EXCEPTION;
135 : }
136 0 : Vector<String> parts = stringToVector(*liter, '\t');
137 0 : for (Vector<String>::iterator viter = parts.begin(); viter != parts.end(); viter++) {
138 0 : viter->trim();
139 0 : }
140 0 : String filename = iter->path().dirName() + "/" + iter->path().baseName();
141 0 : String species = parts[0];
142 0 : String recommended = parts[1];
143 0 : String chemicalName = parts[2];
144 0 : String frequency = parts[3];
145 0 : String qns = parts[4];
146 0 : String intensity = parts[5];
147 0 : String smu2 = parts[6];
148 0 : String logA = parts[7];
149 0 : String eL = parts[8];
150 0 : String eU = parts[9];
151 0 : String lineList = parts[10];
152 0 : if (lineCount == 0) {
153 : // header line look for units
154 0 : Vector<String> freqHeaderParts = stringToVector(frequency, ' ');
155 0 : for (
156 0 : Vector<String>::const_iterator viter=freqHeaderParts.begin();
157 0 : viter!=freqHeaderParts.end(); viter++
158 : ) {
159 : try {
160 0 : Unit freqUnit(*viter);
161 0 : Quantity qq(1, freqUnit);
162 0 : if (qq.isConform(Unit("Hz"))) {
163 0 : _freqUnit = freqUnit.getName();
164 : }
165 0 : }
166 0 : catch (AipsError x) {}
167 0 : }
168 0 : Vector<String> smuHeaderParts = stringToVector(smu2, ' ');
169 0 : for (
170 0 : Vector<String>::const_iterator viter=smuHeaderParts.begin();
171 0 : viter!=smuHeaderParts.end(); viter++
172 : ) {
173 : try {
174 0 : String smu2(*viter);
175 0 : if (smu2.matches(debye2)) {
176 0 : _smu2Unit = "Debye*Debye";
177 :
178 : }
179 0 : }
180 0 : catch (AipsError x) {}
181 0 : }
182 0 : Vector<String> elHeaderParts = stringToVector(eL, ' ');
183 0 : for (
184 0 : Vector<String>::const_iterator viter=elHeaderParts.begin();
185 0 : viter!=elHeaderParts.end(); viter++
186 : ) {
187 : try {
188 0 : String el = *viter;
189 0 : if (el.matches(k)) {
190 0 : _elUnit = "K";
191 : }
192 0 : else if (el.matches(cm1)) {
193 0 : _elUnit = "cm-1";
194 : }
195 0 : }
196 0 : catch (AipsError x) {}
197 0 : }
198 0 : Vector<String> euHeaderParts = stringToVector(eU, ' ');
199 0 : for (
200 0 : Vector<String>::const_iterator viter=euHeaderParts.begin();
201 0 : viter!=euHeaderParts.end(); viter++
202 : ) {
203 : try {
204 0 : String eu = *viter;
205 0 : if (eu.matches(k)) {
206 0 : _euUnit = "K";
207 : }
208 0 : else if (eu.matches(cm1)) {
209 0 : _euUnit = "cm-1";
210 : }
211 0 : }
212 0 : catch (AipsError x) {}
213 0 : }
214 0 : lineCount++;
215 0 : continue;
216 0 : }
217 :
218 0 : if (lineCount % 20000 == 0) {
219 0 : *_log << LogIO::NORMAL << "Parsing line " << lineCount << " of file "
220 0 : << iter->path().baseName() << LogIO::POST;
221 : }
222 0 : Bool recFreq = recommended.length() > 0;
223 0 : if (! frequency.matches(RXdouble) ) {
224 0 : *_log << "File " << filename << ", line " << lineCount
225 : << ": frequency value " << frequency << " is not numeric"
226 0 : << LogIO::EXCEPTION;
227 : }
228 0 : if (! intensity.matches(RXdouble) ) {
229 0 : *_log << "File " << filename << ", line " << lineCount
230 : << ": intensity value " << intensity << " is not numeric"
231 0 : << LogIO::EXCEPTION;
232 : }
233 0 : if (! smu2.matches(RXdouble) ) {
234 0 : *_log << "File " << filename << ", line " << lineCount
235 : << ": S_ij*mu**2 value " << smu2 << " is not numeric"
236 0 : << LogIO::EXCEPTION;
237 : }
238 0 : if (! logA.matches(RXdouble) ) {
239 0 : *_log << "File " << filename << ", line " << lineCount
240 : << ": log(A) value " << logA << " is not numeric"
241 0 : << LogIO::EXCEPTION;
242 : }
243 0 : if (! eL.matches(RXdouble) ) {
244 0 : *_log << "File " << filename << ", line " << lineCount
245 : << ": E_l value " << eL << " is not numeric"
246 0 : << LogIO::EXCEPTION;
247 : }
248 0 : if (! eU.matches(RXdouble) ) {
249 0 : *_log << "File " << filename << ", line " << lineCount
250 : << ": E_u value " << eU << " is not numeric"
251 0 : << LogIO::EXCEPTION;
252 : }
253 0 : _species[transitionIndex] = species;
254 0 : _chemName[transitionIndex] = chemicalName;
255 0 : _qns[transitionIndex] = qns;
256 0 : _lineList[transitionIndex] = lineList;
257 0 : _recommended[transitionIndex] = recFreq;
258 0 : _frequency[transitionIndex] = String::toDouble(frequency);
259 0 : _intensity[transitionIndex] = String::toFloat(intensity);
260 0 : _smu2[transitionIndex] = String::toFloat(smu2);
261 0 : _logA[transitionIndex] = String::toFloat(logA);
262 0 : _eL[transitionIndex] = String::toFloat(eL);
263 0 : _eU[transitionIndex] = String::toFloat(eU);
264 0 : transitionIndex++;
265 0 : lineCount++;
266 0 : }
267 0 : }
268 0 : _species.resize(transitionIndex, true);
269 0 : _chemName.resize(transitionIndex, true);
270 0 : _qns.resize(transitionIndex, true);
271 0 : _lineList.resize(transitionIndex, true);
272 0 : _recommended.resize(transitionIndex, true);
273 0 : _frequency.resize(transitionIndex, true);
274 0 : _intensity.resize(transitionIndex, true);
275 0 : _smu2.resize(transitionIndex, true);
276 0 : _logA.resize(transitionIndex, true);
277 0 : _eL.resize(transitionIndex, true);
278 0 : _eU.resize(transitionIndex, true);
279 0 : }
280 :
281 :
282 0 : SplatalogueTable* ListConverter::_defineTable(const uInt nrows) {
283 0 : String descName = _tableName + "_desc";
284 0 : TableDesc td(descName,TableDesc::Scratch);
285 0 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::SPECIES, SplatalogueTable::SPECIES));
286 0 : td.addColumn (ScalarColumnDesc<Bool>(SplatalogueTable::RECOMMENDED, SplatalogueTable::RECOMMENDED));
287 0 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::CHEMICAL_NAME, SplatalogueTable::CHEMICAL_NAME));
288 0 : td.addColumn (ScalarColumnDesc<Double>(SplatalogueTable::FREQUENCY, SplatalogueTable::FREQUENCY));
289 0 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::QUANTUM_NUMBERS, SplatalogueTable::QUANTUM_NUMBERS));
290 0 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::INTENSITY, SplatalogueTable::INTENSITY));
291 0 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::SMU2, SplatalogueTable::SMU2));
292 0 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::LOGA, SplatalogueTable::LOGA));
293 0 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::EL, SplatalogueTable::EL));
294 0 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::EU, SplatalogueTable::EU));
295 0 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::LINELIST, SplatalogueTable::LINELIST));
296 0 : SetupNewTable tableSetup(_tableName, td, Table::New);
297 0 : StandardStMan stm;
298 0 : tableSetup.bindColumn(SplatalogueTable::SPECIES, stm);
299 0 : tableSetup.bindColumn(SplatalogueTable::RECOMMENDED, stm);
300 0 : tableSetup.bindColumn(SplatalogueTable::CHEMICAL_NAME, stm);
301 0 : tableSetup.bindColumn(SplatalogueTable::FREQUENCY, stm);
302 0 : tableSetup.bindColumn(SplatalogueTable::QUANTUM_NUMBERS, stm);
303 0 : tableSetup.bindColumn(SplatalogueTable::INTENSITY, stm);
304 0 : tableSetup.bindColumn(SplatalogueTable::SMU2, stm);
305 0 : tableSetup.bindColumn(SplatalogueTable::LOGA, stm);
306 0 : tableSetup.bindColumn(SplatalogueTable::EL, stm);
307 0 : tableSetup.bindColumn(SplatalogueTable::EU, stm);
308 0 : tableSetup.bindColumn(SplatalogueTable::LINELIST, stm);
309 0 : if (_freqUnit.empty()) {
310 0 : *_log << LogIO::WARN << "Could not determine frequency unit, assuming GHz" << LogIO::POST;
311 0 : _freqUnit = "GHz";
312 : }
313 0 : return new SplatalogueTable(tableSetup, nrows, _freqUnit, _smu2Unit, _elUnit, _euUnit);
314 0 : }
315 :
316 0 : void ListConverter::_addData(const SplatalogueTable* table) const {
317 0 : ScalarColumn<String> col(*table, SplatalogueTable::SPECIES);
318 0 : col.putColumn(_species);
319 0 : ScalarColumn<Bool> rcol(*table, SplatalogueTable::RECOMMENDED);
320 0 : rcol.putColumn(_recommended);
321 0 : ScalarColumn<String> cncol(*table, SplatalogueTable::CHEMICAL_NAME);
322 0 : cncol.putColumn(_chemName);
323 0 : ScalarColumn<Double> fcol(*table,SplatalogueTable::FREQUENCY);
324 0 : fcol.putColumn(_frequency);
325 0 : ScalarColumn<String> qncol(*table, SplatalogueTable::QUANTUM_NUMBERS);
326 0 : qncol.putColumn(_qns);
327 0 : ScalarColumn<Float> icol(*table, SplatalogueTable::INTENSITY);
328 0 : icol.putColumn(_intensity);
329 0 : ScalarColumn<Float> smu2col(*table, SplatalogueTable::SMU2);
330 0 : smu2col.putColumn(_smu2);
331 0 : ScalarColumn<Float> lacol(*table, SplatalogueTable::LOGA);
332 0 : lacol.putColumn(_logA);
333 0 : ScalarColumn<Float> elcol(*table, SplatalogueTable::EL);
334 0 : elcol.putColumn(_eL);
335 0 : ScalarColumn<Float> eucol(*table, SplatalogueTable::EU);
336 0 : eucol.putColumn(_eU);
337 0 : ScalarColumn<String> llcol(*table, SplatalogueTable::LINELIST);
338 0 : llcol.putColumn(_lineList);
339 0 : }
340 :
341 : } //# NAMESPACE CASA - END
342 :
|
