Line data Source code
1 : //# SelectAverageSpw.cc: Implementation of SelectAverageSpw.h
2 : //# Copyright (C) 1996,1997,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$
27 : //----------------------------------------------------------------------------
28 :
29 :
30 : #include <msvis/MSVis/SelectAverageSpw.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/Slice.h>
34 : #include <casacore/measures/Measures/MeasTable.h>
35 :
36 : #include <stdio.h>
37 : #include <stdlib.h>
38 : #include <errno.h>
39 : #include <iostream>
40 : #include <iomanip>
41 : #include <cmath>
42 :
43 : #define LOG2 0
44 :
45 : using namespace casacore;
46 : namespace casa {
47 :
48 : const Int SelectAverageSpw::maxChan = 100000;
49 :
50 0 : Int SelectAverageSpw::nextSelected(Int sp, Int currId, Matrix<Int>& chanList) {
51 : //cout << "spw=" << sp << " currId=" << currId << endl;
52 0 : Int pick = SelectAverageSpw::maxChan;
53 0 : Int nWin = chanList.shape()(0);
54 0 : for (Int i = 0; i < nWin; i++) {
55 0 : if (chanList(i, 0) == sp) {
56 0 : Int tic = -1;
57 0 : for (Int j = chanList(i, 1); j <= chanList(i, 2);
58 0 : j += chanList(i, 3)) {
59 0 : if (j > currId) {
60 0 : tic = j;
61 0 : break;
62 : }
63 : }
64 : //cout << "tic=" << tic << " pick=" << pick << endl;
65 0 : if (tic > -1 && tic < pick)
66 0 : pick = tic;
67 : }
68 : }
69 0 : if (pick == maxChan)
70 0 : return -1;
71 : else
72 0 : return pick;
73 : }
74 :
75 0 : Int SelectAverageSpw::selectAverageChan(MS* pMS, const Matrix<Int>& cList,
76 : Vector<SAS>& spw, const Int& aveChan) {
77 :
78 0 : uInt nAveChan = 0;
79 :
80 0 : Matrix<Int> chanList;
81 0 : chanList.resize(0, 0);
82 0 : chanList = cList;
83 :
84 0 : Int nWin = chanList.nrow();
85 : //cout << "number of selected spw:chan =" << nWin << endl;
86 :
87 0 : MSSpWindowColumns spwColumn(pMS->spectralWindow());
88 0 : MSDataDescColumns descColumn(pMS->dataDescription());
89 :
90 : //when spw='', the chanList is empty, nWin=0,
91 : //so make up chanList here
92 0 : if (nWin < 1) {
93 0 : Vector<Int> spwNumChan = spwColumn.numChan().getColumn();
94 0 : Int spwNum = spwNumChan.shape()(0);
95 0 : chanList.resize(spwNum, 4);
96 0 : for (Int w = 0; w < spwNum; w++) {
97 0 : chanList(w, 0) = w;
98 0 : chanList(w, 1) = 0;
99 0 : chanList(w, 2) = spwNumChan(w) - 1;
100 0 : chanList(w, 3) = 1;
101 : }
102 0 : nWin = spwNum;
103 : }
104 : //cout << "number of actual spw:chan =" << nWin << endl;
105 :
106 : //MSSelection produced chanList behaves as following:
107 : //the stride for spw='0:50' is 0
108 : //the stride for spw='0:50~50' is 1
109 : //so force the stride=chanList(i, 3)=1
110 0 : for (Int i = 0; i < nWin; i++) {
111 0 : if (chanList(i, 2) < chanList(i, 1)) {
112 0 : chanList(i, 2) = chanList(i, 1);
113 0 : chanList(i, 3) = 1;
114 : }
115 0 : if (chanList(i, 3) < 1)
116 0 : chanList(i, 3) = 1;
117 : }
118 : //cout << "validated chanList=" << chanList;
119 :
120 : //from spw:chan list extract non-repeat spw ids
121 0 : Vector<Int> spwids(nWin);
122 0 : spwids(0) = chanList(0, 0);
123 0 : Int sCnt = 1;
124 0 : for (Int w = 0; w < nWin; w++) {
125 0 : Bool listed = false;
126 0 : for (Int i = 0; i < sCnt; i++) {
127 0 : if (spwids(i) == chanList(w, 0)) {
128 0 : listed = true;
129 0 : break;
130 : }
131 : }
132 0 : if (!listed) {
133 0 : spwids(sCnt++) = chanList(w, 0);
134 : }
135 : }
136 0 : spwids.resize(sCnt, true);
137 : //cout << "actual spws=" << spwids << endl;
138 :
139 0 : Vector<Int> spwDesc = descColumn.spectralWindowId().getColumn();
140 :
141 : //fill attributes for each of the selected spws
142 0 : spw.resize(sCnt);
143 0 : for (Int s = 0; s < sCnt; s++) {
144 :
145 0 : spw(s).spwid = spwids(s);
146 :
147 0 : for (uInt t = 0; t < spwDesc.nelements(); t++) {
148 0 : if (spwDesc(t) == spw(s).spwid) {
149 0 : spw(s).desc = t;
150 0 : break;
151 : }
152 : }
153 :
154 0 : spw(s).measFreqRef = spwColumn.measFreqRef()(spw(s).spwid);
155 0 : Double rfreq = spwColumn.refFrequency()(spw(s).spwid);
156 0 : spw(s).rFreq = rfreq / 1000.;
157 :
158 0 : Vector<Int> chids(SelectAverageSpw::maxChan);
159 0 : Int cCnt = 0;
160 0 : Int ch = -1;
161 0 : while(
162 0 : (ch = SelectAverageSpw::nextSelected(spw(s).spwid, ch, chanList)) > -1
163 0 : && cCnt < SelectAverageSpw::maxChan){
164 0 : chids(cCnt++) = ch;
165 : }
166 0 : chids.resize(cCnt, true);
167 : //cout << "chids=" << chids << endl;
168 0 : spw(s).chans.resize(cCnt);
169 0 : spw(s).chans.assign(chids);
170 :
171 : }
172 :
173 : //fill in the averaged attributes
174 : //cout << "aveChan=" << aveChan << endl;
175 0 : Int ssChan = 0;
176 0 : for (Int s = 0; s < sCnt; s++) {
177 :
178 0 : Int cCnt = Int(ceil(Float(spw(s).chans.size()) / aveChan));
179 : //cout << "cCnt=" << cCnt << endl;
180 :
181 0 : spw(s).sxsChans.resize(cCnt);
182 0 : spw(s).aveChans.resize(cCnt);
183 0 : spw(s).aveFreqs.resize(cCnt);
184 0 : spw(s).aveChanNames.resize(cCnt);
185 :
186 0 : Int aCnt = spw(s).chans.size();
187 0 : for (Int i = 0; i < cCnt; i++) {
188 0 : Int sumChan = 0;
189 : Int j;
190 0 : for (j = 0; j < aveChan && (i * aveChan + j) < aCnt; j++) {
191 0 : sumChan += spw(s).chans(i * aveChan + j);
192 : }
193 0 : spw(s).sxsChans(i) = ssChan++;
194 0 : spw(s).aveChans(i) = sumChan / j;
195 : }
196 :
197 : Vector<Double> chanFrequencies =
198 0 : spwColumn.chanFreq()(spw(s).spwid);
199 : //cout << "chanfreq=" << chanFrequencies << endl;
200 :
201 0 : for (Int i = 0; i < cCnt; i++) {
202 0 : Double freq = 0.;
203 : Int j;
204 0 : for (j = 0; j < aveChan && (i * aveChan + j) < aCnt; j++) {
205 0 : freq += chanFrequencies(spw(s).chans(i * aveChan + j));
206 : }
207 0 : spw(s).aveFreqs(i) = freq / j / 1.e9;
208 : }
209 :
210 0 : for (Int i = 0; i < cCnt; i++) {
211 0 : String nm("");
212 : Int j;
213 0 : Int start = spw(s).chans(i * aveChan);
214 0 : Int end = start;
215 0 : Int step = 1;
216 0 : Int k = 1;
217 0 : for (j = 1; j < aveChan && (i * aveChan + j) < aCnt; j++) {
218 0 : if (k++ == 1) {
219 0 : end = spw(s).chans(i * aveChan + j);
220 0 : step = end - start;
221 : }
222 : else {
223 0 : Int del = spw(s).chans(i * aveChan + j);
224 0 : if (del - end == step) {
225 0 : end = del;
226 : }
227 : else {
228 0 : nm += String("[") + String::toString(start) +
229 0 : String(":") + String::toString(end) +
230 0 : String(":") + String::toString(step) + String("]");
231 0 : start = del;
232 0 : end = start;
233 0 : step = 1;
234 0 : k = 1;
235 : }
236 : }
237 : }
238 0 : if (k > 0) {
239 0 : nm += String("[") + String::toString(start) +
240 0 : String(":") + String::toString(end) +
241 0 : String(":") + String::toString(step) + String("]");
242 : }
243 0 : spw(s).aveChanNames(i) = nm;
244 : }
245 :
246 0 : Matrix<Int> maps(aCnt, 5);
247 0 : Int rCnt = 0;
248 0 : for (Int i = 0; i < cCnt; i++) {
249 : //cout << "++++maps=" << maps;
250 : Int j;
251 0 : Int start = spw(s).chans(i * aveChan);
252 0 : Int end = start;
253 0 : Int step = 1;
254 0 : Int k = 1;
255 0 : for (j = 1; j < aveChan && (i * aveChan + j) < aCnt; j++) {
256 0 : if (k++ == 1) {
257 0 : end = spw(s).chans(i * aveChan + j);
258 0 : step = end - start;
259 : }
260 : else {
261 0 : Int del = spw(s).chans(i * aveChan + j);
262 0 : if (del - end == step) {
263 0 : end = del;
264 : }
265 : else {
266 : //cout << "rCnt=" << rCnt << " start=" << start
267 : // << " end=" << end << " step=" << step << endl;
268 0 : maps(rCnt, 0) = start;
269 0 : maps(rCnt, 1) = end;
270 0 : maps(rCnt, 2) = step;
271 0 : maps(rCnt, 3) = i;
272 0 : maps(rCnt, 4) = spw(s).spwid;
273 : //cout << "----maps=" << maps;
274 0 : start = del;
275 0 : end = start;
276 0 : step = 1;
277 0 : k = 1;
278 0 : rCnt++;
279 : }
280 : }
281 : }
282 0 : if (k > 0) {
283 : //cout << "rCnt=" << rCnt << " start=" << start
284 : // << " end=" << end << " step=" << step << endl;
285 0 : maps(rCnt, 0) = start;
286 0 : maps(rCnt, 1) = end;
287 0 : maps(rCnt, 2) = step;
288 0 : maps(rCnt, 3) = i;
289 0 : maps(rCnt, 4) = spw(s).spwid;
290 : //cout << "----maps=" << maps;
291 0 : rCnt++;
292 : }
293 : }
294 0 : spw(s).aveChanMaps.resize(0, 0);
295 0 : spw(s).aveChanMaps = maps(Slice(0, rCnt), Slice(0, 5));
296 :
297 0 : if (spw(s).aveChans.size() > nAveChan) {
298 0 : nAveChan = spw(s).aveChans.size();
299 : }
300 :
301 : }
302 0 : return nAveChan;
303 :
304 : }
305 :
306 0 : void SelectAverageSpw::chanMap(Matrix<Int>& cmap, const Vector<SAS>& spw) {
307 :
308 0 : Int sCnt = spw.size();
309 0 : Int cCnt = 0;
310 0 : for (Int s = 0; s < sCnt; s++) {
311 0 : cCnt += spw(s).aveChanMaps.nrow();
312 : }
313 :
314 0 : cmap.resize(cCnt, 5);
315 0 : Int k = 0;
316 0 : for (Int s = 0; s < sCnt; s++) {
317 0 : for (uInt t = 0; t < spw(s).aveChanMaps.nrow(); t++) {
318 0 : cmap(k, 0) = spw(s).aveChanMaps(t, 0);
319 0 : cmap(k, 1) = spw(s).aveChanMaps(t, 1);
320 0 : cmap(k, 2) = spw(s).aveChanMaps(t, 2);
321 0 : cmap(k, 3) = spw(s).aveChanMaps(t, 3);
322 0 : cmap(k, 4) = spw(s).aveChanMaps(t, 4);
323 0 : k++;
324 : }
325 : }
326 0 : return;
327 : }
328 :
329 0 : void SelectAverageSpw::averageVelocity(Bool &sorryVel,
330 : MS* pMS, Vector<SAS>& spw,
331 : Vector<Double>& velo, const Int& spwidx, const Int& field,
332 : const String& restfreq, const String& frame,
333 : const String& doppler) {
334 :
335 0 : MSDerivedValues msdv;
336 0 : msdv.setMeasurementSet(*pMS);
337 0 : msdv.setVelocityReference(MDoppler::RADIO);
338 : //msdv.setFrequencyReference(MFrequency::LSRK);
339 :
340 : //fill in the averaged velocities
341 : //cout << "field=" << field << endl;
342 :
343 0 : Int sCnt = spw.size();
344 0 : if (sCnt < 1)
345 0 : return;
346 :
347 0 : Double cspeed = (QC::c( )).getValue() / 1000.;
348 :
349 0 : String itsRestFreq = restfreq;
350 0 : String itsFrame = frame;
351 0 : String itsDoppler = doppler;
352 :
353 : //cout << "itsRestFreq=" << itsRestFreq
354 : // << " itsFrame=" << itsFrame
355 : // << " itsDoppler=" << itsDoppler << endl;
356 :
357 0 : MFrequency::Types freqtp = MFrequency::LSRK;
358 0 : MDoppler::Types doptp = MDoppler::RADIO;
359 0 : MFrequency trans;
360 0 : Quantity qt;
361 :
362 0 : if (restfreq != "") {
363 0 : if (!MFrequency::getType(freqtp, itsFrame))
364 0 : freqtp = MFrequency::LSRK;
365 0 : if (!MDoppler::getType(doptp, itsDoppler))
366 0 : doptp = MDoppler::RADIO;
367 :
368 0 : msdv.setVelocityReference(doptp);
369 0 : msdv.setFrequencyReference(freqtp);
370 :
371 0 : if (MeasTable::Line(trans, itsRestFreq)) {
372 0 : qt = trans.get("GHz");
373 0 : msdv.setRestFrequency(qt);
374 : // LogIO os(LogOrigin("SelectAverageSpw","averageVelocity"));
375 : // os << LogIO::NORMAL << "setRestFrequency: "+ String::toString(qt.getValue()) + " " + qt.getUnit()
376 : // << LogIO::POST;
377 : // SLog::slog()->out("setRestFrequency: " + String::toString(qt.getValue()) +
378 : //" " + qt.getUnit(),
379 : //"", "SAS", LogMessage::NORMAL5);
380 : }
381 : else {
382 :
383 0 : Double fr = 1;
384 0 : String unit = "";
385 0 : String dc = downcase(itsRestFreq);
386 0 : if (dc.contains("ghz")) {
387 0 : dc = dc.before("ghz");
388 0 : fr = atof(dc.chars());
389 0 : unit = "GHz";
390 : }
391 0 : else if (dc.contains("mhz")) {
392 0 : dc = dc.before("mhz");
393 0 : fr = atof(dc.chars());
394 0 : unit = "MHz";
395 : }
396 : //else if (dc.contains("khz")) {
397 : // dc = dc.before("khz");
398 : // fr = atof(dc.chars());
399 : // unit = "KHz";
400 : //}
401 0 : else if (dc.contains("hz")) {
402 0 : dc = dc.before("hz");
403 0 : fr = atof(dc.chars());
404 0 : unit = "Hz";
405 : }
406 :
407 : //cout << "itsRestFreq=" << fr << " " << unit << endl;
408 0 : qt = Quantity(fr, unit);
409 0 : msdv.setRestFrequency(qt);
410 : //SLog::slog()->out("setRestFrequency: " + String::toString(qt.getValue()) +
411 : // " " + qt.getUnit(),
412 : // "", "SAS", LogMessage::NORMAL5);
413 : }
414 :
415 0 : msdv.setFieldCenter(field);
416 :
417 0 : MSSpWindowColumns spwColumn(pMS->spectralWindow());
418 : //Int freqRef = spwColumn.measFreqRef()(spw(spwidx).spwid);
419 0 : Int cCnt = spw(spwidx).aveFreqs.size();
420 0 : for (Int k = 0; k < cCnt; k++) {
421 0 : Double tmp = msdv.toVelocity(
422 0 : Quantity(spw(spwidx).aveFreqs(k), "GHz")).get("km/s").getValue();
423 0 : velo(k) = (tmp < 0) ? max(tmp, -cspeed) : min(tmp, cspeed);
424 : }
425 : }
426 : else {
427 :
428 0 : msdv.setFieldCenter(field);
429 0 : Bool hasRestFreq = false;
430 0 : Int cCnt = spw(spwidx).aveFreqs.size();
431 : //cout << "spwidx=" << spwidx << " cCnt=" << cCnt << endl;
432 0 : velo = 0.;
433 :
434 0 : MSSpWindowColumns spwColumn(pMS->spectralWindow());
435 0 : Int freqRef = spwColumn.measFreqRef()(spw(spwidx).spwid);
436 :
437 0 : hasRestFreq = msdv.setRestFrequency(field, spw(spwidx).spwid);
438 0 : if (hasRestFreq) {
439 0 : msdv.setFrequencyReference(MFrequency::castType((uInt)freqRef));
440 : }
441 : else {
442 : //cout << "sorryVel=" << sorryVel << endl;
443 0 : msdv.setFrequencyReference(MFrequency::LSRK );
444 0 : msdv.setRestFrequency(Quantity(1.420, "GHz"));
445 0 : if (sorryVel) {
446 0 : LogIO os(LogOrigin("SelectAverageSpw","averageVelocity"));
447 : os << LogIO::NORMAL << "No rest frequency found in the MS, can't caluclate velocity properly."
448 0 : << LogIO::POST;
449 : //SLog::slog()->out(String("No rest frequency found in the MS,")
450 : //+ String(" use LSRK 1.420 Ghz for velocity calculation"),
451 : //+ String(" Can't calculate velocity properly."),
452 : //"", "SAS", LogMessage::NORMAL, false);
453 0 : sorryVel = false;
454 : }
455 : }
456 :
457 0 : for (Int k = 0; k < cCnt; k++) {
458 : //velo(k) =
459 : // msdv.toVelocity(
460 : // Quantity(spw(spwidx).aveFreqs(k), "GHz")
461 : // ).get("km/s").getValue();
462 : //cout << " k=" << k << " field=" << field
463 : // << " reffreq=" << freqRef
464 : // << " freq=" << spw(spwidx).aveFreqs(k)
465 : // << endl;
466 0 : Double tmp = msdv.toVelocity(
467 0 : Quantity(spw(spwidx).aveFreqs(k), "GHz")).get("km/s").getValue();
468 0 : velo(k) = (tmp < 0) ? max(tmp, -cspeed) : min(tmp, cspeed);
469 : }
470 : }
471 0 : return;
472 : }
473 :
474 0 : void SelectAverageSpw::showSASC(const Vector<SAS>& spw) {
475 0 : Int sCnt = spw.size();
476 0 : for (Int i = 0; i < sCnt; i++) {
477 0 : cout << "spwid=" << spw(i).spwid
478 0 : << " desc=" << spw(i).desc
479 0 : << " rFreq=" << std::setprecision(16) << spw(i).rFreq
480 0 : << " measFreqRef=" << spw(i).measFreqRef
481 0 : << " chans=" << spw(i).chans
482 0 : << " aveChans" << spw(i).aveChans
483 0 : << " sxsChans" << spw(i).sxsChans
484 0 : << " aveFreqs" << spw(i).aveFreqs
485 0 : << " aveChanNames" << spw(i).aveChanNames
486 0 : << " aveChanMaps" << spw(i).aveChanMaps
487 0 : << endl;
488 : }
489 0 : return;
490 : }
491 :
492 0 : Int SelectAverageSpw::descBySpw(const Int& spid, const Vector<SAS>& spw) {
493 0 : Int descid = -1;
494 0 : for (uInt i = 0; i < spw.size(); i++) {
495 0 : if (spw(i).spwid == spid) {
496 0 : descid = spw(i).desc;
497 0 : break;
498 : }
499 : }
500 0 : return descid;
501 : }
502 :
503 0 : Int SelectAverageSpw::spwByDesc(const Int& desc, const Vector<SAS>& spw) {
504 0 : Int spid = -1;
505 0 : for (uInt i = 0; i < spw.size(); i++) {
506 0 : if (spw(i).desc == desc) {
507 0 : spid = spw(i).spwid;
508 0 : break;
509 : }
510 : }
511 0 : return spid;
512 : }
513 :
514 0 : Int SelectAverageSpw::spwIndexByDesc(
515 : const Int& desc, const Vector<SAS>& spw) {
516 0 : Int idx = -1;
517 0 : for (uInt i = 0; i < spw.size(); i++) {
518 0 : if (spw(i).desc == desc) {
519 0 : idx = i;
520 0 : break;
521 : }
522 : }
523 0 : return idx;
524 : }
525 :
526 0 : Int SelectAverageSpw::spwIndexBySpw(
527 : const Int& spid, const Vector<SAS>& spw) {
528 0 : Int idx = -1;
529 0 : for (uInt i = 0; i < spw.size(); i++) {
530 0 : if (spw(i).spwid == spid) {
531 0 : idx = i;
532 0 : break;
533 : }
534 : }
535 0 : return idx;
536 : }
537 :
538 :
539 : }
540 :
|
