Line data Source code
1 : //# ATCAFiller.cc: ATCA filler - reads rpfits, creates MeasurementSet
2 : //# Copyright (C) 2004
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 :
28 : #include <atnf/atca/ATCAFiller.h>
29 : #include <atnf/atca/ATAtmosphere.h>
30 : #include <casacore/casa/Arrays/Cube.h>
31 : #include <casacore/casa/Utilities/GenSort.h>
32 : #include <casacore/scimath/Mathematics/FFTServer.h>
33 : #include <casacore/casa/OS/DirectoryIterator.h>
34 : #include <casacore/casa/OS/RegularFile.h>
35 : #include <casacore/tables/Tables.h>
36 : #include <atnf/rpfits/RPFITS.h>
37 : #include <casacore/ms/MeasurementSets/MSTileLayout.h>
38 :
39 : using namespace casacore;
40 : namespace casa {
41 :
42 0 : Int myround(Double x) { return Int(floor(x+0.5));}
43 :
44 1 : ATCAFiller::ATCAFiller():
45 : appendMode_p(false),
46 : storedHeader_p(false),
47 : skipScan_p(false),
48 : skipData_p(false),
49 : firstHeader_p(false),
50 : listHeader_p(false),
51 : fileSize_p(0),
52 : birdie_p(false),
53 : reweight_p(false),
54 : noxycorr_p(false),
55 : obsType_p(0),
56 : hires_p(false),
57 : init_p(false),
58 : lastUT_p(0),
59 : bandWidth1_p(0),
60 : numChan1_p(0),
61 : shadow_p(0),
62 : autoFlag_p(true),
63 : flagScanType_p(false),
64 1 : flagCount_p(NFLAG,0)
65 1 : {}
66 :
67 1 : ATCAFiller::~ATCAFiller()
68 1 : {}
69 :
70 1 : Bool ATCAFiller::open(const String& msName, const Vector<String>& rpfitsFiles,
71 : const Vector<String>& options, Int opcor)
72 :
73 : {
74 : Int opt;
75 :
76 3 : LogOrigin orig("ATCAFiller", "open()", WHERE);
77 1 : os_p = LogIO(orig);
78 1 : rpfitsFiles_p = Directory::shellExpand(rpfitsFiles, false);
79 1 : GenSort<String>::sort(rpfitsFiles_p);
80 1 : if (rpfitsFiles_p.nelements() > 0) {
81 1 : os_p << LogIO::NORMAL << "Expanded file names are : " << endl;
82 2 : for (uInt i=0; i<rpfitsFiles_p.nelements(); i++) {
83 1 : if (rpfitsFiles_p(i).empty()) {
84 0 : os_p << "Input file number " << i+1 << " is empty" << LogIO::EXCEPTION;
85 : }
86 1 : os_p << rpfitsFiles_p(i) << endl;
87 : }
88 : } else {
89 0 : os_p << "Input file names vector is empty" << LogIO::EXCEPTION;
90 : }
91 2 : Vector<String> opts;
92 1 : if (options.nelements()==1) {
93 1 : opts=stringToVector(options(0),std::regex(" |,"));
94 : } else {
95 0 : opts=options;
96 : }
97 1 : opcor_p=opcor;
98 :
99 1 : Bool compress=false;
100 : // cerr<<"options="<<opts<<endl;
101 2 : for (opt=0; opt<Int(opts.nelements()); opt++) {
102 1 : if (downcase(opts(opt)) == "birdie") {
103 0 : birdie_p = true;
104 : }
105 1 : if (downcase(opts(opt)) == "reweight") {
106 0 : reweight_p = true;
107 : }
108 1 : if (downcase(opts(opt)) == "noxycorr") {
109 0 : noxycorr_p = true;
110 : }
111 1 : if (downcase(opts(opt)) == "compress") {
112 0 : compress = true;
113 : }
114 1 : if (downcase(opts(opt)) == "fastmosaic") {
115 0 : obsType_p=1;
116 : }
117 1 : if (downcase(opts(opt)) == "noautoflag") {
118 0 : autoFlag_p=false;
119 : }
120 1 : if (downcase(opts(opt)) == "hires") {
121 0 : hires_p=true;
122 : }
123 1 : if (downcase(opts(opt)) == "noac") {
124 1 : noac_p = true;
125 : }
126 : }
127 : // Check if this is CABB data or old ATCA data by checking the
128 : // INSTRUMENT keyword (ATCA, ATCABB)
129 1 : cabb_p = checkCABB(rpfitsFiles_p(0));
130 1 : reweight_p = reweight_p && !cabb_p;
131 1 : atms_p = makeTable(msName,compress,cabb_p);
132 1 : msc_p = new MSColumns(atms_p);
133 1 : msScanInfo_p = atms_p.keywordSet().asTable("ATCA_SCAN_INFO");
134 :
135 1 : String comp;
136 1 : if (compress) comp="Comp";
137 1 : dataAccessor_p = TiledDataStManAccessor(atms_p,"TiledData"+comp);
138 1 : sigmaAccessor_p = TiledDataStManAccessor(atms_p,"TiledSigma");
139 1 : flagAccessor_p = TiledDataStManAccessor(atms_p,"TiledFlag");
140 1 : flagCatAccessor_p = TiledDataStManAccessor(atms_p,"TiledFlagCategory");
141 :
142 1 : prev_fieldId_p = -1;
143 1 : lastWeatherUT_p=0;
144 1 : errCount_p=0;
145 1 : init();
146 1 : init_p=true;
147 2 : return true;
148 : }
149 :
150 :
151 1 : TableDesc ATCAFiller::atcaTableDesc(Bool compress)
152 : {
153 1 : TableDesc td = MS::requiredTableDesc();
154 : // add the data column
155 1 : MS::addColumnToDesc(td,MS::DATA,2);
156 : // and its unit
157 1 : td.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().define("UNIT","Jy");
158 1 : if (compress) {
159 0 : String col = MS::columnName(MS::DATA);
160 0 : td.addColumn(ArrayColumnDesc<Int>(col+"_COMPRESSED",
161 0 : "observed data compressed",2));
162 0 : td.addColumn(ScalarColumnDesc<Float>(col+"_SCALE"));
163 0 : td.addColumn(ScalarColumnDesc<Float>(col+"_OFFSET"));
164 : }
165 :
166 : // add ATCA specific columns here..
167 1 : MS::addColumnToDesc(td,MS::PULSAR_BIN);
168 :
169 2 : td.addColumn(ScalarColumnDesc<Int>("ATCA_SYSCAL_ID_ANT1",
170 1 : "Index into SysCal table for Antenna1"));
171 2 : td.addColumn(ScalarColumnDesc<Int>("ATCA_SYSCAL_ID_ANT2",
172 1 : "Index into SysCal table for Antenna2"));
173 : // the tiled column indices
174 2 : td.addColumn(ScalarColumnDesc<Int>("DATA_HYPERCUBE_ID",
175 1 : "Index for Data Tiling"));
176 2 : td.addColumn(ScalarColumnDesc<Int>("SIGMA_HYPERCUBE_ID",
177 1 : "Index for Sigma Tiling"));
178 2 : td.addColumn(ScalarColumnDesc<Int>("FLAG_HYPERCUBE_ID",
179 1 : "Index for Flag Tiling"));
180 2 : td.addColumn(ScalarColumnDesc<Int>("FLAG_CATEGORY_HYPERCUBE_ID",
181 1 : "Index for Flag Category Tiling"));
182 :
183 1 : return td;
184 : }
185 :
186 1 : MeasurementSet ATCAFiller::makeTable(const String& tableName, Bool compress,
187 : Bool cabb)
188 : {
189 : // make the MeasurementSet Table
190 2 : TableDesc atDesc = atcaTableDesc(compress);
191 :
192 2 : String colData = MS::columnName(MS::DATA);
193 2 : String comp1,comp2;
194 1 : if (compress) { comp1="Comp"; comp2="_COMPRESSED";}
195 :
196 : // define tiled hypercube for the data
197 2 : Vector<String> coordColNames(0); //# don't use coord columns
198 1 : atDesc.defineHypercolumn("TiledData"+comp1,3,
199 2 : stringToVector(colData+comp2),
200 : coordColNames,
201 2 : stringToVector("DATA_HYPERCUBE_ID"));
202 1 : atDesc.defineHypercolumn("TiledSigma",2,
203 2 : stringToVector(MS::columnName(MS::SIGMA)),
204 : coordColNames,
205 2 : stringToVector("SIGMA_HYPERCUBE_ID"));
206 1 : atDesc.defineHypercolumn("TiledFlag",3,
207 2 : stringToVector(MS::columnName(MS::FLAG)),
208 : coordColNames,
209 2 : stringToVector("FLAG_HYPERCUBE_ID"));
210 1 : atDesc.defineHypercolumn("TiledFlagCategory",4,
211 2 : stringToVector(MS::columnName(MS::FLAG_CATEGORY)),
212 : coordColNames,
213 2 : stringToVector("FLAG_CATEGORY_HYPERCUBE_ID"));
214 :
215 2 : SetupNewTable newtab(tableName, atDesc, Table::New);
216 :
217 : // Set the default Storage Manager to be the incremental one
218 2 : IncrementalStMan incrStMan ("IncrementalData");
219 1 : newtab.bindAll(incrStMan, true);
220 : // Make an exception for fast varying data
221 2 : StandardStMan stStMan ("StandardData");
222 1 : newtab.bindColumn(MS::columnName(MS::UVW),stStMan);
223 1 : newtab.bindColumn(MS::columnName(MS::ANTENNA2),stStMan);
224 1 : newtab.bindColumn("ATCA_SYSCAL_ID_ANT2",stStMan);
225 1 : if (compress) {
226 0 : newtab.bindColumn(colData+"_SCALE",stStMan);
227 0 : newtab.bindColumn(colData+"_OFFSET",stStMan);
228 : }
229 :
230 :
231 2 : TiledDataStMan tiledStMan1("TiledData"+comp1);
232 : // Bind the DATA & FLAG column to the tiled stman
233 1 : newtab.bindColumn(MS::columnName(MS::DATA)+comp2,tiledStMan1);
234 1 : newtab.bindColumn("DATA_HYPERCUBE_ID",tiledStMan1);
235 1 : if (compress) {;
236 0 : CompressComplex ccData(colData,colData+"_COMPRESSED",
237 0 : colData+"_SCALE",colData+"_OFFSET", true);
238 0 : newtab.bindColumn(MS::columnName(MS::DATA),ccData);
239 : }
240 :
241 2 : TiledDataStMan tiledStMan2("TiledSigma");
242 : // Bind the SIGMA column to its own tiled stman
243 1 : newtab.bindColumn(MS::columnName(MS::SIGMA),tiledStMan2);
244 1 : newtab.bindColumn("SIGMA_HYPERCUBE_ID",tiledStMan2);
245 :
246 2 : TiledDataStMan tiledStMan3("TiledFlag");
247 : // Bind the FLAG column to its own tiled stman
248 1 : newtab.bindColumn(MS::columnName(MS::FLAG),tiledStMan3);
249 1 : newtab.bindColumn("FLAG_HYPERCUBE_ID",tiledStMan3);
250 :
251 2 : TiledDataStMan tiledStMan3c("TiledFlagCategory");
252 : // Bind the FLAG_CATEGORY column to its own tiled stman
253 1 : newtab.bindColumn(MS::columnName(MS::FLAG_CATEGORY),tiledStMan3c);
254 1 : newtab.bindColumn("FLAG_CATEGORY_HYPERCUBE_ID",tiledStMan3c);
255 :
256 :
257 : // create the table (with 0 rows)
258 1 : MeasurementSet ms(newtab, 0);
259 : // create the subtables
260 1 : makeSubTables(ms, Table::New, cabb);
261 :
262 2 : return ms;
263 : }
264 :
265 1 : void ATCAFiller::makeSubTables(MS& ms, Table::TableOption option, Bool cabb)
266 : {
267 : // This routine is modeled after MS::makeDummySubtables
268 : // we make new tables with 0 rows
269 1 : Int nrow=0;
270 :
271 : // Set up the subtables for the ATCA MS
272 :
273 : // Antenna
274 2 : TableDesc antTD=MSAntenna::requiredTableDesc();
275 1 : MSAntenna::addColumnToDesc(antTD, MSAntenna::PHASED_ARRAY_ID);
276 1 : MSAntenna::addColumnToDesc(antTD, MSAntenna::ORBIT_ID);
277 2 : SetupNewTable antennaSetup(ms.antennaTableName(),antTD,option);
278 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::ANTENNA),
279 2 : Table(antennaSetup,nrow));
280 :
281 : // Data descr
282 2 : TableDesc datadescTD=MSDataDescription::requiredTableDesc();
283 2 : SetupNewTable datadescSetup(ms.dataDescriptionTableName(),datadescTD,option);
284 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::DATA_DESCRIPTION),
285 2 : Table(datadescSetup,nrow));
286 :
287 : // Feed
288 2 : TableDesc feedTD=MSFeed::requiredTableDesc();
289 1 : MSFeed::addColumnToDesc(feedTD, MSFeed::PHASED_FEED_ID);
290 2 : SetupNewTable feedSetup(ms.feedTableName(),feedTD,option);
291 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FEED),
292 2 : Table(feedSetup,nrow));
293 :
294 : // Field
295 2 : TableDesc fieldTD=MSField::requiredTableDesc();
296 2 : SetupNewTable fieldSetup(ms.fieldTableName(),fieldTD,option);
297 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FIELD),
298 2 : Table(fieldSetup,nrow));
299 :
300 : // Flag_cmd
301 2 : TableDesc flagcmdTD=MSFlagCmd::requiredTableDesc();
302 2 : SetupNewTable flagcmdSetup(ms.flagCmdTableName(),flagcmdTD,option);
303 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FLAG_CMD),
304 2 : Table(flagcmdSetup,nrow));
305 :
306 :
307 : // Observation
308 2 : TableDesc obsTD=MSObservation::requiredTableDesc();
309 2 : SetupNewTable observationSetup(ms.observationTableName(),obsTD,option);
310 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::OBSERVATION),
311 2 : Table(observationSetup,nrow));
312 :
313 : // History
314 2 : TableDesc historyTD=MSHistory::requiredTableDesc();
315 2 : SetupNewTable historySetup(ms.historyTableName(),historyTD,option);
316 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::HISTORY),
317 2 : Table(historySetup,nrow));
318 :
319 : // Pointing
320 2 : TableDesc pointingTD=MSPointing::requiredTableDesc();
321 1 : MSPointing::addColumnToDesc(pointingTD, MSPointing::POINTING_OFFSET);
322 2 : SetupNewTable pointingSetup(ms.pointingTableName(),pointingTD,option);
323 : // Pointing table can be large, set some sensible defaults for storageMgrs
324 2 : IncrementalStMan ismPointing ("ISMPointing");
325 2 : StandardStMan ssmPointing("SSMPointing",32768);
326 1 : pointingSetup.bindAll(ismPointing,true);
327 1 : pointingSetup.bindColumn(MSPointing::columnName(MSPointing::ANTENNA_ID),
328 : ssmPointing);
329 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::POINTING),
330 2 : Table(pointingSetup,nrow));
331 :
332 : // Polarization
333 2 : TableDesc polTD=MSPolarization::requiredTableDesc();
334 2 : SetupNewTable polSetup(ms.polarizationTableName(),polTD,option);
335 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::POLARIZATION),
336 2 : Table(polSetup,nrow));
337 :
338 : // Processor
339 2 : TableDesc processorTD=MSProcessor::requiredTableDesc();
340 2 : SetupNewTable processorSetup(ms.processorTableName(),processorTD,option);
341 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::PROCESSOR),
342 2 : Table(processorSetup,nrow));
343 :
344 : // Source
345 2 : TableDesc sourceTD=MSSource::requiredTableDesc();
346 1 : MSSource::addColumnToDesc(sourceTD, MSSource::SYSVEL);
347 1 : MSSource::addColumnToDesc(sourceTD, MSSource::NUM_LINES);
348 1 : MSSource::addColumnToDesc(sourceTD, MSSource::TRANSITION);
349 1 : MSSource::addColumnToDesc(sourceTD, MSSource::REST_FREQUENCY);
350 :
351 2 : SetupNewTable sourceSetup(ms.sourceTableName(),sourceTD,option);
352 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SOURCE),
353 2 : Table(sourceSetup,nrow));
354 :
355 :
356 : // SpectralWindow
357 2 : TableDesc spwTD=MSSpectralWindow::requiredTableDesc();
358 1 : if (!cabb) spwTD.addColumn(ScalarColumnDesc<Int>("ATCA_SAMPLER_BITS",
359 0 : "Number of bits used for sampling"));
360 2 : SetupNewTable spectralWindowSetup(ms.spectralWindowTableName(),spwTD,option);
361 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SPECTRAL_WINDOW),
362 2 : Table(spectralWindowSetup,nrow));
363 :
364 : // State
365 2 : TableDesc stateTD=MSState::requiredTableDesc();
366 2 : SetupNewTable stateSetup(ms.stateTableName(),stateTD,option);
367 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::STATE),
368 2 : Table(stateSetup,nrow));
369 :
370 : // SysCal
371 2 : TableDesc sysCalTD=MSSysCal::requiredTableDesc();
372 2 : Vector<String> cols(4);
373 1 : if (!cabb) {
374 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
375 : ("ATCA_SAMP_STATS_NEG",
376 : "Sampler statistics, negative level",
377 0 : IPosition(1,2),ColumnDesc::Direct));
378 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
379 : ("ATCA_SAMP_STATS_ZERO",
380 : "Sampler statistics, zero level",
381 0 : IPosition(1,2),ColumnDesc::Direct));
382 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
383 : ("ATCA_SAMP_STATS_POS",
384 : "Sampler statistics, positive level",
385 0 : IPosition(1,2),ColumnDesc::Direct));
386 0 : cols(1)="ATCA_SAMP_STATS_NEG";
387 0 : cols(2)="ATCA_SAMP_STATS_ZERO";
388 0 : cols(3)="ATCA_SAMP_STATS_POS";
389 : } else {
390 2 : sysCalTD.addColumn(ArrayColumnDesc<Float>
391 : ("ATCA_GTP",
392 : "Noise Cal On+Off Autocorrelation",
393 3 : IPosition(1,2),ColumnDesc::Direct));
394 2 : sysCalTD.addColumn(ArrayColumnDesc<Float>
395 : ("ATCA_SDO",
396 : "Noise Cal On-Off Autocorrelation",
397 3 : IPosition(1,2),ColumnDesc::Direct));
398 2 : sysCalTD.addColumn(ArrayColumnDesc<Float>
399 : ("ATCA_CALJY",
400 : "Calibration factor",
401 3 : IPosition(1,2),ColumnDesc::Direct));
402 :
403 1 : cols(1)="ATCA_GTP";
404 1 : cols(2)="ATCA_SDO";
405 1 : cols(3)="ATCA_CALJY";
406 : }
407 2 : sysCalTD.addColumn(ScalarColumnDesc<Float>
408 : ("ATCA_XY_AMPLITUDE",
409 1 : "Noise source cross correlation amplitude"));
410 2 : sysCalTD.addColumn(ScalarColumnDesc<Float>
411 : ("ATCA_TRACK_ERR_MAX",
412 1 : "Max tracking error over non blanked cycle"));
413 3 : TableQuantumDesc tqd1(sysCalTD,"ATCA_TRACK_ERR_MAX",Unit("arcsec"));
414 1 : tqd1.write(sysCalTD);
415 2 : sysCalTD.addColumn(ScalarColumnDesc<Float>
416 : ("ATCA_TRACK_ERR_RMS",
417 1 : "RMS tracking error over non blanked cycle"));
418 3 : TableQuantumDesc tqd2(sysCalTD,"ATCA_TRACK_ERR_RMS",Unit("arcsec"));
419 1 : tqd2.write(sysCalTD);
420 :
421 2 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TSYS,
422 2 : IPosition(1, 2), ColumnDesc::Direct);
423 2 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TCAL,
424 2 : IPosition(1, 2), ColumnDesc::Direct);
425 2 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TRX,
426 2 : IPosition(1, 2),ColumnDesc::Direct);
427 1 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TSYS_FLAG);
428 1 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::PHASE_DIFF);
429 1 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::PHASE_DIFF_FLAG);
430 1 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TCAL_FLAG);
431 1 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TRX_FLAG);
432 1 : cols(0) = MSSysCal::columnName(MSSysCal::TSYS);
433 1 : sysCalTD.defineHypercolumn("TiledSysCal",2,cols);
434 2 : SetupNewTable sysCalSetup(ms.sysCalTableName(),sysCalTD,option);
435 2 : IncrementalStMan incStMan;
436 1 : sysCalSetup.bindAll(incStMan);
437 3 : TiledColumnStMan tiledStManSysCal("TiledSysCal",IPosition(2,2,1024));
438 :
439 5 : for (uInt i=0; i<cols.nelements(); i++) {
440 4 : sysCalSetup.bindColumn(cols(i),tiledStManSysCal);
441 : }
442 2 : ms.rwKeywordSet().defineTable(ms.keywordName(MS::SYSCAL),
443 2 : Table(sysCalSetup,nrow));
444 :
445 : // Weather
446 2 : TableDesc weatherTD=MSWeather::requiredTableDesc();
447 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::PRESSURE);
448 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::REL_HUMIDITY);
449 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::TEMPERATURE);
450 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_DIRECTION);
451 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_SPEED);
452 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::PRESSURE_FLAG);
453 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::REL_HUMIDITY_FLAG);
454 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::TEMPERATURE_FLAG);
455 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_DIRECTION_FLAG);
456 1 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_SPEED_FLAG);
457 2 : weatherTD.addColumn(ScalarColumnDesc<Float>
458 : ("ATCA_RAIN_GAUGE",
459 1 : "Rain since 10am local time"));
460 3 : TableQuantumDesc tqds0(weatherTD,"ATCA_RAIN_GAUGE",Unit("mm"));
461 1 : tqds0.write(weatherTD);
462 2 : weatherTD.addColumn(ScalarColumnDesc<Float>
463 : ("ATCA_SEEMON_PHASE",
464 1 : "Seeing monitor raw phase at 22GHz"));
465 3 : TableQuantumDesc tqds1(weatherTD,"ATCA_SEEMON_PHASE",Unit("rad"));
466 1 : tqds1.write(weatherTD);
467 2 : weatherTD.addColumn(ScalarColumnDesc<Float>
468 : ("ATCA_SEEMON_RMS",
469 1 : "Seeing monitor RMS phase"));
470 3 : TableQuantumDesc tqds2(weatherTD,"ATCA_SEEMON_RMS",Unit("mm"));
471 1 : tqds2.write(weatherTD);
472 2 : weatherTD.addColumn(ScalarColumnDesc<Bool>
473 : ("ATCA_SEEMON_FLAG",
474 1 : "Seeing monitor flag"));
475 2 : SetupNewTable weatherSetup(ms.weatherTableName(),weatherTD,option);
476 2 : ms.rwKeywordSet().defineTable(ms.keywordName(MS::WEATHER),
477 2 : Table(weatherSetup,nrow));
478 :
479 : // ATCA_SCAN_INFO
480 2 : TableDesc atsiTD;
481 2 : atsiTD.addColumn(ScalarColumnDesc<Int>("ANTENNA_ID",
482 1 : "Antenna Id"));
483 2 : atsiTD.addColumn(ScalarColumnDesc<Int>("SCAN_ID",
484 1 : "Scan number from main table"));
485 2 : atsiTD.addColumn(ScalarColumnDesc<Int>("SPECTRAL_WINDOW_ID",
486 1 : "Spectral window Id"));
487 2 : atsiTD.addColumn(ArrayColumnDesc<Int>("FINE_ATTENUATOR",
488 : "Fine Attenuator setting A,B",
489 3 : IPosition(1,2),ColumnDesc::Direct));
490 2 : atsiTD.addColumn(ArrayColumnDesc<Int>("COARSE_ATTENUATOR",
491 : "COARSE Attenuator setting A,B",
492 3 : IPosition(1,2),ColumnDesc::Direct));
493 2 : atsiTD.addColumn(ArrayColumnDesc<Int>("MM_ATTENUATOR",
494 : "mm Attenuator setting A,B",
495 3 : IPosition(1,2),ColumnDesc::Direct));
496 2 : atsiTD.addColumn(ArrayColumnDesc<Float>("SUBREFLECTOR",
497 : "Subreflector position(center,edge/tilt)",
498 3 : IPosition(1,2),ColumnDesc::Direct));
499 3 : TableQuantumDesc tqd3(atsiTD,"SUBREFLECTOR",Unit("m"));
500 1 : tqd3.write(atsiTD);
501 2 : atsiTD.addColumn(ScalarColumnDesc<String>("CORR_CONFIG",
502 1 : "Correlator configuration"));
503 2 : atsiTD.addColumn(ScalarColumnDesc<String>("SCAN_TYPE",
504 1 : "Scan type"));
505 2 : atsiTD.addColumn(ScalarColumnDesc<String>("COORD_TYPE",
506 1 : "CAOBS Coordinate type"));
507 2 : atsiTD.addColumn(ScalarColumnDesc<String>("POINTING_INFO",
508 1 : "Pointing info - details of last point scan"));
509 2 : atsiTD.addColumn(ScalarColumnDesc<Bool>("LINE_MODE",
510 1 : "Line Mode: true=spectrum, false=mfs"));
511 2 : atsiTD.addColumn(ScalarColumnDesc<Int>("CACAL_CNT",
512 1 : "Online calibration counter"));
513 :
514 2 : SetupNewTable scanInfoSetup(ms.tableName()+"/ATCA_SCAN_INFO",atsiTD,option);
515 2 : IncrementalStMan incSIStMan("ISMScanInfo");
516 2 : StandardStMan stdSIStMan("SSMScanInfo",32768);
517 1 : scanInfoSetup.bindAll(incSIStMan);
518 1 : scanInfoSetup.bindColumn("ANTENNA_ID",stdSIStMan);
519 2 : ms.rwKeywordSet().defineTable("ATCA_SCAN_INFO",
520 2 : Table(scanInfoSetup,nrow));
521 :
522 : // update the references to the subtable keywords
523 1 : ms.initRefs();
524 1 : }
525 :
526 :
527 1 : void ATCAFiller::init()
528 : {
529 : //Initialize selection
530 1 : scanRange(0, 0);
531 1 : freqRange(0.0,1.e30);
532 2 : Vector<String> fieldNames(0);
533 1 : fields(fieldNames);
534 :
535 : // extra ID columns in main table
536 1 : colSysCalIdAnt1.attach(atms_p,"ATCA_SYSCAL_ID_ANT1");
537 1 : colSysCalIdAnt2.attach(atms_p,"ATCA_SYSCAL_ID_ANT2");
538 :
539 : // extra spectralWindow table columns
540 1 : if (!cabb_p) colSamplerBits.attach(atms_p.spectralWindow(),"ATCA_SAMPLER_BITS");
541 :
542 : // extra syscal table columns
543 1 : if (cabb_p) {
544 1 : colGTP.attach(atms_p.sysCal(),"ATCA_GTP");
545 1 : colSDO.attach(atms_p.sysCal(),"ATCA_SDO");
546 1 : colCalJy.attach(atms_p.sysCal(),"ATCA_CALJY");
547 : } else {
548 0 : colSamplerStatsNeg.attach(atms_p.sysCal(),"ATCA_SAMP_STATS_NEG");
549 0 : colSamplerStatsZero.attach(atms_p.sysCal(),"ATCA_SAMP_STATS_ZERO");
550 0 : colSamplerStatsPos.attach(atms_p.sysCal(),"ATCA_SAMP_STATS_POS");
551 : }
552 : // cParAngle.attach(sysCalTable,"ParalAngle");
553 1 : colXYAmplitude.attach(atms_p.sysCal(),"ATCA_XY_AMPLITUDE");
554 1 : colTrackErrMax.attach(atms_p.sysCal(),"ATCA_TRACK_ERR_MAX");
555 1 : colTrackErrRMS.attach(atms_p.sysCal(),"ATCA_TRACK_ERR_RMS");
556 :
557 : // ScanInfo table columns
558 1 : Table sit(atms_p.keywordSet().asTable("ATCA_SCAN_INFO"));
559 1 : colScanInfoAntId.attach(sit,"ANTENNA_ID");
560 1 : colScanInfoScanId.attach(sit,"SCAN_ID");
561 1 : colScanInfoSpWId.attach(sit,"SPECTRAL_WINDOW_ID");
562 1 : colScanInfoCacal.attach(sit,"CACAL_CNT");
563 1 : colScanInfoFine.attach(sit,"FINE_ATTENUATOR");
564 1 : colScanInfoCoarse.attach(sit,"COARSE_ATTENUATOR");
565 1 : colScanInfommAtt.attach(sit,"MM_ATTENUATOR");
566 1 : colScanInfoSubreflector.attach(sit,"SUBREFLECTOR");
567 1 : colScanInfoCorrConfig.attach(sit,"CORR_CONFIG");
568 1 : colScanInfoScanType.attach(sit,"SCAN_TYPE");
569 1 : colScanInfoCoordType.attach(sit,"COORD_TYPE");
570 1 : colScanInfoPointInfo.attach(sit,"POINTING_INFO");
571 1 : colScanInfoLineMode.attach(sit,"LINE_MODE");
572 :
573 : // Extra weather table columns
574 1 : colWeatherRainGauge.attach(atms_p.weather(),"ATCA_RAIN_GAUGE");
575 1 : colWeatherSeeMonPhase.attach(atms_p.weather(),"ATCA_SEEMON_PHASE");
576 1 : colWeatherSeeMonRMS.attach(atms_p.weather(),"ATCA_SEEMON_RMS");
577 1 : colWeatherSeeMonFlag.attach(atms_p.weather(),"ATCA_SEEMON_FLAG");
578 :
579 1 : nScan_p=nSpW_p=nField_p=scanNo_p=0;
580 1 : gotAN_p=false;
581 1 : }
582 :
583 : // list the current scan and return the current day in mjd
584 1 : void ATCAFiller::listScan(Double & mjd, Int scan, Double ut)
585 : {
586 : // Convert observation date to mjd
587 : Int day,month,year;
588 1 : sscanf(names_.datobs,"%4d-%2d-%2d",&year,&month,&day);
589 1 : MVTime mjd_date(year,month,(Double)day);
590 1 : mjd=mjd_date.second();
591 1 : mjd_date=MVTime((mjd_date.second()+ut)/C::day);
592 1 : os_p << LogIO::NORMAL << "Scan # : "<< scan << endl;
593 1 : os_p << LogIO::NORMAL << "Object : "<< String(names_.object,16) << endl;
594 1 : os_p << LogIO::NORMAL << "Date : "<< mjd_date.string(MVTime::YMD)
595 1 : << LogIO::POST;
596 1 : }
597 :
598 :
599 1 : Bool ATCAFiller::fill() {
600 1 : if (!init_p) return false;
601 1 : if (!appendMode_p) {
602 1 : firstHeader_p = true;
603 1 : skipScan_p=false;
604 1 : skipData_p=false;
605 :
606 1 : nScan_p=1; // we've seen the 1st header
607 1 : scanNo_p=-1; // make zero based for storage in MS
608 :
609 : Int fileno;
610 1 : Int offset=0;
611 :
612 2 : for (fileno=0; fileno<Int(rpfitsFiles_p.nelements())-offset; fileno++) {
613 1 : listHeader_p = true;
614 1 : currentFile_p = rpfitsFiles_p(fileno);
615 1 : fill1(currentFile_p);
616 : }
617 :
618 1 : os_p << LogIO::DEBUGGING << "FillFeed" << LogIO::POST;
619 1 : fillFeedTable();
620 :
621 1 : fillObservationTable();
622 :
623 1 : fillMeasureReferences();
624 1 : os_p << LogIO::DEBUGGING << "#spectral windows " << nSpW_p << LogIO::POST;
625 :
626 :
627 : } else { // appendMode
628 0 : if (!eof_p) {
629 0 : fill1(currentFile_p);
630 : } else {
631 :
632 0 : RegularFile rFile(currentFile_p);
633 0 : uInt newSize = rFile.size();
634 0 : os_p << LogIO::NORMAL << "new file size " << newSize << LogIO::POST;
635 :
636 0 : if ( newSize != fileSize_p) {
637 0 : fill1(currentFile_p);
638 : } else {
639 :
640 0 : Int jstat=1; // close file
641 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w,
642 : &flg, &bin, &if_no, &sourceno);
643 :
644 0 : os_p << LogIO::NORMAL << "Look for next file ..." << LogIO::POST;
645 :
646 0 : Vector<String> elts=stringToVector(currentFile_p,std::regex("/"));
647 0 : String rpfitsDir_p = "";
648 0 : if (elts.nelements()>1) {
649 0 : Int m = elts.nelements() - 1;
650 0 : for (Int n=0; n < m; n++) {
651 0 : rpfitsDir_p = rpfitsDir_p + elts(n) + "/";
652 : }
653 : }
654 0 : os_p << LogIO::DEBUGGING << "RPFITSDIR : " << rpfitsDir_p << LogIO::POST;
655 :
656 0 : Directory dir(rpfitsDir_p);
657 0 : Regex regexp(".*\\.[cxv]+[0-9]+");
658 :
659 0 : DirectoryIterator dirIter(dir, regexp);
660 :
661 0 : String entry;
662 0 : Bool found = false;
663 0 : while (!dirIter.pastEnd()) {
664 0 : entry = rpfitsDir_p+dirIter.name();
665 0 : os_p << LogIO::DEBUGGING << " file: "<< entry << LogIO::POST;
666 0 : if (found) break;
667 0 : if (entry == currentFile_p) found = true;
668 0 : dirIter++;
669 : }
670 :
671 0 : os_p << LogIO::NORMAL << " new file: "<< entry << LogIO::POST;
672 :
673 0 : if (entry == currentFile_p) {
674 0 : os_p << LogIO::NORMAL << " No new file..."<< LogIO::POST;
675 : } else {
676 0 : String oldstr = String(currentFile_p.at(1, currentFile_p.length()));
677 0 : os_p << LogIO::DEBUGGING << " oldstr... "<< oldstr << LogIO::POST;
678 0 : String newstr = String(entry.at(1, entry.length()));
679 0 : os_p << LogIO::DEBUGGING << " newstr... "<< newstr << LogIO::POST;
680 :
681 0 : if (oldstr.after(Regex(".*\\.")) != newstr.after(Regex(".*\\."))) {
682 0 : os_p << LogIO::NORMAL << " Project changed..."<< LogIO::POST;
683 : } else {
684 0 : currentFile_p = entry;
685 0 : listHeader_p = true;
686 0 : fill1(currentFile_p);
687 : }
688 : }
689 : }
690 : }
691 : }
692 : // flag the last integration for shadowing
693 1 : if (shadow_p>0) shadow(0,true);
694 :
695 : // flush the table and unlock it
696 : // NOTE: this still keeps a reference to the table which causes problems
697 : // when the ms do calls close (it wants to write to the table)
698 1 : flush();
699 1 : unlock();
700 :
701 1 : return true;
702 : }
703 :
704 1 : Bool ATCAFiller::checkCABB(const String& rpfitsFile)
705 : {
706 : Int jstat;
707 2 : Regex trailing(" *$"); // trailing blanks
708 :
709 2 : String file = rpfitsFile;
710 1 : os_p << LogIO::NORMAL <<"Checking header of file "<<file<< LogIO::POST;
711 1 : strcpy(names_.file,file.chars());
712 1 : jstat=-2; // open rpfits file and read first header
713 1 : param_.ncard=-1; // collect cards into card array
714 1 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
715 : &if_no, &sourceno);
716 1 : if (jstat==-1) {
717 0 : os_p << LogIO::SEVERE << " Error opening RPFits file: "<<file
718 0 : << LogIO::POST;
719 0 : return false;
720 : }
721 : // read INSTRUMENT keyword
722 2 : String instrument = String(names_.instrument,16).before(trailing);
723 : // close the file
724 1 : jstat=1;
725 1 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
726 : &if_no, &sourceno);
727 1 : if (instrument=="ATCABB") {
728 1 : os_p<< LogIO::NORMAL<<"CABB data detected"<<LogIO::POST;
729 1 : return true;
730 : }
731 0 : else return false;
732 : }
733 :
734 1 : Bool ATCAFiller::fill1(const String& rpfitsFile)
735 : {
736 : Int jstat;
737 2 : Regex trailing(" *$"); // trailing blanks
738 :
739 2 : String file = rpfitsFile;
740 1 : if (listHeader_p == true) {
741 1 : os_p << LogIO::NORMAL <<"Reading file "<<file<< LogIO::POST;
742 1 : strcpy(names_.file,file.chars());
743 1 : jstat=-2; // open rpfits file and read first header
744 1 : param_.ncard=-1; // collect cards into card array
745 1 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
746 : &if_no, &sourceno);
747 1 : if (jstat==-1) {
748 0 : os_p << LogIO::SEVERE << " Error opening RPFits file: "<<file
749 0 : << LogIO::POST;
750 0 : return false;
751 : }
752 : }
753 : // Otherwise we enter fill1 for the second time with the same file name
754 :
755 : // prepare to read data blocks, rpfitsin will tell us if
756 : // there is a header to be read instead
757 1 : jstat=0;
758 :
759 1 : Double lastUT=0,lastUT2=0;
760 1 : Int lastScanNo=-1,lastIF=-1;
761 1 : flagCount_p=0;
762 1502 : while (jstat==0) {
763 1501 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
764 : &if_no, &sourceno);
765 :
766 1501 : switch(jstat) {
767 0 : case -1: // read failed
768 0 : os_p << LogIO::WARN <<
769 0 : "rpfitsin: read failed, retrying"<< LogIO::POST;
770 0 : jstat=0;
771 0 : break;
772 1482 : case 0: // found data: store header (incl Tables) and data
773 1482 : if (listHeader_p) { listScan(mjd0_p,nScan_p,ut); listHeader_p=false;}
774 :
775 : // do we want this scan?
776 1482 : if (skipScan_p || nScan_p<firstScan_p ||
777 1482 : (lastScan_p>0 && nScan_p>lastScan_p)) {
778 : //if (!skipScan_p) cerr<<"ss1=T"<<endl;
779 0 : skipScan_p=true;
780 : } else {
781 : // if_no and sourceno are not (always) set properly
782 : // for syscal records
783 1482 : if (baseline==-1) {if_no=1; sourceno=1;}
784 :
785 : // check if if_no valid
786 1482 : if (!if_.if_found|| if_no>if_.n_if) {
787 : //if (!skipScan_p)cerr<<"ss2=T"<<endl;
788 0 : skipScan_p=true;
789 : // assume not wanted or ^ garbled data
790 : } else {
791 : // check if we want this field
792 1482 : if (lastUT!=ut) {
793 138 : String field=String(names_.object,16).before(trailing);
794 69 : if (fieldSelection_p.nelements()>0 && fieldSelection_p(0).length()>0 &&
795 0 : !(anyEQ(fieldSelection_p,field))) {
796 : //cerr<< "Field:"<<field<<"-"<<fieldSelection_p.nelements()<<endl;
797 : //if (!skipScan_p) cerr<<"ss3=T"<<endl;
798 0 : skipScan_p=true;
799 : }
800 : }
801 1482 : if (!skipScan_p) {
802 1482 : if_no--; // go from 1-based to 0-based indexing (f2c)
803 1482 : sourceno--; // same here
804 : // skip unselected or garbled data
805 1482 : skipData_p=Bool(if_no<0 || !selected(if_no) || sourceno<0);
806 : // reject auto correlations?
807 1482 : if (noac_p) skipData_p |= (Int(baseline)%257==0);
808 : //cerr<<"SkipData="<<skipData_p<< " ifno="<<if_no<< " sourceno="<<sourceno<<endl;
809 : }
810 : }
811 : }
812 1482 : if (!skipScan_p && !skipData_p && firstHeader_p && anten_.nant>0) {
813 1 : nAnt_p=anten_.nant;
814 1 : Int NChan=if_.if_nfreq[if_no];
815 1 : Int NPol=if_.if_nstok[if_no];
816 1 : firstHeader_p=false;
817 1 : os_p << LogIO::NORMAL << " First data/header has NAnt="<<nAnt_p<<
818 1 : ", NChan="<<NChan<<", NPol="<<NPol<< LogIO::POST;
819 : }
820 1482 : if (!skipScan_p && !skipData_p) {
821 912 : if (anten_.nant>0 && anten_.nant!=nAnt_p) {
822 0 : os_p << LogIO::WARN << "#antennas changed from "<< nAnt_p <<
823 0 : " to "<<anten_.nant<<", skipping scan"<< LogIO::POST;
824 : //if (!skipScan_p)cerr<<"ss4=T"<<endl;
825 0 : skipScan_p=true;
826 : }
827 : }
828 1482 : if (!skipScan_p && !skipData_p && !storedHeader_p) {
829 1 : storeHeader();
830 1 : scanNo_p++;
831 1 : storedHeader_p=true;
832 : }
833 1482 : if (!skipScan_p) {
834 1482 : if (baseline==-1) {
835 : // we want at least some of the syscal data for this scan
836 57 : storeSysCal();
837 : }
838 1425 : else if (!skipData_p) {
839 855 : if (lastUT2!=ut || lastIF!=if_no){
840 57 : checkSpW(if_no); // we need to check these every integration to
841 57 : if (lastUT2!=ut) checkField();// cope with freq. switching and mosaicing
842 57 : lastUT2=ut;
843 57 : lastIF=if_no;
844 : }
845 855 : storeData();
846 : }
847 : }
848 1482 : skipData_p=false;
849 1482 : break;
850 0 : case 1:
851 0 : if (skipScan_p || scanNo_p==lastScanNo) {
852 : //cerr<<"ScanNo_p="<<scanNo_p<<" lastScanNo="<<lastScanNo<<" skipScan_p="<<skipScan_p<<endl;
853 0 : os_p << LogIO::NORMAL << "Scan "<<nScan_p <<" skipped"<<LogIO::POST;
854 : } else {
855 0 : os_p << LogIO::NORMAL << "Scan "<<nScan_p <<" stored"<< LogIO::POST;
856 : }
857 0 : lastScanNo=scanNo_p;
858 0 : nScan_p++;
859 0 : errCount_p=0;
860 0 : os_p << LogIO::DEBUGGING << "Read new header "<< LogIO::POST;
861 0 : jstat=-1;
862 0 : param_.ncard=-1;
863 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w,
864 : &flg, &bin, &if_no, &sourceno);
865 0 : jstat=0;
866 :
867 0 : listScan(mjd0_p,nScan_p,ut);
868 0 : storedHeader_p=false;
869 0 : skipScan_p=false; //cerr<<"ss=F"<<endl;
870 0 : eof_p = false;
871 0 : if (!appendMode_p) break;
872 : // suppress this break to return after each end of scan
873 : // (actually each beginning of a new scan) in on line mode
874 :
875 : case 3:
876 1 : if (jstat == 3) { // because of the break suppression in case 1
877 1 : if (!skipScan_p && !skipData_p && !storedHeader_p) {
878 0 : storeHeader(True); // store sole header at end of file
879 0 : scanNo_p++; // to capture last commands & log messages
880 0 : storedHeader_p=true;
881 : }
882 1 : os_p << LogIO::NORMAL << "End of File" << LogIO::POST;
883 1 : eof_p = true;
884 1 : nScan_p++; // increment for next scan
885 : // print flagging stats for last scan
886 1 : if (flagCount_p(COUNT)>0) {
887 2 : Vector<Float> perc(NFLAG);
888 7 : for (Int i=0; i<NFLAG; i++) perc(i)=0.1*((1000*flagCount_p(i))/flagCount_p(COUNT));
889 1 : os_p<< LogIO::NORMAL <<"Number of rows selected : "<<flagCount_p(COUNT)<<endl;
890 1 : os_p<< LogIO::NORMAL <<"Flagged : "<<perc(FLAG)<<"%"<<endl;
891 1 : os_p<< LogIO::NORMAL <<" Antenna off source : "<<perc(ONLINE)<<"%"<<endl;
892 1 : os_p<< LogIO::NORMAL <<" ScanType (Point/Paddle): "<<perc(SCANTYPE)<<"%"<<endl;
893 1 : if (!cabb_p) os_p<< LogIO::NORMAL <<" Bad Sampler stats : "<<perc(SYSCAL)<<"%"<<endl;
894 1 : if (shadow_p>0) os_p<< LogIO::NORMAL <<" Antenna shadowed : "<<perc(SHADOW)<<"%"<<LogIO::POST;
895 : }
896 : }
897 :
898 1 : if (appendMode_p) {
899 : // after bumping one time at the end of the file or after the
900 : // first scan
901 :
902 0 : RegularFile rFile(file);
903 0 : fileSize_p = rFile.size();
904 :
905 0 : os_p << LogIO::NORMAL << "old file size " << fileSize_p
906 0 : << " Waiting ..." << LogIO::POST;
907 :
908 0 : flush();
909 0 : unlock();
910 :
911 : } else {
912 1 : jstat=1; // close file
913 1 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w,
914 : &flg, &bin, &if_no, &sourceno);
915 : }
916 1 : jstat=1; // exit loop
917 1 : break;
918 0 : case 4:
919 0 : os_p << LogIO::WARN
920 : << "rpfitsin: found FG Table, ignoring it"
921 0 : << LogIO::POST;
922 0 : jstat=0;
923 0 : break;
924 18 : case 5:
925 : // jstat == 5 is usually record padding at the end of an integration.
926 : // just continue on loading..
927 18 : jstat=0;
928 18 : break;
929 0 : default:
930 0 : os_p << LogIO::WARN << "unknown rpfitsin return value: "
931 0 : <<jstat<< LogIO::POST;
932 0 : jstat=0;
933 : }
934 1501 : lastUT=ut;
935 : }
936 :
937 : // uncomment this if things take unexpectedly long
938 : // dataAccessor_p.showCacheStatistics(cout);
939 : // sigmaAccessor_p.showCacheStatistics(cout);
940 1 : return true;
941 : }
942 :
943 :
944 :
945 :
946 1 : void ATCAFiller::fillMeasureReferences()
947 : {
948 1 : String key("MEASURE_REFERENCE");
949 1 : msc_p->time().rwKeywordSet().define(key,"UTC");
950 1 : msc_p->uvw().rwKeywordSet().define(key,"J2000");
951 1 : msc_p->antenna().position().rwKeywordSet().define(key,"ITRF");
952 1 : msc_p->feed().time().rwKeywordSet().define(key,"UTC");
953 1 : msc_p->field().time().rwKeywordSet().define(key,"UTC");
954 1 : msc_p->field().delayDir().rwKeywordSet().define(key,"J2000");
955 1 : msc_p->field().phaseDir().rwKeywordSet().define(key,"J2000");
956 1 : msc_p->field().referenceDir().rwKeywordSet().define(key,"J2000");
957 1 : msc_p->source().time().rwKeywordSet().define(key,"UTC");
958 1 : msc_p->source().direction().rwKeywordSet().define(key,"J2000");
959 1 : msc_p->history().time().rwKeywordSet().define(key,"UTC");
960 1 : msc_p->spectralWindow().chanFreq().rwKeywordSet().define(key,"TOPO");
961 1 : msc_p->spectralWindow().refFrequency().rwKeywordSet().define(key,"TOPO");
962 :
963 : MFrequency::Types tp;
964 1 : MFrequency::getType(tp, msc_p->spectralWindow().refFrequency().keywordSet().asString("MEASURE_REFERENCE"));
965 1 : Int meas_freq_ref = tp;
966 1 : msc_p->spectralWindow().measFreqRef().fillColumn(meas_freq_ref);
967 1 : msc_p->sysCal().time().rwKeywordSet().define(key,"UTC");
968 1 : msc_p->weather().time().rwKeywordSet().define(key,"UTC");
969 1 : msc_p->pointing().direction().rwKeywordSet().define(key,"J2000");
970 1 : msc_p->pointing().pointingOffset().rwKeywordSet().define(key,"AZEL");
971 1 : }
972 :
973 2 : ATCAFiller & ATCAFiller::scanRange(Int first, Int last)
974 : {
975 2 : firstScan_p=first;
976 2 : lastScan_p=max(0,last);
977 2 : return *this;
978 : }
979 :
980 2 : ATCAFiller & ATCAFiller::freqRange(Double low, Double high)
981 : {
982 2 : lowFreq_p=low;
983 2 : highFreq_p=high;
984 2 : return *this;
985 : }
986 :
987 1 : ATCAFiller & ATCAFiller::freqSel(const Vector<Int>& spws)
988 : {
989 1 : spws_p=spws;
990 1 : return *this;
991 : }
992 :
993 2 : ATCAFiller & ATCAFiller::fields(const Vector<String>& fieldList)
994 : {
995 2 : fieldSelection_p=fieldList;
996 2 : return *this;
997 : }
998 :
999 0 : ATCAFiller & ATCAFiller::bandwidth1(Int bandwidth1)
1000 : {
1001 0 : bandWidth1_p=0;
1002 0 : for (Int bw=2; bw<=256; bw*=2) {
1003 0 : if(bw==bandwidth1) {
1004 0 : bandWidth1_p=bw;
1005 0 : break;
1006 : }
1007 : }
1008 0 : return *this;
1009 : }
1010 :
1011 0 : ATCAFiller & ATCAFiller::numChan1(Int numchan1)
1012 : {
1013 0 : numChan1_p=0;
1014 0 : for (Int nchan=32; nchan<=16384; nchan*=2) {
1015 0 : if((numchan1==nchan) || (numchan1==(nchan+1))) {
1016 0 : numChan1_p=nchan+1;
1017 0 : break;
1018 : }
1019 : }
1020 0 : return *this;
1021 : }
1022 :
1023 0 : ATCAFiller & ATCAFiller::shadow(Float diam)
1024 : {
1025 0 : shadow_p=0;
1026 0 : if (diam>0) shadow_p=diam;
1027 0 : return *this;
1028 : }
1029 :
1030 1 : ATCAFiller & ATCAFiller::edge(float edge)
1031 : {
1032 1 : edge_p=edge;
1033 1 : return *this;
1034 : }
1035 :
1036 1 : void ATCAFiller::storeHeader(bool last) {
1037 : // Bool skip=false;
1038 2 : Regex trailing(" *$"); // trailing blanks
1039 :
1040 : // First check antenna Table and store any new stations encountered
1041 1 : if (!gotAN_p && anten_.nant>0) {
1042 1 : if (anten_.an_found) {
1043 1 : gotAN_p=true;
1044 :
1045 7 : for (Int i=0; i<nAnt_p; i++) {
1046 6 : Int ant=Int(anten_.ant_num[i]); // 1-based !
1047 6 : if (ant!=i+1) {
1048 0 : os_p << LogIO::SEVERE
1049 : << "AN table corrupt, will try next one " << ant <<":"<<i<<":"<<nAnt_p
1050 0 : << LogIO::POST;
1051 0 : gotAN_p=false;
1052 0 : break;
1053 : }
1054 6 : atms_p.antenna().addRow();
1055 6 : Int n=atms_p.antenna().nrow()-1;
1056 :
1057 12 : String instrument = String(names_.instrument,16).before(trailing);
1058 12 : String station(&names_.sta[i*8],8);
1059 6 : os_p << LogIO::NORMAL <<"Antenna : "<< station << " ";
1060 12 : Vector<Double> xyz(3);
1061 6 : xyz(0)=doubles_.x[i];xyz(1)=doubles_.y[i];xyz(2)=doubles_.z[i];
1062 6 : os_p << " position:"<<xyz << endl;
1063 6 : if (instrument=="ATCA" || cabb_p) {
1064 12 : String antName;
1065 : // construct antenna name
1066 12 : ostringstream ostr; ostr << "CA0" <<i+1;
1067 6 : String str = ostr.str();
1068 6 : msc_p->antenna().name().put(n,str);
1069 6 : msc_p->antenna().station().put(n,atcaPosToStation(xyz));
1070 : } else {
1071 0 : msc_p->antenna().station().put(n,station.before(trailing));
1072 0 : msc_p->antenna().name().put(n,station.before(trailing));
1073 : }
1074 6 : msc_p->antenna().orbitId().put(n,-1);
1075 6 : msc_p->antenna().phasedArrayId().put(n,-1);
1076 6 : msc_p->antenna().dishDiameter().put(n,22.0);
1077 6 : msc_p->antenna().type().put(n, "GROUND-BASED");
1078 :
1079 6 : if (anten_.ant_mount[i]==0) {
1080 6 : msc_p->antenna().mount().put(n,"alt-az");
1081 : }
1082 : else {
1083 0 : msc_p->antenna().mount().put(n,"bizarre");
1084 : }
1085 6 : msc_p->antenna().position().put(n,xyz);
1086 12 : Vector<Double> offset(3); offset=0.;
1087 : // todo: figure out coordinate system of offset
1088 6 : offset(0)=doubles_.axis_offset[i];
1089 6 : msc_p->antenna().offset().put(n,offset);
1090 : // todo: values below may have to be stored elsewhere
1091 : //Vector<String> ft(2);
1092 : //ft(0)=String(&names_.feed_type[i*4],2);
1093 : //ft(1)=String(&names_.feed_type[i*4+2],2);
1094 : // cFeedType.put(n,ft);
1095 : //Vector<Double> pa(2);
1096 : //pa(0)=doubles_.feed_pa[i*2];pa(1)=doubles_.feed_pa[i*2+1];
1097 : // cFeedPA.put(n,pa);
1098 : //Matrix<Double> feedcal(1,1);//feedcal :split out on if and pol?
1099 : //feedcal(0,0)=0.0;
1100 : // cFeedCal.put(n,feedcal);
1101 : }
1102 1 : os_p.post();
1103 : } else {
1104 0 : os_p << LogIO::SEVERE <<
1105 0 : "No AN Table found before start of data!"<< LogIO::POST;
1106 : }
1107 : }
1108 :
1109 : // Check if current if_no already in SpW Table, add it if not.
1110 1 : if (!last) checkSpW(if_no);
1111 :
1112 : // Check if current Observation info already stored, add it if not
1113 1 : checkObservation();
1114 :
1115 : // Store the ATCA header cards
1116 1 : storeATCAHeader(last);
1117 :
1118 : // Check if we've seen current source before, if not add to table
1119 1 : if (!last) checkField();
1120 1 : }
1121 :
1122 1 : void ATCAFiller::storeATCAHeader(Bool last) {
1123 1 : uInt ncard = abs(param_.ncard);
1124 : //cout<<" #cards = "<<ncard<<endl;
1125 1 : String cards = String(names_.card,ncard*80);
1126 1 : const Int nTypes = 10;
1127 1 : Block<String> types(nTypes);
1128 1 : types[0]= "OBSLOG";
1129 1 : types[1] = "ATTEN";
1130 1 : types[2] = "SUBREFL";
1131 1 : types[3] = "CORR_CFG";
1132 1 : types[4] = "SCANTYPE";
1133 1 : types[5] = "COORDTYP";
1134 1 : types[6] = "LINEMODE";
1135 1 : types[7] = "CACALCNT";
1136 1 : types[8] = "POINTCOR";
1137 1 : types[9] = "POINTINF";
1138 1 : String obsLog = "";
1139 1 : String config = "none";
1140 1 : const Int nIF = if_.n_if;
1141 1 : const Int nAnt = anten_.nant;
1142 1 : Cube<Int> fine(2,nIF,nAnt,0),coarse(2,nIF,6,0);
1143 1 : Matrix<Int> mmAtt(nIF,nAnt,0);
1144 1 : Vector<Float> subrPos(nAnt,0),subrTilt(nAnt,0);
1145 1 : Matrix<Float> pointCorr(2,nAnt,0);
1146 1 : newPointingCorr_p=false;
1147 1 : flagScanType_p=false;
1148 1 : String scanType,coordType,pointInfo;
1149 1 : Vector<Bool> lineMode(nIF);
1150 1 : Int cacalCnt=0;
1151 1 : const Regex trailing(" *$"); // trailing blanks
1152 1 : Bool foundAny = false;
1153 1 : String::size_type pos=cards.find("EPHEM12"); // last 'standard' card
1154 : //cout << "pos="<<pos<<" String::npos=="<<String::npos<<endl;
1155 1 : if (pos==String::npos) return;
1156 1 : uInt iFirst=pos/80+1;
1157 43 : for (uInt i=iFirst; i<ncard; i++) {
1158 : // extract card
1159 84 : String card = cards.substr(i*80,80);
1160 : //cout << "card = "<< card<<endl;
1161 : // read antenna number (if present)
1162 42 : Int ant = 0;
1163 42 : pos=card.find("CA0");
1164 42 : if (pos!=String::npos) istringstream(card.substr(pos+2,2))>> ant;
1165 42 : ant-=1; // zero based indexing
1166 42 : if (ant>=nAnt) ant = -1;
1167 462 : for (Int j=0; j<nTypes; j++) {
1168 420 : if (card.find(types[j])==0) {
1169 20 : foundAny=true;
1170 : //cout << "Found card :"<<types[j]<<", ant="<<ant<<endl;
1171 20 : switch (j) {
1172 3 : case 0: obsLog+=card.substr(8,72).before(trailing);
1173 3 : obsLog+="\n";
1174 3 : break;
1175 6 : case 1: { //ATTEN
1176 6 : String::size_type pos=card.find("FINE=");
1177 6 : if (pos!=String::npos && ant>=0) {
1178 0 : for (Int k=0; k<nIF; k++) {
1179 0 : istringstream(card.substr(pos+5+k*2,1))>>fine(0,k,ant);
1180 0 : istringstream(card.substr(pos+6+k*2,1))>>fine(1,k,ant);
1181 : }
1182 : }
1183 6 : pos=card.find("COARSE=");
1184 6 : if (pos!=String::npos && ant>=0) {
1185 0 : for (Int k=0; k<nIF; k++) {
1186 0 : istringstream(card.substr(pos+7+k*2,1))>>coarse(0,k,ant);
1187 0 : istringstream(card.substr(pos+8+k*2,1))>>coarse(1,k,ant);
1188 : }
1189 : }
1190 6 : pos=card.find("MM=");
1191 6 : if (pos!=String::npos && ant>=0) {
1192 0 : istringstream(card.substr(pos+3,2))>>mmAtt(0,ant);
1193 0 : istringstream(card.substr(pos+6,2))>>mmAtt(1,ant);
1194 : }
1195 : }
1196 6 : break;
1197 6 : case 2: { //SUBREFL
1198 6 : String::size_type pos=card.find("POS=");
1199 6 : if (pos!=String::npos && ant>=0) {
1200 6 : istringstream(card.substr(pos+4,6))>>subrPos(ant);
1201 : }
1202 6 : pos=card.find("TILT=");
1203 6 : if (pos!=String::npos && ant>=0) {
1204 6 : istringstream(card.substr(pos+5,6))>>subrTilt(ant);
1205 : }
1206 : }
1207 6 : break;
1208 1 : case 3: {// CORR_CFG
1209 1 : String::size_type pos=card.find("=");
1210 1 : if (pos!=String::npos) {
1211 1 : config=card.substr(pos+3,80-pos-3).before(trailing);
1212 : }
1213 : }
1214 1 : break;
1215 1 : case 4: {// SCANTYPE
1216 1 : String::size_type pos=card.find("=");
1217 1 : if (pos!=String::npos) {
1218 1 : scanType=card.substr(pos+2,80-pos-3).before(trailing);
1219 : }
1220 2 : if (scanType=="PADDLE"||scanType=="POINT"||scanType=="XPOINT"
1221 2 : ||scanType=="EARLY") {
1222 0 : flagScanType_p=true;
1223 : }
1224 : }
1225 1 : break;
1226 1 : case 5: {// COORDTYP
1227 1 : String::size_type pos=card.find("=");
1228 1 : if (pos!=String::npos) {
1229 1 : coordType=card.substr(pos+2,80-pos-3).before(trailing);
1230 : }
1231 : }
1232 1 : break;
1233 1 : case 6: {// LINEMODE
1234 1 : String::size_type pos=card.find("=");
1235 1 : if (pos!=String::npos) {
1236 1 : lineMode(0)=(card[pos+2]=='T');
1237 1 : if (nIF>1) lineMode(1)=(card[pos+4]=='T');
1238 : }
1239 : }
1240 1 : break;
1241 1 : case 7: {// CACALCNT
1242 1 : String::size_type pos=card.find("=");
1243 1 : if (pos!=String::npos) {
1244 1 : istringstream(card.substr(pos+1,7))>>cacalCnt;
1245 : }
1246 : }
1247 1 : break;
1248 0 : case 8: {// POINTCOR
1249 0 : String::size_type pos=card.find("Az=");
1250 0 : if (pos!=String::npos && ant>=0) {
1251 0 : istringstream(card.substr(pos+3,6))>>pointCorr(0,ant);
1252 0 : newPointingCorr_p=true;
1253 : }
1254 0 : pos=card.find("El=");
1255 0 : if (pos!=String::npos && ant>=0) {
1256 0 : istringstream(card.substr(pos+3,6))>>pointCorr(1,ant);
1257 0 : newPointingCorr_p=true;
1258 : }
1259 : }
1260 0 : break;
1261 0 : case 9: {// POINTINF
1262 0 : pointInfo = card.substr(9).before(trailing);
1263 : }
1264 0 : break;
1265 0 : default:
1266 0 : cerr << "Missing SCAN_INFO card type in switch statement"<<endl;
1267 : }
1268 : //cout <<" Match: cardname = "<<types[j]<<" : "<<card<<endl;
1269 : }
1270 : }
1271 : }
1272 1 : if (!foundAny) return;
1273 1 : Vector<String> tmp;
1274 1 : if (msc_p->observation().log().isDefined(obsId_p)) {
1275 0 : tmp=msc_p->observation().log()(obsId_p);
1276 : }
1277 1 : if (tmp.nelements()==0) tmp.resize(1);
1278 1 : Int index=tmp.nelements()-1;
1279 1 : tmp(index)+=obsLog;
1280 1 : msc_p->observation().log().put(obsId_p,tmp);
1281 1 : if ((nAnt*nIF)==0 || last) return;
1282 :
1283 : // find out spectral window index of IFs
1284 2 : Vector<Int> spwId(nIF,-1);
1285 1 : if (if_no>=0) spwId(if_no)=spWId_p;
1286 1 : if (nIF>1) {
1287 4 : for (Int ifNum=0; ifNum<nIF; ifNum++) {
1288 3 : if (ifNum!=if_no) {
1289 2 : if (selected(ifNum)) spwId(ifNum)=checkSpW(ifNum);
1290 : }
1291 : }
1292 : }
1293 : // reset to original spW
1294 1 : if (if_no>0) checkSpW(if_no);
1295 :
1296 1 : Int row=msScanInfo_p.nrow();
1297 1 : msScanInfo_p.addRow(nAnt*nIF);
1298 1 : colScanInfoScanId.put(row,scanNo_p+1);
1299 1 : colScanInfoCacal.put(row,cacalCnt);
1300 1 : colScanInfoCorrConfig.put(row,config);
1301 1 : colScanInfoScanType.put(row,scanType);
1302 1 : colScanInfoCoordType.put(row,coordType);
1303 1 : colScanInfoPointInfo.put(row,pointInfo);
1304 1 : if (newPointingCorr_p) {
1305 0 : pointingCorr_p.reference(pointCorr);
1306 : } else {
1307 1 : pointingCorr_p.resize(0,0);
1308 : }
1309 :
1310 2 : Vector<Int> f(2),c(2),m(2);
1311 2 : Vector<Float> sr(2);
1312 4 : for (Int i=0; i<nIF; i++) {
1313 3 : colScanInfoSpWId.put(row,spwId(i));
1314 3 : colScanInfoLineMode.put(row,lineMode(i));
1315 21 : for (Int ant=0; ant<nAnt; ant++) {
1316 18 : colScanInfoAntId.put(row,ant);
1317 18 : if (!cabb_p) {
1318 0 : f(0)=fine(0,i,ant); f(1)=fine(1,i,ant);
1319 0 : c(0)=coarse(0,i,ant); c(1)=coarse(1,i,ant);
1320 0 : colScanInfoFine.put(row,f);
1321 0 : colScanInfoCoarse.put(row,c);
1322 : }
1323 18 : m(0)=mmAtt(0,ant); m(1)=mmAtt(1,ant);
1324 18 : colScanInfommAtt.put(row,m);
1325 18 : sr(0)=subrPos(ant)/1000.0; sr(1)=subrTilt(ant)/1000.; // convert to meter
1326 18 : colScanInfoSubreflector.put(row,sr);
1327 18 : row++;
1328 : }
1329 : }
1330 : }
1331 :
1332 6 : String ATCAFiller::atcaPosToStation(Vector<Double>& xyz) {
1333 6 : String station("NONE");
1334 : // determine station from xyz position
1335 : // Use W106 as reference
1336 6 : Double x106 = -4751615.0l;
1337 6 : Double y106 = 2791719.246l;
1338 6 : Double z106 = -3200483.747l;
1339 6 : Double incr = 6000.0l/392;
1340 6 : Bool north = (xyz(2)-z106)>1.0;
1341 6 : Bool east = (xyz(0)-x106)<1.0;
1342 6 : Int n = Int(floor(0.5l+
1343 6 : sqrt((xyz(0)-x106)*(xyz(0)-x106)+
1344 6 : (xyz(1)-y106)*(xyz(1)-y106)+
1345 6 : (xyz(2)-z106)*(xyz(2)-z106))/incr));
1346 6 : Bool invalid = (n>392);
1347 6 : if (!invalid) {
1348 6 : if (north) {
1349 0 : ostringstream ostr; ostr << "N"<<n;
1350 0 : station = ostr.str();
1351 : } else {
1352 6 : if (east) n = 106 -n; else n += 106;
1353 6 : ostringstream ostr; ostr << "W"<<n;
1354 6 : station = ostr.str();
1355 : }
1356 : }
1357 6 : return station;
1358 : }
1359 :
1360 117 : Int ATCAFiller::checkSpW(Int ifNumber,Bool log) {
1361 : // Check if current if_no already in SpW Table
1362 : // Add it if not, set SpWId_p index for our SpW Table
1363 : // NOTE we should use ddId instead of SpWid everywhere, for now
1364 : // we create one SpWId per ddId so they always match, this may result
1365 : // in duplicate spw rows if only the pol info changed..
1366 117 : Regex trailing(" *$"); // trailing blanks
1367 117 : if (if_.if_found) {
1368 117 : spWId_p=-1;
1369 234 : for (Int i=0; i<nSpW_p; i++) {
1370 231 : Double freq = msc_p->spectralWindow().refFrequency()(i);
1371 231 : Double bw = msc_p->spectralWindow().totalBandwidth()(i);
1372 231 : Int nchan = msc_p->spectralWindow().numChan()(i);
1373 231 : Int polId = msc_p->dataDescription().polarizationId()(i);
1374 231 : Int npol = msc_p->polarization().numCorr()(polId);
1375 : // compare freq and bw, cope with case where birdie option has reduced bw
1376 231 : if (doubles_.if_freq[ifNumber]==freq &&
1377 114 : doubles_.if_bw[ifNumber]>=bw && doubles_.if_bw[ifNumber]<2*bw &&
1378 114 : (if_.if_nfreq[ifNumber] == nchan || (nchan<33 &&
1379 0 : if_.if_nfreq[ifNumber]==33)) &&
1380 114 : if_.if_nstok[ifNumber]==npol) {
1381 114 : spWId_p=i;
1382 114 : break;
1383 : }
1384 : }
1385 117 : if (spWId_p<0) { // i.e. not found
1386 3 : spWId_p=nSpW_p++;
1387 3 : atms_p.spectralWindow().addRow();
1388 3 : Int n=atms_p.spectralWindow().nrow()-1;
1389 3 : Double refFreq=doubles_.if_freq[ifNumber];
1390 3 : msc_p->spectralWindow().refFrequency().put(n,refFreq);
1391 : // no doppler tracking
1392 3 : Int nchan = if_.if_nfreq[ifNumber];
1393 3 : Int npol =if_.if_nstok[ifNumber];
1394 :
1395 3 : if (log) os_p << LogIO::NORMAL<<
1396 : "IF "<< ifNumber+1 <<
1397 : " : Frequency = "<< refFreq/1.e6 << " MHz" <<
1398 3 : ", bandwidth = " << doubles_.if_bw[ifNumber]/1.e6 << "MHz" <<
1399 3 : ", #channels = "<< nchan << LogIO::POST;
1400 :
1401 3 : Double refChan=doubles_.if_ref[ifNumber]-1; // make zero based
1402 3 : Double chanSpac = doubles_.if_bw[ifNumber]/max(1,nchan-1);
1403 3 : Int inversion=if_.if_invert[ifNumber];
1404 3 : msc_p->spectralWindow().netSideband().put(n,inversion);
1405 3 : chanSpac*=inversion;
1406 3 : Double chanBW = abs(chanSpac);
1407 3 : if (!cabb_p && nchan==33) chanBW=chanBW*1.6; // roughly
1408 3 : if (!cabb_p && nchan==65) chanBW=chanBW*1.3; // guess
1409 :
1410 : // do birdie check here - reduce number of channels we will store
1411 : // for 33 channel data (drop half the channels + edge) and fix
1412 : // spw description
1413 3 : if (birdie_p && nchan == 33 && !cabb_p) {
1414 0 : Int bchan = birdChan(refFreq/1.e9, (Int)refChan, chanSpac/1.e9);
1415 0 : Int edge = 3 + (bchan+1)%2;
1416 0 : nchan = (nchan - 2*edge + 1)/2;
1417 0 : chanSpac*=2;
1418 0 : refChan = (refChan - edge)/2;
1419 : }
1420 6 : Vector<Double> chanFreq(nchan), resolution(nchan),width(nchan);
1421 4134 : for (Int ichan=0; ichan<nchan; ichan++)
1422 4131 : chanFreq(ichan)=refFreq+(ichan-refChan)*chanSpac;
1423 3 : msc_p->spectralWindow().chanFreq().put(n,chanFreq);
1424 3 : resolution=chanBW;
1425 3 : width=chanSpac;
1426 3 : msc_p->spectralWindow().resolution().put(n,resolution);
1427 3 : msc_p->spectralWindow().chanWidth().put(n, width);
1428 3 : msc_p->spectralWindow().effectiveBW().put(n, resolution);
1429 3 : msc_p->spectralWindow().totalBandwidth().put(n,
1430 3 : (nchan<33 ? abs(nchan*chanSpac) : doubles_.if_bw[ifNumber]));
1431 3 : msc_p->spectralWindow().numChan().put(n,nchan);
1432 :
1433 6 : Vector<Int> corrType(npol);
1434 6 : Matrix<Int> corrProduct(2,npol); corrProduct=0;
1435 3 : if (log) os_p << LogIO::NORMAL << " : Polarizations ";
1436 15 : for (Int i=0; i<npol; i++) {
1437 36 : corrType(i) = Stokes::type
1438 12 : (String(&names_.if_cstok[2*(i+ifNumber*MaxNPol)],2).before(trailing));
1439 : }
1440 6 : Vector<Int> tmp(npol); tmp=corrType;
1441 : // Sort the polarizations to standard order
1442 3 : GenSort<Int>::sort(corrType);
1443 3 : if (corrIndex_p.nrow()==0) {
1444 1 : corrIndex_p.resize(64,4);
1445 : }
1446 3 : if (Int(corrIndex_p.nrow())<=spWId_p) {
1447 0 : corrIndex_p.resize(nSpW_p,4,true);
1448 : }
1449 : // Get the sort indices to rearrange the data to standard order
1450 15 : for (Int i=0;i<npol;i++) {
1451 60 : for (Int j=0;j<npol;j++) {
1452 48 : if (corrType(j)==tmp(i)) corrIndex_p(spWId_p,i)=j;
1453 : }
1454 : }
1455 :
1456 : // Figure out the correlation products from the polarizations
1457 3 : corrProduct.resize(2,npol); corrProduct=0;
1458 15 : for (Int i=0; i<npol; i++) {
1459 12 : Stokes::StokesTypes s=Stokes::type(corrType(i));
1460 24 : Fallible<Int> receptor=Stokes::receptor1(s);
1461 12 : if (receptor.isValid()) corrProduct(0,i)=receptor;
1462 12 : receptor=Stokes::receptor2(s);
1463 12 : if (receptor.isValid()) corrProduct(1,i)=receptor;
1464 12 : if (i>0 && log) os_p << " , ";
1465 12 : if (log) os_p << Stokes::name(s)<< " - " << corrType(i);
1466 : }
1467 3 : if (log) os_p << LogIO::POST;
1468 :
1469 : // try to find matching pol row
1470 3 : Int nPolRow = atms_p.polarization().nrow();
1471 3 : Int polRow = -1;
1472 3 : for (Int i = 0; i< nPolRow; i++) {
1473 2 : if (msc_p->polarization().numCorr()(i)== Int(if_.if_nstok[ifNumber])) {
1474 2 : if (allEQ(msc_p->polarization().corrType()(i),corrType)) {
1475 2 : polRow=i;
1476 2 : break;
1477 : }
1478 : }
1479 : }
1480 : // add new pol id if needed
1481 3 : if (polRow==-1) {
1482 1 : atms_p.polarization().addRow();
1483 1 : polRow = nPolRow;
1484 1 : msc_p->polarization().numCorr().put(polRow,Int(if_.if_nstok[ifNumber]));
1485 1 : msc_p->polarization().corrType().put(polRow,corrType);
1486 1 : msc_p->polarization().corrProduct().put(polRow,corrProduct);
1487 1 : msc_p->polarization().flagRow().put(polRow, false);
1488 : }
1489 3 : atms_p.dataDescription().addRow();
1490 3 : msc_p->dataDescription().spectralWindowId().put(n, spWId_p);
1491 3 : msc_p->dataDescription().polarizationId().put(n, polRow);
1492 3 : msc_p->dataDescription().flagRow().put(n, false);
1493 :
1494 3 : msc_p->spectralWindow().ifConvChain().put(n,Int(if_.if_chain[ifNumber])-1);
1495 3 : if (!cabb_p) colSamplerBits.put(n,Int(if_.if_sampl[ifNumber]));
1496 :
1497 : // set up the TiledStorageManagers
1498 6 : Record values1, values2, values3, values3c;
1499 3 : values1.define("DATA_HYPERCUBE_ID",spWId_p);
1500 3 : values2.define("SIGMA_HYPERCUBE_ID",spWId_p);
1501 3 : values3.define("FLAG_HYPERCUBE_ID",spWId_p);
1502 3 : values3c.define("FLAG_CATEGORY_HYPERCUBE_ID",spWId_p);
1503 :
1504 6 : Record values4; values4.define("MODEL_HYPERCUBE_ID",spWId_p);
1505 6 : Record values5; values5.define("CORRECTED_HYPERCUBE_ID",spWId_p);
1506 : //Record values6; values6.define("IMAGING_WT_HYPERCUBE_ID",spWId_p);
1507 :
1508 3 : Int nChan=msc_p->spectralWindow().numChan()(spWId_p);
1509 3 : Int nCorr=msc_p->polarization().numCorr()(polRow);
1510 3 : Int nCat=3;
1511 :
1512 : // Choose an appropriate tileshape
1513 6 : IPosition dataShape(2,nCorr,nChan);
1514 3 : IPosition tileShape = MSTileLayout::tileShape(dataShape,obsType_p,"ATCA");
1515 3 : dataAccessor_p.addHypercube(IPosition(3,nCorr,nChan,0),
1516 : tileShape,values1);
1517 3 : sigmaAccessor_p.addHypercube(IPosition(2,nCorr,0),
1518 6 : IPosition(2,tileShape(0),tileShape(2)),
1519 : values2);
1520 3 : flagAccessor_p.addHypercube(IPosition(3,nCorr,nChan,0),
1521 : tileShape,values3);
1522 3 : flagCatAccessor_p.addHypercube(IPosition(4,nCorr,nChan,nCat,0),
1523 6 : IPosition(4,tileShape(0),tileShape(1),
1524 3 : 1,tileShape(2)),values3c);
1525 :
1526 : }
1527 : }
1528 234 : return spWId_p;
1529 : }
1530 :
1531 :
1532 1 : void ATCAFiller::checkObservation() {
1533 2 : const Regex trailing(" *$"); // trailing blanks
1534 : // Check if current observer already in OBSERVATION table
1535 2 : String observer;
1536 1 : obsId_p=-1;
1537 1 : for (uInt i=0; i<atms_p.observation().nrow(); i++) {
1538 0 : msc_p->observation().observer().get(i,observer);
1539 0 : if (String(names_.rp_observer,16).before(trailing)==observer) {
1540 0 : obsId_p=i;
1541 0 : break;
1542 : }
1543 : }
1544 1 : if (obsId_p<0) {
1545 : // not found, so add it
1546 1 : atms_p.observation().addRow();
1547 1 : obsId_p=atms_p.observation().nrow()-1;
1548 2 : msc_p->observation().observer().put(obsId_p,
1549 2 : String(names_.rp_observer,16).before(trailing));
1550 : // decode project from rpfits file name..
1551 2 : String project=rpfitsFiles_p(0).after(Regex(".*\\."));
1552 1 : if (project.contains(";")) project=project.before(";");
1553 1 : msc_p->observation().project().put(obsId_p,project);
1554 : }
1555 1 : }
1556 :
1557 77 : void ATCAFiller::checkField() {
1558 : // Check if current source already in Field Table
1559 : // Add it if not, set fieldId_p index for our Field Table
1560 : // For now, we have 1 source per field, we may handle mosaicing differently
1561 : // at some point.
1562 154 : const Regex trailing(" *$"); // trailing blanks
1563 : //0.1 arcsec tolerance for positional coincidence
1564 77 : const Double PosTol=0.1*C::arcsec;
1565 :
1566 77 : if (su_.su_found) {
1567 77 : fieldId_p=-1;
1568 77 : Bool seenSource = false;
1569 : Int numPol;
1570 154 : String name;
1571 154 : String su_name=String(&names_.su_name[sourceno*16],16).before(trailing);
1572 77 : for (Int i=0; i<nField_p; i++) {
1573 76 : msc_p->field().name().get(i,name);
1574 76 : msc_p->field().numPoly().get(i,numPol);
1575 :
1576 76 : if (su_name==name) {
1577 76 : IPosition shape(2, 2, numPol+1);
1578 76 : Matrix<Double> phaseDir(shape);
1579 76 : msc_p->field().phaseDir().get(i,phaseDir);
1580 76 : if (abs(phaseDir(0, 0)-doubles_.su_ra[sourceno])<PosTol) {
1581 76 : if (abs(phaseDir(1, 0)-doubles_.su_dec[sourceno])<PosTol) {
1582 76 : fieldId_p=i; // found it!
1583 : // now check if we've seen this field for this spectral window
1584 152 : Vector<Int> sourceIdCol = msc_p->source().sourceId().getColumn();
1585 152 : Vector<Int> spwIdCol = msc_p->source().spectralWindowId().getColumn();
1586 228 : Vector<Int> spwids = spwIdCol(sourceIdCol==i).getCompressedArray();
1587 76 : seenSource = (spwids.nelements()>0) && anyEQ(spwids,spWId_p);
1588 76 : break;
1589 : }
1590 : }
1591 : }
1592 : }
1593 :
1594 77 : if (fieldId_p<0 || !seenSource) { // i.e. not found, or not at this spwid
1595 6 : String src = String(&names_.su_name[sourceno*16],16);
1596 :
1597 3 : nsources_p++;
1598 3 : sources_p = src;
1599 :
1600 3 : Double epoch=mjd0_p+Double(proper_.pm_epoch);
1601 3 : Int numPol = 0;
1602 3 : if (abs(proper_.pm_ra)+abs(proper_.pm_dec) > 0) {
1603 0 : numPol = 1;
1604 : }
1605 6 : IPosition shape(2, 2, numPol+1);
1606 6 : Matrix<Double> dir(shape);
1607 : // convert proper motions from s & " per year to rad/s
1608 3 : const Double arcsecPerYear=C::arcsec/(365.25*C::day);
1609 3 : dir(0, 0)=doubles_.su_ra[sourceno];
1610 3 : dir(1, 0)=doubles_.su_dec[sourceno];
1611 3 : if (numPol>0) {
1612 0 : dir(0, 1)=proper_.pm_ra*15.*arcsecPerYear; // (15"/s)
1613 0 : dir(1, 1)=proper_.pm_dec*arcsecPerYear;
1614 : }
1615 :
1616 3 : if (fieldId_p<0) {
1617 1 : os_p << LogIO::DEBUGGING << "Found field:" << src << LogIO::POST;
1618 1 : if (numPol>0) os_p << LogIO::DEBUGGING << "Field:" << src << " has proper motion parameters"<<LogIO::POST;
1619 1 : fieldId_p=nField_p++;
1620 1 : atms_p.field().addRow();
1621 1 : Int nf=atms_p.field().nrow()-1;
1622 : // for now we have 1 field/source
1623 1 : msc_p->field().sourceId().put(nf,fieldId_p);
1624 1 : msc_p->field().name().put(nf,src.before(trailing));
1625 :
1626 1 : msc_p->field().phaseDir().put(nf,dir);
1627 1 : msc_p->field().delayDir().put(nf,dir);
1628 1 : msc_p->field().referenceDir().put(nf, dir);
1629 1 : msc_p->field().numPoly().put(nf, numPol);
1630 1 : msc_p->field().time().put(nf,epoch);
1631 3 : msc_p->field().code().put(nf,String(&names_.su_cal[sourceno*16],
1632 2 : 16).before(trailing));
1633 : }
1634 3 : dir(0, 0)=doubles_.su_pra[sourceno];
1635 3 : dir(1, 0)=doubles_.su_pdec[sourceno];
1636 6 : Vector<Double> srcdir(2);
1637 3 : srcdir(0)=doubles_.su_ra[sourceno];
1638 3 : srcdir(1)=doubles_.su_dec[sourceno];
1639 3 : atms_p.source().addRow();
1640 3 : Int ns=atms_p.source().nrow()-1;
1641 3 : msc_p->source().sourceId().put(ns,fieldId_p);
1642 3 : msc_p->source().name().put(ns,src.before(trailing));
1643 3 : msc_p->source().direction().put(ns,srcdir);
1644 6 : Vector<Double> rate(2);
1645 3 : rate(0)=proper_.pm_ra*15.*arcsecPerYear; // (15"/s)
1646 3 : rate(1)=proper_.pm_dec*arcsecPerYear;
1647 3 : msc_p->source().properMotion().put(ns,rate);
1648 3 : msc_p->source().time().put(ns,epoch);
1649 3 : msc_p->source().interval().put(ns,DBL_MAX);
1650 : // assume source is at infinity, specify as zero
1651 : //Vector<Double> pos(3); pos=0.;
1652 : //msc_p->source().position().put(ns,pos);
1653 :
1654 : //suRow["VelRefFrame"].put(Int(spect_.ivelref));
1655 : // need to specify reference frame for velocity as well
1656 : // vel1 may be velocity at channel 1 instead of ref freq..
1657 :
1658 3 : msc_p->source().spectralWindowId().put(ns,spWId_p);
1659 6 : Vector<Double> sysv(1),restFreq(1);
1660 3 : Vector<String> transition(1);
1661 3 : sysv = Double(doubles_.vel1);
1662 3 : msc_p->source().sysvel().put(ns, sysv);
1663 3 : msc_p->source().numLines().put(ns,1);
1664 3 : restFreq(0) = Double(doubles_.rfreq);
1665 3 : if (restFreq(0)<1.0) {
1666 : // fill in an appropriate default
1667 3 : Double freq = msc_p->spectralWindow().refFrequency()(spWId_p);
1668 3 : Double bw = msc_p->spectralWindow().totalBandwidth()(spWId_p);
1669 3 : if ( freq>1200.e6 && freq< 1450.e6 && bw <64.e6) {
1670 0 : restFreq(0)= 1420.4e6;
1671 0 : transition(0)="HI";
1672 : } else {
1673 3 : restFreq(0) = freq;
1674 3 : transition(0) = "";
1675 : }
1676 : }
1677 3 : msc_p->source().transition().put(ns,transition);
1678 3 : msc_p->source().restFrequency().put(ns,restFreq);
1679 :
1680 : // dummy fill
1681 3 : msc_p->source().calibrationGroup().put(ns,-1);
1682 : }
1683 :
1684 77 : Bool fillPointingTable=False; //Imaging fails with pointing table filled
1685 77 : if (fillPointingTable &&
1686 0 : (fieldId_p != prev_fieldId_p || newPointingCorr_p)) {
1687 0 : prev_fieldId_p = fieldId_p;
1688 0 : Double epoch=mjd0_p+Double(proper_.pm_epoch);
1689 0 : Int np=atms_p.pointing().nrow();
1690 0 : IPosition shape(2, 2, 1);
1691 0 : Matrix<Double> pointingDir(shape,0.0);
1692 0 : pointingDir(0, 0)=doubles_.su_pra[sourceno];
1693 0 : pointingDir(1, 0)=doubles_.su_pdec[sourceno];
1694 0 : Matrix<Double> pointingOffset(shape,0.0);
1695 0 : for (Int i=0; i<nAnt_p; i++) {
1696 0 : atms_p.pointing().addRow();
1697 0 : msc_p->pointing().antennaId().put(np+i, i+1);
1698 0 : if (newPointingCorr_p) {
1699 0 : pointingOffset(0, 0)=pointingCorr_p(0,i)*C::arcsec;
1700 0 : pointingOffset(1, 0)=pointingCorr_p(0,i)*C::arcsec;
1701 : }
1702 0 : if (i==0) { // ISM storage
1703 0 : msc_p->pointing().time().put(np,epoch);
1704 0 : msc_p->pointing().interval().put(np,DBL_MAX);
1705 0 : msc_p->pointing().numPoly().put(np, 0);
1706 0 : msc_p->pointing().direction().put(np,pointingDir);
1707 0 : msc_p->pointing().pointingOffset().put(np,pointingOffset);
1708 0 : msc_p->pointing().tracking().put(np,True);
1709 : } else {
1710 0 : if (newPointingCorr_p) {
1711 0 : msc_p->pointing().pointingOffset().put(np,pointingOffset);
1712 : }
1713 : }
1714 : }
1715 : }
1716 : }
1717 77 : }
1718 :
1719 :
1720 57 : void ATCAFiller::storeSysCal()
1721 : {
1722 : // Conversion factor for atmospheric pressure
1723 114 : Vector<String> const pressure_unit = msc_p->weather().pressure().keywordSet().asArrayString("QuantumUnits");
1724 57 : auto const pressure_conversion = Quantum<Double>(1.0, Unit("mbar")).getValue(Unit(pressure_unit[0]));
1725 :
1726 : // RPFITS SysCal table layout:
1727 : // sc_.sc_ant = 7 (1-6 is antenna 1-6 syscal data, 7th has ant=0 weather data)
1728 : // sc_cal(q,if,ant) (in practice sc_.sc_if is always 1 since ~1999)
1729 : // q=0 : ant
1730 : // 1 : if
1731 : // 2 : XYPhase (deg)
1732 : // 3 : tsysX (sqrt(tsys*10))
1733 : // 4 : tsysY
1734 : // 5-7: samp-stats X
1735 : // 8-10:samp-stats Y
1736 : // 11 : Parallactic Angle (deg)
1737 : // 12 : Flag
1738 : // 13 : XYAmp (Jy)
1739 : // 15: Tracking error Max (arcsec)
1740 : // 16: Tracking error RMS (arcsec)
1741 57 : Int last_if_no=-2; // if_no can come out to be -1 if syscal rec is blank
1742 57 : Bool skip=false;
1743 57 : sourceno=sc_.sc_srcno-1; // 0-based source number for this syscal record
1744 57 : Int scq=sc_.sc_q, scif=sc_.sc_if;
1745 114 : for (Int i=0; i<scif; i++) {
1746 456 : for (Int ant=0; ant<sc_.sc_ant; ant++) {
1747 399 : if (Int(sc_.sc_cal[scq*(i+scif*ant)+0])==0) {
1748 : // special syscal record with antenna==0
1749 : // field 7-12 contain weather data
1750 57 : if (sc_.sc_ut!=lastWeatherUT_p) {
1751 19 : lastWeatherUT_p = sc_.sc_ut;
1752 19 : Int nAnt = max(1,sc_.sc_ant-1);
1753 19 : Int row=atms_p.weather().nrow();
1754 19 : atms_p.weather().addRow(nAnt);
1755 19 : Double time=mjd0_p+Double(sc_.sc_ut);
1756 133 : for (Int iAnt=0; iAnt<nAnt; iAnt++,row++) {
1757 114 : msc_p->weather().antennaId().put(row,iAnt);
1758 114 : msc_p->weather().time().put(row,time);
1759 114 : msc_p->weather().interval().put(row,Double(param_.intbase));
1760 114 : msc_p->weather().temperature().put(row,
1761 114 : Double(sc_.sc_cal[scq*(i+scif*ant)+1])+273.15); // C to K
1762 114 : msc_p->weather().pressure().put(row,
1763 114 : Double(sc_.sc_cal[scq*(i+scif*ant)+2])*pressure_conversion); // mBar to Pa/hPa
1764 114 : msc_p->weather().relHumidity().put(row,
1765 114 : Double(sc_.sc_cal[scq*(i+scif*ant)+3]));
1766 114 : msc_p->weather().windSpeed().put(row,
1767 114 : Double(sc_.sc_cal[scq*(i+scif*ant)+4])/3.6); // km/s to m/s
1768 114 : msc_p->weather().windDirection().put(row,
1769 114 : Double(sc_.sc_cal[scq*(i+scif*ant)+5])*C::pi/180.0); // deg to rad
1770 114 : msc_p->weather().temperatureFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1771 114 : msc_p->weather().pressureFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1772 114 : msc_p->weather().relHumidityFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1773 114 : msc_p->weather().windSpeedFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1774 114 : msc_p->weather().windDirectionFlag().put(row,Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1775 114 : colWeatherRainGauge.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+7]));
1776 114 : colWeatherSeeMonPhase.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+8]));
1777 114 : colWeatherSeeMonRMS.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+9])/1000);
1778 114 : colWeatherSeeMonFlag.put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+10]));
1779 : }
1780 : }
1781 57 : continue;
1782 : }
1783 342 : if_no=Int(sc_.sc_cal[scq*(i+scif*ant)+1])-1; // make 0-based
1784 342 : if (if_no!=last_if_no) { // check if we want this one
1785 57 : last_if_no=if_no;
1786 57 : skip=Bool(if_no<0 || if_no>=if_.n_if || !selected(if_no));
1787 57 : if (!skip) {
1788 57 : checkSpW(if_no); // sets spWId_p corresponding to this if_no
1789 57 : checkField(); // sets fieldId_p
1790 : }
1791 : }
1792 342 : if (!skip) {
1793 342 : atms_p.sysCal().addRow();
1794 342 : Int n=atms_p.sysCal().nrow()-1;
1795 342 : Double time=mjd0_p+Double(sc_.sc_ut);
1796 342 : msc_p->sysCal().time().put(n,time);
1797 342 : msc_p->sysCal().antennaId()
1798 342 : .put(n,Int(sc_.sc_cal[scq*(i+scif*ant)+0]-1)); // make 0-based
1799 342 : msc_p->sysCal().feedId().put(n,0);
1800 342 : msc_p->sysCal().interval().put(n,Double(param_.intbase));
1801 342 : msc_p->sysCal().spectralWindowId().put(n,spWId_p);
1802 :
1803 342 : msc_p->sysCal().phaseDiff().put(n,sc_.sc_cal[scq*(i+scif*ant)+2]);
1804 :
1805 684 : Vector<Float> tSys(2);
1806 : // convert from sqrt(tsys*10) to true tsys
1807 342 : tSys(0)=square(sc_.sc_cal[scq*(i+scif*ant)+3])/10.;
1808 342 : tSys(1)=square(sc_.sc_cal[scq*(i+scif*ant)+4])/10.;
1809 :
1810 342 : msc_p->sysCal().tsys().put(n,tSys);
1811 :
1812 684 : Vector<Float> a(2),b(2),c(2);
1813 1026 : for (Int j=0; j<2; j++) {
1814 684 : a(j)=sc_.sc_cal[scq*(i+scif*ant)+5+j*3];
1815 684 : b(j)=sc_.sc_cal[scq*(i+scif*ant)+6+j*3];
1816 684 : c(j)=sc_.sc_cal[scq*(i+scif*ant)+7+j*3];
1817 : }
1818 342 : if (cabb_p) {
1819 342 : colGTP.put(n,a);
1820 342 : colSDO.put(n,b);
1821 342 : colCalJy.put(n,c);
1822 : } else {
1823 0 : colSamplerStatsNeg.put(n,a);
1824 0 : colSamplerStatsZero.put(n,b);
1825 0 : colSamplerStatsPos.put(n,c);
1826 : }
1827 : //cParAngle.put(n,sc_.sc_cal[scq*(i+scif*ant)+11]);
1828 342 : msc_p->sysCal().phaseDiffFlag().
1829 342 : put(n,(sc_.sc_cal[scq*(i+scif*ant)+12]!=0));
1830 342 : msc_p->sysCal().tsysFlag().
1831 342 : put(n,(sc_.sc_cal[scq*(i+scif*ant)+12]!=0));
1832 342 : colXYAmplitude.put(n,sc_.sc_cal[scq*(i+scif*ant)+13]);
1833 342 : msc_p->sysCal().tcalFlag().put(n,true);
1834 342 : msc_p->sysCal().trxFlag().put(n,true);
1835 342 : colTrackErrMax.put(n,Float(sc_.sc_cal[scq*(i+scif*ant)+14]));
1836 342 : colTrackErrRMS.put(n,Float(sc_.sc_cal[scq*(i+scif*ant)+15]));
1837 : }
1838 : }
1839 : }
1840 57 : }
1841 :
1842 :
1843 0 : void ATCAFiller::reweight()
1844 : {
1845 0 : Int npol =if_.if_nstok[if_no];
1846 0 : Int nfreq=33;
1847 0 : Int n=2*nfreq-2;
1848 0 : FFTServer<Float, Complex> server;
1849 :
1850 0 : Vector<Complex> cscr(33);
1851 0 : Vector<Float> rscr(64);
1852 :
1853 : static Float wts[64]={ 1.000000,
1854 : 1.028847, 1.0526778, 1.0711329, 1.083963, 1.0909837, 1.0921189,
1855 : 1.0873495, 1.0767593, 1.0605253, 1.0389025, 1.0122122,
1856 : 0.9808576, 0.9453095, 0.9060848, 0.8637353, 0.8188493,
1857 : 0.7720414, 0.7239398, 0.6751758, 0.6263670, 0.5781130,
1858 : 0.5310048, 0.4856412, 0.4426092, 0.4025391, 0.3661798,
1859 : 0.3346139, 0.3096083, 0.2951138, 0.3024814, 0.6209788, 1.000000,
1860 : 0.6209788, 0.3024814, 0.2951138, 0.3096083, 0.3346139,
1861 : 0.3661798, 0.4025391, 0.4426092, 0.4856412, 0.5310048,
1862 : 0.5781130, 0.6263670, 0.6751758, 0.7239398, 0.7720414,
1863 : 0.8188493, 0.8637353, 0.9060848, 0.9453095, 0.9808576,
1864 : 1.0122122, 1.0389025, 1.0605253, 1.0767593, 1.0873495, 1.0921189,
1865 : 1.0909837, 1.083963, 1.0711329, 1.0526778, 1.028847 };
1866 :
1867 0 : for (Int p=0; p<npol; p++) {
1868 0 : for (Int i=0; i<nfreq; i++)
1869 0 : cscr(i) = Complex(vis[0+2*(p+i*npol)],vis[1+2*(p+i*npol)]);
1870 0 : server.fft0(rscr, cscr);
1871 0 : for (Int i=0; i<n; i++) rscr(i) = rscr(i)*wts[i];
1872 0 : server.fft0(cscr, rscr);
1873 0 : for (Int i=0; i<nfreq; i++) {
1874 0 : vis[0+2*(p+i*npol)] = real(cscr(i));
1875 0 : vis[1+2*(p+i*npol)] = imag(cscr(i));
1876 : }
1877 : }
1878 0 : }
1879 :
1880 0 : Int ATCAFiller::birdChan(Double refFreq, Int refChan, Double chanSpac)
1881 : {
1882 0 : Double flo = 0.128* myround(refFreq/0.128);
1883 0 : Int chan = refChan + myround((flo - refFreq)/chanSpac);
1884 0 : if(chan <= 0) {
1885 : // os_p << LogIO::NORMAL << "CHAN " << chan << " ";
1886 0 : chan = chan + myround(0.128/fabs(chanSpac));
1887 : }
1888 : //os_p << LogIO::NORMAL << "birdie channel = " << chan << " refFreq="<<refFreq << " flo="<<flo <<"\n";
1889 0 : return chan;
1890 : }
1891 :
1892 855 : void ATCAFiller::storeData()
1893 : {
1894 855 : const double MJD01Jul2004 = 2453187.5; // 12mm receiver xy inversion end
1895 855 : const double MJD18Oct2007 = 2454390.5; // 3mm CA02 xyphase valid
1896 :
1897 855 : atms_p.addRow();
1898 :
1899 1710 : const RecordFieldId rfid1("DATA_HYPERCUBE_ID");
1900 1710 : const RecordFieldId rfid2("SIGMA_HYPERCUBE_ID");
1901 1710 : const RecordFieldId rfid3("FLAG_HYPERCUBE_ID");
1902 1710 : const RecordFieldId rfid3c("FLAG_CATEGORY_HYPERCUBE_ID");
1903 1710 : Record values1, values2, values3, values3c;
1904 855 : values1.define(rfid1,spWId_p);
1905 855 : values2.define(rfid2,spWId_p);
1906 855 : values3.define(rfid3,spWId_p);
1907 855 : values3c.define(rfid3c,spWId_p);
1908 :
1909 855 : dataAccessor_p.extendHypercube(1,values1);
1910 855 : sigmaAccessor_p.extendHypercube(1,values2);
1911 855 : flagAccessor_p.extendHypercube(1,values3);
1912 855 : flagCatAccessor_p.extendHypercube(1,values3c);
1913 :
1914 1710 : Record values4, values5, values6;
1915 :
1916 855 : Int n=atms_p.nrow()-1;
1917 855 : if (n==0) gotSysCalId_p=false;
1918 :
1919 855 : Int npol =if_.if_nstok[if_no];
1920 855 : Int nfreq=if_.if_nfreq[if_no];
1921 1710 : Regex trailing(" *$"); // trailing blanks
1922 1710 : String instrument = String(names_.instrument,16).before(trailing);
1923 855 : Bool atca = (instrument=="ATCA");
1924 855 : Bool atlba = (instrument=="ATLBA");
1925 :
1926 : // fill in the easy items
1927 : // make antenna numbers 0-based
1928 855 : Int ant1=Int(baseline)/256-1, ant2=Int(baseline)%256-1;
1929 855 : msc_p->antenna1().put(n,ant1);
1930 855 : msc_p->antenna2().put(n,ant2);
1931 855 : Bool flagData=flg;
1932 855 : if (flagData) { flagCount_p(ONLINE)++;}
1933 855 : if (autoFlag_p && flagScanType_p) { flagData=true; flagCount_p(SCANTYPE)++;}
1934 855 : Double exposure=Double(param_.intbase);
1935 : // averaging = number of integration periods averaged by correlator: 1,2 or 3
1936 855 : Int averaging = 1;
1937 855 : if (param_.intime>0) averaging = Int(exposure/param_.intime)+1;
1938 855 : Double interval = averaging*Double(param_.intime);
1939 : // check for old data with intbase set to zero
1940 855 : Double blank=51.e-3; // standard 51 ms blank time at end of integration
1941 855 : if (exposure<0.001) exposure = interval-blank;
1942 855 : if (interval==0.0) interval = exposure;
1943 :
1944 : // Is binning mode active?
1945 855 : Int nBin = Int(floor(param_.intime/exposure+0.01));
1946 855 : if (nBin<4) nBin=1; // short exposure due to long hold period..
1947 855 : msc_p->dataDescId().put(n, spWId_p);
1948 855 : if (ut!=lastUT_p) { // the ISM columns that don't change within integration
1949 19 : msc_p->arrayId().put(n,0);
1950 19 : msc_p->exposure().put(n,exposure);
1951 19 : msc_p->feed1().put(n,0);
1952 19 : msc_p->feed2().put(n,0);
1953 19 : msc_p->fieldId().put(n,fieldId_p);
1954 19 : msc_p->interval().put(n,interval);
1955 19 : msc_p->observationId().put(n,obsId_p);
1956 19 : msc_p->processorId().put(n,-1);
1957 19 : msc_p->scanNumber().put(n,scanNo_p);
1958 19 : msc_p->stateId().put(n,-1);
1959 19 : Double mjd=mjd0_p+ut;
1960 : // try to figure out what the midpoint of the integration interval was
1961 19 : if (nBin>1 || !atca) {
1962 19 : msc_p->time().put(n,mjd);
1963 : } else {
1964 0 : msc_p->time().put(n,floor((mjd+exposure/2+blank*(averaging+1)/2
1965 0 : -interval/2)*1000+0.5)/1000);
1966 : }
1967 : // time centroid is the centroid of the exposure window
1968 : // it is the time used for the uvw & phase calculation
1969 : // [note that cacor corrects uvw for change in centroid from caobs value
1970 : // due to blank, hold and averaging, but not binning]
1971 19 : msc_p->timeCentroid().put(n,mjd);
1972 : }
1973 855 : Int pulsarBin = max(0,bin-1);// make zero-based (but catch unset value)
1974 855 : msc_p->pulsarBin().put(n,pulsarBin);
1975 855 : if (hires_p && nBin>1) {
1976 : // in hires mode we adjust timeCentroid to match the bin centers, but
1977 : // uvw's still refer to center of interval (now given by time column)
1978 : // Time column will still be in time order, but TimeCentroid is not.
1979 0 : Double midTime = msc_p->time()(n);
1980 0 : msc_p->timeCentroid().put(n,midTime+(pulsarBin-nBin/2+0.5)*exposure/averaging);
1981 : }
1982 1710 : Vector<Double> uvw(3); uvw(0)=u; uvw(1)=v; uvw(2)=w;
1983 855 : msc_p->uvw().put(n,uvw);
1984 :
1985 : // use exposure time as weight
1986 1710 : Vector<Float> Weight(npol); Weight=exposure;
1987 855 : msc_p->weight().put(n,Weight);
1988 : //# Vector<Float> weightSpectrum(nfreq); weightSpectrum=exposure;
1989 : //# msc_p->weightSpectrum().put(n,weightSpectrum);
1990 :
1991 : // Find the indices into the sysCal table for this row
1992 855 : if (!gotSysCalId_p || ut!=lastUT_p || spWId_p!=lastSpWId_p) {
1993 855 : lastUT_p=ut;
1994 855 : lastSpWId_p=spWId_p;
1995 : // we need to get the new syscal info
1996 : // search backwards from last syscal row
1997 855 : Bool done=false;
1998 855 : Int nsc=atms_p.sysCal().nrow();
1999 : //# os_p << LogIO::NORMAL << " Current number of syscal rows="<<nsc<<LogIO::POST;
2000 855 : sysCalId_p.resize(nAnt_p); sysCalId_p=-1;
2001 : // search backwards, as it should be a recent entry
2002 5985 : for (Int i=nsc-1; (i>=0 && !done); i--) {
2003 10260 : if (nearAbs(msc_p->sysCal().time()(i),mjd0_p+ut,0.5) &&
2004 5130 : msc_p->sysCal().spectralWindowId()(i)==spWId_p) {
2005 5130 : Int ant=msc_p->sysCal().antennaId()(i);
2006 5130 : if (ant>=0 && ant<nAnt_p) { //valid antenna
2007 5130 : sysCalId_p(ant)=i;
2008 5130 : done=(!(anyEQ(sysCalId_p,-1)));
2009 : }
2010 : }
2011 : }
2012 855 : if (!done && atca) {
2013 0 : errCount_p++;
2014 0 : if (errCount_p<3) os_p << LogIO::WARN <<"Warning: some syscal info is missing"<< LogIO::POST;
2015 : }
2016 : }
2017 : // set a sysCalId in the main table to point to
2018 : // the sysCal subtable row - missing points get a -1
2019 855 : colSysCalIdAnt1.put(n,sysCalId_p(ant1));
2020 855 : colSysCalIdAnt2.put(n,sysCalId_p(ant2));
2021 :
2022 : // Check for bad sampler stats now that we've found the syscal data
2023 855 : if (!flagData && autoFlag_p && (!atlba && !cabb_p && samplerFlag(n))) {
2024 0 : flagData=true; flagCount_p(SYSCAL)++;
2025 : }
2026 855 : flagCount_p(COUNT)++;
2027 855 : if (flagData) flagCount_p(FLAG)++;
2028 855 : msc_p->flagRow().put(n,flagData);
2029 :
2030 :
2031 : // flags and data
2032 855 : const Int nCat = 3; // three initial categories
2033 : // define the categories
2034 1710 : Vector<String> cat(nCat);
2035 855 : cat(0)="FLAG_CMD";
2036 855 : cat(1)="ORIGINAL";
2037 855 : cat(2)="USER";
2038 855 : msc_p->flagCategory().rwKeywordSet().define("CATEGORY",cat);
2039 :
2040 : // Gibbs reweighting
2041 855 : if (nfreq == 33) {
2042 285 : if (reweight_p) reweight();
2043 : }
2044 :
2045 : // do birdie check here - reduce number of channels we will store
2046 : // for 33 channel data (drop half the channels + edge)
2047 855 : Double refFreq=doubles_.if_freq[if_no];
2048 855 : Double refChan=doubles_.if_ref[if_no]-1; // make zero based
2049 855 : Double chanSpac = doubles_.if_bw[if_no]/max(1,nfreq-1);
2050 855 : Int inversion=if_.if_invert[if_no];
2051 855 : chanSpac*=inversion;
2052 855 : Int bchan = -1;
2053 855 : Int edge = 0;
2054 855 : if (birdie_p && atca) {
2055 0 : bchan = birdChan(refFreq/1.e9, (Int)refChan, chanSpac/1.e9);
2056 0 : if (nfreq == 33) { // ATCA continuum mode - flag edge + every other channel
2057 0 : edge = 3 + (bchan+1)%2;
2058 0 : nfreq = (nfreq - 2*edge + 1)/2;
2059 0 : chanSpac*=2;
2060 0 : refChan = (refChan - edge)/2;
2061 0 : bchan = (bchan -edge)/2;
2062 : }
2063 : }
2064 855 : Bool band12mm = (refFreq>13.e9 && refFreq<28.e9);
2065 855 : Bool band3mm = (refFreq>75.e9);
2066 :
2067 1710 : Matrix<Complex> VIS(npol,nfreq);
2068 1710 : Cube<Bool> flagCat(npol,nfreq,nCat,false);
2069 1710 : Matrix<Bool> flags= flagCat.xyPlane(1); // references flagCat's storage
2070 :
2071 :
2072 855 : if (birdie_p) {
2073 0 : if (bchan >= 0 && bchan < nfreq) {
2074 0 : for (Int i=0; i<npol; i++) {
2075 0 : flags(i, bchan) = true;
2076 : }
2077 : }
2078 : }
2079 :
2080 : // Some corrections as specified by Bob Sault(2000/09/07):
2081 : //Conjugate all data if the IF_INVERT switch is negative
2082 : //Negate the XY and YX correlations (both real and imag parts) for all
2083 : // except 12mm observations.
2084 : //The xyphase measurement should also be multiplied by the IF_INVERT switch
2085 : //sign. The xyphase measurement is the phase of a correlation. So to
2086 : //correct the data with the xyphase measurement, you will want to multiply
2087 : //by gains whose phase has the opposite sign to the measurements.
2088 : //Note any Gibbs reweighting needs to be done before any xyphase correction
2089 : //is done.
2090 :
2091 :
2092 : // correct for inversion - conjugate data; skip birdie channels
2093 : // resort pols
2094 1178190 : for (Int j=0; j<nfreq; j++) {
2095 1177335 : Int k = (edge>0 ? (2*j+edge) : j);
2096 5886675 : for (Int i=0; i<npol; i++) {
2097 4709340 : Int ipol = corrIndex_p(spWId_p,i);
2098 4709340 : VIS(ipol,j)=Complex(vis[0+2*(i+k*npol)],inversion*vis[1+2*(i+k*npol)]);
2099 : }
2100 : }
2101 : // Flag NaNs
2102 2565 : Matrix<Bool> nanFlags = isNaN(real(VIS));
2103 855 : nanFlags |= isNaN(imag(VIS));
2104 855 : flags |= nanFlags;
2105 :
2106 : // Flag CABB birdies and RFI
2107 855 : if (cabb_p) rfiFlag(flags);
2108 :
2109 : // correct for xy phase - multiply y correlation with opposite of xyphase
2110 855 : Int id1 = sysCalId_p(ant1);
2111 855 : Complex gain1 = (id1!=-1 && !msc_p->sysCal().phaseDiffFlag()(id1) ?
2112 540 : exp(Complex(0,1)*msc_p->sysCal().phaseDiff()(id1)*inversion):
2113 855 : Complex(1,0));
2114 : // negate xyphase for second antenna since the correlation product uses conjugate for ant2
2115 855 : Int id2 = sysCalId_p(ant2);
2116 855 : Complex gain2 = (id2!=-1 && !msc_p->sysCal().phaseDiffFlag()(id2) ?
2117 540 : exp(-Complex(0,1)*msc_p->sysCal().phaseDiff()(id2)*inversion):
2118 855 : Complex(1,0));
2119 855 : Int polId = msc_p->dataDescription().polarizationId()(spWId_p);
2120 1710 : Vector<Int> corrType = msc_p->polarization().corrType()(polId);
2121 : // 3mm receiver on CA02 got valid xyphase on 18Oct2007
2122 855 : if (noxycorr_p || (band3mm && mjd0_p<MJD18Oct2007)) { gain1=gain2=1; }
2123 855 : if (band3mm && mjd0_p>=MJD18Oct2007) {
2124 0 : if (ant1!=1) gain1 = 1;
2125 0 : if (ant2!=1) gain2 = 1;
2126 : }
2127 : // "-" because XY and YX corr need to be negated, except for 12mm until 2004
2128 855 : if (!band12mm || mjd0_p>MJD01Jul2004) { gain1=-gain1; gain2=-gain2;}
2129 4275 : for (Int i=0; i<npol; i++) {
2130 3420 : switch (corrType(i)) {
2131 855 : case Stokes::XX :
2132 855 : break;
2133 855 : case Stokes::XY :
2134 : {
2135 1710 : Vector<Complex> tmp(VIS.row(i));
2136 855 : tmp*=gain2;
2137 : }
2138 855 : break;
2139 855 : case Stokes::YX :
2140 : {
2141 1710 : Vector<Complex> tmp(VIS.row(i));
2142 855 : tmp*=gain1;
2143 : }
2144 855 : break;
2145 855 : case Stokes::YY :
2146 : {
2147 855 : Vector<Complex> tmp(VIS.row(i));
2148 855 : tmp*=(gain1*gain2);
2149 : }
2150 855 : break;
2151 0 : default:
2152 0 : break;
2153 : }
2154 : }
2155 :
2156 855 : msc_p->data().put(n,VIS);
2157 : // CASA wants all flags set if flagRow is set
2158 855 : if (flagData) flags=flagData;
2159 855 : msc_p->flag().put(n,flags);
2160 855 : msc_p->flagCategory().put(n,flagCat);
2161 :
2162 : // now calculate the sigma for this row
2163 : //
2164 1710 : Vector<Double> chnbw;
2165 855 : msc_p->spectralWindow().resolution().get(spWId_p,chnbw);
2166 1710 : Vector<Float> tsys1,tsys2;
2167 : //# os_p << LogIO::NORMAL << "Looking up syscal info at rows: "<<sysCalId_p(ant1)<<" and "<<
2168 : //# sysCalId_p(ant2)<< LogIO::POST;
2169 :
2170 855 : if (sysCalId_p(ant1)!=-1 && !msc_p->sysCal().tsysFlag()(sysCalId_p(ant1)))
2171 540 : msc_p->sysCal().tsys().get(sysCalId_p(ant1),tsys1);
2172 315 : else { tsys1.resize(nAnt_p); tsys1=50.;}
2173 855 : if (sysCalId_p(ant2)!=-1 && !msc_p->sysCal().tsysFlag()(sysCalId_p(ant2)))
2174 540 : msc_p->sysCal().tsys().get(sysCalId_p(ant2),tsys2);
2175 315 : else { tsys2.resize(nAnt_p); tsys2=50.;}
2176 1710 : Vector<Float> sigma(npol); sigma=0.;
2177 1710 : Matrix<Int> corrProduct(2,npol);
2178 855 : msc_p->polarization().corrProduct().get(polId,corrProduct);
2179 :
2180 : // sigma=sqrt(Tx1*Tx2)/sqrt(chnbw*intTime)*JyPerK;
2181 855 : Float JyPerK=13.; // guess for ATCA dishes at 3-20cm
2182 855 : Float factor=sqrt(chnbw(0)*exposure)/JyPerK;
2183 4275 : for (Int pol=0; pol<npol; pol++) {
2184 3420 : Float tsysAv=tsys1(corrProduct(0,pol))*tsys2(corrProduct(1,pol));
2185 3420 : if (tsysAv>=0 && factor>0 ) sigma(pol)=sqrt(tsysAv)/factor;
2186 : }
2187 855 : msc_p->sigma().put(n,sigma);
2188 :
2189 : // determine shadowing & flag data with shadowed antennas
2190 855 : if (shadow_p>0) shadow(n);
2191 :
2192 855 : }
2193 :
2194 855 : void ATCAFiller::rfiFlag(Matrix<Bool> & flags) {
2195 : // CABB birdies
2196 : // 1 MHz mode
2197 855 : const int nb1=11;
2198 : static const int b1[nb1] = {640,256,768,1408,1280,1920,1792,1176,156,128,1152};
2199 : // 64 MHz mode
2200 855 : const int nb2=3;
2201 : static const int b2[nb2] = {8,16,24};
2202 :
2203 : // Get details of spectrum
2204 855 : Int nfreq=if_.if_nfreq[if_no];
2205 855 : Double bw=doubles_.if_bw[if_no];
2206 855 : Int chn=0;
2207 : // Flag birdies and edge channels
2208 855 : if (bw>2.e9) {
2209 570 : chn = nfreq*edge_p/200;
2210 570 : if (nfreq==2049) {
2211 : // CABB 2049 chan * 1 MHz continuum mode
2212 3420 : for (Int i=0; i<nb1; i++) {
2213 3135 : flags.column(b1[i])=true;
2214 : }
2215 285 : } else if (nfreq==33) {
2216 : // CABB 33 chan * 64 MHz continuum mode
2217 1140 : for (Int i=0; i<nb2; i++) {
2218 855 : flags.column(b2[i])=true;
2219 : }
2220 : }
2221 285 : } else if (bw<1.e9 && nfreq>=2049) {
2222 : // CABB zoom mode, 2049 or more channels (combined zooms)
2223 285 : chn = 2049*edge_p/200;
2224 : }
2225 47310 : for (Int i=0; i<chn; i++) flags.column(i)=true;
2226 47310 : for (Int i=nfreq-chn; i<nfreq; i++) flags.column(i)=true;
2227 855 : }
2228 :
2229 1 : void ATCAFiller::fillFeedTable() {
2230 : // At present there is always only a single feed per antenna (at any
2231 : // given spectralwindow).
2232 1 : Int nAnt=atms_p.antenna().nrow();
2233 : // Int nSpW=atms_p.spectralWindow().nrow();
2234 : // Only X and Y receptors available
2235 2 : Vector<String> rec_type(2); rec_type(0)="X"; rec_type(1)="Y";
2236 2 : Matrix<Complex> polResponse(2,2);
2237 1 : polResponse=0.; polResponse(0,0)=polResponse(1,1)=1.;
2238 2 : Matrix<Double> offset(2,2); offset=0.;
2239 2 : Vector<Double> position(3); position=0.;
2240 : // X feed is at 45 degrees, Y at 135 degrees for all except 7mm
2241 2 : Vector<Double> receptorAngle(2);
2242 1 : receptorAngle(0)=45*C::degree;
2243 1 : receptorAngle(1)=135*C::degree;
2244 : // Single entry, so we use the first spectral window
2245 1 : Double f = msc_p->spectralWindow().refFrequency()(0);
2246 1 : if (f>30e9 && f<50e9) receptorAngle+=90*C::degree;
2247 :
2248 : // fill the feed table
2249 1 : Int row=-1;
2250 : // in principle we should have a separate entry for each SpectralWindow
2251 : // but at present all this is 'dummy' filled
2252 : // for (Int spw=0; spw<nSpW; spw++) {
2253 7 : for (Int ant=0; ant<nAnt; ant++) {
2254 6 : atms_p.feed().addRow(); row++;
2255 6 : msc_p->feed().antennaId().put(row,ant);
2256 6 : msc_p->feed().beamId().put(row,-1);
2257 6 : msc_p->feed().feedId().put(row,0);
2258 6 : msc_p->feed().interval().put(row,DBL_MAX);
2259 6 : msc_p->feed().phasedFeedId().put(row,-1);
2260 6 : msc_p->feed().spectralWindowId().put(row,-1); //spw);
2261 6 : msc_p->feed().time().put(row,0.);
2262 6 : msc_p->feed().numReceptors().put(row,2);
2263 6 : msc_p->feed().beamOffset().put(row,offset);
2264 6 : msc_p->feed().polarizationType().put(row,rec_type);
2265 6 : msc_p->feed().polResponse().put(row,polResponse);
2266 6 : msc_p->feed().position().put(row,position);
2267 6 : msc_p->feed().receptorAngle().put(row,receptorAngle);
2268 : }
2269 1 : }
2270 :
2271 1 : void ATCAFiller::fillObservationTable()
2272 : {
2273 2 : Vector<Double> tim = msc_p->time().getColumn();
2274 2 : Vector<Int> obsid = msc_p->observationId().getColumn();
2275 : Int startInd;
2276 : Int endInd;
2277 :
2278 1 : Int nObs=atms_p.observation().nrow();
2279 2 : Vector<Double> vt(2);
2280 2 : for (Int obs=0; obs<nObs; obs++) {
2281 : // fill time range
2282 1 : startInd = 0;
2283 1 : endInd = atms_p.nrow()-1;
2284 1 : for (uInt i=0; i<atms_p.nrow(); i++) {
2285 1 : if (obsid(i) == obs) { startInd = i; break; }
2286 : }
2287 856 : for (uInt i=startInd; i<atms_p.nrow(); i++) {
2288 855 : if (obsid(i) > obs) { endInd = i-1; break; }
2289 : }
2290 1 : vt(0) = tim(startInd);
2291 1 : vt(1) = tim(endInd);
2292 1 : msc_p->observation().timeRange().put(obs, vt);
2293 1 : msc_p->observation().releaseDate().put(obs,vt(1)+18*30.5*86400);
2294 : // telescope name
2295 1 : msc_p->observation().telescopeName().put(obs, "ATCA");
2296 : }
2297 1 : }
2298 :
2299 1 : void ATCAFiller::flush()
2300 : {
2301 1 : atms_p.flush();
2302 1 : atms_p.antenna().flush();
2303 1 : atms_p.dataDescription().flush();
2304 1 : atms_p.feed().flush();
2305 1 : atms_p.field().flush();
2306 1 : atms_p.observation().flush();
2307 1 : atms_p.history().flush();
2308 1 : atms_p.pointing().flush();
2309 1 : atms_p.polarization().flush();
2310 1 : atms_p.source().flush();
2311 1 : atms_p.spectralWindow().flush();
2312 1 : atms_p.sysCal().flush();
2313 1 : atms_p.weather().flush();
2314 1 : }
2315 :
2316 1 : void ATCAFiller::unlock()
2317 : {
2318 1 : atms_p.unlock();
2319 1 : atms_p.antenna().unlock();
2320 1 : atms_p.dataDescription().unlock();
2321 1 : atms_p.feed().unlock();
2322 1 : atms_p.field().unlock();
2323 1 : atms_p.observation().unlock();
2324 1 : atms_p.history().unlock();
2325 1 : atms_p.pointing().unlock();
2326 1 : atms_p.polarization().unlock();
2327 1 : atms_p.source().unlock();
2328 1 : atms_p.spectralWindow().unlock();
2329 1 : atms_p.sysCal().unlock();
2330 1 : atms_p.weather().unlock();
2331 1 : }
2332 :
2333 1541 : Bool ATCAFiller::selected(Int ifNum)
2334 : {
2335 1541 : Bool select=true;
2336 : //Check for rotten eggs
2337 1541 : if (if_.if_nfreq[ifNum]==0)
2338 0 : select=false;
2339 : // check if we want this frequency
2340 1541 : if (spws_p.nelements()>0 && spws_p(0)>=0 && !anyEQ(spws_p,ifNum)) {
2341 0 : select=false;
2342 : } else {
2343 : // check if we want this frequency
2344 1541 : if (lowFreq_p>0 && (doubles_.if_freq[ifNum]-lowFreq_p)<0)
2345 0 : select=false;
2346 : else {
2347 1541 : if (highFreq_p>0 && (highFreq_p-doubles_.if_freq[ifNum])<0)
2348 0 : select=false;
2349 : else {
2350 : // check if we want this bandwidth on the first IF
2351 1541 : if (bandWidth1_p>0 &&
2352 0 : (bandWidth1_p != myround(doubles_.if_bw[ifNum]/1.e6)))
2353 0 : select=false;
2354 : else {
2355 : // check if we want this number of channels on the first IF
2356 1541 : if (numChan1_p>0 && numChan1_p != if_.if_nfreq[ifNum])
2357 0 : select=false;
2358 : }
2359 : }
2360 : }
2361 : }
2362 1541 : return select;
2363 : }
2364 :
2365 0 : void ATCAFiller::shadow(Int row, Bool last)
2366 : {
2367 : // 1. collect together all rows with the same time
2368 : // 2. determine which antennas are shadowed
2369 : // 3. flag all baselines involving a shadowed antenna
2370 :
2371 0 : if (last || msc_p->time()(row)!=prevTime_p) {
2372 : // flag previous rows if needed
2373 0 : if (nRowCache_p>0) {
2374 0 : Vector<Bool> flag(nAnt_p); flag=false;
2375 0 : Vector<Double> uvw;
2376 0 : for (Int i=0; i< nRowCache_p; i++) {
2377 0 : Int k=rowCache_p[i];
2378 : // check for shadowing: projected baseline < shadow diameter
2379 0 : uvw = msc_p->uvw()(k);
2380 0 : double uvd2=uvw(0)*uvw(0)+uvw(1)*uvw(1);
2381 0 : if (uvd2>0) {
2382 0 : Bool shadowed = (uvd2<shadow_p*shadow_p);
2383 0 : if (shadowed) {
2384 : // w term decides which antenna is being shadowed
2385 0 : if (uvw(2)>0) {
2386 0 : flag(msc_p->antenna2()(k))=true;
2387 : // cout << " Shadowed: antenna2="<<msc_p->antenna2()(k)<<" for row "
2388 : // <<k<<" time="<< msc_p->timeMeas()(k) <<endl;
2389 : } else {
2390 0 : flag(msc_p->antenna1()(k))=true;
2391 : // cout << " Shadowed: antenna1="<<msc_p->antenna1()(k)<<" for row "
2392 : // <<k<<" time="<<msc_p->timeMeas()(k) <<endl;
2393 : }
2394 : }
2395 : }
2396 : }
2397 : // now flag rows with shadowed antennas
2398 0 : for (Int i=0; i< nRowCache_p; i++) {
2399 0 : Int k=rowCache_p[i];
2400 0 : Int ant1 = msc_p->antenna1()(k);
2401 0 : Int ant2 = msc_p->antenna2()(k);
2402 0 : if (flag(ant1) || flag(ant2)) {
2403 0 : flagCount_p(SHADOW)++;
2404 0 : msc_p->flagRow().put(k,true);
2405 : }
2406 : }
2407 : }
2408 : // reinitialize
2409 0 : if (!last) {
2410 0 : nRowCache_p=0;
2411 0 : prevTime_p=msc_p->time()(row);
2412 : }
2413 : }
2414 0 : if (!last) {
2415 0 : if (Int(rowCache_p.nelements())<=nRowCache_p) {
2416 0 : rowCache_p.resize(2*(nRowCache_p+1),true);
2417 : }
2418 0 : rowCache_p[nRowCache_p++]=row;
2419 : }
2420 0 : }
2421 :
2422 0 : Bool ATCAFiller::samplerFlag(Int row, Double posNegTolerance,
2423 : Double zeroTolerance) {
2424 : // check the sampler stats for the current row, return false if data needs to
2425 : // be flagged
2426 :
2427 0 : const Float posNegRef = 17.3;
2428 0 : const Float zeroRef =50.0;
2429 0 : Vector<Int> rows(2);
2430 0 : rows(0) = colSysCalIdAnt1(row);
2431 0 : rows(1) = colSysCalIdAnt2(row);
2432 0 : Bool flag = rows(0)<0 || rows(1)<0;
2433 0 : for (Int i=0; (!flag && i<2); i++) {
2434 0 : Vector<Float> neg = colSamplerStatsNeg(rows(i));
2435 0 : Vector<Float> pos = colSamplerStatsPos(rows(i));
2436 0 : Vector<Float> zero = colSamplerStatsZero(rows(i));
2437 0 : flag |=(abs(neg(0)-posNegRef)>posNegTolerance) ||
2438 0 : (abs(neg(1)-posNegRef)>posNegTolerance) ||
2439 0 : (abs(pos(0)-posNegRef)>posNegTolerance) ||
2440 0 : (abs(pos(1)-posNegRef)>posNegTolerance) ||
2441 0 : (abs(zero(0)-zeroRef)>zeroTolerance) ||
2442 0 : (abs(zero(1)-zeroRef)>zeroTolerance);
2443 : }
2444 0 : return flag;
2445 : }
2446 :
2447 0 : Vector<Double> ATCAFiller::opacities(Vector<Double> fGHz, Float tempK,Float humi,Float press,
2448 : Float height) {
2449 : // use (former) ASAP/Miriad code to calculate zenith opacities at a range of frequencies
2450 : // given the surface weather conditions
2451 :
2452 0 : ATAtmosphere atm = ATAtmosphere(tempK,press*100.,humi/100.);
2453 0 : atm.setObservatoryElevation(height);
2454 0 : return atm.zenithOpacities(fGHz*1e9);
2455 : }
2456 :
2457 : } //# end casa namespace
|
