Line data Source code
1 : //# MSTransformRegridder.cc: This file contains the implementation of the MSTransformRegridder class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <mstransform/MSTransform/MSTransformRegridder.h>
24 :
25 : using namespace casacore;
26 : namespace casa { //# NAMESPACE CASA - BEGIN
27 :
28 : /////////////////////////////////////////////
29 : /// MSTransformRegridder implementation ///
30 : /////////////////////////////////////////////
31 :
32 : // -----------------------------------------------------------------------
33 : //
34 : // -----------------------------------------------------------------------
35 0 : MSTransformRegridder::MSTransformRegridder()
36 : {
37 0 : return;
38 : }
39 :
40 : // -----------------------------------------------------------------------
41 : //
42 : // -----------------------------------------------------------------------
43 0 : MSTransformRegridder::~MSTransformRegridder()
44 : {
45 0 : return;
46 0 : }
47 :
48 :
49 : // -----------------------------------------------------------------------
50 : // Make one spectral window from all SPWs given by the SPW Ids vector
51 : // -----------------------------------------------------------------------
52 0 : Bool MSTransformRegridder::combineSpws( LogIO& os,
53 : String msName,
54 : const Vector<Int>& spwids,
55 : Vector<Double>& newCHAN_FREQ,
56 : Vector<Double>& newCHAN_WIDTH,
57 : vector<vector<Int> >& averageWhichChan,
58 : vector<vector<Int> >& averageWhichSPW,
59 : vector<vector<Double> >& averageChanFrac,
60 : Bool verbose)
61 : {
62 0 : MeasurementSet ms_p(msName, Table::Old);
63 0 : Bool result = false;
64 0 : result = combineSpwsCore(os,ms_p,spwids,newCHAN_FREQ,newCHAN_WIDTH,
65 : averageWhichChan, averageWhichSPW, averageChanFrac, verbose);
66 0 : return result;
67 : }
68 :
69 : // -----------------------------------------------------------------------
70 : // Make one spectral window from all SPWs given by the SPW Ids vector
71 : // -----------------------------------------------------------------------
72 0 : Bool MSTransformRegridder::combineSpwsCore( LogIO& os,
73 : MeasurementSet& ms_p,
74 : const Vector<Int>& spwids,
75 : Vector<Double>& newCHAN_FREQ,
76 : Vector<Double>& newCHAN_WIDTH,
77 : vector<vector<Int> >& averageWhichChan,
78 : vector<vector<Int> >& averageWhichSPW,
79 : vector<vector<Double> >& averageChanFrac,
80 : Bool verbose)
81 : {
82 : // Analyze SPW Ids
83 0 : if (spwids.nelements() == 0)
84 : {
85 0 : os << LogIO::WARN << LogOrigin("MSTransformRegridder", __FUNCTION__)
86 0 : << "No SPWs selected for combination ..." << LogIO::POST;
87 0 : return true;
88 : }
89 :
90 : // Find all existing spws,
91 0 : MSSpectralWindow spwtable = ms_p.spectralWindow();
92 0 : Int origNumSPWs = spwtable.nrow();
93 :
94 0 : vector<Int> spwsToCombine;
95 0 : Vector<Bool> includeIt(origNumSPWs, false);
96 :
97 : // jagonzal: This covers for the case when we want to combine all the input SPWs
98 0 : if (spwids(0) == -1)
99 : {
100 0 : for (Int i = 0; i < origNumSPWs; i++)
101 : {
102 0 : spwsToCombine.push_back(i);
103 0 : includeIt(i) = true;
104 : }
105 : }
106 : // jagonzal: Nominal case when we want to combine a sub-set of the input SPWs
107 : else
108 : {
109 0 : for (uInt i = 0; i < spwids.nelements(); i++)
110 : {
111 0 : if (spwids(i) < origNumSPWs && spwids(i) >= 0)
112 : {
113 0 : spwsToCombine.push_back(spwids(i));
114 0 : includeIt(spwids(i)) = true;
115 : }
116 : else
117 : {
118 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
119 : << "Invalid SPW ID selected for combination "
120 0 : << spwids(i) << "valid range is 0 - "
121 0 : << origNumSPWs - 1 << ")" << LogIO::POST;
122 0 : return false;
123 : }
124 : }
125 : }
126 :
127 : // jagonzal: Marginal case when there is no actual SPW combination
128 0 : if (spwsToCombine.size() <= 1)
129 : {
130 0 : if (verbose)
131 : {
132 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
133 : << "Less than two SPWs selected. No combination necessary."
134 0 : << LogIO::POST;
135 : }
136 0 : return true;
137 : }
138 :
139 : // Sort the SPW Ids
140 0 : std::sort(spwsToCombine.begin(), spwsToCombine.end());
141 :
142 0 : uInt nSpwsToCombine = spwsToCombine.size();
143 :
144 : // Prepare access to the SPW table
145 0 : MSSpWindowColumns SPWColrs(spwtable);
146 0 : ScalarColumn<Int> numChanColr = SPWColrs.numChan();
147 0 : ArrayColumn<Double> chanFreqColr = SPWColrs.chanFreq();
148 0 : ArrayColumn<Double> chanWidthColr = SPWColrs.chanWidth();
149 0 : ScalarColumn<Int> measFreqRefColr = SPWColrs.measFreqRef();
150 0 : ArrayColumn<Double> effectiveBWColr = SPWColrs.effectiveBW();
151 0 : ScalarColumn<Double> refFrequencyColr = SPWColrs.refFrequency();
152 0 : ArrayColumn<Double> resolutionColr = SPWColrs.resolution();
153 0 : ScalarColumn<Double> totalBandwidthColr = SPWColrs.totalBandwidth();
154 :
155 : // Create a list of the SPW Ids sorted by first (lowest) channel frequency
156 0 : vector<Int> spwsSorted(nSpwsToCombine);
157 0 : Vector<Bool> isDescending(origNumSPWs, false);
158 0 : Bool negChanWidthWarned = false;
159 :
160 :
161 0 : Double* firstFreq = new Double[nSpwsToCombine];
162 0 : uInt count = 0;
163 0 : for (uInt i = 0; (Int) i < origNumSPWs; i++)
164 : {
165 0 : if (includeIt(i))
166 : {
167 0 : Vector<Double> CHAN_FREQ(chanFreqColr(i));
168 :
169 : // If frequencies are ascending, take the first channel, otherwise the last
170 0 : uInt nCh = CHAN_FREQ.nelements();
171 0 : if (CHAN_FREQ(0) <= CHAN_FREQ(nCh - 1))
172 : {
173 0 : firstFreq[count] = CHAN_FREQ(0);
174 : }
175 : else
176 : {
177 0 : firstFreq[count] = CHAN_FREQ(nCh - 1);
178 0 : isDescending(i) = true;
179 : }
180 0 : count++;
181 : }
182 : }
183 :
184 0 : Sort sort;
185 0 : sort.sortKey(firstFreq, TpDouble); // define sort key
186 0 : Vector<uInt> inx(nSpwsToCombine);
187 0 : sort.sort(inx, nSpwsToCombine);
188 0 : for (uInt i = 0; i < nSpwsToCombine; i++)
189 : {
190 0 : spwsSorted[i] = spwsToCombine[inx(i)];
191 : }
192 0 : delete[] firstFreq;
193 :
194 :
195 0 : Int id0 = spwsSorted[0];
196 :
197 0 : uInt newNUM_CHAN = numChanColr(id0);
198 0 : newCHAN_FREQ.assign(chanFreqColr(id0));
199 0 : newCHAN_WIDTH.assign(chanWidthColr(id0));
200 0 : Bool newIsDescending = isDescending(id0);
201 : {
202 0 : Bool negativeWidths = false;
203 0 : for (uInt i = 0; i < newNUM_CHAN; i++)
204 : {
205 0 : if (newCHAN_WIDTH(i) < 0.)
206 : {
207 0 : negativeWidths = true;
208 0 : newCHAN_WIDTH(i) = -newCHAN_WIDTH(i);
209 : }
210 : }
211 0 : if (negativeWidths && verbose)
212 : {
213 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
214 : << " *** Encountered negative channel widths in SPECTRAL_WINDOW table."
215 0 : << LogIO::POST;
216 0 : negChanWidthWarned = true;
217 : }
218 : }
219 :
220 : // Need to reverse the order for processing
221 0 : if (newIsDescending)
222 : {
223 0 : Vector<Double> tempF, tempW;
224 0 : tempF.assign(newCHAN_FREQ);
225 0 : tempW.assign(newCHAN_WIDTH);
226 0 : for (uInt i = 0; i < newNUM_CHAN; i++)
227 : {
228 0 : newCHAN_FREQ(i) = tempF(newNUM_CHAN - 1 - i);
229 0 : newCHAN_WIDTH(i) = tempW(newNUM_CHAN - 1 - i);
230 : }
231 : }
232 :
233 0 : Vector<Double> newEFFECTIVE_BW(effectiveBWColr(id0));
234 0 : Int newMEAS_FREQ_REF = measFreqRefColr(id0);
235 0 : Vector<Double> newRESOLUTION(resolutionColr(id0));
236 :
237 0 : vector<Int> averageN; // For each new channel store the number of old channels to average over
238 : //vector<vector<Int> > averageWhichSPW; // For each new channel store the (old) SPWs to average over
239 : //vector<vector<Int> > averageWhichChan; // For each new channel store the channel numbers to av. over
240 : //vector<vector<Double> > averageChanFrac; // For each new channel store the channel fraction for each old channel
241 :
242 : // Initialize the averaging vectors
243 0 : for (uInt i = 0; i < newNUM_CHAN; i++)
244 : {
245 0 : averageN.push_back(1);
246 0 : vector<Int> tv; // Just a temporary auxiliary vector
247 0 : tv.push_back(id0);
248 0 : averageWhichSPW.push_back(tv);
249 0 : if (newIsDescending)
250 : {
251 0 : tv[0] = newNUM_CHAN - 1 - i;
252 : }
253 : else
254 : {
255 0 : tv[0] = i;
256 : }
257 0 : averageWhichChan.push_back(tv);
258 0 : vector<Double> tvd; // another one
259 0 : tvd.push_back(1.);
260 0 : averageChanFrac.push_back(tvd);
261 : }
262 :
263 0 : if (verbose)
264 : {
265 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
266 : << "Input SPWs sorted by first (lowest) channel frequency:"
267 0 : << LogIO::POST;
268 :
269 0 : ostringstream oss; // needed for iomanip functions
270 0 : oss << " SPW " << std::setw(3) << id0 << ": " << std::setw(5)
271 0 : << newNUM_CHAN << " channels, first channel = "
272 0 : << std::setprecision(9) << std::setw(14) << std::scientific
273 0 : << newCHAN_FREQ(0) << " Hz";
274 0 : if (newNUM_CHAN > 1)
275 : {
276 : oss << ", last channel = " << std::setprecision(9)
277 0 : << std::setw(14) << std::scientific << newCHAN_FREQ(newNUM_CHAN - 1) << " Hz";
278 : }
279 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
280 0 : << oss.str() << LogIO::POST;
281 : }
282 :
283 : // Loop over remaining given SPWs
284 0 : for (uInt i = 1; i < nSpwsToCombine; i++)
285 : {
286 0 : Int idi = spwsSorted[i];
287 :
288 0 : uInt newNUM_CHANi = numChanColr(idi);
289 0 : Vector<Double> newCHAN_FREQi(chanFreqColr(idi));
290 0 : Vector<Double> newCHAN_WIDTHi(chanWidthColr(idi));
291 0 : Bool newIsDescendingI = isDescending(idi);
292 : {
293 0 : Bool negativeWidths = false;
294 0 : for (uInt ii = 0; ii < newNUM_CHANi; ii++)
295 : {
296 0 : if (newCHAN_WIDTHi(ii) < 0.)
297 : {
298 0 : negativeWidths = true;
299 0 : newCHAN_WIDTHi(ii) = -newCHAN_WIDTHi(ii);
300 : }
301 : }
302 0 : if (negativeWidths && !negChanWidthWarned && verbose)
303 : {
304 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
305 : << " *** Encountered negative channel widths in SPECTRAL_WINDOW table."
306 0 : << LogIO::POST;
307 0 : negChanWidthWarned = true;
308 : }
309 : }
310 :
311 : // need to reverse the order for processing
312 0 : if (newIsDescendingI)
313 : {
314 0 : Vector<Double> tempF, tempW;
315 0 : tempF.assign(newCHAN_FREQi);
316 0 : tempW.assign(newCHAN_WIDTHi);
317 0 : for (uInt ii = 0; ii < newNUM_CHANi; ii++)
318 : {
319 0 : newCHAN_FREQi(ii) = tempF(newNUM_CHANi - 1 - ii);
320 0 : newCHAN_WIDTHi(ii) = tempW(newNUM_CHANi - 1 - ii);
321 : }
322 : }
323 :
324 0 : Vector<Double> newEFFECTIVE_BWi(effectiveBWColr(idi));
325 0 : Int newMEAS_FREQ_REFi = measFreqRefColr(idi);
326 0 : Vector<Double> newRESOLUTIONi(resolutionColr(idi));
327 :
328 0 : if (verbose)
329 : {
330 0 : ostringstream oss;
331 0 : oss << " SPW " << std::setw(3) << idi << ": " << std::setw(5)
332 0 : << newNUM_CHANi << " channels, first channel = "
333 0 : << std::setprecision(9) << std::setw(14)
334 0 : << std::scientific << newCHAN_FREQi(0) << " Hz";
335 :
336 0 : if (newNUM_CHANi > 1)
337 : {
338 : oss << ", last channel = " << std::setprecision(9)
339 0 : << std::setw(14) << std::scientific
340 0 : << newCHAN_FREQi(newNUM_CHANi - 1) << " Hz";
341 : }
342 0 : os << LogIO::NORMAL << oss.str() << LogIO::POST;
343 : }
344 :
345 0 : vector<Double> mergedChanFreq;
346 0 : vector<Double> mergedChanWidth;
347 0 : vector<Double> mergedEffBW;
348 0 : vector<Double> mergedRes;
349 0 : vector<Int> mergedAverageN;
350 0 : vector<vector<Int> > mergedAverageWhichSPW;
351 0 : vector<vector<Int> > mergedAverageWhichChan;
352 0 : vector<vector<Double> > mergedAverageChanFrac;
353 :
354 : // check for compatibility
355 0 : if (newMEAS_FREQ_REFi != newMEAS_FREQ_REF)
356 : {
357 0 : os << LogIO::WARN << LogOrigin("MSTransformRegridder", __FUNCTION__)
358 : << "SPW " << idi
359 : << " cannot be combined with SPW " << id0
360 0 : << ". Non-matching ref. frame." << LogIO::POST;
361 0 : return false;
362 : }
363 :
364 : // Append or prepend spw to new spw overlap at all?
365 0 : if (newCHAN_FREQ(newNUM_CHAN - 1) + newCHAN_WIDTH(newNUM_CHAN - 1)
366 0 : / 2. < newCHAN_FREQi(0) - newCHAN_WIDTHi(0) / 2.)
367 : {
368 : // No overlap, and need to append
369 0 : for (uInt j = 0; j < newNUM_CHAN; j++)
370 : {
371 0 : mergedChanFreq.push_back(newCHAN_FREQ(j));
372 0 : mergedChanWidth.push_back(newCHAN_WIDTH(j));
373 0 : mergedEffBW.push_back(newEFFECTIVE_BW(j));
374 0 : mergedRes.push_back(newRESOLUTION(j));
375 0 : mergedAverageN.push_back(averageN[j]);
376 0 : mergedAverageWhichSPW.push_back(averageWhichSPW[j]);
377 0 : mergedAverageWhichChan.push_back(averageWhichChan[j]);
378 0 : mergedAverageChanFrac.push_back(averageChanFrac[j]);
379 : }
380 :
381 0 : vector<Int> tv;
382 0 : tv.push_back(idi); // Origin is SPW Id-i
383 0 : vector<Int> tv2;
384 0 : tv2.push_back(0);
385 0 : vector<Double> tvd;
386 0 : tvd.push_back(1.); // Fraction is 1.
387 0 : for (uInt j = 0; j < newNUM_CHANi; j++)
388 : {
389 0 : mergedChanFreq.push_back(newCHAN_FREQi(j));
390 0 : mergedChanWidth.push_back(newCHAN_WIDTHi(j));
391 0 : mergedEffBW.push_back(newEFFECTIVE_BWi(j));
392 0 : mergedRes.push_back(newRESOLUTIONi(j));
393 0 : mergedAverageN.push_back(1); // so far only one channel
394 0 : mergedAverageWhichSPW.push_back(tv);
395 0 : if (newIsDescendingI)
396 : {
397 0 : tv2[0] = newNUM_CHANi - 1 - j;
398 : }
399 : else
400 : {
401 0 : tv2[0] = j;
402 : }
403 0 : mergedAverageWhichChan.push_back(tv2); // channel number is j
404 0 : mergedAverageChanFrac.push_back(tvd);
405 : }
406 : }
407 0 : else if (newCHAN_FREQ(0) - newCHAN_WIDTH(0) / 2. > newCHAN_FREQi(
408 0 : newNUM_CHANi - 1) + newCHAN_WIDTHi(newNUM_CHANi - 1) / 2.)
409 : {
410 : // No overlap, need to prepend
411 0 : vector<Int> tv;
412 0 : tv.push_back(idi); // origin is SPW Id-i
413 0 : vector<Int> tv2;
414 0 : tv2.push_back(0);
415 0 : vector<Double> tvd;
416 0 : tvd.push_back(1.); // fraction is 1.
417 0 : for (uInt j = 0; j < newNUM_CHANi; j++)
418 : {
419 0 : mergedChanFreq.push_back(newCHAN_FREQi(j));
420 0 : mergedChanWidth.push_back(newCHAN_WIDTHi(j));
421 0 : mergedEffBW.push_back(newEFFECTIVE_BWi(j));
422 0 : mergedRes.push_back(newRESOLUTIONi(j));
423 0 : mergedAverageN.push_back(1); // so far only one channel
424 0 : mergedAverageWhichSPW.push_back(tv);
425 0 : if (newIsDescendingI)
426 : {
427 0 : tv2[0] = newNUM_CHANi - 1 - j;
428 : }
429 : else
430 : {
431 0 : tv2[0] = j;
432 : }
433 0 : mergedAverageWhichChan.push_back(tv2); // channel number is j
434 0 : mergedAverageChanFrac.push_back(tvd);
435 : }
436 :
437 0 : for (uInt j = 0; j < newNUM_CHAN; j++)
438 : {
439 0 : mergedChanFreq.push_back(newCHAN_FREQ(j));
440 0 : mergedChanWidth.push_back(newCHAN_WIDTH(j));
441 0 : mergedEffBW.push_back(newEFFECTIVE_BW(j));
442 0 : mergedRes.push_back(newRESOLUTION(j));
443 0 : mergedAverageN.push_back(averageN[j]);
444 0 : mergedAverageWhichSPW.push_back(averageWhichSPW[j]);
445 0 : mergedAverageWhichChan.push_back(averageWhichChan[j]);
446 0 : mergedAverageChanFrac.push_back(averageChanFrac[j]);
447 : }
448 : }
449 : // there is overlap
450 : else
451 : {
452 0 : Int id0StartChan = 0;
453 : // SPW Id-i starts before SPW Id-0
454 0 : if (newCHAN_FREQi(0) - newCHAN_WIDTHi(0) / 2. < newCHAN_FREQ(
455 0 : newNUM_CHAN - 1) - newCHAN_WIDTH(newNUM_CHAN - 1) / 2.)
456 : {
457 :
458 :
459 : // Some utilities for the averaging info
460 0 : vector<Int> tv; // temporary vector
461 0 : tv.push_back(idi); // origin is spw idi
462 0 : vector<Int> tv2;
463 0 : tv2.push_back(0);
464 0 : vector<Double> tvd;
465 0 : tvd.push_back(1.); // fraction is 1.
466 :
467 : // Find the first overlapping channel and prepend non-overlapping channels
468 0 : Double ubound0 = newCHAN_FREQ(0) + newCHAN_WIDTH(0) / 2.;
469 0 : Double lbound0 = newCHAN_FREQ(0) - newCHAN_WIDTH(0) / 2.;
470 0 : Double uboundk = 0.;
471 0 : Double lboundk = 0.;
472 : uInt k;
473 0 : for (k = 0; k < newNUM_CHANi; k++)
474 : {
475 0 : uboundk = newCHAN_FREQi(k) + newCHAN_WIDTHi(k) / 2.;
476 0 : lboundk = newCHAN_FREQi(k) - newCHAN_WIDTHi(k) / 2.;
477 0 : if (lbound0 < uboundk)
478 : {
479 0 : break;
480 : }
481 0 : mergedChanFreq.push_back(newCHAN_FREQi(k));
482 0 : mergedChanWidth.push_back(newCHAN_WIDTHi(k));
483 0 : mergedEffBW.push_back(newEFFECTIVE_BWi(k));
484 0 : mergedRes.push_back(newRESOLUTIONi(k));
485 0 : mergedAverageN.push_back(1); // So far only one channel
486 0 : mergedAverageWhichSPW.push_back(tv);
487 0 : if (newIsDescendingI)
488 : {
489 0 : tv2[0] = newNUM_CHANi - 1 - k;
490 : }
491 : else
492 : {
493 0 : tv2[0] = k;
494 : }
495 0 : mergedAverageWhichChan.push_back(tv2); // Channel number is k
496 0 : mergedAverageChanFrac.push_back(tvd);
497 : }
498 :
499 : // k-th is the one, actual overlap, need to merge channel k with channel 0
500 0 : if (lbound0 < uboundk && lboundk < lbound0)
501 : {
502 0 : Double newWidth = ubound0 - lboundk;
503 0 : Double newCenter = lboundk + newWidth / 2.;
504 0 : mergedChanFreq.push_back(newCenter);
505 0 : mergedChanWidth.push_back(newWidth);
506 0 : mergedEffBW.push_back(newWidth);
507 0 : mergedRes.push_back(newWidth);
508 0 : mergedAverageN.push_back(averageN[0] + 1); // one more channel contributes
509 :
510 : // Channel k is from SPW Id-i
511 0 : if (newIsDescendingI)
512 : {
513 0 : tv2[0] = newNUM_CHANi - 1 - k;
514 : }
515 : else
516 : {
517 0 : tv2[0] = k;
518 : }
519 :
520 0 : for (int j = 0; j < averageN[0]; j++)
521 : {
522 0 : tv.push_back(averageWhichSPW[0][j]); // additional contributors
523 0 : tv2.push_back(averageWhichChan[0][j]); // channel 0 from SPW Id-0
524 0 : tvd.push_back(averageChanFrac[0][j]);
525 : }
526 0 : mergedAverageWhichSPW.push_back(tv);
527 0 : mergedAverageWhichChan.push_back(tv2);
528 0 : mergedAverageChanFrac.push_back(tvd);
529 0 : id0StartChan = 1;
530 : }
531 : }
532 :
533 : // Now move along SPW id0 and merge until end of id0 is reached, then just copy
534 0 : for (uInt j = id0StartChan; j < newNUM_CHAN; j++)
535 : {
536 0 : mergedChanFreq.push_back(newCHAN_FREQ(j));
537 0 : mergedChanWidth.push_back(newCHAN_WIDTH(j));
538 0 : mergedEffBW.push_back(newEFFECTIVE_BW(j));
539 0 : mergedRes.push_back(newRESOLUTION(j));
540 :
541 0 : for (uInt k = 0; k < newNUM_CHANi; k++)
542 : {
543 0 : Double overlap_frac = 0.;
544 :
545 : // Does channel j in SPW Id-0 overlap with channel k in SPW Id-i?
546 0 : Double uboundj = newCHAN_FREQ(j) + newCHAN_WIDTH(j)/ 2.;
547 0 : Double uboundk = newCHAN_FREQi(k) + newCHAN_WIDTHi(k)/ 2.;
548 0 : Double lboundj = newCHAN_FREQ(j) - newCHAN_WIDTH(j)/ 2.;
549 0 : Double lboundk = newCHAN_FREQi(k) - newCHAN_WIDTHi(k)/ 2.;
550 :
551 : // Determine fraction
552 :
553 : // Chan k is completely covered by chan j
554 0 : if (lboundj <= lboundk && uboundk <= uboundj)
555 : {
556 0 : overlap_frac = 1.;
557 : }
558 : // chan j is completely covered by chan k
559 0 : else if (lboundk <= lboundj && uboundj <= uboundk)
560 : {
561 0 : overlap_frac = newCHAN_WIDTH(j) / newCHAN_WIDTHi(k);
562 : }
563 : // lower end of k is overlapping with j
564 0 : else if (lboundj < lboundk && lboundk < uboundj && uboundj < uboundk)
565 : {
566 0 : overlap_frac = (uboundj - lboundk) / newCHAN_WIDTHi(k);
567 : }
568 : // upper end of k is overlapping with j
569 0 : else if (lboundk < lboundj && lboundj < uboundk && lboundj < uboundk)
570 : {
571 0 : overlap_frac = (uboundk - lboundj) / newCHAN_WIDTHi(k);
572 : }
573 :
574 : // Update averaging info
575 0 : if (overlap_frac > 0.)
576 : {
577 0 : averageN[j] += 1;
578 0 : averageWhichSPW[j].push_back(idi);
579 0 : if (newIsDescendingI)
580 : {
581 0 : averageWhichChan[j].push_back(newNUM_CHANi - 1 - k);
582 : }
583 : else
584 : {
585 0 : averageWhichChan[j].push_back(k);
586 : }
587 0 : averageChanFrac[j].push_back(overlap_frac);
588 : }
589 : } // End loop over SPW Id-i
590 :
591 :
592 : // Append this channel with updated averaging info
593 0 : mergedAverageN.push_back(averageN[j]);
594 0 : mergedAverageWhichSPW.push_back(averageWhichSPW[j]);
595 0 : mergedAverageWhichChan.push_back(averageWhichChan[j]);
596 0 : mergedAverageChanFrac.push_back(averageChanFrac[j]);
597 :
598 : } // End loop over SPW Id-0
599 :
600 : // SPW Id-i still continues, find the last overlapping channel
601 0 : if (newCHAN_FREQ(newNUM_CHAN - 1) + newCHAN_WIDTH(
602 0 : newNUM_CHAN - 1) / 2. < newCHAN_FREQi(newNUM_CHANi - 1)
603 0 : + newCHAN_WIDTHi(newNUM_CHANi - 1) / 2.)
604 : {
605 0 : Int j = newNUM_CHAN - 1;
606 0 : Double uboundj = newCHAN_FREQ(j) + newCHAN_WIDTH(j) / 2.;
607 0 : Double lboundj = newCHAN_FREQ(j) - newCHAN_WIDTH(j) / 2.;
608 0 : Double uboundk = 0;
609 0 : Double lboundk = 0;
610 : Int k;
611 0 : for (k = newNUM_CHANi - 1; k >= 0; k--)
612 : {
613 0 : uboundk = newCHAN_FREQi(k) + newCHAN_WIDTHi(k) / 2.;
614 0 : lboundk = newCHAN_FREQi(k) - newCHAN_WIDTHi(k) / 2.;
615 0 : if (lboundk <= uboundj)
616 : {
617 0 : break;
618 : }
619 : }
620 :
621 : // k-th is the one
622 :
623 : // actual overlap
624 0 : if (lboundk < uboundj && uboundj < uboundk)
625 : {
626 0 : Double overlap_frac = (uboundj - lboundk) / newCHAN_WIDTHi(k);
627 :
628 : // Merge channel k completely with channel j
629 0 : if (overlap_frac > 0.01)
630 : {
631 0 : Double newWidth = uboundk - lboundj;
632 0 : Double newCenter = (lboundj + uboundk) / 2.;
633 0 : mergedChanFreq[j] = newCenter;
634 0 : mergedChanWidth[j] = newWidth;
635 0 : mergedEffBW[j] = newWidth;
636 0 : mergedRes[j] = newWidth;
637 0 : mergedAverageChanFrac[j][mergedAverageN[j] - 1] = 1.;
638 : }
639 : // Create separate, (slightly) more narrow channel
640 : else
641 : {
642 0 : Double newWidth = uboundk - uboundj;
643 0 : Double newCenter = (uboundj + uboundk) / 2.;
644 0 : vector<Int> tv;
645 0 : tv.push_back(idi); // Origin is SPW Id-i
646 0 : vector<Int> tv2;
647 0 : tv2.push_back(0);
648 0 : vector<Double> tvd;
649 0 : tvd.push_back(1.); // Fraction is 1.
650 0 : mergedChanFreq.push_back(newCenter);
651 0 : mergedChanWidth.push_back(newWidth);
652 0 : mergedEffBW.push_back(newWidth);
653 0 : mergedRes.push_back(newWidth);
654 0 : mergedAverageN.push_back(1); // So far only one channel
655 0 : mergedAverageWhichSPW.push_back(tv);
656 0 : if (newIsDescendingI)
657 : {
658 0 : tv2[0] = newNUM_CHANi - 1 - k;
659 : }
660 : else
661 : {
662 0 : tv2[0] = k;
663 : }
664 0 : mergedAverageWhichChan.push_back(tv2); // Channel number is k
665 0 : mergedAverageChanFrac.push_back(tvd);
666 : }
667 :
668 0 : k++; // Start appending remaining channels after k
669 : }
670 :
671 : // Append the remaining channels
672 0 : vector<Int> tv;
673 0 : tv.push_back(idi); // Origin is SPW Id-i
674 0 : vector<Int> tv2;
675 0 : tv2.push_back(0);
676 0 : vector<Double> tvd;
677 0 : tvd.push_back(1.); // Fraction is 1.
678 :
679 0 : for (uInt m = k; m < newNUM_CHANi; m++)
680 : {
681 0 : mergedChanFreq.push_back(newCHAN_FREQi(m));
682 0 : mergedChanWidth.push_back(newCHAN_WIDTHi(m));
683 0 : mergedEffBW.push_back(newEFFECTIVE_BWi(m));
684 0 : mergedRes.push_back(newRESOLUTIONi(m));
685 0 : mergedAverageN.push_back(1); // So far only one channel
686 0 : mergedAverageWhichSPW.push_back(tv);
687 :
688 0 : if (newIsDescendingI)
689 : {
690 0 : tv2[0] = newNUM_CHANi - 1 - m;
691 : }
692 : else
693 : {
694 0 : tv2[0] = m;
695 : }
696 0 : mergedAverageWhichChan.push_back(tv2); // Channel number is m
697 0 : mergedAverageChanFrac.push_back(tvd);
698 : }
699 : } // End if SPW Id-i still continues
700 : } // End if there is overlap
701 :
702 :
703 0 : newNUM_CHAN = mergedChanFreq.size();
704 0 : newCHAN_FREQ.assign(Vector<Double> (mergedChanFreq));
705 0 : newCHAN_WIDTH.assign(Vector<Double> (mergedChanWidth));
706 0 : newEFFECTIVE_BW.assign(Vector<Double> (mergedEffBW));
707 0 : newRESOLUTION.assign(Vector<Double> (mergedRes));
708 0 : averageN = mergedAverageN;
709 0 : averageWhichSPW.assign(mergedAverageWhichSPW.begin(), mergedAverageWhichSPW.end());
710 0 : averageChanFrac.assign(mergedAverageChanFrac.begin(), mergedAverageChanFrac.end());
711 0 : averageWhichChan.assign(mergedAverageWhichChan.begin(), mergedAverageWhichChan.end());
712 : }
713 :
714 0 : return true;
715 : }
716 :
717 : // -----------------------------------------------------------------------
718 : // A wrapper for regridChanBounds() which takes the user interface type re-gridding parameters
719 : // The ready-made grid is returned in newCHAN_FREQ and newCHAN_WIDTH
720 : // -----------------------------------------------------------------------
721 0 : Bool MSTransformRegridder::calcChanFreqs( LogIO& os,
722 : Vector<Double>& newCHAN_FREQ,
723 : Vector<Double>& newCHAN_WIDTH,
724 : Double& weightScale,
725 : const Vector<Double>& oldCHAN_FREQ,
726 : const Vector<Double>& oldCHAN_WIDTH,
727 : const MDirection phaseCenter,
728 : const MFrequency::Types theOldRefFrame,
729 : const MEpoch theObsTime,
730 : const MPosition mObsPos,
731 : const String& mode,
732 : const int nchan,
733 : const String& start,
734 : const String& width,
735 : const String& restfreq,
736 : const String& outframe,
737 : const String& veltype,
738 : const Bool verbose,
739 : const MRadialVelocity mRV)
740 : {
741 0 : Vector<Double> newChanLoBound;
742 0 : Vector<Double> newChanHiBound;
743 0 : String t_phasec;
744 :
745 0 : String t_mode;
746 0 : String t_outframe;
747 0 : String t_regridQuantity;
748 : Double t_restfreq;
749 0 : String t_regridInterpMeth;
750 : Double t_cstart;
751 : Double t_bandwidth;
752 : Double t_cwidth;
753 : Bool t_centerIsStart;
754 : Bool t_startIsEnd;
755 : Int t_nchan;
756 : Int t_width;
757 : Int t_start;
758 :
759 0 : if (!convertGridPars(os, mode, nchan, start,width,
760 : "linear", // A dummy value in this context
761 : restfreq, outframe,veltype,
762 : ////// output ////
763 : t_mode, t_outframe, t_regridQuantity, t_restfreq,
764 : t_regridInterpMeth, t_cstart, t_bandwidth, t_cwidth,
765 : t_centerIsStart, t_startIsEnd, t_nchan, t_width, t_start))
766 : {
767 : // An error occurred
768 0 : return false;
769 : }
770 :
771 : // Reference frame transformation
772 0 : Bool needTransform = true;
773 0 : Bool doRadVelCorr = false;
774 : MFrequency::Types theFrame;
775 0 : String oframe = outframe;
776 0 : oframe.upcase();
777 :
778 : // No output reference frame given
779 0 : if (outframe == "")
780 : {
781 : // Keep the reference frame as is
782 0 : theFrame = theOldRefFrame;
783 0 : needTransform = false;
784 : }
785 : // GEO transformation + radial velocity correction
786 0 : else if (oframe == "SOURCE")
787 : {
788 0 : theFrame = MFrequency::GEO;
789 0 : doRadVelCorr = true;
790 : }
791 0 : else if (!MFrequency::getType(theFrame, outframe))
792 : {
793 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
794 0 : << "Parameter \"outframe\" value " << outframe << " is invalid." << LogIO::POST;
795 0 : return false;
796 : }
797 0 : else if (theFrame == theOldRefFrame)
798 : {
799 0 : needTransform = false;
800 : }
801 :
802 0 : uInt oldNUM_CHAN = oldCHAN_FREQ.size();
803 0 : if (oldNUM_CHAN == 0)
804 : {
805 0 : newCHAN_FREQ.resize(0);
806 0 : newCHAN_WIDTH.resize(0);
807 0 : return true;
808 : }
809 :
810 0 : if (oldNUM_CHAN != oldCHAN_WIDTH.size())
811 : {
812 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
813 : << "Internal error: inconsistent dimensions of input channel frequency and width arrays."
814 0 : << LogIO::POST;
815 0 : return false;
816 : }
817 :
818 0 : Vector<Double> absOldCHAN_WIDTH;
819 0 : absOldCHAN_WIDTH.assign(oldCHAN_WIDTH);
820 0 : Bool negativeWidths = false;
821 0 : for (uInt i = 0; i < oldCHAN_WIDTH.size(); i++)
822 : {
823 0 : if (oldCHAN_WIDTH(i) < 0.)
824 : {
825 0 : negativeWidths = true;
826 0 : absOldCHAN_WIDTH(i) = -oldCHAN_WIDTH(i);
827 : }
828 : }
829 :
830 0 : if (negativeWidths && verbose)
831 : {
832 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
833 : << " *** Encountered negative channel widths in input spectral window."
834 0 : << LogIO::POST;
835 : }
836 :
837 0 : Vector<Double> transNewXin;
838 0 : Vector<Double> transCHAN_WIDTH(oldNUM_CHAN);
839 :
840 0 : if (needTransform)
841 : {
842 0 : transNewXin.resize(oldNUM_CHAN);
843 :
844 : // Set up conversion
845 0 : Unit unit(String("Hz"));
846 : MFrequency::Ref fromFrame = MFrequency::Ref(theOldRefFrame,
847 0 : MeasFrame(phaseCenter, mObsPos, theObsTime));
848 : MFrequency::Ref toFrame = MFrequency::Ref(theFrame,
849 0 : MeasFrame(phaseCenter, mObsPos, theObsTime));
850 0 : MFrequency::Convert freqTrans(unit, fromFrame, toFrame);
851 :
852 0 : MDoppler radVelCorr;
853 0 : Bool radVelSignificant = false;
854 :
855 : // Prepare correction for radial velocity if requested and possible
856 0 : if (doRadVelCorr)
857 : {
858 0 : Quantity mrv = mRV.get("m/s");
859 0 : radVelCorr = MDoppler(-mrv); // NOTE: Opposite sign to achieve correction!
860 0 : Double mRVVal = mrv.getValue();
861 0 : if (fabs(mRVVal) > 1E-6)
862 : {
863 0 : radVelSignificant = true;
864 : }
865 :
866 0 : if (verbose)
867 : {
868 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__)
869 : << "Note: The given additional radial velocity of "
870 : << mRVVal << " m/s will be taken into account."
871 0 : << LogIO::POST;
872 : }
873 : }
874 :
875 0 : for (uInt i = 0; i < oldNUM_CHAN; i++)
876 : {
877 0 : transNewXin[i] = freqTrans(oldCHAN_FREQ[i]).get(unit).getValue();
878 :
879 : // eliminate possible offsets
880 0 : transCHAN_WIDTH[i] = freqTrans(oldCHAN_FREQ[i] + absOldCHAN_WIDTH[i] / 2.).
881 0 : get(unit).getValue() - freqTrans(oldCHAN_FREQ[i] - absOldCHAN_WIDTH[i] / 2.).
882 0 : get(unit).getValue();
883 : }
884 :
885 : // Correct in addition for radial velocity
886 0 : if (radVelSignificant)
887 : {
888 0 : transNewXin = radVelCorr.shiftFrequency(transNewXin);
889 :
890 : //shiftFrequency is a scaling, so channel widths scale as well
891 0 : transCHAN_WIDTH = radVelCorr.shiftFrequency(transCHAN_WIDTH);
892 : }
893 :
894 : }
895 : else
896 : {
897 : // Just copy
898 0 : transNewXin.assign(oldCHAN_FREQ);
899 0 : transCHAN_WIDTH.assign(absOldCHAN_WIDTH);
900 : }
901 :
902 : // Calculate new grid
903 :
904 0 : String message;
905 :
906 0 : if (!regridChanBounds( newChanLoBound, newChanHiBound, t_cstart,
907 : t_bandwidth, t_cwidth, t_restfreq, t_regridQuantity, transNewXin,
908 : transCHAN_WIDTH, message, t_centerIsStart, t_startIsEnd, t_nchan,
909 : t_width, t_start))
910 : {
911 : // There was an error
912 0 : os << LogIO::WARN << LogOrigin("MSTransformRegridder", __FUNCTION__) << message << LogIO::POST;
913 0 : return false;
914 : }
915 :
916 0 : if (verbose)
917 : {
918 0 : os << LogIO::NORMAL << LogOrigin("MSTransformRegridder", __FUNCTION__) << message << LogIO::POST;
919 : }
920 :
921 : // We have a useful set of channel boundaries
922 0 : uInt newNUM_CHAN = newChanLoBound.size();
923 :
924 : // Complete the calculation of the new channel centers and widths
925 : // from newNUM_CHAN, newChanLoBound, and newChanHiBound
926 0 : newCHAN_FREQ.resize(newNUM_CHAN);
927 0 : newCHAN_WIDTH.resize(newNUM_CHAN);
928 0 : for (uInt i = 0; i < newNUM_CHAN; i++)
929 : {
930 0 : newCHAN_FREQ[i] = (newChanLoBound[i] + newChanHiBound[i]) / 2.;
931 0 : newCHAN_WIDTH[i] = newChanHiBound[i] - newChanLoBound[i];
932 : }
933 :
934 0 : weightScale = newCHAN_WIDTH[0]/transCHAN_WIDTH[0];
935 :
936 0 : return true;
937 : }
938 :
939 : // -----------------------------------------------------------------------
940 : // Helper function for handling the re-gridding parameter user input
941 : // -----------------------------------------------------------------------
942 0 : Bool MSTransformRegridder::convertGridPars( LogIO& os,
943 : const String& mode,
944 : const int nchan,
945 : const String& start,
946 : const String& width,
947 : const String& interp,
948 : const String& restfreq,
949 : const String& outframe,
950 : const String& veltype,
951 : String& t_mode,
952 : String& t_outframe,
953 : String& t_regridQuantity,
954 : Double& t_restfreq,
955 : String& t_regridInterpMeth,
956 : Double& t_cstart,
957 : Double& t_bandwidth,
958 : Double& t_cwidth,
959 : Bool& t_centerIsStart,
960 : Bool& t_startIsEnd,
961 : Int& t_nchan,
962 : Int& t_width,
963 : Int& t_start)
964 : {
965 0 : Bool rstat(false);
966 :
967 : try {
968 :
969 0 : casacore::QuantumHolder qh;
970 0 : String error;
971 :
972 0 : t_mode = mode;
973 0 : t_restfreq = 0.;
974 0 : if (!restfreq.empty() && !(restfreq == "[]"))
975 : {
976 0 : if (qh.fromString(error, restfreq))
977 : {
978 0 : t_restfreq = qh.asQuantity().getValue("Hz");
979 : }
980 : else
981 : {
982 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
983 0 : << "restfreq: " << error << LogIO::POST;
984 0 : return false;
985 : }
986 : }
987 :
988 : // Determine grid
989 0 : t_cstart = -9e99; // Default value indicating that the original start of the SPW should be used
990 0 : t_bandwidth = -1.; // Default value indicating that the original width of the SPW should be used
991 0 : t_cwidth = -1.; // Default value indicating that the original channel width of the SPW should be used
992 0 : t_nchan = -1;
993 0 : t_width = 0;
994 0 : t_start = -1;
995 0 : t_startIsEnd = false; // false means that start specifies the lower end in frequency (default)
996 : // true means that start specifies the upper end in frequency
997 :
998 : // start was set
999 0 : if (!start.empty() && !(start == "[]"))
1000 : {
1001 0 : if (t_mode == "channel")
1002 : {
1003 0 : t_start = atoi(start.c_str());
1004 : }
1005 0 : if (t_mode == "channel_b")
1006 : {
1007 0 : t_cstart = Double(atoi(start.c_str()));
1008 : }
1009 0 : else if (t_mode == "frequency")
1010 : {
1011 0 : if (qh.fromString(error, start))
1012 : {
1013 0 : t_cstart = qh.asQuantity().getValue("Hz");
1014 : }
1015 : else
1016 : {
1017 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1018 0 : << "start: " << error << LogIO::POST;
1019 0 : return false;
1020 : }
1021 : }
1022 0 : else if (t_mode == "velocity")
1023 : {
1024 0 : if (qh.fromString(error, start))
1025 : {
1026 0 : t_cstart = qh.asQuantity().getValue("m/s");
1027 : }
1028 : else
1029 : {
1030 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1031 0 : << "start: " << error << LogIO::POST;
1032 0 : return false;
1033 : }
1034 : }
1035 : }
1036 :
1037 : // channel width was set
1038 0 : if (!width.empty() && !(width == "[]"))
1039 : {
1040 0 : if (t_mode == "channel")
1041 : {
1042 0 : Int w = atoi(width.c_str());
1043 0 : t_width = abs(w);
1044 0 : if (w < 0)
1045 : {
1046 0 : t_startIsEnd = true;
1047 : }
1048 : }
1049 0 : else if (t_mode == "channel_b")
1050 : {
1051 0 : Double w = atoi(width.c_str());
1052 0 : t_cwidth = abs(w);
1053 0 : if (w < 0)
1054 : {
1055 0 : t_startIsEnd = true;
1056 : }
1057 : }
1058 0 : else if (t_mode == "frequency")
1059 : {
1060 0 : if (qh.fromString(error, width))
1061 : {
1062 0 : Double w = qh.asQuantity().getValue("Hz");
1063 0 : t_cwidth = abs(w);
1064 0 : if (w < 0)
1065 : {
1066 0 : t_startIsEnd = true;
1067 : }
1068 : }
1069 : else
1070 : {
1071 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1072 0 : << "width: " << error << LogIO::POST;
1073 0 : return false;
1074 : }
1075 : }
1076 0 : else if (t_mode == "velocity")
1077 : {
1078 0 : if (qh.fromString(error, width))
1079 : {
1080 0 : Double w = qh.asQuantity().getValue("m/s");
1081 0 : t_cwidth = abs(w);
1082 0 : if (w >= 0)
1083 : {
1084 0 : t_startIsEnd = true;
1085 : }
1086 : }
1087 : else
1088 : {
1089 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1090 0 : << "width: " << error << LogIO::POST;
1091 0 : return false;
1092 : }
1093 : }
1094 : }
1095 : // width was not set
1096 : else
1097 : {
1098 : // For the velocity mode the default t_startIsEnd is true if the sign of width is not known
1099 0 : if (t_mode == "velocity")
1100 : {
1101 0 : t_startIsEnd = true;
1102 : }
1103 : }
1104 :
1105 : // Number of output channels was set
1106 0 : if (nchan > 0)
1107 : {
1108 0 : if (t_mode == "channel_b")
1109 : {
1110 0 : if (t_cwidth > 0)
1111 : {
1112 0 : t_bandwidth = Double(nchan * t_cwidth);
1113 : }
1114 : else
1115 : {
1116 0 : t_bandwidth = Double(nchan);
1117 : }
1118 : }
1119 : else
1120 : {
1121 0 : t_nchan = nchan;
1122 : }
1123 : }
1124 :
1125 0 : if (t_mode == "channel")
1126 : {
1127 0 : t_regridQuantity = "freq";
1128 : }
1129 0 : else if (t_mode == "channel_b")
1130 : {
1131 0 : t_regridQuantity = "chan";
1132 : }
1133 0 : else if (t_mode == "frequency")
1134 : {
1135 0 : t_regridQuantity = "freq";
1136 : }
1137 0 : else if (t_mode == "velocity")
1138 : {
1139 0 : if (t_restfreq == 0.)
1140 : {
1141 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1142 : << "Need to set restfreq in velocity mode."
1143 0 : << LogIO::POST;
1144 0 : return false;
1145 : }
1146 :
1147 0 : t_regridQuantity = "vrad";
1148 0 : if (veltype == "optical")
1149 : {
1150 0 : t_regridQuantity = "vopt";
1151 : }
1152 0 : else if (veltype != "radio")
1153 : {
1154 0 : os << LogIO::WARN << LogOrigin("MSTransformRegridder", __FUNCTION__)
1155 : << "Invalid velocity type " << veltype
1156 0 : << ", setting type to \"radio\"" << LogIO::POST;
1157 : }
1158 :
1159 : }
1160 : else
1161 : {
1162 0 : os << LogIO::WARN << LogOrigin("MSTransformRegridder", __FUNCTION__)
1163 0 : << "Invalid mode " << t_mode << LogIO::POST;
1164 0 : return false;
1165 : }
1166 :
1167 0 : t_outframe = outframe;
1168 0 : t_regridInterpMeth = interp;
1169 0 : t_centerIsStart = true;
1170 :
1171 : // End prepare re-gridding parameters
1172 0 : rstat = true;
1173 :
1174 : }
1175 0 : catch (AipsError x)
1176 : {
1177 0 : os << LogIO::SEVERE << LogOrigin("MSTransformRegridder", __FUNCTION__)
1178 0 : << "Exception Reported: " << x.getMesg() << LogIO::POST;
1179 0 : rstat = false;
1180 : }
1181 :
1182 0 : return rstat;
1183 : }
1184 :
1185 : // -----------------------------------------------------------------------
1186 : // Calculate the final new channel boundaries from the re-regridding parameters and
1187 : // the old channel boundaries (already transformed to the desired reference frame).
1188 : // Returns false if input parameters were invalid and no useful boundaries could be created
1189 : // -----------------------------------------------------------------------
1190 0 : Bool MSTransformRegridder::regridChanBounds( Vector<Double>& newChanLoBound,
1191 : Vector<Double>& newChanHiBound,
1192 : const Double regridCenterC,
1193 : const Double regridBandwidth,
1194 : const Double regridChanWidthC,
1195 : const Double regridVeloRestfrq,
1196 : const String regridQuant,
1197 : const Vector<Double>& transNewXinC,
1198 : const Vector<Double>& transCHAN_WIDTHC,
1199 : String& message,
1200 : const Bool centerIsStartC,
1201 : const Bool startIsEndC,
1202 : const Int nchanC,
1203 : const Int width,
1204 : const Int startC)
1205 : {
1206 0 : ostringstream oss;
1207 :
1208 0 : Vector<Double> transNewXin(transNewXinC);
1209 0 : Vector<Double> transCHAN_WIDTH(transCHAN_WIDTHC);
1210 0 : Bool centerIsStart = centerIsStartC;
1211 0 : Bool startIsEnd = startIsEndC;
1212 0 : Double regridChanWidth = regridChanWidthC;
1213 0 : Double regridCenter = regridCenterC;
1214 0 : Int nchan = nchanC;
1215 0 : Int start = startC;
1216 :
1217 0 : Int oldNUM_CHAN = transNewXin.size();
1218 :
1219 : // Detect spectral windows defined with descending frequency
1220 0 : Bool isDescending = false;
1221 0 : for (uInt i = 1; i < transNewXin.size(); i++)
1222 : {
1223 0 : if (transNewXin(i) < transNewXin(i - 1))
1224 : {
1225 0 : isDescending = true;
1226 : }
1227 : // i.e. descending was detected but now we encounter ascending
1228 0 : else if (isDescending)
1229 : {
1230 0 : oss << "Channel frequencies are neither in ascending nor in descending order. Cannot process.";
1231 0 : message = oss.str();
1232 0 : return false;
1233 : }
1234 : }
1235 :
1236 : // Need to reverse the order for processing and later reverse the result
1237 0 : if (isDescending)
1238 : {
1239 0 : uInt n = transNewXin.size();
1240 0 : Vector<Double> tempF, tempW;
1241 0 : tempF.assign(transNewXin);
1242 0 : tempW.assign(transCHAN_WIDTH);
1243 0 : for (uInt i = 0; i < n; i++)
1244 : {
1245 0 : transNewXin(i) = tempF(n - 1 - i);
1246 0 : transCHAN_WIDTH(i) = tempW(n - 1 - i);
1247 : }
1248 :
1249 : // Also need to adjust the start values
1250 0 : if (startC >= 0)
1251 : {
1252 0 : start = n - 1 - startC;
1253 0 : if (centerIsStartC)
1254 : {
1255 0 : startIsEnd = !startIsEnd;
1256 : }
1257 : }
1258 : }
1259 :
1260 : // Verify regridCenter, regridBandwidth, and regridChanWidth
1261 : // Note: these are in the units corresponding to regridQuant!
1262 :
1263 0 : if (regridQuant == "chan")
1264 : {
1265 : // Channel numbers ...
1266 0 : Int regridCenterChan = -1;
1267 0 : Int regridBandwidthChan = -1;
1268 0 : Int regridChanWidthChan = -1;
1269 :
1270 : // Not set
1271 0 : if (regridCenter < -1E30)
1272 : {
1273 : // Find channel center closest to center of bandwidth
1274 0 : lDouble BWCenterF = (transNewXin[0] + transNewXin[oldNUM_CHAN - 1]) / 2.;
1275 0 : for (Int i = 0; i < oldNUM_CHAN; i++)
1276 : {
1277 0 : if (transNewXin[i] >= BWCenterF)
1278 : {
1279 0 : regridCenterChan = i;
1280 0 : break;
1281 : }
1282 : }
1283 0 : centerIsStart = false;
1284 : }
1285 : // Valid input
1286 0 : else if (0. <= regridCenter && regridCenter < Double(oldNUM_CHAN))
1287 : {
1288 0 : regridCenterChan = (Int) floor(regridCenter);
1289 : }
1290 : // Invalid
1291 : else
1292 : {
1293 0 : if (centerIsStart)
1294 : {
1295 0 : oss << "SPW start ";
1296 : }
1297 : else
1298 : {
1299 0 : oss << "SPW center ";
1300 : }
1301 0 : oss << regridCenter << " outside valid range which is " << 0 << " - " << oldNUM_CHAN - 1 << ".";
1302 0 : message = oss.str();
1303 0 : return false;
1304 : }
1305 :
1306 : // Not set or nchan set
1307 0 : if (regridBandwidth <= 0. || nchan > 0)
1308 : {
1309 0 : if (nchan > 0)
1310 : {
1311 0 : regridBandwidthChan = nchan;
1312 : }
1313 : else
1314 : {
1315 0 : regridBandwidthChan = oldNUM_CHAN;
1316 : }
1317 : }
1318 : else
1319 : {
1320 0 : regridBandwidthChan = (Int) floor(regridBandwidth);
1321 : }
1322 :
1323 0 : if (centerIsStart)
1324 : {
1325 0 : if (startIsEnd)
1326 : {
1327 0 : regridCenterChan = regridCenterChan - regridBandwidthChan / 2;
1328 : }
1329 : else
1330 : {
1331 0 : regridCenterChan = regridCenterChan + regridBandwidthChan / 2;
1332 : }
1333 0 : centerIsStart = false;
1334 : }
1335 :
1336 : // Center too close to lower edge
1337 0 : if (regridCenterChan - regridBandwidthChan / 2 < 0)
1338 : {
1339 0 : regridBandwidthChan = 2 * regridCenterChan + 1;
1340 0 : oss << " *** Requested output SPW width too large." << endl;
1341 : }
1342 : // Center too close to upper edge
1343 0 : if (oldNUM_CHAN < regridCenterChan + regridBandwidthChan / 2)
1344 : {
1345 0 : regridBandwidthChan = 2 * (oldNUM_CHAN - regridCenterChan);
1346 0 : oss << " *** Requested output SPW width too large." << endl;
1347 : }
1348 :
1349 0 : if (regridChanWidth < 1.)
1350 : {
1351 0 : regridChanWidthChan = 1;
1352 : }
1353 0 : else if (regridChanWidth > Double(regridBandwidthChan))
1354 : {
1355 0 : regridChanWidthChan = regridBandwidthChan; // i.e. SPW = a single channel
1356 0 : oss << " *** Requested output channel width too large. Adjusted to maximum possible value." << endl;
1357 : }
1358 : // Valid input
1359 : else
1360 : {
1361 0 : regridChanWidthChan = (Int) floor(regridChanWidth);
1362 0 : if (nchan > 0)
1363 : {
1364 0 : regridBandwidthChan = nchan * regridChanWidthChan;
1365 : }
1366 : }
1367 :
1368 0 : if (regridBandwidthChan != floor(regridBandwidth))
1369 : {
1370 0 : oss << " *** Output SPW width set to " << regridBandwidthChan << " original channels" << endl;
1371 0 : oss << " in an attempt to keep center of output SPW close to center of requested SPW." << endl;
1372 : }
1373 :
1374 : // Calculate newChanLoBound and newChanHiBound from
1375 : // regridCenterChan, regridBandwidthChan, and regridChanWidthChan
1376 0 : Int bwLowerEndChan = regridCenterChan - regridBandwidthChan / 2;
1377 0 : Int bwUpperEndChan = bwLowerEndChan + regridBandwidthChan - 1;
1378 0 : Int numNewChanDown = 0;
1379 0 : Int numNewChanUp = 0;
1380 :
1381 : // Only one new channel
1382 0 : if (regridChanWidthChan == regridBandwidthChan)
1383 : {
1384 0 : newChanLoBound.resize(1);
1385 0 : newChanHiBound.resize(1);
1386 0 : newChanLoBound[0] = transNewXin[bwLowerEndChan] - transCHAN_WIDTH[bwLowerEndChan] / 2.;
1387 0 : newChanHiBound[0] = transNewXin[bwUpperEndChan] + transCHAN_WIDTH[bwUpperEndChan] / 2.;
1388 0 : numNewChanUp = 1;
1389 : }
1390 : // Have more than one new channel
1391 : else
1392 : {
1393 : // Need to accommodate the possibility that the original channels are not contiguous!
1394 :
1395 : // The numbers of the Channels from which the lower bounds will be taken for the new channels
1396 0 : vector<Int> loNCBup;
1397 : // Starting from the central channel going up
1398 0 : vector<Int> hiNCBup; // The numbers of the Channels from which the high
1399 : // Bounds will be taken for the new channels starting from the central channel going up
1400 0 : vector<Int> loNCBdown; // The numbers of the Channels from which the
1401 : // Lower bounds will be taken for the new channels starting from the central channel going down
1402 0 : vector<Int> hiNCBdown; // The numbers of the Channels from which the high bounds will be taken
1403 : // for the new channels starting from the central channel going down.
1404 : // Want to keep the center of the center channel at the center of the new center
1405 : // channel if the bandwidth is an odd multiple of the new channel width
1406 : // otherwise the center channel is the lower edge of the new center channel
1407 : Int startChan;
1408 0 : lDouble tnumChan = regridBandwidthChan / regridChanWidthChan;
1409 0 : if ((Int) tnumChan % 2 != 0)
1410 : {
1411 : // Odd multiple
1412 0 : startChan = regridCenterChan - regridChanWidthChan / 2;
1413 : }
1414 : else
1415 : {
1416 0 : startChan = regridCenterChan;
1417 : }
1418 :
1419 : // Upper half
1420 0 : for (Int i = startChan; i <= bwUpperEndChan; i += regridChanWidthChan)
1421 : {
1422 0 : loNCBup.push_back(i);
1423 0 : if (i + regridChanWidthChan - 1 <= bwUpperEndChan)
1424 : {
1425 : // Can go one more normal step up
1426 0 : hiNCBup.push_back(i + regridChanWidthChan - 1);
1427 : }
1428 : else
1429 : {
1430 : // Create narrower channels at the edges if necessary
1431 0 : oss << " *** Last channel at upper edge of new SPW made only "
1432 0 : << bwUpperEndChan - i + 1
1433 0 : << " original channels wide to fit given total bandwidth."
1434 0 : << endl;
1435 0 : hiNCBup.push_back(bwUpperEndChan);
1436 : }
1437 : }
1438 :
1439 : // Lower half
1440 0 : for (Int i = startChan - 1; i >= bwLowerEndChan; i -= regridChanWidthChan)
1441 : {
1442 0 : hiNCBdown.push_back(i);
1443 0 : if (i - regridChanWidthChan + 1 >= bwLowerEndChan)
1444 : {
1445 : // Can go one more normal step down
1446 0 : loNCBdown.push_back(i - regridChanWidthChan + 1);
1447 : }
1448 : else
1449 : {
1450 : // Create narrower channels at the edges if necessary
1451 0 : oss << " *** First channel at lower edge of new SPW made only "
1452 0 : << i - bwLowerEndChan + 1
1453 0 : << " original channels wide to fit given total bandwidth."
1454 0 : << endl;
1455 0 : loNCBdown.push_back(bwLowerEndChan);
1456 : }
1457 : }
1458 :
1459 : // The number of channels below the central one
1460 0 : numNewChanDown = loNCBdown.size();
1461 :
1462 : // The number of channels above and including the central one
1463 0 : numNewChanUp = loNCBup.size();
1464 :
1465 0 : newChanLoBound.resize(numNewChanDown + numNewChanUp);
1466 0 : newChanHiBound.resize(numNewChanDown + numNewChanUp);
1467 0 : for (Int i = 0; i < numNewChanDown; i++)
1468 : {
1469 0 : Int k = numNewChanDown - i - 1; // Need to assign in reverse
1470 0 : newChanLoBound[i] = transNewXin[loNCBdown[k]] - transCHAN_WIDTH[loNCBdown[k]] / 2.;
1471 0 : newChanHiBound[i] = transNewXin[hiNCBdown[k]] + transCHAN_WIDTH[hiNCBdown[k]] / 2.;
1472 : }
1473 :
1474 0 : for (Int i = 0; i < numNewChanUp; i++)
1475 : {
1476 0 : newChanLoBound[i + numNewChanDown] = transNewXin[loNCBup[i]] - transCHAN_WIDTH[loNCBup[i]] / 2.;
1477 0 : newChanHiBound[i + numNewChanDown] = transNewXin[hiNCBup[i]] + transCHAN_WIDTH[hiNCBup[i]] / 2.;
1478 : }
1479 : } // end if
1480 :
1481 0 : oss << " New channels defined based on original channels" << endl
1482 0 : << " Central channel contains original channel "
1483 0 : << regridCenterChan << endl << " Channel width = "
1484 0 : << regridChanWidthChan << " original channels" << endl
1485 0 : << " Total width of SPW = " << regridBandwidthChan
1486 0 : << " original channels == " << numNewChanDown + numNewChanUp
1487 0 : << " new channels" << endl;
1488 :
1489 0 : uInt nc = newChanLoBound.size();
1490 0 : oss << " Total width of SPW (in output frame) = "
1491 0 : << newChanHiBound[nc- 1] - newChanLoBound[0] << " Hz" << endl;
1492 0 : oss << " Lower edge = " << std::setprecision(9) << newChanLoBound[0] << " Hz,"
1493 0 : << " upper edge = " << std::setprecision(9) << newChanHiBound[nc - 1] << " Hz" << endl;
1494 :
1495 0 : if (isDescending)
1496 : {
1497 0 : Vector<Double> tempL, tempU;
1498 0 : tempL.assign(newChanLoBound);
1499 0 : tempU.assign(newChanHiBound);
1500 0 : for (uInt i = 0; i < nc; i++)
1501 : {
1502 0 : newChanLoBound(i) = tempL(nc - 1 - i);
1503 0 : newChanHiBound(i) = tempU(nc - 1 - i);
1504 : }
1505 : }
1506 :
1507 0 : message = oss.str();
1508 :
1509 0 : return true;
1510 : }
1511 : // We operate on real numbers /////////////////
1512 : else
1513 : {
1514 : // First transform them to frequencies
1515 0 : lDouble regridCenterF = -1.; // Initialize as "not set"
1516 0 : lDouble regridBandwidthF = -1.;
1517 0 : lDouble regridChanWidthF = -1.;
1518 :
1519 0 : if (regridQuant == "vrad")
1520 : {
1521 : // radio velocity need restfrq
1522 0 : if (regridVeloRestfrq < -1E30) // means "not set"
1523 : {
1524 0 : oss << "Parameter \"restfreq\" needs to be set if regrid_quantity==vrad. Cannot proceed with regridSpw ...";
1525 0 : message = oss.str();
1526 0 : return false;
1527 : }
1528 0 : else if (regridVeloRestfrq < 0. || regridVeloRestfrq > 1E30)
1529 : {
1530 0 : oss << "Parameter \"restfreq\" value " << regridVeloRestfrq << " is invalid.";
1531 0 : message = oss.str();
1532 0 : return false;
1533 : }
1534 :
1535 : lDouble regridCenterVel;
1536 0 : if (regridCenter > -C::c)
1537 : {
1538 : // (We deal with invalid values later)
1539 0 : if (centerIsStart)
1540 : {
1541 : Double tcWidth;
1542 0 : if (regridChanWidth > 0.)
1543 : {
1544 0 : tcWidth = regridChanWidth;
1545 : }
1546 : else
1547 : {
1548 0 : tcWidth = vrad( transNewXin[0] - transCHAN_WIDTH[0] / 2.,
1549 0 : regridVeloRestfrq) - vrad(
1550 0 : transNewXin[0] + transCHAN_WIDTH[0] / 2.,
1551 : regridVeloRestfrq);
1552 : }
1553 :
1554 : // start is the center of the last channel (in freq)
1555 0 : if (startIsEnd)
1556 : {
1557 0 : regridCenter -= tcWidth / 2.;
1558 : }
1559 : // start is the center of the first channel (in freq)
1560 : else
1561 : {
1562 0 : regridCenter += tcWidth / 2.;
1563 : }
1564 : }
1565 :
1566 0 : regridCenterF = freq_from_vrad(regridCenter, regridVeloRestfrq);
1567 :
1568 0 : regridCenterVel = regridCenter;
1569 : }
1570 : // center was not specified
1571 : else
1572 : {
1573 0 : regridCenterF = (transNewXin[0] + transNewXin[oldNUM_CHAN - 1])/ 2.;
1574 0 : regridCenterVel = vrad(regridCenterF, regridVeloRestfrq);
1575 0 : centerIsStart = false;
1576 : }
1577 0 : if (nchan > 0)
1578 : {
1579 0 : if (regridChanWidth > 0.)
1580 : {
1581 0 : lDouble chanUpperEdgeF = freq_from_vrad(regridCenterVel - regridChanWidth / 2.,regridVeloRestfrq);
1582 0 : regridChanWidthF = 2. * (chanUpperEdgeF - regridCenterF);
1583 : }
1584 : // take channel width from first channel
1585 : else
1586 : {
1587 0 : regridChanWidthF = transCHAN_WIDTH[0];
1588 : }
1589 :
1590 0 : regridBandwidthF = nchan * regridChanWidthF;
1591 : // Can convert start to center
1592 0 : if (centerIsStart)
1593 : {
1594 0 : if (startIsEnd)
1595 : {
1596 0 : regridCenterF = regridCenterF - regridBandwidthF / 2.;
1597 : }
1598 : else
1599 : {
1600 0 : regridCenterF = regridCenterF + regridBandwidthF / 2.;
1601 : }
1602 0 : centerIsStart = false;
1603 : }
1604 : }
1605 0 : else if (regridBandwidth > 0.)
1606 : {
1607 : // Can convert start to center
1608 0 : if (centerIsStart)
1609 : {
1610 0 : if (startIsEnd)
1611 : {
1612 0 : regridCenterVel = regridCenter + regridBandwidth / 2.;
1613 : }
1614 : else
1615 : {
1616 0 : regridCenterVel = regridCenter - regridBandwidth / 2.;
1617 : }
1618 0 : regridCenterF = freq_from_vrad(regridCenterVel,regridVeloRestfrq);
1619 0 : centerIsStart = false;
1620 : }
1621 0 : lDouble bwUpperEndF = freq_from_vrad(regridCenterVel - regridBandwidth / 2.,regridVeloRestfrq);
1622 0 : regridBandwidthF = 2. * (bwUpperEndF - regridCenterF);
1623 : }
1624 :
1625 0 : if (regridChanWidth > 0. && regridChanWidthF < 0.)
1626 : {
1627 0 : lDouble chanUpperEdgeF = freq_from_vrad(regridCenterVel - regridChanWidth / 2.,regridVeloRestfrq);
1628 0 : regridChanWidthF = 2. * (chanUpperEdgeF - freq_from_vrad(regridCenterVel, regridVeloRestfrq));
1629 : }
1630 : }
1631 0 : else if (regridQuant == "vopt")
1632 : {
1633 : // Optical velocity need restfrq
1634 0 : if (regridVeloRestfrq < -1E30) // means "not set"
1635 : {
1636 0 : oss << "Parameter \"restfreq\" needs to be set if regrid_quantity==vopt. Cannot proceed with regridSpw ...";
1637 0 : message = oss.str();
1638 0 : return false;
1639 : }
1640 0 : else if (regridVeloRestfrq <= 0. || regridVeloRestfrq > 1E30)
1641 : {
1642 0 : oss << "Parameter \"restfreq\" value " << regridVeloRestfrq << " is invalid.";
1643 0 : message = oss.str();
1644 0 : return false;
1645 : }
1646 :
1647 : lDouble regridCenterVel;
1648 0 : if (regridCenter > -C::c)
1649 : {
1650 0 : if (centerIsStart)
1651 : {
1652 : Double tcWidth;
1653 0 : if (regridChanWidth > 0.)
1654 : {
1655 0 : tcWidth = regridChanWidth;
1656 : }
1657 : else
1658 : {
1659 0 : tcWidth = vopt( transNewXin[0] - transCHAN_WIDTH[0] / 2.,
1660 0 : regridVeloRestfrq) - vopt(
1661 0 : transNewXin[0] + transCHAN_WIDTH[0] / 2.,
1662 : regridVeloRestfrq);
1663 : }
1664 :
1665 : // start is the center of the last channel (in freq)
1666 0 : if (startIsEnd)
1667 : {
1668 0 : regridCenter -= tcWidth / 2.;
1669 : }
1670 : // start is the center of the first channel (in freq)
1671 : else
1672 : {
1673 0 : regridCenter += tcWidth / 2.;
1674 : }
1675 : }
1676 :
1677 : // (We deal with invalid values later)
1678 0 : regridCenterF = freq_from_vopt(regridCenter, regridVeloRestfrq);
1679 0 : regridCenterVel = regridCenter;
1680 : }
1681 : // Center was not specified
1682 : else
1683 : {
1684 0 : regridCenterF = (transNewXin[0] - transCHAN_WIDTH[0]
1685 0 : + transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1]) / 2.;
1686 0 : regridCenterVel = vopt(regridCenterF, regridVeloRestfrq);
1687 0 : centerIsStart = false;
1688 : }
1689 :
1690 0 : if (nchan > 0)
1691 : {
1692 : lDouble cw;
1693 0 : lDouble divbytwo = 0.5;
1694 0 : if (centerIsStart)
1695 : {
1696 0 : divbytwo = 1.;
1697 : }
1698 0 : if (regridChanWidth > 0.)
1699 : {
1700 0 : cw = regridChanWidth;
1701 : }
1702 : // Determine channel width from first channel
1703 : else
1704 : {
1705 0 : lDouble upEdge = vopt(transNewXin[0] - transCHAN_WIDTH[0],regridVeloRestfrq);
1706 0 : lDouble loEdge = vopt(transNewXin[0] + transCHAN_WIDTH[0],regridVeloRestfrq);
1707 0 : cw = abs(upEdge - loEdge);
1708 : }
1709 0 : lDouble bwUpperEndF = 0.;
1710 :
1711 : // Start is end in velocity
1712 0 : if (centerIsStart && !startIsEnd)
1713 : {
1714 0 : bwUpperEndF = freq_from_vopt(regridCenterVel - (lDouble) nchan * cw * divbytwo,regridVeloRestfrq);
1715 : }
1716 : else
1717 : {
1718 0 : bwUpperEndF = freq_from_vopt(regridCenterVel + (lDouble) nchan * cw * divbytwo,regridVeloRestfrq);
1719 : }
1720 :
1721 0 : regridBandwidthF = abs(bwUpperEndF - regridCenterF) / divbytwo;
1722 :
1723 : // Can convert start to center
1724 0 : if (centerIsStart)
1725 : {
1726 0 : if (startIsEnd)
1727 : {
1728 0 : regridCenterVel = regridCenterVel + (lDouble) nchan * cw / 2.;
1729 : }
1730 : else
1731 : {
1732 0 : regridCenterVel = regridCenterVel - (lDouble) nchan * cw / 2.;
1733 : }
1734 :
1735 0 : regridCenterF = freq_from_vopt(regridCenterVel,regridVeloRestfrq);
1736 0 : centerIsStart = false;
1737 : }
1738 0 : nchan = 0; // indicate that nchan should not be used in the following
1739 : }
1740 0 : else if (regridBandwidth > 0.)
1741 : {
1742 : // can convert start to center
1743 0 : if (centerIsStart)
1744 : {
1745 0 : if (startIsEnd)
1746 : {
1747 0 : regridCenterVel = regridCenter + regridBandwidth / 2.;
1748 : }
1749 : else
1750 : {
1751 0 : regridCenterVel = regridCenter - regridBandwidth / 2.;
1752 : }
1753 0 : regridCenterF = freq_from_vopt(regridCenterVel, regridVeloRestfrq);
1754 0 : centerIsStart = false;
1755 : }
1756 :
1757 0 : lDouble bwUpperEndF = freq_from_vopt(regridCenterVel - regridBandwidth / 2.,regridVeloRestfrq);
1758 0 : regridBandwidthF = 2. * (bwUpperEndF - regridCenterF);
1759 : }
1760 :
1761 0 : if (regridChanWidth > 0. && regridChanWidthF < 0.)
1762 : {
1763 0 : lDouble chanUpperEdgeF = freq_from_vopt(regridCenterVel - regridChanWidth / 2.,regridVeloRestfrq);
1764 0 : regridChanWidthF = 2. * (chanUpperEdgeF - freq_from_vopt(regridCenterVel, regridVeloRestfrq));
1765 : }
1766 : }
1767 0 : else if (regridQuant == "freq")
1768 : {
1769 : // width parameter overrides regridChanWidth
1770 0 : if (width > 0)
1771 : {
1772 0 : regridChanWidth = width * transCHAN_WIDTH[0];
1773 : }
1774 :
1775 0 : if (start >= 0)
1776 : {
1777 0 : Int firstChan = start;
1778 0 : if (start >= (Int) transNewXin.size())
1779 : {
1780 0 : oss << " *** Parameter start exceeds total number of channels which is "
1781 0 : << transNewXin.size() << ". Set to 0." << endl;
1782 0 : firstChan = 0;
1783 0 : startIsEnd = false;
1784 : }
1785 :
1786 0 : if (startIsEnd)
1787 : {
1788 0 : regridCenter = transNewXin[firstChan] + transCHAN_WIDTH[firstChan] / 2.;
1789 : }
1790 : else
1791 : {
1792 0 : regridCenter = transNewXin[firstChan] - transCHAN_WIDTH[firstChan] / 2.;
1793 : }
1794 0 : centerIsStart = true;
1795 : }
1796 : else
1797 : {
1798 : // start is the center of the first channel
1799 0 : if (centerIsStart)
1800 : {
1801 : Double tcWidth;
1802 0 : if (regridChanWidth > 0.)
1803 : {
1804 0 : tcWidth = regridChanWidth;
1805 : }
1806 : else
1807 : {
1808 0 : tcWidth = transCHAN_WIDTH[0];
1809 : }
1810 :
1811 0 : if (startIsEnd)
1812 : {
1813 0 : regridCenter += tcWidth / 2.;
1814 : }
1815 : else
1816 : {
1817 0 : regridCenter -= tcWidth / 2.;
1818 : }
1819 : }
1820 : }
1821 0 : regridCenterF = regridCenter;
1822 0 : regridBandwidthF = regridBandwidth;
1823 0 : regridChanWidthF = regridChanWidth;
1824 : }
1825 0 : else if (regridQuant == "wave") {
1826 : // wavelength ...
1827 : lDouble regridCenterWav;
1828 0 : if (regridCenter > 0.)
1829 : {
1830 0 : if (centerIsStart)
1831 : {
1832 : Double tcWidth;
1833 0 : if (regridChanWidth > 0.)
1834 : {
1835 0 : tcWidth = regridChanWidth;
1836 : }
1837 : else
1838 : {
1839 0 : tcWidth = lambda(transNewXin[0] - transCHAN_WIDTH[0] / 2.) -
1840 0 : lambda(transNewXin[0] + transCHAN_WIDTH[0]/ 2.);
1841 : }
1842 :
1843 : // start is the center of the last channel (in freq)
1844 0 : if (startIsEnd)
1845 : {
1846 0 : regridCenter -= tcWidth / 2.;
1847 : }
1848 : // start is the center of the first channel (in freq)
1849 : else
1850 : {
1851 0 : regridCenter += tcWidth / 2.;
1852 : }
1853 : }
1854 0 : regridCenterF = freq_from_lambda(regridCenter);
1855 0 : regridCenterWav = regridCenter;
1856 : }
1857 : // center was not specified
1858 : else
1859 : {
1860 0 : regridCenterF = (transNewXin[0] + transNewXin[oldNUM_CHAN - 1])/ 2.;
1861 0 : regridCenterWav = lambda(regridCenterF);
1862 0 : centerIsStart = false;
1863 : }
1864 0 : if (nchan > 0)
1865 : {
1866 : lDouble cw;
1867 0 : lDouble divbytwo = 0.5;
1868 0 : if (centerIsStart)
1869 : {
1870 0 : divbytwo = 1.;
1871 : }
1872 :
1873 0 : if (regridChanWidth > 0.)
1874 : {
1875 0 : cw = regridChanWidth;
1876 : }
1877 : // Determine channel width from first channel
1878 : else
1879 : {
1880 0 : lDouble upEdge = lambda(transNewXin[0] - transCHAN_WIDTH[0]);
1881 0 : lDouble loEdge = lambda(transNewXin[0] + transCHAN_WIDTH[0]);
1882 0 : cw = abs(upEdge - loEdge);
1883 : }
1884 :
1885 0 : lDouble bwUpperEndF = 0.;
1886 0 : if (centerIsStart && !startIsEnd)
1887 : {
1888 0 : bwUpperEndF = freq_from_lambda(regridCenterWav - (lDouble) nchan * cw * divbytwo);
1889 : }
1890 : else
1891 : {
1892 0 : bwUpperEndF = freq_from_lambda(regridCenterWav + (lDouble) nchan * cw * divbytwo);
1893 : }
1894 :
1895 0 : regridBandwidthF = (bwUpperEndF - regridCenterF) / divbytwo;
1896 :
1897 : // Can convert start to center
1898 0 : if (centerIsStart)
1899 : {
1900 0 : if (startIsEnd)
1901 : {
1902 0 : regridCenterWav = regridCenterWav + (lDouble) nchan* cw / 2.;
1903 : }
1904 : else
1905 : {
1906 0 : regridCenterWav = regridCenterWav - (lDouble) nchan* cw / 2.;
1907 : }
1908 0 : regridCenterF = freq_from_lambda(regridCenterWav);
1909 0 : centerIsStart = false;
1910 : }
1911 0 : nchan = 0; // indicate that nchan should not be used in the following
1912 :
1913 : }
1914 0 : else if (regridBandwidth > 0. && regridBandwidth / 2.< regridCenterWav)
1915 : {
1916 : // Can convert start to center
1917 0 : if (centerIsStart)
1918 : {
1919 0 : if (startIsEnd)
1920 : {
1921 0 : regridCenterWav = regridCenter + regridBandwidth / 2.;
1922 : }
1923 : else
1924 : {
1925 0 : regridCenterWav = regridCenter - regridBandwidth / 2.;
1926 : }
1927 0 : regridCenterF = freq_from_lambda(regridCenterWav);
1928 0 : centerIsStart = false;
1929 : }
1930 0 : lDouble bwUpperEndF = lambda(regridCenterWav - regridBandwidth / 2.);
1931 0 : regridBandwidthF = 2. * (bwUpperEndF - regridCenterF);
1932 : }
1933 :
1934 0 : if (regridChanWidth > 0. && regridChanWidth / 2. < regridCenterWav)
1935 : {
1936 0 : lDouble chanUpperEdgeF = lambda(regridCenterWav - regridChanWidth / 2.);
1937 0 : regridChanWidthF = 2. * (chanUpperEdgeF - regridCenterF);
1938 : }
1939 : }
1940 : else
1941 : {
1942 0 : oss << "Invalid value " << regridQuant << " for parameter \"mode\".";
1943 0 : message = oss.str();
1944 0 : return false;
1945 : }
1946 : // (transformation of regrid parameters to frequencies completed)
1947 :
1948 : // Then determine the actually possible parameters
1949 : lDouble theRegridCenterF;
1950 : lDouble theRegridBWF;
1951 : lDouble theCentralChanWidthF;
1952 :
1953 : // For vrad and vopt also need to keep this adjusted value
1954 0 : lDouble theChanWidthX = -1.;
1955 :
1956 0 : if (regridCenterF < 0.) // means "not set"
1957 : {
1958 : // Keep regrid center as it is in the data
1959 0 : theRegridCenterF = (transNewXin[0] - transCHAN_WIDTH[0] / 2.
1960 0 : + transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.) / 2.;
1961 0 : centerIsStart = false;
1962 : }
1963 : // regridCenterF was set
1964 : else
1965 : {
1966 : // Keep center in limits
1967 0 : theRegridCenterF = regridCenterF;
1968 0 : if ((theRegridCenterF - (transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.)) > 1.) // 1 Hz tolerance
1969 : {
1970 0 : oss << "*** Requested center of SPW " << theRegridCenterF
1971 0 : << " Hz is too large by "
1972 0 : << theRegridCenterF - transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.
1973 0 : << " Hz\n";
1974 0 : theRegridCenterF = transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.;
1975 0 : oss << "*** Reset to maximum possible value " << theRegridCenterF << " Hz";
1976 : }
1977 0 : else if (theRegridCenterF < (transNewXin[0] - transCHAN_WIDTH[0] / 2.))
1978 : {
1979 0 : Double diff = (transNewXin[0] - transCHAN_WIDTH[0] / 2.) - theRegridCenterF;
1980 :
1981 : // Cope with numerical accuracy problems
1982 0 : if (diff > 1.)
1983 : {
1984 0 : oss << "*** Requested center of SPW (" << theRegridCenterF
1985 0 : << " Hz) is smaller than minimum possible value";
1986 0 : oss << " by " << diff << " Hz";
1987 : }
1988 :
1989 0 : theRegridCenterF = transNewXin[0] - transCHAN_WIDTH[0] / 2.;
1990 0 : if (diff > 1.)
1991 : {
1992 0 : oss << "\n*** Reset to minimum possible value " << theRegridCenterF << " Hz";
1993 : }
1994 : }
1995 : }
1996 :
1997 0 : if (regridBandwidthF <= 0. || nchan != 0) // "not set" or use nchan instead
1998 : {
1999 : // Keep bandwidth as is
2000 0 : theRegridBWF = transNewXin[oldNUM_CHAN - 1] - transNewXin[0]
2001 0 : + transCHAN_WIDTH[0] / 2. + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.;
2002 :
2003 : // Use nchan parameter if available
2004 0 : if (nchan != 0)
2005 : {
2006 0 : if (nchan < 0)
2007 : {
2008 0 : if (regridQuant == "freq" || regridQuant == "vrad") // i.e. equidistant in freq
2009 : {
2010 : // Define via width of first channel to avoid numerical problems
2011 :
2012 : // channel width not set
2013 0 : if (regridChanWidthF <= 0.)
2014 : {
2015 0 : theRegridBWF = transCHAN_WIDTH[0] *
2016 0 : floor((theRegridBWF + transCHAN_WIDTH[0] * 0.01)/ transCHAN_WIDTH[0]);
2017 : }
2018 : else
2019 : {
2020 0 : theRegridBWF = regridChanWidthF *
2021 0 : floor((theRegridBWF + regridChanWidthF * 0.01)/ regridChanWidthF);
2022 : }
2023 : }
2024 : }
2025 : // Channel width not set
2026 0 : else if (regridChanWidthF <= 0.)
2027 : {
2028 0 : theRegridBWF = transCHAN_WIDTH[0] * nchan;
2029 : }
2030 : else
2031 : {
2032 0 : theRegridBWF = regridChanWidthF * nchan;
2033 : }
2034 :
2035 : // Center was not set by user but calculated
2036 0 : if (regridCenterF <= 0. || regridCenter < -C::c)
2037 : {
2038 : // Need to update
2039 0 : theRegridCenterF = transNewXin[0] - transCHAN_WIDTH[0] / 2. + theRegridBWF / 2.;
2040 0 : centerIsStart = false;
2041 : }
2042 : // Center but not nchan was set by user
2043 0 : else if (nchan < 0)
2044 : {
2045 : // Verify that the bandwidth is correct
2046 0 : if (centerIsStart)
2047 : {
2048 0 : if (startIsEnd)
2049 : {
2050 0 : theRegridBWF = theRegridCenterF - transNewXin[0] + transCHAN_WIDTH[0] / 2.;
2051 : }
2052 : // start is start
2053 : else
2054 : {
2055 0 : theRegridBWF = transNewXin[oldNUM_CHAN - 1] +
2056 0 : transCHAN_WIDTH[oldNUM_CHAN - 1] / 2. -
2057 : theRegridCenterF;
2058 : }
2059 0 : if (regridQuant == "freq" || regridQuant == "vrad") // i.e. equidistant in freq
2060 : {
2061 : // Define via width of first channel to avoid numerical problems
2062 0 : if (regridChanWidthF <= 0.) { // channel width not set
2063 0 : theRegridBWF = transCHAN_WIDTH[0]
2064 0 : * floor((theRegridBWF + transCHAN_WIDTH[0]* 0.01) / transCHAN_WIDTH[0]);
2065 : }
2066 : else
2067 : {
2068 0 : theRegridBWF = regridChanWidthF
2069 0 : * floor((theRegridBWF+ regridChanWidthF* 0.01)/ regridChanWidthF);
2070 : }
2071 : }
2072 : }
2073 : // center is center
2074 : else {
2075 0 : theRegridBWF = 2. * min((Double) (theRegridCenterF - transNewXin[0]- transCHAN_WIDTH[0]),
2076 0 : (Double) (transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] - theRegridCenterF));
2077 : }
2078 : }
2079 : }
2080 :
2081 : // Now can convert start to center
2082 0 : if (centerIsStart)
2083 : {
2084 0 : if (startIsEnd)
2085 : {
2086 0 : theRegridCenterF = theRegridCenterF - theRegridBWF / 2.;
2087 : }
2088 : else
2089 : {
2090 0 : theRegridCenterF = theRegridCenterF + theRegridBWF / 2.;
2091 : }
2092 0 : centerIsStart = false;
2093 : }
2094 : }
2095 : // regridBandwidthF was set
2096 : else {
2097 : // Determine actually possible bandwidth
2098 : // width will be truncated to the maximum width possible
2099 : // symmetrically around the value given by "regrid_center"
2100 0 : theRegridBWF = regridBandwidthF;
2101 :
2102 : // Now can convert start to center
2103 0 : if (centerIsStart)
2104 : {
2105 0 : if (startIsEnd)
2106 : {
2107 0 : theRegridCenterF = theRegridCenterF - theRegridBWF / 2.;
2108 : }
2109 : else
2110 : {
2111 0 : theRegridCenterF = theRegridCenterF + theRegridBWF / 2.;
2112 : }
2113 0 : centerIsStart = false;
2114 : }
2115 :
2116 0 : Double rangeTol = 1.; // Hz
2117 0 : if ((regridQuant == "vopt" || regridQuant == "wave")) // i.e. if the center is the center w.r.t. wavelength
2118 : {
2119 0 : rangeTol = transCHAN_WIDTH[0];
2120 : }
2121 :
2122 0 : if ((theRegridCenterF + theRegridBWF / 2.) - (transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.) > rangeTol)
2123 : {
2124 : oss << " *** Input spectral window exceeds upper end of original window. "
2125 0 : "Adjusting to max. possible value." << endl;
2126 0 : theRegridBWF = min( (Double) fabs(transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.- theRegridCenterF),
2127 0 : (Double) fabs(theRegridCenterF - transNewXin[0]+ transCHAN_WIDTH[0] / 2.)) * 2.;
2128 :
2129 0 : if (theRegridBWF < transCHAN_WIDTH[0])
2130 : {
2131 0 : theRegridCenterF = ( transNewXin[0]
2132 0 : + transCHAN_WIDTH[oldNUM_CHAN - 1]
2133 0 : + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.
2134 0 : - transCHAN_WIDTH[0] / 2.) / 2.;
2135 0 : theRegridBWF = transCHAN_WIDTH[oldNUM_CHAN - 1]
2136 0 : - transNewXin[0] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2. + transCHAN_WIDTH[0] / 2.;
2137 : }
2138 : }
2139 0 : if ((theRegridCenterF - theRegridBWF / 2.) - (transNewXin[0] - transCHAN_WIDTH[0] / 2.) < -rangeTol)
2140 : {
2141 : oss << " *** Input spectral window exceeds lower end of original window. "
2142 0 : " Adjusting to max. possible value."<< endl;
2143 :
2144 0 : theRegridBWF = min( (Double) fabs(transNewXin[oldNUM_CHAN - 1] + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2. - theRegridCenterF),
2145 0 : (Double) fabs(theRegridCenterF - transNewXin[0] + transCHAN_WIDTH[0] / 2.)) * 2.;
2146 :
2147 0 : if (theRegridBWF < transCHAN_WIDTH[0])
2148 : {
2149 0 : theRegridCenterF = (transNewXin[0]
2150 0 : + transCHAN_WIDTH[oldNUM_CHAN - 1]
2151 0 : + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.
2152 0 : - transCHAN_WIDTH[0] / 2.) / 2.;
2153 0 : theRegridBWF = transCHAN_WIDTH[oldNUM_CHAN - 1]
2154 0 : - transNewXin[0]
2155 0 : + transCHAN_WIDTH[oldNUM_CHAN - 1] / 2.
2156 0 : + transCHAN_WIDTH[0] / 2.;
2157 : }
2158 : }
2159 : }
2160 :
2161 0 : if (regridChanWidthF <= 0.) // "not set"
2162 : {
2163 0 : if (nchan != 0 || centerIsStartC) // use first channel
2164 : {
2165 0 : theCentralChanWidthF = transCHAN_WIDTH[0];
2166 : }
2167 : else
2168 : {
2169 : // keep channel width similar to the old one
2170 0 : theCentralChanWidthF = transCHAN_WIDTH[oldNUM_CHAN / 2]; // use channel width from
2171 : // near central channel
2172 : }
2173 : }
2174 : // regridChanWidthF was set
2175 : else
2176 : {
2177 : // Keep in limits
2178 0 : theCentralChanWidthF = regridChanWidthF;
2179 :
2180 : // Too large => make a single channel
2181 0 : if (theCentralChanWidthF > theRegridBWF)
2182 : {
2183 0 : theCentralChanWidthF = theRegridBWF;
2184 : oss
2185 0 : << " *** Requested new channel width exceeds defined SPW width."
2186 0 : << endl
2187 0 : << " Creating a single channel with the defined SPW width."
2188 0 : << endl;
2189 : }
2190 : // Check if too small
2191 0 : else if (theCentralChanWidthF < transCHAN_WIDTH[0])
2192 : {
2193 : // Determine smallest channel width
2194 0 : lDouble smallestChanWidth = 1E30;
2195 0 : Int ii = 0;
2196 0 : for (Int i = 0; i < oldNUM_CHAN; i++)
2197 : {
2198 0 : if (transCHAN_WIDTH[i] < smallestChanWidth)
2199 : {
2200 0 : smallestChanWidth = transCHAN_WIDTH[i];
2201 0 : ii = i;
2202 : }
2203 : }
2204 :
2205 : // 1 Hz tolerance to cope with numerical accuracy problems
2206 0 : if (theCentralChanWidthF < smallestChanWidth - 1.)
2207 : {
2208 0 : oss << " *** Requested new channel width (" << theCentralChanWidthF << " Hz) "
2209 0 : << "is smaller than smallest original channel width" << endl;
2210 0 : oss << " which is " << smallestChanWidth << " Hz" << endl;
2211 :
2212 0 : if (regridQuant == "vrad")
2213 : {
2214 0 : oss << " or "
2215 0 : << (vrad(transNewXin[ii],regridVeloRestfrq)
2216 0 : - vrad(transNewXin[ii] + transCHAN_WIDTH[ii] / 2.,regridVeloRestfrq)) * 2. << " m/s";
2217 : }
2218 :
2219 0 : if (regridQuant == "vopt")
2220 : {
2221 0 : oss << " or "
2222 0 : << (vopt(transNewXin[ii],regridVeloRestfrq)
2223 0 : - vopt(transNewXin[ii] + transCHAN_WIDTH[ii] / 2.,regridVeloRestfrq)) * 2. << " m/s";
2224 : }
2225 :
2226 0 : message = oss.str();
2227 0 : return false;
2228 :
2229 : }
2230 : // input channel width was OK, memorize
2231 : else
2232 : {
2233 0 : theChanWidthX = regridChanWidth;
2234 : }
2235 : }
2236 : }
2237 :
2238 0 : oss << " Channels equidistant in " << regridQuant << endl
2239 0 : << " Central frequency (in output frame) = " << theRegridCenterF << " Hz";
2240 :
2241 0 : if (regridQuant == "vrad")
2242 : {
2243 0 : oss << " == " << vrad(theRegridCenterF, regridVeloRestfrq) << " m/s radio velocity";
2244 : }
2245 0 : else if (regridQuant == "vopt")
2246 : {
2247 0 : oss << " == " << vopt(theRegridCenterF, regridVeloRestfrq) << " m/s optical velocity";
2248 : }
2249 0 : else if (regridQuant == "wave")
2250 : {
2251 0 : oss << " == " << lambda(theRegridCenterF) << " m wavelength";
2252 : }
2253 0 : oss << endl;
2254 :
2255 0 : if (isDescending)
2256 : {
2257 0 : oss << " Channel central frequency is decreasing with increasing channel number." << endl;
2258 : }
2259 :
2260 0 : oss << " Width of central channel (in output frame) = " << theCentralChanWidthF << " Hz";
2261 :
2262 0 : if (regridQuant == "vrad")
2263 : {
2264 0 : oss << " == " << vrad(theRegridCenterF - theCentralChanWidthF,regridVeloRestfrq)
2265 0 : - vrad(theRegridCenterF,regridVeloRestfrq) << " m/s radio velocity";
2266 : }
2267 0 : else if (regridQuant == "vopt")
2268 : {
2269 0 : oss << " == " << vopt(theRegridCenterF - theCentralChanWidthF,regridVeloRestfrq)
2270 0 : - vopt(theRegridCenterF,regridVeloRestfrq) << " m/s optical velocity";
2271 : }
2272 0 : else if (regridQuant == "wave")
2273 : {
2274 0 : oss << " == " << lambda(theRegridCenterF - theCentralChanWidthF)
2275 0 : - lambda(theRegridCenterF) << " m wavelength";
2276 : }
2277 0 : oss << endl;
2278 :
2279 : // Now calculate newChanLoBound, and newChanHiBound from theRegridCenterF, theRegridBWF, theCentralChanWidthF
2280 0 : vector<lDouble> loFBup; // The lower bounds for the new channels starting from the central channel going up
2281 0 : vector<lDouble> hiFBup; // The lower bounds for the new channels starting from the central channel going up
2282 0 : vector<lDouble> loFBdown; // The lower bounds for the new channels starting from the central channel going down
2283 0 : vector<lDouble> hiFBdown; // The lower bounds for the new channels starting from the central channel going down
2284 :
2285 0 : lDouble edgeTolerance = theCentralChanWidthF * 0.01; // Needed to avoid numerical accuracy problems
2286 :
2287 : // Re-gridding in radio velocity ...
2288 0 : if (regridQuant == "vrad")
2289 : {
2290 : // Create freq boundaries equidistant and contiguous in radio velocity
2291 0 : lDouble upperEndF = theRegridCenterF + theRegridBWF / 2.;
2292 0 : lDouble lowerEndF = theRegridCenterF - theRegridBWF / 2.;
2293 0 : lDouble upperEndV = vrad(upperEndF, regridVeloRestfrq);
2294 0 : lDouble lowerEndV = vrad(lowerEndF, regridVeloRestfrq);
2295 : lDouble velLo;
2296 : lDouble velHi;
2297 :
2298 : // Want to keep the center of the center channel at the center
2299 : // of the new center channel if the bandwidth is an odd multiple
2300 : // of the new channel width, otherwise the center channel is the
2301 : // lower edge of the new center channel
2302 0 : lDouble tnumChan = floor((theRegridBWF + edgeTolerance) / theCentralChanWidthF);
2303 :
2304 0 : if ((Int) tnumChan % 2 != 0)
2305 : {
2306 : // Odd multiple
2307 0 : loFBup.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2308 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2309 0 : loFBdown.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2310 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2311 : }
2312 : else
2313 : {
2314 0 : loFBup.push_back(theRegridCenterF);
2315 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF);
2316 0 : loFBdown.push_back(theRegridCenterF);
2317 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF);
2318 : }
2319 :
2320 : // Cannot use original channel width in velocity units
2321 0 : if (theChanWidthX < 0)
2322 : {
2323 : // Need to calculate back from central channel width in Hz
2324 0 : theChanWidthX = vrad(loFBup[0], regridVeloRestfrq) - vrad(hiFBup[0], regridVeloRestfrq);
2325 : }
2326 :
2327 : // calc velocity corresponding to the upper end (in freq) of the
2328 : // last added channel which is the lower end of the next channel
2329 0 : velLo = vrad(hiFBup[0], regridVeloRestfrq);
2330 :
2331 : // calc velocity corresponding to the upper end (in freq) of the next channel
2332 0 : velHi = velLo - theChanWidthX; // vrad goes down as freq goes up!
2333 :
2334 : // (Preventing accuracy problems)
2335 0 : while (upperEndV - theChanWidthX / 10. < velHi)
2336 : {
2337 : // calc frequency of the upper end (in freq) of the next channel
2338 0 : lDouble freqHi = freq_from_vrad(velHi, regridVeloRestfrq);
2339 :
2340 : // End of bandwidth not yet reached
2341 0 : if (freqHi <= upperEndF + edgeTolerance)
2342 : {
2343 0 : loFBup.push_back(hiFBup.back());
2344 0 : hiFBup.push_back(freqHi);
2345 : }
2346 0 : else if (freqHi < upperEndF + edgeTolerance)
2347 : {
2348 0 : loFBup.push_back(hiFBup.back());
2349 0 : hiFBup.push_back(upperEndF);
2350 0 : break;
2351 : }
2352 : else
2353 : {
2354 0 : break;
2355 : }
2356 :
2357 : // calc velocity corresponding to the upper end (in freq) of the added channel
2358 0 : velLo = vrad(hiFBup.back(), regridVeloRestfrq);
2359 :
2360 : // calc velocity corresponding to the upper end (in freq) of the next channel
2361 0 : velHi = velLo - theChanWidthX; // vrad goes down as freq goes up
2362 : }
2363 :
2364 : // calc velocity corresponding to the lower end (in freq) of the
2365 : // Last added channel which is the upper end of the next channel
2366 0 : velHi = vrad(loFBdown[0], regridVeloRestfrq);
2367 :
2368 : // calc velocity corresponding to the lower end (in freq) of the next channel
2369 0 : velLo = velHi + theChanWidthX; // vrad goes up as freq goes down!
2370 :
2371 : // (Preventing accuracy problems)
2372 0 : while (velLo < lowerEndV + theChanWidthX / 10.)
2373 : {
2374 : // calc frequency of the lower end (in freq) of the next channel
2375 0 : lDouble freqLo = freq_from_vrad(velLo, regridVeloRestfrq);
2376 :
2377 : // End of bandwidth not yet reached
2378 0 : if (freqLo >= lowerEndF - edgeTolerance)
2379 : {
2380 0 : hiFBdown.push_back(loFBdown.back());
2381 0 : loFBdown.push_back(freqLo);
2382 : }
2383 0 : else if (freqLo > lowerEndF - edgeTolerance)
2384 : {
2385 0 : hiFBdown.push_back(loFBdown.back());
2386 0 : loFBdown.push_back(lowerEndF);
2387 0 : break;
2388 : }
2389 : else
2390 : {
2391 0 : break;
2392 : }
2393 :
2394 : // calc velocity corresponding to the upper end of the next channel
2395 0 : velHi = vrad(loFBdown.back(), regridVeloRestfrq);
2396 :
2397 : // calc velocity corresponding to the lower end (in freq) of the next channel
2398 0 : velLo = velHi + theChanWidthX; // vrad goes up as freq goes down
2399 : }
2400 : }
2401 : // Regridding in optical velocity ...
2402 0 : else if (regridQuant == "vopt")
2403 : {
2404 : // Create freq boundaries equidistant and contiguous in optical velocity
2405 0 : lDouble upperEndF = theRegridCenterF + theRegridBWF / 2.;
2406 0 : lDouble lowerEndF = theRegridCenterF - theRegridBWF / 2.;
2407 0 : lDouble upperEndV = vopt(upperEndF, regridVeloRestfrq);
2408 0 : lDouble lowerEndV = vopt(lowerEndF, regridVeloRestfrq);
2409 : lDouble velLo;
2410 : lDouble velHi;
2411 :
2412 : // Want to keep the center of the center channel at the center
2413 : // of the new center channel if the bandwidth is an odd multiple
2414 : // of the new channel width, otherwise the center channel is the
2415 : // lower edge of the new center channel
2416 :
2417 : // Enlarged edge tolerance since channels non-equidistant in freq
2418 0 : lDouble tnumChan = floor((theRegridBWF + edgeTolerance) / theCentralChanWidthF);
2419 :
2420 : // Odd multiple
2421 0 : if ((Int) tnumChan % 2 != 0)
2422 : {
2423 :
2424 0 : loFBup.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2425 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2426 0 : loFBdown.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2427 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2428 : }
2429 : else
2430 : {
2431 0 : loFBup.push_back(theRegridCenterF);
2432 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF);
2433 0 : loFBdown.push_back(theRegridCenterF);
2434 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF);
2435 : }
2436 :
2437 : // Cannot use original channel width in velocity units
2438 0 : if (theChanWidthX < 0)
2439 : {
2440 : // Need to calculate back from central channel width in Hz
2441 0 : theChanWidthX = vopt(loFBup[0], regridVeloRestfrq) - vopt(hiFBup[0], regridVeloRestfrq);
2442 : }
2443 :
2444 : // calc velocity corresponding to the upper end (in freq) of the
2445 : // last added channel which is the lower end of the next channel
2446 0 : velLo = vopt(hiFBup[0], regridVeloRestfrq);
2447 :
2448 : // calc velocity corresponding to the upper end (in freq) of the next channel
2449 0 : velHi = velLo - theChanWidthX; // vopt goes down as freq goes up!
2450 :
2451 : // (Preventing accuracy problems)
2452 0 : while (upperEndV - velHi < theChanWidthX / 10.)
2453 : {
2454 : // calc frequency of the upper end (in freq) of the next channel
2455 0 : lDouble freqHi = freq_from_vopt(velHi, regridVeloRestfrq);
2456 :
2457 : // End of bandwidth not yet reached
2458 0 : if (freqHi <= upperEndF + edgeTolerance)
2459 : {
2460 0 : loFBup.push_back(hiFBup.back());
2461 0 : hiFBup.push_back(freqHi);
2462 : }
2463 0 : else if (freqHi < upperEndF + edgeTolerance)
2464 : {
2465 0 : loFBup.push_back(hiFBup.back());
2466 0 : hiFBup.push_back(upperEndF);
2467 0 : break;
2468 : }
2469 : else
2470 : {
2471 0 : break;
2472 : }
2473 :
2474 : // calc velocity corresponding to the upper end (in freq) of the added channel
2475 0 : velLo = vopt(hiFBup.back(), regridVeloRestfrq);
2476 : // calc velocity corresponding to the upper end (in freq) of the next channel
2477 0 : velHi = velLo - theChanWidthX; // vopt goes down as freq goes up
2478 : }
2479 :
2480 : // calc velocity corresponding to the lower end (in freq) of the
2481 : // last added channel which is the upper end of the next channel
2482 0 : velHi = vopt(loFBdown[0], regridVeloRestfrq);
2483 :
2484 : // calc velocity corresponding to the lower end (in freq) of the next channel
2485 0 : velLo = velHi + theChanWidthX; // vopt goes up as freq goes down!
2486 :
2487 : // (Preventing accuracy problems)
2488 0 : while (velLo - lowerEndV < theChanWidthX / 10.)
2489 : {
2490 : // calc frequency of the lower end (in freq) of the next channel
2491 0 : lDouble freqLo = freq_from_vopt(velLo, regridVeloRestfrq);
2492 :
2493 : // End of bandwidth not yet reached
2494 0 : if (freqLo >= lowerEndF - edgeTolerance)
2495 : {
2496 0 : hiFBdown.push_back(loFBdown.back());
2497 0 : loFBdown.push_back(freqLo);
2498 : }
2499 0 : else if (freqLo > lowerEndF - edgeTolerance)
2500 : {
2501 0 : hiFBdown.push_back(loFBdown.back());
2502 0 : loFBdown.push_back(lowerEndF);
2503 0 : break;
2504 : }
2505 : else
2506 : {
2507 0 : break;
2508 : }
2509 :
2510 : // calc velocity corresponding to the upper end of the next channel
2511 0 : velHi = vopt(loFBdown.back(), regridVeloRestfrq);
2512 : // calc velocity corresponding to the lower end (in freq) of the next channel
2513 0 : velLo = velHi + theChanWidthX; // vopt goes up as freq goes down
2514 : }
2515 : }
2516 : // Re-gridding in frequency ...
2517 0 : else if (regridQuant == "freq")
2518 : {
2519 : // Create freq boundaries equidistant and contiguous in frequency
2520 0 : lDouble upperEndF = theRegridCenterF + theRegridBWF / 2.;
2521 0 : lDouble lowerEndF = theRegridCenterF - theRegridBWF / 2.;
2522 :
2523 : // Want to keep the center of the center channel at the center
2524 : // of the new center channel if the bandwidth is an odd multiple
2525 : // of the new channel width, otherwise the center channel is the
2526 : // lower edge of the new center channel
2527 0 : lDouble tnumChan = floor((theRegridBWF + edgeTolerance) / theCentralChanWidthF);
2528 :
2529 : // Odd multiple
2530 0 : if ((Int) tnumChan % 2 != 0)
2531 : {
2532 0 : loFBup.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2533 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2534 0 : loFBdown.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2535 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2536 : }
2537 : else
2538 : {
2539 0 : loFBup.push_back(theRegridCenterF);
2540 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF);
2541 0 : loFBdown.push_back(theRegridCenterF);
2542 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF);
2543 : }
2544 :
2545 0 : while (hiFBup.back() < upperEndF + edgeTolerance)
2546 : {
2547 : // calc frequency of the upper end of the next channel
2548 0 : lDouble freqHi = hiFBup.back() + theCentralChanWidthF;
2549 :
2550 : // End of bandwidth not yet reached
2551 0 : if (freqHi <= upperEndF + edgeTolerance)
2552 : {
2553 0 : loFBup.push_back(hiFBup.back());
2554 0 : hiFBup.push_back(freqHi);
2555 : }
2556 : else
2557 : {
2558 0 : break;
2559 : }
2560 : }
2561 :
2562 0 : while (loFBdown.back() > lowerEndF - edgeTolerance)
2563 : {
2564 : // calc frequency of the lower end of the next channel
2565 0 : lDouble freqLo = loFBdown.back() - theCentralChanWidthF;
2566 :
2567 : // End of bandwidth not yet reached
2568 0 : if (freqLo >= lowerEndF - edgeTolerance)
2569 : {
2570 0 : hiFBdown.push_back(loFBdown.back());
2571 0 : loFBdown.push_back(freqLo);
2572 : }
2573 : else
2574 : {
2575 0 : break;
2576 : }
2577 : }
2578 : }
2579 : // Re-gridding in wavelength ...
2580 0 : else if (regridQuant == "wave")
2581 : {
2582 : // Create freq boundaries equidistant and contiguous in wavelength
2583 0 : lDouble upperEndF = theRegridCenterF + theRegridBWF / 2.;
2584 0 : lDouble lowerEndF = theRegridCenterF - theRegridBWF / 2.;
2585 0 : lDouble upperEndL = lambda(upperEndF);
2586 0 : lDouble lowerEndL = lambda(lowerEndF);
2587 : lDouble lambdaLo;
2588 : lDouble lambdaHi;
2589 :
2590 : // Want to keep the center of the center channel at the center
2591 : // of the new center channel if the bandwidth is an odd multiple
2592 : // of the new channel width, otherwise the center channel is the
2593 : // lower edge of the new center channel
2594 0 : lDouble tnumChan = floor((theRegridBWF + edgeTolerance) / theCentralChanWidthF);
2595 :
2596 : // Odd multiple
2597 0 : if ((Int) tnumChan % 2 != 0)
2598 : {
2599 0 : loFBup.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2600 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2601 0 : loFBdown.push_back(theRegridCenterF - theCentralChanWidthF / 2.);
2602 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF / 2.);
2603 : }
2604 : else
2605 : {
2606 0 : loFBup.push_back(theRegridCenterF);
2607 0 : hiFBup.push_back(theRegridCenterF + theCentralChanWidthF);
2608 0 : loFBdown.push_back(theRegridCenterF);
2609 0 : hiFBdown.push_back(theRegridCenterF + theCentralChanWidthF);
2610 : }
2611 :
2612 : // Cannot use original channel width in wavelength units
2613 0 : if (theChanWidthX < 0)
2614 : {
2615 : // Need to calculate back from central channel width in Hz
2616 0 : theChanWidthX = lambda(loFBup[0]) - lambda(hiFBup[0]);
2617 : }
2618 :
2619 : // calc wavelength corresponding to the upper end (in freq) of the
2620 : // last added channel which is the lower end of the next channel
2621 0 : lambdaLo = lambda(hiFBup[0]);
2622 :
2623 : // calc wavelength corresponding to the upper end (in freq) of the next channel
2624 0 : lambdaHi = lambdaLo - theChanWidthX; // lambda goes down as freq goes up!
2625 :
2626 : // (Preventing accuracy problems)
2627 0 : while (upperEndL - lambdaHi < theChanWidthX / 10.)
2628 : {
2629 : // calc frequency of the upper end (in freq) of the next channel
2630 0 : lDouble freqHi = freq_from_lambda(lambdaHi);
2631 :
2632 : // End of bandwidth not yet reached
2633 0 : if (freqHi <= upperEndF + edgeTolerance)
2634 : {
2635 0 : loFBup.push_back(hiFBup.back());
2636 0 : hiFBup.push_back(freqHi);
2637 : }
2638 0 : else if (freqHi < upperEndF + edgeTolerance)
2639 : {
2640 0 : loFBup.push_back(hiFBup.back());
2641 0 : hiFBup.push_back(upperEndF);
2642 0 : break;
2643 : }
2644 : else
2645 : {
2646 0 : break;
2647 : }
2648 :
2649 : // calc wavelength corresponding to the upper end (in freq) of the added channel
2650 0 : lambdaLo = lambda(hiFBup.back());
2651 : // calc wavelength corresponding to the upper end (in freq) of the next channel
2652 0 : lambdaHi = lambdaLo - theChanWidthX; // lambda goes down as freq goes up
2653 : }
2654 :
2655 : // calc wavelength corresponding to the lower end (in freq) of the
2656 : // last added channel which is the upper end of the next channel
2657 0 : lambdaHi = lambda(loFBdown[0]);
2658 :
2659 : // calc wavelength corresponding to the lower end (in freq) of the next channel
2660 0 : lambdaLo = lambdaHi + theChanWidthX; // lambda goes up as freq goes down!
2661 :
2662 :
2663 : // (Preventing accuracy problems)
2664 0 : while (lambdaLo - lowerEndL < theChanWidthX / 10.)
2665 : {
2666 : // calc frequency of the lower end (in freq) of the next channel
2667 0 : lDouble freqLo = freq_from_lambda(lambdaLo);
2668 :
2669 : // End of bandwidth not yet reached
2670 0 : if (freqLo >= lowerEndF - edgeTolerance)
2671 : {
2672 0 : hiFBdown.push_back(loFBdown.back());
2673 0 : loFBdown.push_back(freqLo);
2674 : }
2675 0 : else if (freqLo > lowerEndF - edgeTolerance)
2676 : {
2677 0 : hiFBdown.push_back(loFBdown.back());
2678 0 : loFBdown.push_back(lowerEndF);
2679 0 : break;
2680 : }
2681 : else
2682 : {
2683 0 : break;
2684 : }
2685 :
2686 : // calc wavelength corresponding to the upper end of the next channel
2687 0 : lambdaHi = lambda(loFBdown.back());
2688 : // calc wavelength corresponding to the lower end (in freq) of the next channel
2689 0 : lambdaLo = lambdaHi + theChanWidthX; // wavelength goes up as freq goes down
2690 : }
2691 :
2692 : }
2693 : // should not get here
2694 : else
2695 : {
2696 0 : oss << "Invalid value " << regridQuant << " for parameter \"mode\".";
2697 0 : message = oss.str();
2698 0 : return false;
2699 : }
2700 :
2701 0 : Int numNewChanDown = loFBdown.size();
2702 0 : Int numNewChanUp = loFBup.size();
2703 :
2704 : // central channel contained in both vectors
2705 0 : newChanLoBound.resize(numNewChanDown + numNewChanUp - 1);
2706 :
2707 0 : newChanHiBound.resize(numNewChanDown + numNewChanUp - 1);
2708 0 : for (Int i = 0; i < numNewChanDown; i++)
2709 : {
2710 0 : Int k = numNewChanDown - i - 1; // Need to assign in reverse
2711 0 : newChanLoBound[i] = loFBdown[k];
2712 0 : newChanHiBound[i] = hiFBdown[k];
2713 : }
2714 0 : for (Int i = 1; i < numNewChanUp; i++) // Start at 1 to omit the central channel here
2715 : {
2716 0 : newChanLoBound[i + numNewChanDown - 1] = loFBup[i];
2717 0 : newChanHiBound[i + numNewChanDown - 1] = hiFBup[i];
2718 : }
2719 :
2720 0 : uInt nc = newChanLoBound.size();
2721 0 : oss << " Number of channels = " << nc << endl;
2722 0 : oss << " Total width of SPW (in output frame) = " << newChanHiBound[nc - 1] - newChanLoBound[0] << " Hz" << endl;
2723 0 : oss << " Lower edge = " << newChanLoBound[0] << " Hz,"
2724 0 : << " upper edge = " << newChanHiBound[nc - 1] << " Hz" << endl;
2725 :
2726 : // Original SPW was in reverse order, need to restore that
2727 0 : if (isDescending)
2728 : {
2729 0 : Vector<Double> tempL, tempU;
2730 0 : tempL.assign(newChanLoBound);
2731 0 : tempU.assign(newChanHiBound);
2732 0 : for (uInt i = 0; i < nc; i++)
2733 : {
2734 0 : newChanLoBound(i) = tempL(nc - 1 - i);
2735 0 : newChanHiBound(i) = tempU(nc - 1 - i);
2736 : }
2737 : }
2738 :
2739 0 : message = oss.str();
2740 :
2741 0 : return true;
2742 :
2743 : } // end if (regridQuant== ...
2744 : }
2745 :
2746 : } //# NAMESPACE CASA - END
|
