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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
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 0 : ATCAFiller::ATCAFiller():
45 0 : appendMode_p(false),
46 0 : storedHeader_p(false),
47 0 : skipScan_p(false),
48 0 : skipData_p(false),
49 0 : firstHeader_p(false),
50 0 : listHeader_p(false),
51 0 : fileSize_p(0),
52 0 : birdie_p(false),
53 0 : reweight_p(false),
54 0 : noxycorr_p(false),
55 0 : obsType_p(0),
56 0 : hires_p(false),
57 0 : init_p(false),
58 0 : lastUT_p(0),
59 0 : bandWidth1_p(0),
60 0 : numChan1_p(0),
61 0 : shadow_p(0),
62 0 : autoFlag_p(true),
63 0 : flagScanType_p(false),
64 0 : flagCount_p(NFLAG,0)
65 0 : {}
66 :
67 0 : ATCAFiller::~ATCAFiller()
68 0 : {}
69 :
70 0 : Bool ATCAFiller::open(const String& msName, const Vector<String>& rpfitsFiles,
71 : const Vector<String>& options, Int opcor)
72 :
73 : {
74 : Int opt;
75 :
76 0 : LogOrigin orig("ATCAFiller", "open()", WHERE);
77 0 : os_p = LogIO(orig);
78 0 : rpfitsFiles_p = Directory::shellExpand(rpfitsFiles, false);
79 0 : GenSort<String>::sort(rpfitsFiles_p);
80 0 : if (rpfitsFiles_p.nelements() > 0) {
81 0 : os_p << LogIO::NORMAL << "Expanded file names are : " << endl;
82 0 : for (uInt i=0; i<rpfitsFiles_p.nelements(); i++) {
83 0 : if (rpfitsFiles_p(i).empty()) {
84 0 : os_p << "Input file number " << i+1 << " is empty" << LogIO::EXCEPTION;
85 : }
86 0 : os_p << rpfitsFiles_p(i) << endl;
87 : }
88 : } else {
89 0 : os_p << "Input file names vector is empty" << LogIO::EXCEPTION;
90 : }
91 0 : Vector<String> opts;
92 0 : if (options.nelements()==1) {
93 0 : opts=stringToVector(options(0),std::regex(" |,"));
94 : } else {
95 0 : opts=options;
96 : }
97 0 : opcor_p=opcor;
98 :
99 0 : Bool compress=false;
100 : // cerr<<"options="<<opts<<endl;
101 0 : for (opt=0; opt<Int(opts.nelements()); opt++) {
102 0 : if (downcase(opts(opt)) == "birdie") {
103 0 : birdie_p = true;
104 : }
105 0 : if (downcase(opts(opt)) == "reweight") {
106 0 : reweight_p = true;
107 : }
108 0 : if (downcase(opts(opt)) == "noxycorr") {
109 0 : noxycorr_p = true;
110 : }
111 0 : if (downcase(opts(opt)) == "compress") {
112 0 : compress = true;
113 : }
114 0 : if (downcase(opts(opt)) == "fastmosaic") {
115 0 : obsType_p=1;
116 : }
117 0 : if (downcase(opts(opt)) == "noautoflag") {
118 0 : autoFlag_p=false;
119 : }
120 0 : if (downcase(opts(opt)) == "hires") {
121 0 : hires_p=true;
122 : }
123 0 : if (downcase(opts(opt)) == "noac") {
124 0 : 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 0 : cabb_p = checkCABB(rpfitsFiles_p(0));
130 0 : reweight_p = reweight_p && !cabb_p;
131 0 : atms_p = makeTable(msName,compress,cabb_p);
132 0 : msc_p = new MSColumns(atms_p);
133 0 : msScanInfo_p = atms_p.keywordSet().asTable("ATCA_SCAN_INFO");
134 :
135 0 : String comp;
136 0 : if (compress) comp="Comp";
137 0 : dataAccessor_p = TiledDataStManAccessor(atms_p,"TiledData"+comp);
138 0 : sigmaAccessor_p = TiledDataStManAccessor(atms_p,"TiledSigma");
139 0 : flagAccessor_p = TiledDataStManAccessor(atms_p,"TiledFlag");
140 0 : flagCatAccessor_p = TiledDataStManAccessor(atms_p,"TiledFlagCategory");
141 :
142 0 : prev_fieldId_p = -1;
143 0 : lastWeatherUT_p=0;
144 0 : errCount_p=0;
145 0 : init();
146 0 : init_p=true;
147 0 : return true;
148 0 : }
149 :
150 :
151 0 : TableDesc ATCAFiller::atcaTableDesc(Bool compress)
152 : {
153 0 : TableDesc td = MS::requiredTableDesc();
154 : // add the data column
155 0 : MS::addColumnToDesc(td,MS::DATA,2);
156 : // and its unit
157 0 : td.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().define("UNIT","Jy");
158 0 : if (compress) {
159 0 : String col = MS::columnName(MS::DATA);
160 0 : td.addColumn(ArrayColumnDesc<Int>(col+"_COMPRESSED",
161 : "observed data compressed",2));
162 0 : td.addColumn(ScalarColumnDesc<Float>(col+"_SCALE"));
163 0 : td.addColumn(ScalarColumnDesc<Float>(col+"_OFFSET"));
164 0 : }
165 :
166 : // add ATCA specific columns here..
167 0 : MS::addColumnToDesc(td,MS::PULSAR_BIN);
168 :
169 0 : td.addColumn(ScalarColumnDesc<Int>("ATCA_SYSCAL_ID_ANT1",
170 : "Index into SysCal table for Antenna1"));
171 0 : td.addColumn(ScalarColumnDesc<Int>("ATCA_SYSCAL_ID_ANT2",
172 : "Index into SysCal table for Antenna2"));
173 : // the tiled column indices
174 0 : td.addColumn(ScalarColumnDesc<Int>("DATA_HYPERCUBE_ID",
175 : "Index for Data Tiling"));
176 0 : td.addColumn(ScalarColumnDesc<Int>("SIGMA_HYPERCUBE_ID",
177 : "Index for Sigma Tiling"));
178 0 : td.addColumn(ScalarColumnDesc<Int>("FLAG_HYPERCUBE_ID",
179 : "Index for Flag Tiling"));
180 0 : td.addColumn(ScalarColumnDesc<Int>("FLAG_CATEGORY_HYPERCUBE_ID",
181 : "Index for Flag Category Tiling"));
182 :
183 0 : return td;
184 0 : }
185 :
186 0 : MeasurementSet ATCAFiller::makeTable(const String& tableName, Bool compress,
187 : Bool cabb)
188 : {
189 : // make the MeasurementSet Table
190 0 : TableDesc atDesc = atcaTableDesc(compress);
191 :
192 0 : String colData = MS::columnName(MS::DATA);
193 0 : String comp1,comp2;
194 0 : if (compress) { comp1="Comp"; comp2="_COMPRESSED";}
195 :
196 : // define tiled hypercube for the data
197 0 : Vector<String> coordColNames(0); //# don't use coord columns
198 0 : atDesc.defineHypercolumn("TiledData"+comp1,3,
199 0 : stringToVector(colData+comp2),
200 : coordColNames,
201 0 : stringToVector("DATA_HYPERCUBE_ID"));
202 0 : atDesc.defineHypercolumn("TiledSigma",2,
203 0 : stringToVector(MS::columnName(MS::SIGMA)),
204 : coordColNames,
205 0 : stringToVector("SIGMA_HYPERCUBE_ID"));
206 0 : atDesc.defineHypercolumn("TiledFlag",3,
207 0 : stringToVector(MS::columnName(MS::FLAG)),
208 : coordColNames,
209 0 : stringToVector("FLAG_HYPERCUBE_ID"));
210 0 : atDesc.defineHypercolumn("TiledFlagCategory",4,
211 0 : stringToVector(MS::columnName(MS::FLAG_CATEGORY)),
212 : coordColNames,
213 0 : stringToVector("FLAG_CATEGORY_HYPERCUBE_ID"));
214 :
215 0 : SetupNewTable newtab(tableName, atDesc, Table::New);
216 :
217 : // Set the default Storage Manager to be the incremental one
218 0 : IncrementalStMan incrStMan ("IncrementalData");
219 0 : newtab.bindAll(incrStMan, true);
220 : // Make an exception for fast varying data
221 0 : StandardStMan stStMan ("StandardData");
222 0 : newtab.bindColumn(MS::columnName(MS::UVW),stStMan);
223 0 : newtab.bindColumn(MS::columnName(MS::ANTENNA2),stStMan);
224 0 : newtab.bindColumn("ATCA_SYSCAL_ID_ANT2",stStMan);
225 0 : if (compress) {
226 0 : newtab.bindColumn(colData+"_SCALE",stStMan);
227 0 : newtab.bindColumn(colData+"_OFFSET",stStMan);
228 : }
229 :
230 :
231 0 : TiledDataStMan tiledStMan1("TiledData"+comp1);
232 : // Bind the DATA & FLAG column to the tiled stman
233 0 : newtab.bindColumn(MS::columnName(MS::DATA)+comp2,tiledStMan1);
234 0 : newtab.bindColumn("DATA_HYPERCUBE_ID",tiledStMan1);
235 0 : 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 0 : }
240 :
241 0 : TiledDataStMan tiledStMan2("TiledSigma");
242 : // Bind the SIGMA column to its own tiled stman
243 0 : newtab.bindColumn(MS::columnName(MS::SIGMA),tiledStMan2);
244 0 : newtab.bindColumn("SIGMA_HYPERCUBE_ID",tiledStMan2);
245 :
246 0 : TiledDataStMan tiledStMan3("TiledFlag");
247 : // Bind the FLAG column to its own tiled stman
248 0 : newtab.bindColumn(MS::columnName(MS::FLAG),tiledStMan3);
249 0 : newtab.bindColumn("FLAG_HYPERCUBE_ID",tiledStMan3);
250 :
251 0 : TiledDataStMan tiledStMan3c("TiledFlagCategory");
252 : // Bind the FLAG_CATEGORY column to its own tiled stman
253 0 : newtab.bindColumn(MS::columnName(MS::FLAG_CATEGORY),tiledStMan3c);
254 0 : newtab.bindColumn("FLAG_CATEGORY_HYPERCUBE_ID",tiledStMan3c);
255 :
256 :
257 : // create the table (with 0 rows)
258 0 : MeasurementSet ms(newtab, 0);
259 : // create the subtables
260 0 : makeSubTables(ms, Table::New, cabb);
261 :
262 0 : return ms;
263 0 : }
264 :
265 0 : 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 0 : Int nrow=0;
270 :
271 : // Set up the subtables for the ATCA MS
272 :
273 : // Antenna
274 0 : TableDesc antTD=MSAntenna::requiredTableDesc();
275 0 : MSAntenna::addColumnToDesc(antTD, MSAntenna::PHASED_ARRAY_ID);
276 0 : MSAntenna::addColumnToDesc(antTD, MSAntenna::ORBIT_ID);
277 0 : SetupNewTable antennaSetup(ms.antennaTableName(),antTD,option);
278 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::ANTENNA),
279 0 : Table(antennaSetup,nrow));
280 :
281 : // Data descr
282 0 : TableDesc datadescTD=MSDataDescription::requiredTableDesc();
283 0 : SetupNewTable datadescSetup(ms.dataDescriptionTableName(),datadescTD,option);
284 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::DATA_DESCRIPTION),
285 0 : Table(datadescSetup,nrow));
286 :
287 : // Feed
288 0 : TableDesc feedTD=MSFeed::requiredTableDesc();
289 0 : MSFeed::addColumnToDesc(feedTD, MSFeed::PHASED_FEED_ID);
290 0 : SetupNewTable feedSetup(ms.feedTableName(),feedTD,option);
291 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FEED),
292 0 : Table(feedSetup,nrow));
293 :
294 : // Field
295 0 : TableDesc fieldTD=MSField::requiredTableDesc();
296 0 : SetupNewTable fieldSetup(ms.fieldTableName(),fieldTD,option);
297 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FIELD),
298 0 : Table(fieldSetup,nrow));
299 :
300 : // Flag_cmd
301 0 : TableDesc flagcmdTD=MSFlagCmd::requiredTableDesc();
302 0 : SetupNewTable flagcmdSetup(ms.flagCmdTableName(),flagcmdTD,option);
303 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::FLAG_CMD),
304 0 : Table(flagcmdSetup,nrow));
305 :
306 :
307 : // Observation
308 0 : TableDesc obsTD=MSObservation::requiredTableDesc();
309 0 : SetupNewTable observationSetup(ms.observationTableName(),obsTD,option);
310 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::OBSERVATION),
311 0 : Table(observationSetup,nrow));
312 :
313 : // History
314 0 : TableDesc historyTD=MSHistory::requiredTableDesc();
315 0 : SetupNewTable historySetup(ms.historyTableName(),historyTD,option);
316 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::HISTORY),
317 0 : Table(historySetup,nrow));
318 :
319 : // Pointing
320 0 : TableDesc pointingTD=MSPointing::requiredTableDesc();
321 0 : MSPointing::addColumnToDesc(pointingTD, MSPointing::POINTING_OFFSET);
322 0 : SetupNewTable pointingSetup(ms.pointingTableName(),pointingTD,option);
323 : // Pointing table can be large, set some sensible defaults for storageMgrs
324 0 : IncrementalStMan ismPointing ("ISMPointing");
325 0 : StandardStMan ssmPointing("SSMPointing",32768);
326 0 : pointingSetup.bindAll(ismPointing,true);
327 0 : pointingSetup.bindColumn(MSPointing::columnName(MSPointing::ANTENNA_ID),
328 : ssmPointing);
329 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::POINTING),
330 0 : Table(pointingSetup,nrow));
331 :
332 : // Polarization
333 0 : TableDesc polTD=MSPolarization::requiredTableDesc();
334 0 : SetupNewTable polSetup(ms.polarizationTableName(),polTD,option);
335 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::POLARIZATION),
336 0 : Table(polSetup,nrow));
337 :
338 : // Processor
339 0 : TableDesc processorTD=MSProcessor::requiredTableDesc();
340 0 : SetupNewTable processorSetup(ms.processorTableName(),processorTD,option);
341 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::PROCESSOR),
342 0 : Table(processorSetup,nrow));
343 :
344 : // Source
345 0 : TableDesc sourceTD=MSSource::requiredTableDesc();
346 0 : MSSource::addColumnToDesc(sourceTD, MSSource::SYSVEL);
347 0 : MSSource::addColumnToDesc(sourceTD, MSSource::NUM_LINES);
348 0 : MSSource::addColumnToDesc(sourceTD, MSSource::TRANSITION);
349 0 : MSSource::addColumnToDesc(sourceTD, MSSource::REST_FREQUENCY);
350 :
351 0 : SetupNewTable sourceSetup(ms.sourceTableName(),sourceTD,option);
352 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SOURCE),
353 0 : Table(sourceSetup,nrow));
354 :
355 :
356 : // SpectralWindow
357 0 : TableDesc spwTD=MSSpectralWindow::requiredTableDesc();
358 0 : if (!cabb) spwTD.addColumn(ScalarColumnDesc<Int>("ATCA_SAMPLER_BITS",
359 : "Number of bits used for sampling"));
360 0 : SetupNewTable spectralWindowSetup(ms.spectralWindowTableName(),spwTD,option);
361 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SPECTRAL_WINDOW),
362 0 : Table(spectralWindowSetup,nrow));
363 :
364 : // State
365 0 : TableDesc stateTD=MSState::requiredTableDesc();
366 0 : SetupNewTable stateSetup(ms.stateTableName(),stateTD,option);
367 0 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::STATE),
368 0 : Table(stateSetup,nrow));
369 :
370 : // SysCal
371 0 : TableDesc sysCalTD=MSSysCal::requiredTableDesc();
372 0 : Vector<String> cols(4);
373 0 : 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 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
391 : ("ATCA_GTP",
392 : "Noise Cal On+Off Autocorrelation",
393 0 : IPosition(1,2),ColumnDesc::Direct));
394 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
395 : ("ATCA_SDO",
396 : "Noise Cal On-Off Autocorrelation",
397 0 : IPosition(1,2),ColumnDesc::Direct));
398 0 : sysCalTD.addColumn(ArrayColumnDesc<Float>
399 : ("ATCA_CALJY",
400 : "Calibration factor",
401 0 : IPosition(1,2),ColumnDesc::Direct));
402 :
403 0 : cols(1)="ATCA_GTP";
404 0 : cols(2)="ATCA_SDO";
405 0 : cols(3)="ATCA_CALJY";
406 : }
407 0 : sysCalTD.addColumn(ScalarColumnDesc<Float>
408 : ("ATCA_XY_AMPLITUDE",
409 : "Noise source cross correlation amplitude"));
410 0 : sysCalTD.addColumn(ScalarColumnDesc<Float>
411 : ("ATCA_TRACK_ERR_MAX",
412 : "Max tracking error over non blanked cycle"));
413 0 : TableQuantumDesc tqd1(sysCalTD,"ATCA_TRACK_ERR_MAX",Unit("arcsec"));
414 0 : tqd1.write(sysCalTD);
415 0 : sysCalTD.addColumn(ScalarColumnDesc<Float>
416 : ("ATCA_TRACK_ERR_RMS",
417 : "RMS tracking error over non blanked cycle"));
418 0 : TableQuantumDesc tqd2(sysCalTD,"ATCA_TRACK_ERR_RMS",Unit("arcsec"));
419 0 : tqd2.write(sysCalTD);
420 :
421 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TSYS,
422 0 : IPosition(1, 2), ColumnDesc::Direct);
423 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TCAL,
424 0 : IPosition(1, 2), ColumnDesc::Direct);
425 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TRX,
426 0 : IPosition(1, 2),ColumnDesc::Direct);
427 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TSYS_FLAG);
428 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::PHASE_DIFF);
429 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::PHASE_DIFF_FLAG);
430 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TCAL_FLAG);
431 0 : MSSysCal::addColumnToDesc(sysCalTD, MSSysCal::TRX_FLAG);
432 0 : cols(0) = MSSysCal::columnName(MSSysCal::TSYS);
433 0 : sysCalTD.defineHypercolumn("TiledSysCal",2,cols);
434 0 : SetupNewTable sysCalSetup(ms.sysCalTableName(),sysCalTD,option);
435 0 : IncrementalStMan incStMan;
436 0 : sysCalSetup.bindAll(incStMan);
437 0 : TiledColumnStMan tiledStManSysCal("TiledSysCal",IPosition(2,2,1024));
438 :
439 0 : for (uInt i=0; i<cols.nelements(); i++) {
440 0 : sysCalSetup.bindColumn(cols(i),tiledStManSysCal);
441 : }
442 0 : ms.rwKeywordSet().defineTable(ms.keywordName(MS::SYSCAL),
443 0 : Table(sysCalSetup,nrow));
444 :
445 : // Weather
446 0 : TableDesc weatherTD=MSWeather::requiredTableDesc();
447 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::PRESSURE);
448 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::REL_HUMIDITY);
449 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::TEMPERATURE);
450 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_DIRECTION);
451 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_SPEED);
452 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::PRESSURE_FLAG);
453 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::REL_HUMIDITY_FLAG);
454 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::TEMPERATURE_FLAG);
455 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_DIRECTION_FLAG);
456 0 : MSWeather::addColumnToDesc(weatherTD,MSWeather::WIND_SPEED_FLAG);
457 0 : weatherTD.addColumn(ScalarColumnDesc<Float>
458 : ("ATCA_RAIN_GAUGE",
459 : "Rain since 10am local time"));
460 0 : TableQuantumDesc tqds0(weatherTD,"ATCA_RAIN_GAUGE",Unit("mm"));
461 0 : tqds0.write(weatherTD);
462 0 : weatherTD.addColumn(ScalarColumnDesc<Float>
463 : ("ATCA_SEEMON_PHASE",
464 : "Seeing monitor raw phase at 22GHz"));
465 0 : TableQuantumDesc tqds1(weatherTD,"ATCA_SEEMON_PHASE",Unit("rad"));
466 0 : tqds1.write(weatherTD);
467 0 : weatherTD.addColumn(ScalarColumnDesc<Float>
468 : ("ATCA_SEEMON_RMS",
469 : "Seeing monitor RMS phase"));
470 0 : TableQuantumDesc tqds2(weatherTD,"ATCA_SEEMON_RMS",Unit("mm"));
471 0 : tqds2.write(weatherTD);
472 0 : weatherTD.addColumn(ScalarColumnDesc<Bool>
473 : ("ATCA_SEEMON_FLAG",
474 : "Seeing monitor flag"));
475 0 : SetupNewTable weatherSetup(ms.weatherTableName(),weatherTD,option);
476 0 : ms.rwKeywordSet().defineTable(ms.keywordName(MS::WEATHER),
477 0 : Table(weatherSetup,nrow));
478 :
479 : // ATCA_SCAN_INFO
480 0 : TableDesc atsiTD;
481 0 : atsiTD.addColumn(ScalarColumnDesc<Int>("ANTENNA_ID",
482 : "Antenna Id"));
483 0 : atsiTD.addColumn(ScalarColumnDesc<Int>("SCAN_ID",
484 : "Scan number from main table"));
485 0 : atsiTD.addColumn(ScalarColumnDesc<Int>("SPECTRAL_WINDOW_ID",
486 : "Spectral window Id"));
487 0 : atsiTD.addColumn(ArrayColumnDesc<Int>("FINE_ATTENUATOR",
488 : "Fine Attenuator setting A,B",
489 0 : IPosition(1,2),ColumnDesc::Direct));
490 0 : atsiTD.addColumn(ArrayColumnDesc<Int>("COARSE_ATTENUATOR",
491 : "COARSE Attenuator setting A,B",
492 0 : IPosition(1,2),ColumnDesc::Direct));
493 0 : atsiTD.addColumn(ArrayColumnDesc<Int>("MM_ATTENUATOR",
494 : "mm Attenuator setting A,B",
495 0 : IPosition(1,2),ColumnDesc::Direct));
496 0 : atsiTD.addColumn(ArrayColumnDesc<Float>("SUBREFLECTOR",
497 : "Subreflector position(center,edge/tilt)",
498 0 : IPosition(1,2),ColumnDesc::Direct));
499 0 : TableQuantumDesc tqd3(atsiTD,"SUBREFLECTOR",Unit("m"));
500 0 : tqd3.write(atsiTD);
501 0 : atsiTD.addColumn(ScalarColumnDesc<String>("CORR_CONFIG",
502 : "Correlator configuration"));
503 0 : atsiTD.addColumn(ScalarColumnDesc<String>("SCAN_TYPE",
504 : "Scan type"));
505 0 : atsiTD.addColumn(ScalarColumnDesc<String>("COORD_TYPE",
506 : "CAOBS Coordinate type"));
507 0 : atsiTD.addColumn(ScalarColumnDesc<String>("POINTING_INFO",
508 : "Pointing info - details of last point scan"));
509 0 : atsiTD.addColumn(ScalarColumnDesc<Bool>("LINE_MODE",
510 : "Line Mode: true=spectrum, false=mfs"));
511 0 : atsiTD.addColumn(ScalarColumnDesc<Int>("CACAL_CNT",
512 : "Online calibration counter"));
513 :
514 0 : SetupNewTable scanInfoSetup(ms.tableName()+"/ATCA_SCAN_INFO",atsiTD,option);
515 0 : IncrementalStMan incSIStMan("ISMScanInfo");
516 0 : StandardStMan stdSIStMan("SSMScanInfo",32768);
517 0 : scanInfoSetup.bindAll(incSIStMan);
518 0 : scanInfoSetup.bindColumn("ANTENNA_ID",stdSIStMan);
519 0 : ms.rwKeywordSet().defineTable("ATCA_SCAN_INFO",
520 0 : Table(scanInfoSetup,nrow));
521 :
522 : // update the references to the subtable keywords
523 0 : ms.initRefs();
524 0 : }
525 :
526 :
527 0 : void ATCAFiller::init()
528 : {
529 : //Initialize selection
530 0 : scanRange(0, 0);
531 0 : freqRange(0.0,1.e30);
532 0 : Vector<String> fieldNames(0);
533 0 : fields(fieldNames);
534 :
535 : // extra ID columns in main table
536 0 : colSysCalIdAnt1.attach(atms_p,"ATCA_SYSCAL_ID_ANT1");
537 0 : colSysCalIdAnt2.attach(atms_p,"ATCA_SYSCAL_ID_ANT2");
538 :
539 : // extra spectralWindow table columns
540 0 : if (!cabb_p) colSamplerBits.attach(atms_p.spectralWindow(),"ATCA_SAMPLER_BITS");
541 :
542 : // extra syscal table columns
543 0 : if (cabb_p) {
544 0 : colGTP.attach(atms_p.sysCal(),"ATCA_GTP");
545 0 : colSDO.attach(atms_p.sysCal(),"ATCA_SDO");
546 0 : 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 0 : colXYAmplitude.attach(atms_p.sysCal(),"ATCA_XY_AMPLITUDE");
554 0 : colTrackErrMax.attach(atms_p.sysCal(),"ATCA_TRACK_ERR_MAX");
555 0 : colTrackErrRMS.attach(atms_p.sysCal(),"ATCA_TRACK_ERR_RMS");
556 :
557 : // ScanInfo table columns
558 0 : Table sit(atms_p.keywordSet().asTable("ATCA_SCAN_INFO"));
559 0 : colScanInfoAntId.attach(sit,"ANTENNA_ID");
560 0 : colScanInfoScanId.attach(sit,"SCAN_ID");
561 0 : colScanInfoSpWId.attach(sit,"SPECTRAL_WINDOW_ID");
562 0 : colScanInfoCacal.attach(sit,"CACAL_CNT");
563 0 : colScanInfoFine.attach(sit,"FINE_ATTENUATOR");
564 0 : colScanInfoCoarse.attach(sit,"COARSE_ATTENUATOR");
565 0 : colScanInfommAtt.attach(sit,"MM_ATTENUATOR");
566 0 : colScanInfoSubreflector.attach(sit,"SUBREFLECTOR");
567 0 : colScanInfoCorrConfig.attach(sit,"CORR_CONFIG");
568 0 : colScanInfoScanType.attach(sit,"SCAN_TYPE");
569 0 : colScanInfoCoordType.attach(sit,"COORD_TYPE");
570 0 : colScanInfoPointInfo.attach(sit,"POINTING_INFO");
571 0 : colScanInfoLineMode.attach(sit,"LINE_MODE");
572 :
573 : // Extra weather table columns
574 0 : colWeatherRainGauge.attach(atms_p.weather(),"ATCA_RAIN_GAUGE");
575 0 : colWeatherSeeMonPhase.attach(atms_p.weather(),"ATCA_SEEMON_PHASE");
576 0 : colWeatherSeeMonRMS.attach(atms_p.weather(),"ATCA_SEEMON_RMS");
577 0 : colWeatherSeeMonFlag.attach(atms_p.weather(),"ATCA_SEEMON_FLAG");
578 :
579 0 : nScan_p=nSpW_p=nField_p=scanNo_p=0;
580 0 : gotAN_p=false;
581 0 : }
582 :
583 : // list the current scan and return the current day in mjd
584 0 : void ATCAFiller::listScan(Double & mjd, Int scan, Double ut)
585 : {
586 : // Convert observation date to mjd
587 : Int day,month,year;
588 0 : sscanf(names_.datobs,"%4d-%2d-%2d",&year,&month,&day);
589 0 : MVTime mjd_date(year,month,(Double)day);
590 0 : mjd=mjd_date.second();
591 0 : mjd_date=MVTime((mjd_date.second()+ut)/C::day);
592 0 : os_p << LogIO::NORMAL << "Scan # : "<< scan << endl;
593 0 : os_p << LogIO::NORMAL << "Object : "<< String(names_.object,16) << endl;
594 0 : os_p << LogIO::NORMAL << "Date : "<< mjd_date.string(MVTime::YMD)
595 0 : << LogIO::POST;
596 0 : }
597 :
598 :
599 0 : Bool ATCAFiller::fill() {
600 0 : if (!init_p) return false;
601 0 : if (!appendMode_p) {
602 0 : firstHeader_p = true;
603 0 : skipScan_p=false;
604 0 : skipData_p=false;
605 :
606 0 : nScan_p=1; // we've seen the 1st header
607 0 : scanNo_p=-1; // make zero based for storage in MS
608 :
609 : Int fileno;
610 0 : Int offset=0;
611 :
612 0 : for (fileno=0; fileno<Int(rpfitsFiles_p.nelements())-offset; fileno++) {
613 0 : listHeader_p = true;
614 0 : currentFile_p = rpfitsFiles_p(fileno);
615 0 : fill1(currentFile_p);
616 : }
617 :
618 0 : os_p << LogIO::DEBUGGING << "FillFeed" << LogIO::POST;
619 0 : fillFeedTable();
620 :
621 0 : fillObservationTable();
622 :
623 0 : fillMeasureReferences();
624 0 : 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 0 : }
689 0 : }
690 0 : }
691 : }
692 : // flag the last integration for shadowing
693 0 : 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 0 : flush();
699 0 : unlock();
700 :
701 0 : return true;
702 : }
703 :
704 0 : Bool ATCAFiller::checkCABB(const String& rpfitsFile)
705 : {
706 : Int jstat;
707 0 : Regex trailing(" *$"); // trailing blanks
708 :
709 0 : String file = rpfitsFile;
710 0 : os_p << LogIO::NORMAL <<"Checking header of file "<<file<< LogIO::POST;
711 0 : strcpy(names_.file,file.chars());
712 0 : jstat=-2; // open rpfits file and read first header
713 0 : param_.ncard=-1; // collect cards into card array
714 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
715 : &if_no, &sourceno);
716 0 : 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 0 : String instrument = String(names_.instrument,16).before(trailing);
723 : // close the file
724 0 : jstat=1;
725 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
726 : &if_no, &sourceno);
727 0 : if (instrument=="ATCABB") {
728 0 : os_p<< LogIO::NORMAL<<"CABB data detected"<<LogIO::POST;
729 0 : return true;
730 : }
731 0 : else return false;
732 0 : }
733 :
734 0 : Bool ATCAFiller::fill1(const String& rpfitsFile)
735 : {
736 : Int jstat;
737 0 : Regex trailing(" *$"); // trailing blanks
738 :
739 0 : String file = rpfitsFile;
740 0 : if (listHeader_p == true) {
741 0 : os_p << LogIO::NORMAL <<"Reading file "<<file<< LogIO::POST;
742 0 : strcpy(names_.file,file.chars());
743 0 : jstat=-2; // open rpfits file and read first header
744 0 : param_.ncard=-1; // collect cards into card array
745 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
746 : &if_no, &sourceno);
747 0 : 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 0 : jstat=0;
758 :
759 0 : Double lastUT=0,lastUT2=0;
760 0 : Int lastScanNo=-1,lastIF=-1;
761 0 : flagCount_p=0;
762 0 : while (jstat==0) {
763 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w, &flg, &bin,
764 : &if_no, &sourceno);
765 :
766 0 : 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 0 : case 0: // found data: store header (incl Tables) and data
773 0 : if (listHeader_p) { listScan(mjd0_p,nScan_p,ut); listHeader_p=false;}
774 :
775 : // do we want this scan?
776 0 : if (skipScan_p || nScan_p<firstScan_p ||
777 0 : (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 0 : if (baseline==-1) {if_no=1; sourceno=1;}
784 :
785 : // check if if_no valid
786 0 : 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 0 : if (lastUT!=ut) {
793 0 : String field=String(names_.object,16).before(trailing);
794 0 : 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 0 : }
801 0 : if (!skipScan_p) {
802 0 : if_no--; // go from 1-based to 0-based indexing (f2c)
803 0 : sourceno--; // same here
804 : // skip unselected or garbled data
805 0 : skipData_p=Bool(if_no<0 || !selected(if_no) || sourceno<0);
806 : // reject auto correlations?
807 0 : if (noac_p) skipData_p |= (Int(baseline)%257==0);
808 : //cerr<<"SkipData="<<skipData_p<< " ifno="<<if_no<< " sourceno="<<sourceno<<endl;
809 : }
810 : }
811 : }
812 0 : if (!skipScan_p && !skipData_p && firstHeader_p && anten_.nant>0) {
813 0 : nAnt_p=anten_.nant;
814 0 : Int NChan=if_.if_nfreq[if_no];
815 0 : Int NPol=if_.if_nstok[if_no];
816 0 : firstHeader_p=false;
817 0 : os_p << LogIO::NORMAL << " First data/header has NAnt="<<nAnt_p<<
818 0 : ", NChan="<<NChan<<", NPol="<<NPol<< LogIO::POST;
819 : }
820 0 : if (!skipScan_p && !skipData_p) {
821 0 : 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 0 : if (!skipScan_p && !skipData_p && !storedHeader_p) {
829 0 : storeHeader();
830 0 : scanNo_p++;
831 0 : storedHeader_p=true;
832 : }
833 0 : if (!skipScan_p) {
834 0 : if (baseline==-1) {
835 : // we want at least some of the syscal data for this scan
836 0 : storeSysCal();
837 : }
838 0 : else if (!skipData_p) {
839 0 : if (lastUT2!=ut || lastIF!=if_no){
840 0 : checkSpW(if_no); // we need to check these every integration to
841 0 : if (lastUT2!=ut) checkField();// cope with freq. switching and mosaicing
842 0 : lastUT2=ut;
843 0 : lastIF=if_no;
844 : }
845 0 : storeData();
846 : }
847 : }
848 0 : skipData_p=false;
849 0 : 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 0 : if (jstat == 3) { // because of the break suppression in case 1
877 0 : 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 0 : os_p << LogIO::NORMAL << "End of File" << LogIO::POST;
883 0 : eof_p = true;
884 0 : nScan_p++; // increment for next scan
885 : // print flagging stats for last scan
886 0 : if (flagCount_p(COUNT)>0) {
887 0 : Vector<Float> perc(NFLAG);
888 0 : for (Int i=0; i<NFLAG; i++) perc(i)=0.1*((1000*flagCount_p(i))/flagCount_p(COUNT));
889 0 : os_p<< LogIO::NORMAL <<"Number of rows selected : "<<flagCount_p(COUNT)<<endl;
890 0 : os_p<< LogIO::NORMAL <<"Flagged : "<<perc(FLAG)<<"%"<<endl;
891 0 : os_p<< LogIO::NORMAL <<" Antenna off source : "<<perc(ONLINE)<<"%"<<endl;
892 0 : os_p<< LogIO::NORMAL <<" ScanType (Point/Paddle): "<<perc(SCANTYPE)<<"%"<<endl;
893 0 : if (!cabb_p) os_p<< LogIO::NORMAL <<" Bad Sampler stats : "<<perc(SYSCAL)<<"%"<<endl;
894 0 : if (shadow_p>0) os_p<< LogIO::NORMAL <<" Antenna shadowed : "<<perc(SHADOW)<<"%"<<LogIO::POST;
895 0 : }
896 : }
897 :
898 0 : 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 0 : } else {
912 0 : jstat=1; // close file
913 0 : rpfitsin_(&jstat, vis, weight, &baseline, &ut, &u, &v, &w,
914 : &flg, &bin, &if_no, &sourceno);
915 : }
916 0 : jstat=1; // exit loop
917 0 : 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 0 : case 5:
925 : // jstat == 5 is usually record padding at the end of an integration.
926 : // just continue on loading..
927 0 : jstat=0;
928 0 : 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 0 : lastUT=ut;
935 : }
936 :
937 : // uncomment this if things take unexpectedly long
938 : // dataAccessor_p.showCacheStatistics(cout);
939 : // sigmaAccessor_p.showCacheStatistics(cout);
940 0 : return true;
941 0 : }
942 :
943 :
944 :
945 :
946 0 : void ATCAFiller::fillMeasureReferences()
947 : {
948 0 : String key("MEASURE_REFERENCE");
949 0 : msc_p->time().rwKeywordSet().define(key,"UTC");
950 0 : msc_p->uvw().rwKeywordSet().define(key,"J2000");
951 0 : msc_p->antenna().position().rwKeywordSet().define(key,"ITRF");
952 0 : msc_p->feed().time().rwKeywordSet().define(key,"UTC");
953 0 : msc_p->field().time().rwKeywordSet().define(key,"UTC");
954 0 : msc_p->field().delayDir().rwKeywordSet().define(key,"J2000");
955 0 : msc_p->field().phaseDir().rwKeywordSet().define(key,"J2000");
956 0 : msc_p->field().referenceDir().rwKeywordSet().define(key,"J2000");
957 0 : msc_p->source().time().rwKeywordSet().define(key,"UTC");
958 0 : msc_p->source().direction().rwKeywordSet().define(key,"J2000");
959 0 : msc_p->history().time().rwKeywordSet().define(key,"UTC");
960 0 : msc_p->spectralWindow().chanFreq().rwKeywordSet().define(key,"TOPO");
961 0 : msc_p->spectralWindow().refFrequency().rwKeywordSet().define(key,"TOPO");
962 :
963 : MFrequency::Types tp;
964 0 : MFrequency::getType(tp, msc_p->spectralWindow().refFrequency().keywordSet().asString("MEASURE_REFERENCE"));
965 0 : Int meas_freq_ref = tp;
966 0 : msc_p->spectralWindow().measFreqRef().fillColumn(meas_freq_ref);
967 0 : msc_p->sysCal().time().rwKeywordSet().define(key,"UTC");
968 0 : msc_p->weather().time().rwKeywordSet().define(key,"UTC");
969 0 : msc_p->pointing().direction().rwKeywordSet().define(key,"J2000");
970 0 : msc_p->pointing().pointingOffset().rwKeywordSet().define(key,"AZEL");
971 0 : }
972 :
973 0 : ATCAFiller & ATCAFiller::scanRange(Int first, Int last)
974 : {
975 0 : firstScan_p=first;
976 0 : lastScan_p=max(0,last);
977 0 : return *this;
978 : }
979 :
980 0 : ATCAFiller & ATCAFiller::freqRange(Double low, Double high)
981 : {
982 0 : lowFreq_p=low;
983 0 : highFreq_p=high;
984 0 : return *this;
985 : }
986 :
987 0 : ATCAFiller & ATCAFiller::freqSel(const Vector<Int>& spws)
988 : {
989 0 : spws_p=spws;
990 0 : return *this;
991 : }
992 :
993 0 : ATCAFiller & ATCAFiller::fields(const Vector<String>& fieldList)
994 : {
995 0 : fieldSelection_p=fieldList;
996 0 : 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 0 : ATCAFiller & ATCAFiller::edge(float edge)
1031 : {
1032 0 : edge_p=edge;
1033 0 : return *this;
1034 : }
1035 :
1036 0 : void ATCAFiller::storeHeader(bool last) {
1037 : // Bool skip=false;
1038 0 : Regex trailing(" *$"); // trailing blanks
1039 :
1040 : // First check antenna Table and store any new stations encountered
1041 0 : if (!gotAN_p && anten_.nant>0) {
1042 0 : if (anten_.an_found) {
1043 0 : gotAN_p=true;
1044 :
1045 0 : for (Int i=0; i<nAnt_p; i++) {
1046 0 : Int ant=Int(anten_.ant_num[i]); // 1-based !
1047 0 : 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 0 : atms_p.antenna().addRow();
1055 0 : Int n=atms_p.antenna().nrow()-1;
1056 :
1057 0 : String instrument = String(names_.instrument,16).before(trailing);
1058 0 : String station(&names_.sta[i*8],8);
1059 0 : os_p << LogIO::NORMAL <<"Antenna : "<< station << " ";
1060 0 : Vector<Double> xyz(3);
1061 0 : xyz(0)=doubles_.x[i];xyz(1)=doubles_.y[i];xyz(2)=doubles_.z[i];
1062 0 : os_p << " position:"<<xyz << endl;
1063 0 : if (instrument=="ATCA" || cabb_p) {
1064 0 : String antName;
1065 : // construct antenna name
1066 0 : ostringstream ostr; ostr << "CA0" <<i+1;
1067 0 : String str = ostr.str();
1068 0 : msc_p->antenna().name().put(n,str);
1069 0 : msc_p->antenna().station().put(n,atcaPosToStation(xyz));
1070 0 : } 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 0 : msc_p->antenna().orbitId().put(n,-1);
1075 0 : msc_p->antenna().phasedArrayId().put(n,-1);
1076 0 : msc_p->antenna().dishDiameter().put(n,22.0);
1077 0 : msc_p->antenna().type().put(n, "GROUND-BASED");
1078 :
1079 0 : if (anten_.ant_mount[i]==0) {
1080 0 : msc_p->antenna().mount().put(n,"alt-az");
1081 : }
1082 : else {
1083 0 : msc_p->antenna().mount().put(n,"bizarre");
1084 : }
1085 0 : msc_p->antenna().position().put(n,xyz);
1086 0 : Vector<Double> offset(3); offset=0.;
1087 : // todo: figure out coordinate system of offset
1088 0 : offset(0)=doubles_.axis_offset[i];
1089 0 : 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 0 : }
1102 0 : 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 0 : if (!last) checkSpW(if_no);
1111 :
1112 : // Check if current Observation info already stored, add it if not
1113 0 : checkObservation();
1114 :
1115 : // Store the ATCA header cards
1116 0 : storeATCAHeader(last);
1117 :
1118 : // Check if we've seen current source before, if not add to table
1119 0 : if (!last) checkField();
1120 0 : }
1121 :
1122 0 : void ATCAFiller::storeATCAHeader(Bool last) {
1123 0 : uInt ncard = abs(param_.ncard);
1124 : //cout<<" #cards = "<<ncard<<endl;
1125 0 : String cards = String(names_.card,ncard*80);
1126 0 : const Int nTypes = 10;
1127 0 : Block<String> types(nTypes);
1128 0 : types[0]= "OBSLOG";
1129 0 : types[1] = "ATTEN";
1130 0 : types[2] = "SUBREFL";
1131 0 : types[3] = "CORR_CFG";
1132 0 : types[4] = "SCANTYPE";
1133 0 : types[5] = "COORDTYP";
1134 0 : types[6] = "LINEMODE";
1135 0 : types[7] = "CACALCNT";
1136 0 : types[8] = "POINTCOR";
1137 0 : types[9] = "POINTINF";
1138 0 : String obsLog = "";
1139 0 : String config = "none";
1140 0 : const Int nIF = if_.n_if;
1141 0 : const Int nAnt = anten_.nant;
1142 0 : Cube<Int> fine(2,nIF,nAnt,0),coarse(2,nIF,6,0);
1143 0 : Matrix<Int> mmAtt(nIF,nAnt,0);
1144 0 : Vector<Float> subrPos(nAnt,0),subrTilt(nAnt,0);
1145 0 : Matrix<Float> pointCorr(2,nAnt,0);
1146 0 : newPointingCorr_p=false;
1147 0 : flagScanType_p=false;
1148 0 : String scanType,coordType,pointInfo;
1149 0 : Vector<Bool> lineMode(nIF);
1150 0 : Int cacalCnt=0;
1151 0 : const Regex trailing(" *$"); // trailing blanks
1152 0 : Bool foundAny = false;
1153 0 : String::size_type pos=cards.find("EPHEM12"); // last 'standard' card
1154 : //cout << "pos="<<pos<<" String::npos=="<<String::npos<<endl;
1155 0 : if (pos==String::npos) return;
1156 0 : uInt iFirst=pos/80+1;
1157 0 : for (uInt i=iFirst; i<ncard; i++) {
1158 : // extract card
1159 0 : String card = cards.substr(i*80,80);
1160 : //cout << "card = "<< card<<endl;
1161 : // read antenna number (if present)
1162 0 : Int ant = 0;
1163 0 : pos=card.find("CA0");
1164 0 : if (pos!=String::npos) istringstream(card.substr(pos+2,2))>> ant;
1165 0 : ant-=1; // zero based indexing
1166 0 : if (ant>=nAnt) ant = -1;
1167 0 : for (Int j=0; j<nTypes; j++) {
1168 0 : if (card.find(types[j])==0) {
1169 0 : foundAny=true;
1170 : //cout << "Found card :"<<types[j]<<", ant="<<ant<<endl;
1171 0 : switch (j) {
1172 0 : case 0: obsLog+=card.substr(8,72).before(trailing);
1173 0 : obsLog+="\n";
1174 0 : break;
1175 0 : case 1: { //ATTEN
1176 0 : String::size_type pos=card.find("FINE=");
1177 0 : 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 0 : pos=card.find("COARSE=");
1184 0 : 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 0 : pos=card.find("MM=");
1191 0 : 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 0 : break;
1197 0 : case 2: { //SUBREFL
1198 0 : String::size_type pos=card.find("POS=");
1199 0 : if (pos!=String::npos && ant>=0) {
1200 0 : istringstream(card.substr(pos+4,6))>>subrPos(ant);
1201 : }
1202 0 : pos=card.find("TILT=");
1203 0 : if (pos!=String::npos && ant>=0) {
1204 0 : istringstream(card.substr(pos+5,6))>>subrTilt(ant);
1205 : }
1206 : }
1207 0 : break;
1208 0 : case 3: {// CORR_CFG
1209 0 : String::size_type pos=card.find("=");
1210 0 : if (pos!=String::npos) {
1211 0 : config=card.substr(pos+3,80-pos-3).before(trailing);
1212 : }
1213 : }
1214 0 : break;
1215 0 : case 4: {// SCANTYPE
1216 0 : String::size_type pos=card.find("=");
1217 0 : if (pos!=String::npos) {
1218 0 : scanType=card.substr(pos+2,80-pos-3).before(trailing);
1219 : }
1220 0 : if (scanType=="PADDLE"||scanType=="POINT"||scanType=="XPOINT"
1221 0 : ||scanType=="EARLY") {
1222 0 : flagScanType_p=true;
1223 : }
1224 : }
1225 0 : break;
1226 0 : case 5: {// COORDTYP
1227 0 : String::size_type pos=card.find("=");
1228 0 : if (pos!=String::npos) {
1229 0 : coordType=card.substr(pos+2,80-pos-3).before(trailing);
1230 : }
1231 : }
1232 0 : break;
1233 0 : case 6: {// LINEMODE
1234 0 : String::size_type pos=card.find("=");
1235 0 : if (pos!=String::npos) {
1236 0 : lineMode(0)=(card[pos+2]=='T');
1237 0 : if (nIF>1) lineMode(1)=(card[pos+4]=='T');
1238 : }
1239 : }
1240 0 : break;
1241 0 : case 7: {// CACALCNT
1242 0 : String::size_type pos=card.find("=");
1243 0 : if (pos!=String::npos) {
1244 0 : istringstream(card.substr(pos+1,7))>>cacalCnt;
1245 : }
1246 : }
1247 0 : 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 0 : }
1272 0 : if (!foundAny) return;
1273 0 : Vector<String> tmp;
1274 0 : if (msc_p->observation().log().isDefined(obsId_p)) {
1275 0 : tmp=msc_p->observation().log()(obsId_p);
1276 : }
1277 0 : if (tmp.nelements()==0) tmp.resize(1);
1278 0 : Int index=tmp.nelements()-1;
1279 0 : tmp(index)+=obsLog;
1280 0 : msc_p->observation().log().put(obsId_p,tmp);
1281 0 : if ((nAnt*nIF)==0 || last) return;
1282 :
1283 : // find out spectral window index of IFs
1284 0 : Vector<Int> spwId(nIF,-1);
1285 0 : if (if_no>=0) spwId(if_no)=spWId_p;
1286 0 : if (nIF>1) {
1287 0 : for (Int ifNum=0; ifNum<nIF; ifNum++) {
1288 0 : if (ifNum!=if_no) {
1289 0 : if (selected(ifNum)) spwId(ifNum)=checkSpW(ifNum);
1290 : }
1291 : }
1292 : }
1293 : // reset to original spW
1294 0 : if (if_no>0) checkSpW(if_no);
1295 :
1296 0 : Int row=msScanInfo_p.nrow();
1297 0 : msScanInfo_p.addRow(nAnt*nIF);
1298 0 : colScanInfoScanId.put(row,scanNo_p+1);
1299 0 : colScanInfoCacal.put(row,cacalCnt);
1300 0 : colScanInfoCorrConfig.put(row,config);
1301 0 : colScanInfoScanType.put(row,scanType);
1302 0 : colScanInfoCoordType.put(row,coordType);
1303 0 : colScanInfoPointInfo.put(row,pointInfo);
1304 0 : if (newPointingCorr_p) {
1305 0 : pointingCorr_p.reference(pointCorr);
1306 : } else {
1307 0 : pointingCorr_p.resize(0,0);
1308 : }
1309 :
1310 0 : Vector<Int> f(2),c(2),m(2);
1311 0 : Vector<Float> sr(2);
1312 0 : for (Int i=0; i<nIF; i++) {
1313 0 : colScanInfoSpWId.put(row,spwId(i));
1314 0 : colScanInfoLineMode.put(row,lineMode(i));
1315 0 : for (Int ant=0; ant<nAnt; ant++) {
1316 0 : colScanInfoAntId.put(row,ant);
1317 0 : 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 0 : m(0)=mmAtt(0,ant); m(1)=mmAtt(1,ant);
1324 0 : colScanInfommAtt.put(row,m);
1325 0 : sr(0)=subrPos(ant)/1000.0; sr(1)=subrTilt(ant)/1000.; // convert to meter
1326 0 : colScanInfoSubreflector.put(row,sr);
1327 0 : row++;
1328 : }
1329 : }
1330 0 : }
1331 :
1332 0 : String ATCAFiller::atcaPosToStation(Vector<Double>& xyz) {
1333 0 : String station("NONE");
1334 : // determine station from xyz position
1335 : // Use W106 as reference
1336 0 : Double x106 = -4751615.0l;
1337 0 : Double y106 = 2791719.246l;
1338 0 : Double z106 = -3200483.747l;
1339 0 : Double incr = 6000.0l/392;
1340 0 : Bool north = (xyz(2)-z106)>1.0;
1341 0 : Bool east = (xyz(0)-x106)<1.0;
1342 0 : Int n = Int(floor(0.5l+
1343 0 : sqrt((xyz(0)-x106)*(xyz(0)-x106)+
1344 0 : (xyz(1)-y106)*(xyz(1)-y106)+
1345 0 : (xyz(2)-z106)*(xyz(2)-z106))/incr));
1346 0 : Bool invalid = (n>392);
1347 0 : if (!invalid) {
1348 0 : if (north) {
1349 0 : ostringstream ostr; ostr << "N"<<n;
1350 0 : station = ostr.str();
1351 0 : } else {
1352 0 : if (east) n = 106 -n; else n += 106;
1353 0 : ostringstream ostr; ostr << "W"<<n;
1354 0 : station = ostr.str();
1355 0 : }
1356 : }
1357 0 : return station;
1358 0 : }
1359 :
1360 0 : 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 0 : Regex trailing(" *$"); // trailing blanks
1367 0 : if (if_.if_found) {
1368 0 : spWId_p=-1;
1369 0 : for (Int i=0; i<nSpW_p; i++) {
1370 0 : Double freq = msc_p->spectralWindow().refFrequency()(i);
1371 0 : Double bw = msc_p->spectralWindow().totalBandwidth()(i);
1372 0 : Int nchan = msc_p->spectralWindow().numChan()(i);
1373 0 : Int polId = msc_p->dataDescription().polarizationId()(i);
1374 0 : Int npol = msc_p->polarization().numCorr()(polId);
1375 : // compare freq and bw, cope with case where birdie option has reduced bw
1376 0 : if (doubles_.if_freq[ifNumber]==freq &&
1377 0 : doubles_.if_bw[ifNumber]>=bw && doubles_.if_bw[ifNumber]<2*bw &&
1378 0 : (if_.if_nfreq[ifNumber] == nchan || (nchan<33 &&
1379 0 : if_.if_nfreq[ifNumber]==33)) &&
1380 0 : if_.if_nstok[ifNumber]==npol) {
1381 0 : spWId_p=i;
1382 0 : break;
1383 : }
1384 : }
1385 0 : if (spWId_p<0) { // i.e. not found
1386 0 : spWId_p=nSpW_p++;
1387 0 : atms_p.spectralWindow().addRow();
1388 0 : Int n=atms_p.spectralWindow().nrow()-1;
1389 0 : Double refFreq=doubles_.if_freq[ifNumber];
1390 0 : msc_p->spectralWindow().refFrequency().put(n,refFreq);
1391 : // no doppler tracking
1392 0 : Int nchan = if_.if_nfreq[ifNumber];
1393 0 : Int npol =if_.if_nstok[ifNumber];
1394 :
1395 0 : if (log) os_p << LogIO::NORMAL<<
1396 : "IF "<< ifNumber+1 <<
1397 : " : Frequency = "<< refFreq/1.e6 << " MHz" <<
1398 0 : ", bandwidth = " << doubles_.if_bw[ifNumber]/1.e6 << "MHz" <<
1399 0 : ", #channels = "<< nchan << LogIO::POST;
1400 :
1401 0 : Double refChan=doubles_.if_ref[ifNumber]-1; // make zero based
1402 0 : Double chanSpac = doubles_.if_bw[ifNumber]/max(1,nchan-1);
1403 0 : Int inversion=if_.if_invert[ifNumber];
1404 0 : msc_p->spectralWindow().netSideband().put(n,inversion);
1405 0 : chanSpac*=inversion;
1406 0 : Double chanBW = abs(chanSpac);
1407 0 : if (!cabb_p && nchan==33) chanBW=chanBW*1.6; // roughly
1408 0 : 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 0 : 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 0 : Vector<Double> chanFreq(nchan), resolution(nchan),width(nchan);
1421 0 : for (Int ichan=0; ichan<nchan; ichan++)
1422 0 : chanFreq(ichan)=refFreq+(ichan-refChan)*chanSpac;
1423 0 : msc_p->spectralWindow().chanFreq().put(n,chanFreq);
1424 0 : resolution=chanBW;
1425 0 : width=chanSpac;
1426 0 : msc_p->spectralWindow().resolution().put(n,resolution);
1427 0 : msc_p->spectralWindow().chanWidth().put(n, width);
1428 0 : msc_p->spectralWindow().effectiveBW().put(n, resolution);
1429 0 : msc_p->spectralWindow().totalBandwidth().put(n,
1430 0 : (nchan<33 ? abs(nchan*chanSpac) : doubles_.if_bw[ifNumber]));
1431 0 : msc_p->spectralWindow().numChan().put(n,nchan);
1432 :
1433 0 : Vector<Int> corrType(npol);
1434 0 : Matrix<Int> corrProduct(2,npol); corrProduct=0;
1435 0 : if (log) os_p << LogIO::NORMAL << " : Polarizations ";
1436 0 : for (Int i=0; i<npol; i++) {
1437 0 : corrType(i) = Stokes::type
1438 0 : (String(&names_.if_cstok[2*(i+ifNumber*MaxNPol)],2).before(trailing));
1439 : }
1440 0 : Vector<Int> tmp(npol); tmp=corrType;
1441 : // Sort the polarizations to standard order
1442 0 : GenSort<Int>::sort(corrType);
1443 0 : if (corrIndex_p.nrow()==0) {
1444 0 : corrIndex_p.resize(64,4);
1445 : }
1446 0 : 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 0 : for (Int i=0;i<npol;i++) {
1451 0 : for (Int j=0;j<npol;j++) {
1452 0 : if (corrType(j)==tmp(i)) corrIndex_p(spWId_p,i)=j;
1453 : }
1454 : }
1455 :
1456 : // Figure out the correlation products from the polarizations
1457 0 : corrProduct.resize(2,npol); corrProduct=0;
1458 0 : for (Int i=0; i<npol; i++) {
1459 0 : Stokes::StokesTypes s=Stokes::type(corrType(i));
1460 0 : Fallible<Int> receptor=Stokes::receptor1(s);
1461 0 : if (receptor.isValid()) corrProduct(0,i)=receptor;
1462 0 : receptor=Stokes::receptor2(s);
1463 0 : if (receptor.isValid()) corrProduct(1,i)=receptor;
1464 0 : if (i>0 && log) os_p << " , ";
1465 0 : if (log) os_p << Stokes::name(s)<< " - " << corrType(i);
1466 0 : }
1467 0 : if (log) os_p << LogIO::POST;
1468 :
1469 : // try to find matching pol row
1470 0 : Int nPolRow = atms_p.polarization().nrow();
1471 0 : Int polRow = -1;
1472 0 : for (Int i = 0; i< nPolRow; i++) {
1473 0 : if (msc_p->polarization().numCorr()(i)== Int(if_.if_nstok[ifNumber])) {
1474 0 : if (allEQ(msc_p->polarization().corrType()(i),corrType)) {
1475 0 : polRow=i;
1476 0 : break;
1477 : }
1478 : }
1479 : }
1480 : // add new pol id if needed
1481 0 : if (polRow==-1) {
1482 0 : atms_p.polarization().addRow();
1483 0 : polRow = nPolRow;
1484 0 : msc_p->polarization().numCorr().put(polRow,Int(if_.if_nstok[ifNumber]));
1485 0 : msc_p->polarization().corrType().put(polRow,corrType);
1486 0 : msc_p->polarization().corrProduct().put(polRow,corrProduct);
1487 0 : msc_p->polarization().flagRow().put(polRow, false);
1488 : }
1489 0 : atms_p.dataDescription().addRow();
1490 0 : msc_p->dataDescription().spectralWindowId().put(n, spWId_p);
1491 0 : msc_p->dataDescription().polarizationId().put(n, polRow);
1492 0 : msc_p->dataDescription().flagRow().put(n, false);
1493 :
1494 0 : msc_p->spectralWindow().ifConvChain().put(n,Int(if_.if_chain[ifNumber])-1);
1495 0 : if (!cabb_p) colSamplerBits.put(n,Int(if_.if_sampl[ifNumber]));
1496 :
1497 : // set up the TiledStorageManagers
1498 0 : Record values1, values2, values3, values3c;
1499 0 : values1.define("DATA_HYPERCUBE_ID",spWId_p);
1500 0 : values2.define("SIGMA_HYPERCUBE_ID",spWId_p);
1501 0 : values3.define("FLAG_HYPERCUBE_ID",spWId_p);
1502 0 : values3c.define("FLAG_CATEGORY_HYPERCUBE_ID",spWId_p);
1503 :
1504 0 : Record values4; values4.define("MODEL_HYPERCUBE_ID",spWId_p);
1505 0 : Record values5; values5.define("CORRECTED_HYPERCUBE_ID",spWId_p);
1506 : //Record values6; values6.define("IMAGING_WT_HYPERCUBE_ID",spWId_p);
1507 :
1508 0 : Int nChan=msc_p->spectralWindow().numChan()(spWId_p);
1509 0 : Int nCorr=msc_p->polarization().numCorr()(polRow);
1510 0 : Int nCat=3;
1511 :
1512 : // Choose an appropriate tileshape
1513 0 : IPosition dataShape(2,nCorr,nChan);
1514 0 : IPosition tileShape = MSTileLayout::tileShape(dataShape,obsType_p,"ATCA");
1515 0 : dataAccessor_p.addHypercube(IPosition(3,nCorr,nChan,0),
1516 : tileShape,values1);
1517 0 : sigmaAccessor_p.addHypercube(IPosition(2,nCorr,0),
1518 0 : IPosition(2,tileShape(0),tileShape(2)),
1519 : values2);
1520 0 : flagAccessor_p.addHypercube(IPosition(3,nCorr,nChan,0),
1521 : tileShape,values3);
1522 0 : flagCatAccessor_p.addHypercube(IPosition(4,nCorr,nChan,nCat,0),
1523 0 : IPosition(4,tileShape(0),tileShape(1),
1524 0 : 1,tileShape(2)),values3c);
1525 :
1526 0 : }
1527 : }
1528 0 : return spWId_p;
1529 0 : }
1530 :
1531 :
1532 0 : void ATCAFiller::checkObservation() {
1533 0 : const Regex trailing(" *$"); // trailing blanks
1534 : // Check if current observer already in OBSERVATION table
1535 0 : String observer;
1536 0 : obsId_p=-1;
1537 0 : 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 0 : if (obsId_p<0) {
1545 : // not found, so add it
1546 0 : atms_p.observation().addRow();
1547 0 : obsId_p=atms_p.observation().nrow()-1;
1548 0 : msc_p->observation().observer().put(obsId_p,
1549 0 : String(names_.rp_observer,16).before(trailing));
1550 : // decode project from rpfits file name..
1551 0 : String project=rpfitsFiles_p(0).after(Regex(".*\\."));
1552 0 : if (project.contains(";")) project=project.before(";");
1553 0 : msc_p->observation().project().put(obsId_p,project);
1554 0 : }
1555 0 : }
1556 :
1557 0 : 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 0 : const Regex trailing(" *$"); // trailing blanks
1563 : //0.1 arcsec tolerance for positional coincidence
1564 0 : const Double PosTol=0.1*C::arcsec;
1565 :
1566 0 : if (su_.su_found) {
1567 0 : fieldId_p=-1;
1568 0 : Bool seenSource = false;
1569 : Int numPol;
1570 0 : String name;
1571 0 : String su_name=String(&names_.su_name[sourceno*16],16).before(trailing);
1572 0 : for (Int i=0; i<nField_p; i++) {
1573 0 : msc_p->field().name().get(i,name);
1574 0 : msc_p->field().numPoly().get(i,numPol);
1575 :
1576 0 : if (su_name==name) {
1577 0 : IPosition shape(2, 2, numPol+1);
1578 0 : Matrix<Double> phaseDir(shape);
1579 0 : msc_p->field().phaseDir().get(i,phaseDir);
1580 0 : if (abs(phaseDir(0, 0)-doubles_.su_ra[sourceno])<PosTol) {
1581 0 : if (abs(phaseDir(1, 0)-doubles_.su_dec[sourceno])<PosTol) {
1582 0 : fieldId_p=i; // found it!
1583 : // now check if we've seen this field for this spectral window
1584 0 : Vector<Int> sourceIdCol = msc_p->source().sourceId().getColumn();
1585 0 : Vector<Int> spwIdCol = msc_p->source().spectralWindowId().getColumn();
1586 0 : Vector<Int> spwids = spwIdCol(sourceIdCol==i).getCompressedArray();
1587 0 : seenSource = (spwids.nelements()>0) && anyEQ(spwids,spWId_p);
1588 0 : break;
1589 0 : }
1590 : }
1591 0 : }
1592 : }
1593 :
1594 0 : if (fieldId_p<0 || !seenSource) { // i.e. not found, or not at this spwid
1595 0 : String src = String(&names_.su_name[sourceno*16],16);
1596 :
1597 0 : nsources_p++;
1598 0 : sources_p = src;
1599 :
1600 0 : Double epoch=mjd0_p+Double(proper_.pm_epoch);
1601 0 : Int numPol = 0;
1602 0 : if (abs(proper_.pm_ra)+abs(proper_.pm_dec) > 0) {
1603 0 : numPol = 1;
1604 : }
1605 0 : IPosition shape(2, 2, numPol+1);
1606 0 : Matrix<Double> dir(shape);
1607 : // convert proper motions from s & " per year to rad/s
1608 0 : const Double arcsecPerYear=C::arcsec/(365.25*C::day);
1609 0 : dir(0, 0)=doubles_.su_ra[sourceno];
1610 0 : dir(1, 0)=doubles_.su_dec[sourceno];
1611 0 : 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 0 : if (fieldId_p<0) {
1617 0 : os_p << LogIO::DEBUGGING << "Found field:" << src << LogIO::POST;
1618 0 : if (numPol>0) os_p << LogIO::DEBUGGING << "Field:" << src << " has proper motion parameters"<<LogIO::POST;
1619 0 : fieldId_p=nField_p++;
1620 0 : atms_p.field().addRow();
1621 0 : Int nf=atms_p.field().nrow()-1;
1622 : // for now we have 1 field/source
1623 0 : msc_p->field().sourceId().put(nf,fieldId_p);
1624 0 : msc_p->field().name().put(nf,src.before(trailing));
1625 :
1626 0 : msc_p->field().phaseDir().put(nf,dir);
1627 0 : msc_p->field().delayDir().put(nf,dir);
1628 0 : msc_p->field().referenceDir().put(nf, dir);
1629 0 : msc_p->field().numPoly().put(nf, numPol);
1630 0 : msc_p->field().time().put(nf,epoch);
1631 0 : msc_p->field().code().put(nf,String(&names_.su_cal[sourceno*16],
1632 0 : 16).before(trailing));
1633 : }
1634 0 : dir(0, 0)=doubles_.su_pra[sourceno];
1635 0 : dir(1, 0)=doubles_.su_pdec[sourceno];
1636 0 : Vector<Double> srcdir(2);
1637 0 : srcdir(0)=doubles_.su_ra[sourceno];
1638 0 : srcdir(1)=doubles_.su_dec[sourceno];
1639 0 : atms_p.source().addRow();
1640 0 : Int ns=atms_p.source().nrow()-1;
1641 0 : msc_p->source().sourceId().put(ns,fieldId_p);
1642 0 : msc_p->source().name().put(ns,src.before(trailing));
1643 0 : msc_p->source().direction().put(ns,srcdir);
1644 0 : Vector<Double> rate(2);
1645 0 : rate(0)=proper_.pm_ra*15.*arcsecPerYear; // (15"/s)
1646 0 : rate(1)=proper_.pm_dec*arcsecPerYear;
1647 0 : msc_p->source().properMotion().put(ns,rate);
1648 0 : msc_p->source().time().put(ns,epoch);
1649 0 : 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 0 : msc_p->source().spectralWindowId().put(ns,spWId_p);
1659 0 : Vector<Double> sysv(1),restFreq(1);
1660 0 : Vector<String> transition(1);
1661 0 : sysv = Double(doubles_.vel1);
1662 0 : msc_p->source().sysvel().put(ns, sysv);
1663 0 : msc_p->source().numLines().put(ns,1);
1664 0 : restFreq(0) = Double(doubles_.rfreq);
1665 0 : if (restFreq(0)<1.0) {
1666 : // fill in an appropriate default
1667 0 : Double freq = msc_p->spectralWindow().refFrequency()(spWId_p);
1668 0 : Double bw = msc_p->spectralWindow().totalBandwidth()(spWId_p);
1669 0 : if ( freq>1200.e6 && freq< 1450.e6 && bw <64.e6) {
1670 0 : restFreq(0)= 1420.4e6;
1671 0 : transition(0)="HI";
1672 : } else {
1673 0 : restFreq(0) = freq;
1674 0 : transition(0) = "";
1675 : }
1676 : }
1677 0 : msc_p->source().transition().put(ns,transition);
1678 0 : msc_p->source().restFrequency().put(ns,restFreq);
1679 :
1680 : // dummy fill
1681 0 : msc_p->source().calibrationGroup().put(ns,-1);
1682 0 : }
1683 :
1684 0 : Bool fillPointingTable=False; //Imaging fails with pointing table filled
1685 0 : 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 0 : }
1716 0 : }
1717 0 : }
1718 :
1719 :
1720 0 : void ATCAFiller::storeSysCal()
1721 : {
1722 : // Conversion factor for atmospheric pressure
1723 0 : Vector<String> const pressure_unit = msc_p->weather().pressure().keywordSet().asArrayString("QuantumUnits");
1724 0 : 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 0 : Int last_if_no=-2; // if_no can come out to be -1 if syscal rec is blank
1742 0 : Bool skip=false;
1743 0 : sourceno=sc_.sc_srcno-1; // 0-based source number for this syscal record
1744 0 : Int scq=sc_.sc_q, scif=sc_.sc_if;
1745 0 : for (Int i=0; i<scif; i++) {
1746 0 : for (Int ant=0; ant<sc_.sc_ant; ant++) {
1747 0 : 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 0 : if (sc_.sc_ut!=lastWeatherUT_p) {
1751 0 : lastWeatherUT_p = sc_.sc_ut;
1752 0 : Int nAnt = max(1,sc_.sc_ant-1);
1753 0 : Int row=atms_p.weather().nrow();
1754 0 : atms_p.weather().addRow(nAnt);
1755 0 : Double time=mjd0_p+Double(sc_.sc_ut);
1756 0 : for (Int iAnt=0; iAnt<nAnt; iAnt++,row++) {
1757 0 : msc_p->weather().antennaId().put(row,iAnt);
1758 0 : msc_p->weather().time().put(row,time);
1759 0 : msc_p->weather().interval().put(row,Double(param_.intbase));
1760 0 : msc_p->weather().temperature().put(row,
1761 0 : Double(sc_.sc_cal[scq*(i+scif*ant)+1])+273.15); // C to K
1762 0 : msc_p->weather().pressure().put(row,
1763 0 : Double(sc_.sc_cal[scq*(i+scif*ant)+2])*pressure_conversion); // mBar to Pa/hPa
1764 0 : msc_p->weather().relHumidity().put(row,
1765 0 : Double(sc_.sc_cal[scq*(i+scif*ant)+3]));
1766 0 : msc_p->weather().windSpeed().put(row,
1767 0 : Double(sc_.sc_cal[scq*(i+scif*ant)+4])/3.6); // km/s to m/s
1768 0 : msc_p->weather().windDirection().put(row,
1769 0 : Double(sc_.sc_cal[scq*(i+scif*ant)+5])*C::pi/180.0); // deg to rad
1770 0 : msc_p->weather().temperatureFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1771 0 : msc_p->weather().pressureFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1772 0 : msc_p->weather().relHumidityFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1773 0 : msc_p->weather().windSpeedFlag().put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1774 0 : msc_p->weather().windDirectionFlag().put(row,Bool(sc_.sc_cal[scq*(i+scif*ant)+6]));
1775 0 : colWeatherRainGauge.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+7]));
1776 0 : colWeatherSeeMonPhase.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+8]));
1777 0 : colWeatherSeeMonRMS.put(row, Float(sc_.sc_cal[scq*(i+scif*ant)+9])/1000);
1778 0 : colWeatherSeeMonFlag.put(row, Bool(sc_.sc_cal[scq*(i+scif*ant)+10]));
1779 : }
1780 : }
1781 0 : continue;
1782 0 : }
1783 0 : if_no=Int(sc_.sc_cal[scq*(i+scif*ant)+1])-1; // make 0-based
1784 0 : if (if_no!=last_if_no) { // check if we want this one
1785 0 : last_if_no=if_no;
1786 0 : skip=Bool(if_no<0 || if_no>=if_.n_if || !selected(if_no));
1787 0 : if (!skip) {
1788 0 : checkSpW(if_no); // sets spWId_p corresponding to this if_no
1789 0 : checkField(); // sets fieldId_p
1790 : }
1791 : }
1792 0 : if (!skip) {
1793 0 : atms_p.sysCal().addRow();
1794 0 : Int n=atms_p.sysCal().nrow()-1;
1795 0 : Double time=mjd0_p+Double(sc_.sc_ut);
1796 0 : msc_p->sysCal().time().put(n,time);
1797 0 : msc_p->sysCal().antennaId()
1798 0 : .put(n,Int(sc_.sc_cal[scq*(i+scif*ant)+0]-1)); // make 0-based
1799 0 : msc_p->sysCal().feedId().put(n,0);
1800 0 : msc_p->sysCal().interval().put(n,Double(param_.intbase));
1801 0 : msc_p->sysCal().spectralWindowId().put(n,spWId_p);
1802 :
1803 0 : msc_p->sysCal().phaseDiff().put(n,sc_.sc_cal[scq*(i+scif*ant)+2]);
1804 :
1805 0 : Vector<Float> tSys(2);
1806 : // convert from sqrt(tsys*10) to true tsys
1807 0 : tSys(0)=square(sc_.sc_cal[scq*(i+scif*ant)+3])/10.;
1808 0 : tSys(1)=square(sc_.sc_cal[scq*(i+scif*ant)+4])/10.;
1809 :
1810 0 : msc_p->sysCal().tsys().put(n,tSys);
1811 :
1812 0 : Vector<Float> a(2),b(2),c(2);
1813 0 : for (Int j=0; j<2; j++) {
1814 0 : a(j)=sc_.sc_cal[scq*(i+scif*ant)+5+j*3];
1815 0 : b(j)=sc_.sc_cal[scq*(i+scif*ant)+6+j*3];
1816 0 : c(j)=sc_.sc_cal[scq*(i+scif*ant)+7+j*3];
1817 : }
1818 0 : if (cabb_p) {
1819 0 : colGTP.put(n,a);
1820 0 : colSDO.put(n,b);
1821 0 : 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 0 : msc_p->sysCal().phaseDiffFlag().
1829 0 : put(n,(sc_.sc_cal[scq*(i+scif*ant)+12]!=0));
1830 0 : msc_p->sysCal().tsysFlag().
1831 0 : put(n,(sc_.sc_cal[scq*(i+scif*ant)+12]!=0));
1832 0 : colXYAmplitude.put(n,sc_.sc_cal[scq*(i+scif*ant)+13]);
1833 0 : msc_p->sysCal().tcalFlag().put(n,true);
1834 0 : msc_p->sysCal().trxFlag().put(n,true);
1835 0 : colTrackErrMax.put(n,Float(sc_.sc_cal[scq*(i+scif*ant)+14]));
1836 0 : colTrackErrRMS.put(n,Float(sc_.sc_cal[scq*(i+scif*ant)+15]));
1837 0 : }
1838 : }
1839 : }
1840 0 : }
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 0 : void ATCAFiller::storeData()
1893 : {
1894 0 : const double MJD01Jul2004 = 2453187.5; // 12mm receiver xy inversion end
1895 0 : const double MJD18Oct2007 = 2454390.5; // 3mm CA02 xyphase valid
1896 :
1897 0 : atms_p.addRow();
1898 :
1899 0 : const RecordFieldId rfid1("DATA_HYPERCUBE_ID");
1900 0 : const RecordFieldId rfid2("SIGMA_HYPERCUBE_ID");
1901 0 : const RecordFieldId rfid3("FLAG_HYPERCUBE_ID");
1902 0 : const RecordFieldId rfid3c("FLAG_CATEGORY_HYPERCUBE_ID");
1903 0 : Record values1, values2, values3, values3c;
1904 0 : values1.define(rfid1,spWId_p);
1905 0 : values2.define(rfid2,spWId_p);
1906 0 : values3.define(rfid3,spWId_p);
1907 0 : values3c.define(rfid3c,spWId_p);
1908 :
1909 0 : dataAccessor_p.extendHypercube(1,values1);
1910 0 : sigmaAccessor_p.extendHypercube(1,values2);
1911 0 : flagAccessor_p.extendHypercube(1,values3);
1912 0 : flagCatAccessor_p.extendHypercube(1,values3c);
1913 :
1914 0 : Record values4, values5, values6;
1915 :
1916 0 : Int n=atms_p.nrow()-1;
1917 0 : if (n==0) gotSysCalId_p=false;
1918 :
1919 0 : Int npol =if_.if_nstok[if_no];
1920 0 : Int nfreq=if_.if_nfreq[if_no];
1921 0 : Regex trailing(" *$"); // trailing blanks
1922 0 : String instrument = String(names_.instrument,16).before(trailing);
1923 0 : Bool atca = (instrument=="ATCA");
1924 0 : Bool atlba = (instrument=="ATLBA");
1925 :
1926 : // fill in the easy items
1927 : // make antenna numbers 0-based
1928 0 : Int ant1=Int(baseline)/256-1, ant2=Int(baseline)%256-1;
1929 0 : msc_p->antenna1().put(n,ant1);
1930 0 : msc_p->antenna2().put(n,ant2);
1931 0 : Bool flagData=flg;
1932 0 : if (flagData) { flagCount_p(ONLINE)++;}
1933 0 : if (autoFlag_p && flagScanType_p) { flagData=true; flagCount_p(SCANTYPE)++;}
1934 0 : Double exposure=Double(param_.intbase);
1935 : // averaging = number of integration periods averaged by correlator: 1,2 or 3
1936 0 : Int averaging = 1;
1937 0 : if (param_.intime>0) averaging = Int(exposure/param_.intime)+1;
1938 0 : Double interval = averaging*Double(param_.intime);
1939 : // check for old data with intbase set to zero
1940 0 : Double blank=51.e-3; // standard 51 ms blank time at end of integration
1941 0 : if (exposure<0.001) exposure = interval-blank;
1942 0 : if (interval==0.0) interval = exposure;
1943 :
1944 : // Is binning mode active?
1945 0 : Int nBin = Int(floor(param_.intime/exposure+0.01));
1946 0 : if (nBin<4) nBin=1; // short exposure due to long hold period..
1947 0 : msc_p->dataDescId().put(n, spWId_p);
1948 0 : if (ut!=lastUT_p) { // the ISM columns that don't change within integration
1949 0 : msc_p->arrayId().put(n,0);
1950 0 : msc_p->exposure().put(n,exposure);
1951 0 : msc_p->feed1().put(n,0);
1952 0 : msc_p->feed2().put(n,0);
1953 0 : msc_p->fieldId().put(n,fieldId_p);
1954 0 : msc_p->interval().put(n,interval);
1955 0 : msc_p->observationId().put(n,obsId_p);
1956 0 : msc_p->processorId().put(n,-1);
1957 0 : msc_p->scanNumber().put(n,scanNo_p);
1958 0 : msc_p->stateId().put(n,-1);
1959 0 : Double mjd=mjd0_p+ut;
1960 : // try to figure out what the midpoint of the integration interval was
1961 0 : if (nBin>1 || !atca) {
1962 0 : 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 0 : msc_p->timeCentroid().put(n,mjd);
1972 : }
1973 0 : Int pulsarBin = max(0,bin-1);// make zero-based (but catch unset value)
1974 0 : msc_p->pulsarBin().put(n,pulsarBin);
1975 0 : 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 0 : Vector<Double> uvw(3); uvw(0)=u; uvw(1)=v; uvw(2)=w;
1983 0 : msc_p->uvw().put(n,uvw);
1984 :
1985 : // use exposure time as weight
1986 0 : Vector<Float> Weight(npol); Weight=exposure;
1987 0 : 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 0 : if (!gotSysCalId_p || ut!=lastUT_p || spWId_p!=lastSpWId_p) {
1993 0 : lastUT_p=ut;
1994 0 : lastSpWId_p=spWId_p;
1995 : // we need to get the new syscal info
1996 : // search backwards from last syscal row
1997 0 : Bool done=false;
1998 0 : Int nsc=atms_p.sysCal().nrow();
1999 : //# os_p << LogIO::NORMAL << " Current number of syscal rows="<<nsc<<LogIO::POST;
2000 0 : sysCalId_p.resize(nAnt_p); sysCalId_p=-1;
2001 : // search backwards, as it should be a recent entry
2002 0 : for (Int i=nsc-1; (i>=0 && !done); i--) {
2003 0 : if (nearAbs(msc_p->sysCal().time()(i),mjd0_p+ut,0.5) &&
2004 0 : msc_p->sysCal().spectralWindowId()(i)==spWId_p) {
2005 0 : Int ant=msc_p->sysCal().antennaId()(i);
2006 0 : if (ant>=0 && ant<nAnt_p) { //valid antenna
2007 0 : sysCalId_p(ant)=i;
2008 0 : done=(!(anyEQ(sysCalId_p,-1)));
2009 : }
2010 : }
2011 : }
2012 0 : 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 0 : colSysCalIdAnt1.put(n,sysCalId_p(ant1));
2020 0 : colSysCalIdAnt2.put(n,sysCalId_p(ant2));
2021 :
2022 : // Check for bad sampler stats now that we've found the syscal data
2023 0 : if (!flagData && autoFlag_p && (!atlba && !cabb_p && samplerFlag(n))) {
2024 0 : flagData=true; flagCount_p(SYSCAL)++;
2025 : }
2026 0 : flagCount_p(COUNT)++;
2027 0 : if (flagData) flagCount_p(FLAG)++;
2028 0 : msc_p->flagRow().put(n,flagData);
2029 :
2030 :
2031 : // flags and data
2032 0 : const Int nCat = 3; // three initial categories
2033 : // define the categories
2034 0 : Vector<String> cat(nCat);
2035 0 : cat(0)="FLAG_CMD";
2036 0 : cat(1)="ORIGINAL";
2037 0 : cat(2)="USER";
2038 0 : msc_p->flagCategory().rwKeywordSet().define("CATEGORY",cat);
2039 :
2040 : // Gibbs reweighting
2041 0 : if (nfreq == 33) {
2042 0 : 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 0 : Double refFreq=doubles_.if_freq[if_no];
2048 0 : Double refChan=doubles_.if_ref[if_no]-1; // make zero based
2049 0 : Double chanSpac = doubles_.if_bw[if_no]/max(1,nfreq-1);
2050 0 : Int inversion=if_.if_invert[if_no];
2051 0 : chanSpac*=inversion;
2052 0 : Int bchan = -1;
2053 0 : Int edge = 0;
2054 0 : 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 0 : Bool band12mm = (refFreq>13.e9 && refFreq<28.e9);
2065 0 : Bool band3mm = (refFreq>75.e9);
2066 :
2067 0 : Matrix<Complex> VIS(npol,nfreq);
2068 0 : Cube<Bool> flagCat(npol,nfreq,nCat,false);
2069 0 : Matrix<Bool> flags= flagCat.xyPlane(1); // references flagCat's storage
2070 :
2071 :
2072 0 : 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 0 : for (Int j=0; j<nfreq; j++) {
2095 0 : Int k = (edge>0 ? (2*j+edge) : j);
2096 0 : for (Int i=0; i<npol; i++) {
2097 0 : Int ipol = corrIndex_p(spWId_p,i);
2098 0 : VIS(ipol,j)=Complex(vis[0+2*(i+k*npol)],inversion*vis[1+2*(i+k*npol)]);
2099 : }
2100 : }
2101 : // Flag NaNs
2102 0 : Matrix<Bool> nanFlags = isNaN(real(VIS));
2103 0 : nanFlags |= isNaN(imag(VIS));
2104 0 : flags |= nanFlags;
2105 :
2106 : // Flag CABB birdies and RFI
2107 0 : if (cabb_p) rfiFlag(flags);
2108 :
2109 : // correct for xy phase - multiply y correlation with opposite of xyphase
2110 0 : Int id1 = sysCalId_p(ant1);
2111 0 : Complex gain1 = (id1!=-1 && !msc_p->sysCal().phaseDiffFlag()(id1) ?
2112 0 : exp(Complex(0,1)*msc_p->sysCal().phaseDiff()(id1)*inversion):
2113 0 : Complex(1,0));
2114 : // negate xyphase for second antenna since the correlation product uses conjugate for ant2
2115 0 : Int id2 = sysCalId_p(ant2);
2116 0 : Complex gain2 = (id2!=-1 && !msc_p->sysCal().phaseDiffFlag()(id2) ?
2117 0 : exp(-Complex(0,1)*msc_p->sysCal().phaseDiff()(id2)*inversion):
2118 0 : Complex(1,0));
2119 0 : Int polId = msc_p->dataDescription().polarizationId()(spWId_p);
2120 0 : Vector<Int> corrType = msc_p->polarization().corrType()(polId);
2121 : // 3mm receiver on CA02 got valid xyphase on 18Oct2007
2122 0 : if (noxycorr_p || (band3mm && mjd0_p<MJD18Oct2007)) { gain1=gain2=1; }
2123 0 : 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 0 : if (!band12mm || mjd0_p>MJD01Jul2004) { gain1=-gain1; gain2=-gain2;}
2129 0 : for (Int i=0; i<npol; i++) {
2130 0 : switch (corrType(i)) {
2131 0 : case Stokes::XX :
2132 0 : break;
2133 0 : case Stokes::XY :
2134 : {
2135 0 : Vector<Complex> tmp(VIS.row(i));
2136 0 : tmp*=gain2;
2137 0 : }
2138 0 : break;
2139 0 : case Stokes::YX :
2140 : {
2141 0 : Vector<Complex> tmp(VIS.row(i));
2142 0 : tmp*=gain1;
2143 0 : }
2144 0 : break;
2145 0 : case Stokes::YY :
2146 : {
2147 0 : Vector<Complex> tmp(VIS.row(i));
2148 0 : tmp*=(gain1*gain2);
2149 0 : }
2150 0 : break;
2151 0 : default:
2152 0 : break;
2153 : }
2154 : }
2155 :
2156 0 : msc_p->data().put(n,VIS);
2157 : // CASA wants all flags set if flagRow is set
2158 0 : if (flagData) flags=flagData;
2159 0 : msc_p->flag().put(n,flags);
2160 0 : msc_p->flagCategory().put(n,flagCat);
2161 :
2162 : // now calculate the sigma for this row
2163 : //
2164 0 : Vector<Double> chnbw;
2165 0 : msc_p->spectralWindow().resolution().get(spWId_p,chnbw);
2166 0 : 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 0 : if (sysCalId_p(ant1)!=-1 && !msc_p->sysCal().tsysFlag()(sysCalId_p(ant1)))
2171 0 : msc_p->sysCal().tsys().get(sysCalId_p(ant1),tsys1);
2172 0 : else { tsys1.resize(nAnt_p); tsys1=50.;}
2173 0 : if (sysCalId_p(ant2)!=-1 && !msc_p->sysCal().tsysFlag()(sysCalId_p(ant2)))
2174 0 : msc_p->sysCal().tsys().get(sysCalId_p(ant2),tsys2);
2175 0 : else { tsys2.resize(nAnt_p); tsys2=50.;}
2176 0 : Vector<Float> sigma(npol); sigma=0.;
2177 0 : Matrix<Int> corrProduct(2,npol);
2178 0 : msc_p->polarization().corrProduct().get(polId,corrProduct);
2179 :
2180 : // sigma=sqrt(Tx1*Tx2)/sqrt(chnbw*intTime)*JyPerK;
2181 0 : Float JyPerK=13.; // guess for ATCA dishes at 3-20cm
2182 0 : Float factor=sqrt(chnbw(0)*exposure)/JyPerK;
2183 0 : for (Int pol=0; pol<npol; pol++) {
2184 0 : Float tsysAv=tsys1(corrProduct(0,pol))*tsys2(corrProduct(1,pol));
2185 0 : if (tsysAv>=0 && factor>0 ) sigma(pol)=sqrt(tsysAv)/factor;
2186 : }
2187 0 : msc_p->sigma().put(n,sigma);
2188 :
2189 : // determine shadowing & flag data with shadowed antennas
2190 0 : if (shadow_p>0) shadow(n);
2191 :
2192 0 : }
2193 :
2194 0 : void ATCAFiller::rfiFlag(Matrix<Bool> & flags) {
2195 : // CABB birdies
2196 : // 1 MHz mode
2197 0 : 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 0 : const int nb2=3;
2201 : static const int b2[nb2] = {8,16,24};
2202 :
2203 : // Get details of spectrum
2204 0 : Int nfreq=if_.if_nfreq[if_no];
2205 0 : Double bw=doubles_.if_bw[if_no];
2206 0 : Int chn=0;
2207 : // Flag birdies and edge channels
2208 0 : if (bw>2.e9) {
2209 0 : chn = nfreq*edge_p/200;
2210 0 : if (nfreq==2049) {
2211 : // CABB 2049 chan * 1 MHz continuum mode
2212 0 : for (Int i=0; i<nb1; i++) {
2213 0 : flags.column(b1[i])=true;
2214 : }
2215 0 : } else if (nfreq==33) {
2216 : // CABB 33 chan * 64 MHz continuum mode
2217 0 : for (Int i=0; i<nb2; i++) {
2218 0 : flags.column(b2[i])=true;
2219 : }
2220 : }
2221 0 : } else if (bw<1.e9 && nfreq>=2049) {
2222 : // CABB zoom mode, 2049 or more channels (combined zooms)
2223 0 : chn = 2049*edge_p/200;
2224 : }
2225 0 : for (Int i=0; i<chn; i++) flags.column(i)=true;
2226 0 : for (Int i=nfreq-chn; i<nfreq; i++) flags.column(i)=true;
2227 0 : }
2228 :
2229 0 : void ATCAFiller::fillFeedTable() {
2230 : // At present there is always only a single feed per antenna (at any
2231 : // given spectralwindow).
2232 0 : Int nAnt=atms_p.antenna().nrow();
2233 : // Int nSpW=atms_p.spectralWindow().nrow();
2234 : // Only X and Y receptors available
2235 0 : Vector<String> rec_type(2); rec_type(0)="X"; rec_type(1)="Y";
2236 0 : Matrix<Complex> polResponse(2,2);
2237 0 : polResponse=0.; polResponse(0,0)=polResponse(1,1)=1.;
2238 0 : Matrix<Double> offset(2,2); offset=0.;
2239 0 : Vector<Double> position(3); position=0.;
2240 : // X feed is at 45 degrees, Y at 135 degrees for all except 7mm
2241 0 : Vector<Double> receptorAngle(2);
2242 0 : receptorAngle(0)=45*C::degree;
2243 0 : receptorAngle(1)=135*C::degree;
2244 : // Single entry, so we use the first spectral window
2245 0 : Double f = msc_p->spectralWindow().refFrequency()(0);
2246 0 : if (f>30e9 && f<50e9) receptorAngle+=90*C::degree;
2247 :
2248 : // fill the feed table
2249 0 : 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 0 : for (Int ant=0; ant<nAnt; ant++) {
2254 0 : atms_p.feed().addRow(); row++;
2255 0 : msc_p->feed().antennaId().put(row,ant);
2256 0 : msc_p->feed().beamId().put(row,-1);
2257 0 : msc_p->feed().feedId().put(row,0);
2258 0 : msc_p->feed().interval().put(row,DBL_MAX);
2259 0 : msc_p->feed().phasedFeedId().put(row,-1);
2260 0 : msc_p->feed().spectralWindowId().put(row,-1); //spw);
2261 0 : msc_p->feed().time().put(row,0.);
2262 0 : msc_p->feed().numReceptors().put(row,2);
2263 0 : msc_p->feed().beamOffset().put(row,offset);
2264 0 : msc_p->feed().polarizationType().put(row,rec_type);
2265 0 : msc_p->feed().polResponse().put(row,polResponse);
2266 0 : msc_p->feed().position().put(row,position);
2267 0 : msc_p->feed().receptorAngle().put(row,receptorAngle);
2268 : }
2269 0 : }
2270 :
2271 0 : void ATCAFiller::fillObservationTable()
2272 : {
2273 0 : Vector<Double> tim = msc_p->time().getColumn();
2274 0 : Vector<Int> obsid = msc_p->observationId().getColumn();
2275 : Int startInd;
2276 : Int endInd;
2277 :
2278 0 : Int nObs=atms_p.observation().nrow();
2279 0 : Vector<Double> vt(2);
2280 0 : for (Int obs=0; obs<nObs; obs++) {
2281 : // fill time range
2282 0 : startInd = 0;
2283 0 : endInd = atms_p.nrow()-1;
2284 0 : for (uInt i=0; i<atms_p.nrow(); i++) {
2285 0 : if (obsid(i) == obs) { startInd = i; break; }
2286 : }
2287 0 : for (uInt i=startInd; i<atms_p.nrow(); i++) {
2288 0 : if (obsid(i) > obs) { endInd = i-1; break; }
2289 : }
2290 0 : vt(0) = tim(startInd);
2291 0 : vt(1) = tim(endInd);
2292 0 : msc_p->observation().timeRange().put(obs, vt);
2293 0 : msc_p->observation().releaseDate().put(obs,vt(1)+18*30.5*86400);
2294 : // telescope name
2295 0 : msc_p->observation().telescopeName().put(obs, "ATCA");
2296 : }
2297 0 : }
2298 :
2299 0 : void ATCAFiller::flush()
2300 : {
2301 0 : atms_p.flush();
2302 0 : atms_p.antenna().flush();
2303 0 : atms_p.dataDescription().flush();
2304 0 : atms_p.feed().flush();
2305 0 : atms_p.field().flush();
2306 0 : atms_p.observation().flush();
2307 0 : atms_p.history().flush();
2308 0 : atms_p.pointing().flush();
2309 0 : atms_p.polarization().flush();
2310 0 : atms_p.source().flush();
2311 0 : atms_p.spectralWindow().flush();
2312 0 : atms_p.sysCal().flush();
2313 0 : atms_p.weather().flush();
2314 0 : }
2315 :
2316 0 : void ATCAFiller::unlock()
2317 : {
2318 0 : atms_p.unlock();
2319 0 : atms_p.antenna().unlock();
2320 0 : atms_p.dataDescription().unlock();
2321 0 : atms_p.feed().unlock();
2322 0 : atms_p.field().unlock();
2323 0 : atms_p.observation().unlock();
2324 0 : atms_p.history().unlock();
2325 0 : atms_p.pointing().unlock();
2326 0 : atms_p.polarization().unlock();
2327 0 : atms_p.source().unlock();
2328 0 : atms_p.spectralWindow().unlock();
2329 0 : atms_p.sysCal().unlock();
2330 0 : atms_p.weather().unlock();
2331 0 : }
2332 :
2333 0 : Bool ATCAFiller::selected(Int ifNum)
2334 : {
2335 0 : Bool select=true;
2336 : //Check for rotten eggs
2337 0 : if (if_.if_nfreq[ifNum]==0)
2338 0 : select=false;
2339 : // check if we want this frequency
2340 0 : 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 0 : if (lowFreq_p>0 && (doubles_.if_freq[ifNum]-lowFreq_p)<0)
2345 0 : select=false;
2346 : else {
2347 0 : 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 0 : 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 0 : if (numChan1_p>0 && numChan1_p != if_.if_nfreq[ifNum])
2357 0 : select=false;
2358 : }
2359 : }
2360 : }
2361 : }
2362 0 : 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 0 : }
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 0 : }
2444 0 : return flag;
2445 0 : }
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 0 : }
2456 :
2457 : } //# end casa namespace
|
