Line data Source code
1 : /*
2 : * NRO2MSReader.cc
3 : *
4 : * Created on: May 9, 2016
5 : * Author: wataru kawasaki
6 : */
7 :
8 : #include <singledishfiller/Filler/NRO2MSReader.h>
9 :
10 : #include <iostream>
11 : #include <string>
12 : #include <map>
13 :
14 : #include <casacore/casa/OS/File.h>
15 : #include <casacore/casa/Containers/Record.h>
16 : #include <casacore/casa/Utilities/Regex.h>
17 : #include <casacore/tables/Tables/TableRecord.h>
18 : #include <casacore/casa/Arrays/Vector.h>
19 : #include <casacore/casa/Arrays/ArrayMath.h>
20 : #include <casacore/casa/IO/ArrayIO.h>
21 : #include <casacore/tables/Tables/Table.h>
22 :
23 : using namespace casacore;
24 :
25 : namespace {
26 0 : Double queryAntennaDiameter(String const &name) {
27 0 : String capitalized = name;
28 0 : capitalized.upcase();
29 0 : Double diameter = 0.0;
30 0 : if (capitalized.matches(Regex(".*(DV|DA|PM)[0-9]+$"))) {
31 0 : diameter = 12.0;
32 0 : } else if (capitalized.matches(Regex(".*CM[0-9]+$"))) {
33 0 : diameter = 7.0;
34 0 : } else if (capitalized.contains("GBT")) {
35 0 : diameter = 104.9;
36 0 : } else if (capitalized.contains("MOPRA")) {
37 0 : diameter = 22.0;
38 0 : } else if (capitalized.contains("PKS") || capitalized.contains("PARKS")) {
39 0 : diameter = 64.0;
40 0 : } else if (capitalized.contains("TIDBINBILLA")) {
41 0 : diameter = 70.0;
42 0 : } else if (capitalized.contains("CEDUNA")) {
43 0 : diameter = 30.0;
44 0 : } else if (capitalized.contains("HOBART")) {
45 0 : diameter = 26.0;
46 0 : } else if (capitalized.contains("APEX")) {
47 0 : diameter = 12.0;
48 0 : } else if (capitalized.contains("ASTE")) {
49 0 : diameter = 10.0;
50 0 : } else if (capitalized.contains("NRO")) {
51 0 : diameter = 45.0;
52 : }
53 :
54 0 : return diameter;
55 : }
56 :
57 : template<class T, class U>
58 0 : U getMapValue(std::map<T, U> const mymap, T const key, U const default_value) {
59 0 : auto iter = mymap.find(key);
60 0 : if (iter != mymap.end()) {
61 0 : return iter->second;
62 : } else {
63 0 : return default_value;
64 : }
65 : }
66 :
67 0 : String getIntent(Int srctype) {
68 0 : static std::map<Int, String> intent_map;
69 0 : if (intent_map.size() == 0) {
70 0 : String sep1 = "#";
71 0 : String sep2 = ",";
72 0 : String target = "OBSERVE_TARGET";
73 0 : String atmcal = "CALIBRATE_ATMOSPHERE";
74 0 : String anycal = "CALIBRATE_SOMETHING";
75 0 : String onstr = "ON_SOURCE";
76 0 : String offstr = "OFF_SOURCE";
77 0 : String pswitch = "POSITION_SWITCH";
78 0 : String nod = "NOD";
79 0 : String fswitch = "FREQUENCY_SWITCH";
80 0 : String sigstr = "SIG";
81 0 : String refstr = "REF";
82 0 : String hot = "HOT";
83 0 : String warm = "WARM";
84 0 : String cold = "COLD";
85 0 : String unspecified = "UNSPECIFIED";
86 0 : String ftlow = "LOWER";
87 0 : String fthigh = "HIGHER";
88 0 : intent_map[0] = target + sep1 + onstr + sep2 + pswitch;
89 0 : intent_map[1] = target + sep1 + offstr + sep2 + pswitch;
90 0 : intent_map[2] = target + sep1 + onstr + sep2 + nod;
91 0 : intent_map[3] = target + sep1 + onstr + sep2 + fswitch + sep1 + sigstr;
92 0 : intent_map[4] = target + sep1 + onstr + sep2 + fswitch + sep1 + refstr;
93 0 : intent_map[6] = atmcal + sep1 + offstr + sep2 + unspecified;
94 0 : intent_map[7] = atmcal + sep1 + hot + sep2 + unspecified;
95 0 : intent_map[8] = atmcal + sep1 + warm + sep2 + unspecified;
96 0 : intent_map[9] = atmcal + sep1 + cold + sep2 + unspecified;
97 0 : intent_map[10] = atmcal + sep1 + onstr + sep2 + pswitch;
98 0 : intent_map[11] = atmcal + sep1 + offstr + sep2 + pswitch;
99 0 : intent_map[12] = atmcal + sep1 + onstr + sep2 + nod;
100 0 : intent_map[13] = atmcal + sep1 + onstr + sep2 + fswitch + sep1 + sigstr;
101 0 : intent_map[14] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + refstr;
102 0 : intent_map[20] = target + sep1 + onstr + sep2 + fswitch + sep1 + ftlow;
103 0 : intent_map[21] = target + sep1 + offstr + sep2 + fswitch + sep1 + ftlow;
104 0 : intent_map[26] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + ftlow;
105 0 : intent_map[27] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + ftlow;
106 0 : intent_map[28] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + ftlow;
107 0 : intent_map[29] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + ftlow;
108 0 : intent_map[30] = target + sep1 + onstr + sep2 + fswitch + sep1 + fthigh;
109 0 : intent_map[31] = target + sep1 + offstr + sep2 + fswitch + sep1 + fthigh;
110 0 : intent_map[36] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + fthigh;
111 0 : intent_map[37] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + fthigh;
112 0 : intent_map[38] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + fthigh;
113 0 : intent_map[39] = atmcal + sep1 + offstr + sep2 + fswitch + sep1 + fthigh;
114 0 : intent_map[90] = target + sep1 + onstr + sep2 + unspecified;
115 0 : intent_map[91] = target + sep1 + offstr + sep2 + unspecified;
116 0 : intent_map[92] = anycal + sep1 + offstr + sep2 + unspecified;
117 : }
118 0 : String default_type = "UNKNOWN_INTENT";
119 0 : return getMapValue(intent_map, srctype, default_type);
120 : }
121 :
122 0 : Int getSubscan(Int srctype) {
123 0 : static std::map<Int, Int> subscan_map;
124 0 : if (subscan_map.size() == 0) {
125 0 : subscan_map[0] = 1;
126 0 : subscan_map[1] = 2;
127 0 : subscan_map[2] = 1;
128 0 : subscan_map[3] = 1;
129 0 : subscan_map[4] = 2;
130 0 : subscan_map[6] = 1;
131 0 : subscan_map[7] = 2;
132 0 : subscan_map[8] = 3;
133 0 : subscan_map[9] = 4;
134 0 : subscan_map[10] = 5;
135 0 : subscan_map[11] = 6;
136 0 : subscan_map[12] = 7;
137 0 : subscan_map[13] = 8;
138 0 : subscan_map[14] = 9;
139 0 : subscan_map[20] = 1;
140 0 : subscan_map[21] = 2;
141 0 : subscan_map[26] = 1;
142 0 : subscan_map[27] = 2;
143 0 : subscan_map[28] = 3;
144 0 : subscan_map[29] = 4;
145 0 : subscan_map[30] = 3;
146 0 : subscan_map[31] = 4;
147 0 : subscan_map[36] = 5;
148 0 : subscan_map[37] = 6;
149 0 : subscan_map[38] = 7;
150 0 : subscan_map[39] = 8;
151 0 : subscan_map[90] = 1;
152 0 : subscan_map[91] = 2;
153 0 : subscan_map[92] = 1;
154 : }
155 0 : Int default_subscan = 1;
156 0 : return getMapValue(subscan_map, srctype, default_subscan);
157 : }
158 :
159 : constexpr double kDay2Sec = 86400.0;
160 : // constexpr double kSec2Day = 1.0 / kDay2Sec;
161 :
162 : // CAS-11223
163 : // Time difference between JST and UTC is 9 hours
164 : constexpr double kJSTOffsetHour = 9.0;
165 : constexpr double kJSTOffsetMin = kJSTOffsetHour * 60.0;
166 : constexpr double kJSTOffsetSec = kJSTOffsetMin * 60.0;
167 : //inline double jst2utcHour(double jst_time) {
168 : // return jst_time - kJSTOffsetHour;
169 : //}
170 0 : inline double jst2utcSec(double jst_time) {
171 0 : return jst_time - kJSTOffsetSec;
172 : }
173 : }
174 :
175 : using namespace casacore;
176 : namespace casa { //# NAMESPACE CASA - BEGIN
177 :
178 : using namespace sdfiller;
179 :
180 0 : NRO2MSReader::NRO2MSReader(std::string const &scantable_name) :
181 : ReaderInterface(scantable_name), fp_(NULL), obs_header_(),
182 : beam_id_counter_(0), source_spw_id_counter_(0), spw_id_counter_(0),
183 : time_range_sec_(),
184 0 : get_antenna_row_([&](sdfiller::AntennaRecord &r) {return NRO2MSReader::getAntennaRowImpl(r);}),
185 0 : get_field_row_([&](sdfiller::FieldRecord &r) {return NRO2MSReader::getFieldRowImpl(r);}),
186 0 : get_observation_row_([&](sdfiller::ObservationRecord &r) {return NRO2MSReader::getObservationRowImpl(r);}),
187 0 : get_processor_row_([&](sdfiller::ProcessorRecord &r) {return NRO2MSReader::getProcessorRowImpl(r);}),
188 0 : get_source_row_([&](sdfiller::SourceRecord &r) {return NRO2MSReader::getSourceRowImpl(r);}),
189 0 : get_spw_row_([&](sdfiller::SpectralWindowRecord &r) {return NRO2MSReader::getSpectralWindowRowImpl(r);})
190 : {
191 : // std::cout << "NRO2MSReader::NRO2MSReader" << std::endl;
192 0 : }
193 :
194 0 : NRO2MSReader::~NRO2MSReader() {
195 : // std::cout << "NRO2MSReader::~NRO2MSReader" << std::endl;
196 : // to ensure file pointer is closed
197 0 : finalizeSpecific();
198 0 : }
199 :
200 0 : void NRO2MSReader::checkEndian() {
201 0 : fseek(fp_, 144, SEEK_SET) ;
202 : int tmp ;
203 0 : if (fread(&tmp, 1, sizeof(int), fp_) != sizeof(int)) {
204 0 : return ;
205 : }
206 0 : if ((0 < tmp)&&(tmp <= NRO_ARYMAX)) {
207 0 : same_endian_ = true;
208 : } else {
209 0 : same_endian_ = false;
210 : }
211 : }
212 :
213 0 : void NRO2MSReader::readObsHeader() {
214 0 : fseek(fp_, 0, SEEK_SET);
215 :
216 0 : readHeader(obs_header_.LOFIL0, 8);
217 0 : readHeader(obs_header_.VER0, 8);
218 0 : readHeader(obs_header_.GROUP0, 16);
219 0 : readHeader(obs_header_.PROJ0, 16);
220 0 : readHeader(obs_header_.SCHED0, 24);
221 0 : readHeader(obs_header_.OBSVR0, 40);
222 0 : readHeader(obs_header_.LOSTM0, 16);
223 0 : readHeader(obs_header_.LOETM0, 16);
224 0 : readHeader(obs_header_.ARYNM0);
225 0 : readHeader(obs_header_.NSCAN0);
226 0 : readHeader(obs_header_.TITLE0, 120);
227 0 : readHeader(obs_header_.OBJ0, 16);
228 0 : readHeader(obs_header_.EPOCH0, 8);
229 0 : readHeader(obs_header_.RA00);
230 0 : readHeader(obs_header_.DEC00);
231 0 : readHeader(obs_header_.GL00);
232 0 : readHeader(obs_header_.GB00);
233 0 : readHeader(obs_header_.NCALB0);
234 0 : readHeader(obs_header_.SCNCD0);
235 0 : readHeader(obs_header_.SCMOD0, 120);
236 0 : readHeader<double>(obs_header_.VEL0);
237 0 : readHeader(obs_header_.VREF0, 4);
238 0 : readHeader(obs_header_.VDEF0, 4);
239 0 : readHeader(obs_header_.SWMOD0, 8);
240 0 : readHeader<double>(obs_header_.FRQSW0);
241 0 : readHeader<double>(obs_header_.DBEAM0);
242 0 : readHeader<double>(obs_header_.MLTOF0);
243 0 : readHeader<double>(obs_header_.CMTQ0);
244 0 : readHeader<double>(obs_header_.CMTE0);
245 0 : readHeader<double>(obs_header_.CMTSOM0);
246 0 : readHeader<double>(obs_header_.CMTNODE0);
247 0 : readHeader<double>(obs_header_.CMTI0);
248 0 : readHeader(obs_header_.CMTTMO0, 24);
249 0 : readHeader<double>(obs_header_.SBDX0);
250 0 : readHeader<double>(obs_header_.SBDY0);
251 0 : readHeader<double>(obs_header_.SBDZ10);
252 0 : readHeader<double>(obs_header_.SBDZ20);
253 0 : readHeader<double>(obs_header_.DAZP0);
254 0 : readHeader<double>(obs_header_.DELP0);
255 0 : readHeader<int>(obs_header_.CBIND0);
256 0 : readHeader<int>(obs_header_.NCH0);
257 0 : readHeader<int>(obs_header_.CHRANGE0, 2);
258 0 : readHeader<double>(obs_header_.ALCTM0);
259 0 : readHeader<double>(obs_header_.IPTIM0);
260 0 : readHeader<double>(obs_header_.PA0);
261 0 : readHeader(obs_header_.RX0, 16, NRO_ARYMAX);
262 0 : readHeader<double>(obs_header_.HPBW0, NRO_ARYMAX);
263 0 : readHeader<double>(obs_header_.EFFA0, NRO_ARYMAX);
264 0 : readHeader<double>(obs_header_.EFFB0, NRO_ARYMAX);
265 0 : readHeader<double>(obs_header_.EFFL0, NRO_ARYMAX);
266 0 : readHeader<double>(obs_header_.EFSS0, NRO_ARYMAX);
267 0 : readHeader<double>(obs_header_.GAIN0, NRO_ARYMAX);
268 0 : readHeader(obs_header_.HORN0, 4, NRO_ARYMAX);
269 0 : readHeader(obs_header_.POLTP0, 4, NRO_ARYMAX);
270 0 : readHeader<double>(obs_header_.POLDR0, NRO_ARYMAX);
271 0 : readHeader<double>(obs_header_.POLAN0, NRO_ARYMAX);
272 0 : readHeader<double>(obs_header_.DFRQ0, NRO_ARYMAX);
273 0 : readHeader(obs_header_.SIDBD0, 4, NRO_ARYMAX);
274 0 : readHeader<int>(obs_header_.REFN0, NRO_ARYMAX);
275 0 : readHeader<int>(obs_header_.IPINT0, NRO_ARYMAX);
276 0 : readHeader<int>(obs_header_.MULTN0, NRO_ARYMAX);
277 0 : readHeader<double>(obs_header_.MLTSCF0, NRO_ARYMAX);
278 0 : readHeader(obs_header_.LAGWIN0, 8, NRO_ARYMAX);
279 0 : readHeader<double>(obs_header_.BEBW0, NRO_ARYMAX);
280 0 : readHeader<double>(obs_header_.BERES0, NRO_ARYMAX);
281 0 : readHeader<double>(obs_header_.CHWID0, NRO_ARYMAX);
282 0 : readHeader<int>(obs_header_.ARRY0, NRO_ARYMAX);
283 0 : readHeader<int>(obs_header_.NFCAL0, NRO_ARYMAX);
284 0 : readHeader<double>(obs_header_.F0CAL0, NRO_ARYMAX);
285 0 : for (size_t i = 0; i < NRO_ARYMAX; ++i) {
286 0 : readHeader<double>(obs_header_.FQCAL0[i], 10);
287 : }
288 0 : for (size_t i = 0; i < NRO_ARYMAX; ++i) {
289 0 : readHeader<double>(obs_header_.CHCAL0[i], 10);
290 : }
291 0 : for (size_t i = 0; i < NRO_ARYMAX; ++i) {
292 0 : readHeader<double>(obs_header_.CWCAL0[i], 10);
293 : }
294 0 : readHeader<int>(obs_header_.SCNLEN0);
295 0 : readHeader<int>(obs_header_.SBIND0);
296 0 : readHeader<int>(obs_header_.IBIT0);
297 0 : readHeader(obs_header_.SITE0, 8);
298 0 : readHeader(obs_header_.TRK_TYPE, 8);
299 0 : readHeader(obs_header_.SCAN_COORD, 8);
300 0 : readHeader<int>(obs_header_.NBEAM);
301 0 : readHeader<int>(obs_header_.NPOL);
302 0 : readHeader<int>(obs_header_.NSPWIN);
303 0 : readHeader<int>(obs_header_.CHMAX);
304 0 : readHeader(obs_header_.VERSION, 40);
305 0 : readHeader<int16_t>(obs_header_.ARRYTB, 36);
306 0 : readHeader(obs_header_.POLNAME, 3, 12);
307 : // readHeader(obs_header_.CDMY1, 108);
308 : { // Debug output
309 0 : LogIO os(LogOrigin("NRODataset", "readObsHeader", WHERE) );
310 0 : os << LogIO::DEBUGGING << "NRO ARRAY TABLE from NOSTAR header:" << LogIO::POST;
311 0 : for (int i = 0; i < NRO_ARYMAX; ++i) {
312 0 : os << LogIO::DEBUGGING << "- Array" << i << " : " << obs_header_.ARRYTB[i] << LogIO::POST;
313 : }
314 0 : os << LogIO::DEBUGGING << "POL NAMES:" << LogIO::POST;
315 0 : for (size_t i = 0; i < 12; ++i) {
316 0 : if (obs_header_.POLNAME[i][0] == '\0') continue;
317 0 : os << LogIO::DEBUGGING << i << " : " << obs_header_.POLNAME[i] << LogIO::POST;
318 : }
319 : }
320 0 : }
321 :
322 0 : void NRO2MSReader::readScanData(int const irow, sdfiller::NRODataScanData &data) {
323 0 : if (irow < 0) throw AipsError("Negative row number");
324 0 : size_t offset = len_obs_header_ + static_cast<size_t>(irow) * obs_header_.SCNLEN0;
325 0 : fseek(fp_, offset, SEEK_SET);
326 :
327 0 : readHeader(data.LSFIL0, 4);
328 0 : readHeader(data.ISCN0);
329 0 : readHeader(data.LAVST0, 24);
330 0 : readHeader(data.SCNTP0, 8);
331 0 : readHeader(data.DSCX0);
332 0 : readHeader(data.DSCY0);
333 0 : readHeader(data.SCX0);
334 0 : readHeader(data.SCY0);
335 0 : readHeader(data.PAZ0);
336 0 : readHeader(data.PEL0);
337 0 : readHeader(data.RAZ0);
338 0 : readHeader(data.REL0);
339 0 : readHeader(data.XX0);
340 0 : readHeader(data.YY0);
341 0 : readHeader(data.ARRYT0, 4);
342 0 : readHeader(data.TEMP0);
343 0 : readHeader(data.PATM0);
344 0 : readHeader(data.PH200);
345 0 : readHeader(data.VWIND0);
346 0 : readHeader(data.DWIND0);
347 0 : readHeader(data.TAU0);
348 0 : readHeader(data.TSYS0);
349 0 : readHeader(data.BATM0);
350 0 : readHeader(data.LINE0);
351 0 : readHeader(data.IDMY1, 4);
352 0 : readHeader(data.VRAD0);
353 0 : readHeader(data.FRQ00);
354 0 : readHeader(data.FQTRK0);
355 0 : readHeader(data.FQIF10);
356 0 : readHeader(data.ALCV0);
357 0 : for (size_t i = 0; i < 2; ++i) {
358 0 : readHeader(data.OFFCD0[i], 2);
359 : }
360 0 : readHeader(data.IDMY0);
361 0 : readHeader(data.IDMY2);
362 0 : readHeader(data.DPFRQ0);
363 0 : readHeader(data.ARRYSCN, 8);
364 0 : readHeader(data.CDMY1, 136);
365 0 : readHeader(data.SFCTR0);
366 0 : readHeader(data.ADOFF0);
367 : // KS: use of NCH0 instead of CHMAX (2017/01/25)
368 0 : readHeader(data.LDATA, obs_header_.NCH0 * obs_header_.IBIT0 / 8);
369 0 : }
370 :
371 0 : bool NRO2MSReader::checkScanArray(string const scan_array,
372 : NROArrayData const *header_array) {
373 : // ARRYSCN is a string in the form, ##SS## (beam id in 2 digits, pol name, and spw id in 2 digits)
374 0 : if (scan_array.size() < 6) {
375 : throw AipsError(
376 0 : "Internal Data ERROR: insufficient length of ARRAY information from scan header");
377 : }
378 : // indices in NOSTAR data are 1-base
379 0 : int const sbeam = atoi(scan_array.substr(0, 2).c_str()) -1;
380 0 : string const spol = scan_array.substr(2, 2);
381 0 : int const sspw = atoi(scan_array.substr(4, 2).c_str()) -1;
382 0 : return (sbeam == header_array->getBeamId()
383 0 : && spol == header_array->getPolName()
384 0 : && sspw == header_array->getSpwId());
385 : }
386 :
387 :
388 0 : double NRO2MSReader::getMJD(string const &strtime) {
389 : // TODO: should be checked which time zone the time depends on
390 : // 2008/11/14 Takeshi Nakazato
391 0 : string strYear = strtime.substr(0, 4);
392 0 : string strMonth = strtime.substr(4, 2);
393 0 : string strDay = strtime.substr(6, 2);
394 0 : string strHour = strtime.substr(8, 2);
395 0 : string strMinute = strtime.substr(10, 2);
396 0 : string strSecond = strtime.substr(12, strtime.size() - 12);
397 0 : unsigned int year = atoi(strYear.c_str());
398 0 : unsigned int month = atoi(strMonth.c_str());
399 0 : unsigned int day = atoi(strDay.c_str());
400 0 : unsigned int hour = atoi(strHour.c_str());
401 0 : unsigned int minute = atoi(strMinute.c_str());
402 0 : double second = atof(strSecond.c_str());
403 0 : Time t(year, month, day, hour, minute, second);
404 :
405 0 : return t.modifiedJulianDay() ;
406 : }
407 :
408 0 : double NRO2MSReader::getIntMiddleTimeSec(sdfiller::NRODataScanData const &data) {
409 0 : return getMJD(data.LAVST0) * kDay2Sec + 0.5 * obs_header_.IPTIM0;
410 : }
411 :
412 0 : double NRO2MSReader::getIntStartTimeSec(int const scanno) {
413 0 : if (scanno < 0) throw AipsError("Negative scan number");
414 0 : size_t offset = len_obs_header_ + static_cast<size_t>(scanno) * obs_header_.ARYNM0 * obs_header_.SCNLEN0 + 8;
415 0 : fseek(fp_, offset, SEEK_SET);
416 0 : string time_header;
417 0 : readHeader(time_header, 24);
418 0 : return getMJD(time_header) * kDay2Sec;
419 : }
420 :
421 0 : double NRO2MSReader::getIntEndTimeSec(int const scanno) {
422 0 : double interval = obs_header_.IPTIM0;
423 0 : return getIntStartTimeSec(scanno) + interval;
424 : }
425 :
426 0 : void NRO2MSReader::getFullTimeRange() {
427 0 : if (time_range_sec_.size() != 2) {
428 0 : time_range_sec_.resize(2);
429 : }
430 0 : time_range_sec_[0] = getIntStartTimeSec(0);
431 0 : time_range_sec_[1] = getIntEndTimeSec(obs_header_.NSCAN0 - 1);
432 0 : }
433 :
434 0 : double NRO2MSReader::getMiddleOfTimeRangeSec() {
435 0 : return 0.5 * (time_range_sec_[0] + time_range_sec_[1]);
436 : }
437 :
438 0 : double NRO2MSReader::getRestFrequency(int const spwno) {
439 0 : size_t offset = len_obs_header_ + static_cast<size_t>(getFirstArrayIdWithSpwID(spwno)) * obs_header_.SCNLEN0 + 184;
440 0 : fseek(fp_, offset, SEEK_SET);
441 : double restfreq_header;
442 0 : readHeader(restfreq_header);
443 0 : return restfreq_header;
444 : }
445 :
446 0 : void NRO2MSReader::constructArrayTable() {
447 0 : size_t array_max = NRO_ARYMAX;
448 0 : array_mapper_.resize(array_max);
449 0 : LogIO os(LogOrigin("NRODataset", "constructArrayTable", WHERE) );
450 0 : os << LogIO::DEBUG1 << "NRO ARRAY TABLE:" << LogIO::POST;
451 :
452 0 : for (size_t i = 0; i < array_mapper_.size(); ++i) {
453 0 : array_mapper_[i].set(obs_header_.ARRYTB[i], obs_header_.POLNAME);
454 0 : int beam_id = array_mapper_[i].beam_id;
455 0 : int spw_id = array_mapper_[i].spw_id;
456 0 : int stokes = array_mapper_[i].stokes_type;
457 0 : if (array_mapper_[i].isUsed()) {
458 0 : if (beam_id < 0 || beam_id >= obs_header_.NBEAM
459 0 : || stokes == Stokes::Undefined || spw_id < 0
460 0 : || spw_id >= obs_header_.NSPWIN) {
461 0 : throw AipsError("Internal Data ERROR: inconsistent ARRAY table");
462 : }
463 : os << LogIO::DEBUG1 << "- Array " << i << " : (beam, pol, spw) = ("
464 0 : << beam_id << ", " << array_mapper_[i].pol_name
465 0 : << ", " << spw_id << ")" << LogIO::POST;
466 : }
467 : else {
468 0 : os << LogIO::DEBUG1 << "- Array " << i << " : Unused" << LogIO::POST;
469 : }
470 : }
471 0 : }
472 :
473 0 : string NRO2MSReader::convertVRefName(string const &vref0) {
474 0 : string res;
475 0 : if (vref0 == "LSR") {
476 0 : res = "LSRK";
477 0 : } else if (vref0 == "HEL") {
478 0 : res = "BARY";
479 0 : } else if (vref0 == "GAL") {
480 0 : res = "GALACTO";
481 : } else {
482 0 : res = "Undefined";
483 : }
484 0 : return res;
485 : }
486 :
487 0 : void NRO2MSReader::shiftFrequency(string const &vdef,
488 : double const v,
489 : std::vector<double> &freqs) {
490 0 : double factor = v/2.99792458e8;
491 0 : if (vdef == "RAD") {
492 : //factor = 1.0 / (1.0 + factor);
493 0 : factor = 1.0 - factor;
494 0 : } else if (vdef == "OPT") {
495 : //factor += 1.0;
496 0 : factor = 1.0 / (1.0 + factor);
497 : } else {
498 0 : cout << "vdef=" << vdef << " is not supported." << endl;
499 0 : factor = 1.0;
500 : }
501 0 : for (size_t i = 0; i < 2; ++i) {
502 0 : freqs[i] *= factor;
503 : }
504 0 : }
505 :
506 0 : std::vector<double> NRO2MSReader::getSpectrum(int const irow, sdfiller::NRODataScanData const &data) {
507 : // size of spectrum is not (SCNLEN-HEADER_SIZE)*8/IBIT0
508 : // but obs_header_.NCH0 after binding
509 0 : int const nchan = obs_header_.NCH0;
510 0 : int const chmax_ = (obs_header_.SCNLEN0 - SCNLEN_HEADER_SIZE) * 8 / obs_header_.IBIT0;
511 0 : vector<double> spec( chmax_ ) ; // spectrum "after" binding
512 : // DEBUG
513 : //cout << "NRODataset::getSpectrum() nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
514 :
515 0 : const int bit = obs_header_.IBIT0; // fixed to 12 bit
516 0 : double scale = data.SFCTR0;
517 : // DEBUG
518 : //cout << "NRODataset::getSpectrum() scale = " << scale << endl ;
519 0 : double offset = data.ADOFF0;
520 : // DEBUG
521 : //cout << "NRODataset::getSpectrum() offset = " << offset << endl ;
522 :
523 0 : if ((scale == 0.0) && (offset == 0.0)) {
524 0 : LogIO os( LogOrigin("NRODataset","getSpectrum",WHERE) );
525 0 : os << LogIO::WARN << "zero spectrum for row " << irow << LogIO::POST;
526 0 : if (spec.size() != (unsigned int)nchan) {
527 0 : spec.resize(nchan);
528 : }
529 0 : for (vector<double>::iterator i = spec.begin(); i != spec.end(); ++i) {
530 0 : *i = 0.0;
531 : }
532 0 : return spec;
533 : }
534 :
535 0 : unsigned char *cdata = reinterpret_cast<unsigned char *>(const_cast<char *>(data.LDATA.c_str()));
536 0 : vector<double> mscale;
537 0 : mscale.resize(NRO_ARYMAX);
538 0 : for (size_t i = 0; i < NRO_ARYMAX; ++i) {
539 0 : mscale[i] = obs_header_.MLTSCF0[i];
540 : }
541 0 : string sbeamno = data.ARRYT0.substr(1, data.ARRYT0.size() - 1);
542 0 : int index = atoi(sbeamno.c_str()) - 1;
543 0 : double dscale = mscale[index];
544 :
545 : // char -> int -> double
546 0 : vector<double>::iterator iter = spec.begin();
547 :
548 : static const int shift_right[] = { 4, 0 };
549 : static const int start_pos[] = { 0, 1 };
550 : static const int incr[] = { 0, 3 };
551 0 : int j = 0;
552 0 : for (int i = 0; i < chmax_; ++i) {
553 : // char -> int
554 0 : int ivalue = 0;
555 0 : if (bit == 12) { // 12 bit qunatization
556 0 : const int ialt = i & 1;
557 0 : const int idx = j + start_pos[ialt];
558 0 : const unsigned tmp = unsigned(cdata[idx]) << 8 | cdata[idx + 1];
559 0 : ivalue = int((tmp >> shift_right[ialt]) & 0xFFF);
560 0 : j += incr[ialt];
561 : }
562 :
563 0 : if ((ivalue < 0) || (ivalue > 4096)) {
564 : //cerr << "NRODataset::getSpectrum() ispec[" << i << "] is out of range" << endl ;
565 0 : LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE));
566 0 : os << LogIO::SEVERE << "ivalue for row " << i << " is out of range" << LogIO::EXCEPTION;
567 0 : if (spec.size() != (unsigned int)nchan) {
568 0 : spec.resize(nchan);
569 : }
570 0 : for (vector<double>::iterator i = spec.begin(); i != spec.end(); ++i) {
571 0 : *i = 0.0;
572 : }
573 0 : return spec;
574 : }
575 :
576 : // int -> double
577 0 : *iter = (double)(ivalue * scale + offset) * dscale;
578 : // DEBUG
579 : //cout << "NRODataset::getSpectrum() spec[" << i << "] = " << *iter << endl ;
580 0 : iter++;
581 : }
582 :
583 : // DEBUG
584 : //cout << "NRODataset::getSpectrum() end process" << endl ;
585 0 : return spec;
586 : }
587 :
588 : //Int NRO2MSReader::getPolNo(string const &rx) {
589 : // Int polno = 0;
590 : // // 2013/01/23 TN
591 : // // In NRO 45m telescope, naming convension for dual-polarization
592 : // // receiver is as follows:
593 : // //
594 : // // xxxH for horizontal component,
595 : // // xxxV for vertical component.
596 : // //
597 : // // Exception is H20ch1/ch2.
598 : // // Here, POLNO is assigned as follows:
599 : // //
600 : // // POLNO=0: xxxH or H20ch1
601 : // // 1: xxxV or H20ch2
602 : // //
603 : // // For others, POLNO is always 0.
604 : // string last_letter = rx.substr(rx.size()-1, 1);
605 : // if ((last_letter == "V") || (rx == "H20ch2")) {
606 : // polno = 1;
607 : // }
608 : //
609 : // return polno;
610 : //}
611 :
612 0 : void NRO2MSReader::initializeSpecific() {
613 : // std::cout << "NRO2MSReader::initialize" << std::endl;
614 0 : fp_ = fopen(STRING2CHAR(name_), "rb");
615 0 : if (fp_ == NULL) {
616 0 : throw AipsError("Input data doesn't exist or is invalid");
617 : }
618 0 : checkEndian();
619 :
620 0 : readObsHeader();
621 0 : getFullTimeRange();
622 0 : constructArrayTable();
623 0 : }
624 :
625 0 : void NRO2MSReader::finalizeSpecific() {
626 : // std::cout << "NRO2MSReader::finalize" << std::endl;
627 0 : if (fp_ != NULL) {
628 0 : fclose(fp_);
629 : }
630 0 : fp_ = NULL;
631 0 : }
632 :
633 0 : size_t NRO2MSReader::getNumberOfRows() const {
634 0 : int nrows = obs_header_.ARYNM0 * obs_header_.NSCAN0;
635 0 : return (nrows >= 0) ? nrows : 0;
636 : }
637 :
638 0 : MDirection::Types NRO2MSReader::getDirectionFrame() const {
639 : MDirection::Types res;
640 0 : int scan_coord = obs_header_.SCNCD0;
641 : //std::cout << "********** SCNCD0 = " << scan_coord << std::endl;
642 0 : if (scan_coord == 0) { // RADEC
643 0 : string epoch = obs_header_.EPOCH0;
644 0 : if (epoch == "J2000") {
645 0 : res = MDirection::J2000;
646 0 : } else if (epoch == "B1950") {
647 0 : res = MDirection::B1950;
648 : } else {
649 0 : throw AipsError("Unsupported epoch.");
650 : }
651 0 : } else if (scan_coord == 1) { // LB
652 0 : res = MDirection::GALACTIC;
653 0 : } else if (scan_coord == 2) { // AZEL
654 0 : res = MDirection::AZEL;
655 : } else {
656 0 : throw AipsError("Unsupported epoch.");
657 : }
658 0 : return res;
659 : }
660 :
661 0 : Bool NRO2MSReader::getAntennaRowImpl(AntennaRecord &record) {
662 : // std::cout << "NRO2MSReader::getAntennaRowImpl" << std::endl;
663 :
664 0 : record.station = "";
665 0 : String header_antenna_name = obs_header_.SITE0;
666 0 : ostringstream oss;
667 0 : oss << header_antenna_name << "-BEAM" << beam_id_counter_;
668 0 : record.name = String(oss);
669 0 : record.mount = "ALT-AZ";
670 0 : record.dish_diameter = queryAntennaDiameter(header_antenna_name);
671 0 : record.type = "GROUND-BASED";
672 0 : record.position = MPosition(MVPosition(posx_, posy_, posz_), MPosition::ITRF);
673 :
674 0 : beam_id_counter_++;
675 0 : if (obs_header_.NBEAM <= beam_id_counter_) {
676 0 : get_antenna_row_ = [&](sdfiller::AntennaRecord &r) {return NRO2MSReader::noMoreRowImpl<AntennaRecord>(r);};
677 : }
678 :
679 0 : return true;
680 : }
681 :
682 0 : Bool NRO2MSReader::getObservationRowImpl(ObservationRecord &record) {
683 : // std::cout << "NRO2MSReader::getObservationRowImpl" << std::endl;
684 :
685 0 : if (record.time_range.size() != 2) {
686 0 : record.time_range.resize(2);
687 : }
688 0 : for (int i = 0; i < 2; ++i) {
689 : // 2018/05/30 TN
690 : // time_range should be in sec
691 : // CAS-11223 NRO timestamp is in JST so it should be converted to UTC
692 0 : record.time_range[i] = jst2utcSec(time_range_sec_[i]);
693 : }
694 :
695 0 : record.observer = obs_header_.OBSVR0;
696 0 : record.project = obs_header_.PROJ0;
697 0 : record.telescope_name = obs_header_.SITE0;
698 :
699 : // only one entry so redirect function pointer to noMoreRowImpl
700 0 : get_observation_row_ = [&](sdfiller::ObservationRecord &r) {return NRO2MSReader::noMoreRowImpl<ObservationRecord>(r);};
701 :
702 0 : return true;
703 : }
704 :
705 :
706 0 : Bool NRO2MSReader::getProcessorRowImpl(ProcessorRecord &/*record*/) {
707 : // std::cout << "NRO2MSReader::getProcessorRowImpl" << std::endl;
708 :
709 : // just add empty row once
710 :
711 : // only one entry so redirect function pointer to noMoreRowImpl
712 0 : get_processor_row_ = [&](sdfiller::ProcessorRecord &r) {return NRO2MSReader::noMoreRowImpl<ProcessorRecord>(r);};
713 :
714 0 : return true;
715 : }
716 :
717 0 : Bool NRO2MSReader::getSourceRowImpl(SourceRecord &record) {
718 0 : record.name = obs_header_.OBJ0;
719 0 : record.source_id = 0;
720 0 : record.spw_id = source_spw_id_counter_;
721 0 : record.direction = MDirection(
722 0 : Quantity(obs_header_.RA00, "rad"),
723 0 : Quantity(obs_header_.DEC00, "rad"),
724 0 : getDirectionFrame());
725 0 : record.proper_motion = Vector<Double>(0, 0.0);
726 0 : record.rest_frequency.resize(1);
727 0 : record.rest_frequency[0] = getRestFrequency(source_spw_id_counter_);
728 : //record.transition
729 0 : record.num_lines = 1;
730 0 : record.sysvel.resize(1);
731 0 : record.sysvel[0] = obs_header_.VEL0;
732 0 : double time_sec = getMiddleOfTimeRangeSec();
733 : // 2018/05/30 TN
734 : // CAS-11223
735 0 : record.time = jst2utcSec(time_sec);
736 : // 2018/05/30 TN
737 : // CAS-11442
738 0 : record.interval = time_range_sec_[1] - time_range_sec_[0];
739 :
740 0 : source_spw_id_counter_++;
741 0 : if (obs_header_.NSPWIN <= source_spw_id_counter_) {
742 0 : get_source_row_ = [&](sdfiller::SourceRecord &r) {return NRO2MSReader::noMoreRowImpl<SourceRecord>(r);};
743 : }
744 :
745 0 : return true;
746 : }
747 :
748 0 : Bool NRO2MSReader::getFieldRowImpl(FieldRecord &record) {
749 0 : record.name = obs_header_.OBJ0;
750 0 : record.field_id = 0;
751 : // 2018/05/30 TN
752 : // CAS-11223
753 0 : record.time = jst2utcSec(getMiddleOfTimeRangeSec());
754 0 : record.source_name = obs_header_.OBJ0;
755 0 : record.frame = getDirectionFrame();
756 0 : Matrix<Double> direction0(2, 2, 0.0);
757 0 : Matrix<Double> direction(direction0(IPosition(2, 0, 0), IPosition(2, 1, 0)));
758 0 : direction(0, 0) = obs_header_.RA00;
759 0 : direction(1, 0) = obs_header_.DEC00;
760 0 : record.direction = direction;
761 :
762 : // only one entry so redirect function pointer to noMoreRowImpl
763 0 : get_field_row_ = [&](sdfiller::FieldRecord &r) {return NRO2MSReader::noMoreRowImpl<FieldRecord>(r);};
764 :
765 0 : return true;
766 : }
767 :
768 0 : Bool NRO2MSReader::getSpectralWindowRowImpl(
769 : SpectralWindowRecord &record) {
770 0 : record.spw_id = spw_id_counter_;
771 0 : record.num_chan = obs_header_.NCH0;
772 : MFrequency::Types frame_type;
773 0 : Bool status = MFrequency::getType(frame_type, convertVRefName(obs_header_.VREF0));
774 0 : if (!status) {
775 0 : frame_type = MFrequency::N_Types;
776 : }
777 0 : record.meas_freq_ref = frame_type;
778 :
779 0 : NRODataScanData scan_data;
780 0 : int spw_id_array = getFirstArrayIdWithSpwID(spw_id_counter_);
781 0 : readScanData(spw_id_array, scan_data);
782 0 : double freq_offset = scan_data.FRQ00 - obs_header_.F0CAL0[spw_id_array];
783 0 : std::vector<double> freqs(2, freq_offset);
784 0 : std::vector<double> chcal(2);
785 0 : for (size_t i = 0; i < 2; ++i) {
786 0 : freqs[i] += obs_header_.FQCAL0[spw_id_array][i];
787 0 : chcal[i] = obs_header_.CHCAL0[spw_id_array][i];
788 : }
789 : //-------------(change 2016/9/23)---------
790 0 : shiftFrequency(obs_header_.VDEF0, obs_header_.VEL0, freqs);
791 : //-------------(end change 2016/9/23)-----
792 0 : double band_width = freqs[1] - freqs[0];
793 0 : double chan_width = band_width / (chcal[1] - chcal[0]);
794 0 : if (scan_data.FQIF10 > 0.0) { // USB
795 0 : double tmp = freqs[1];
796 0 : freqs[1] = freqs[0];
797 0 : freqs[0] = tmp;
798 0 : chan_width *= -1.0;
799 : }
800 0 : record.refpix = chcal[0] - 1; // 0-based
801 0 : record.refval = freqs[0];
802 0 : record.increment = chan_width;
803 :
804 0 : spw_id_counter_++;
805 0 : if (obs_header_.NSPWIN <= spw_id_counter_) {
806 0 : get_spw_row_ = [&](sdfiller::SpectralWindowRecord &r) {return NRO2MSReader::noMoreRowImpl<SpectralWindowRecord>(r);};
807 : }
808 :
809 0 : return true;
810 : }
811 :
812 0 : Bool NRO2MSReader::getData(size_t irow, DataRecord &record) {
813 : // std::cout << "NRO2MSReader::getData(irow=" << irow << ")" << std::endl;
814 :
815 0 : if (irow >= getNumberOfRows()) {
816 0 : return false;
817 : }
818 :
819 : // std::cout << "Accessing row " << irow << std::endl;
820 0 : NRODataScanData scan_data;
821 0 : readScanData(irow, scan_data);
822 :
823 : // Verify Array INFO in scan header
824 0 : string array_number = scan_data.ARRYT0.substr(1, scan_data.ARRYT0.size() - 1);
825 0 : size_t const array_id = atoi(array_number.c_str()) - 1; // 1-base -> 0-base
826 0 : if (!checkScanArray(scan_data.ARRYSCN, &array_mapper_[array_id])) {
827 0 : throw AipsError("Internal Data ERROR: inconsistent ARRAY information in scan header");
828 : }
829 :
830 : // 2018/05/30 TN
831 : // CAS-11223
832 0 : record.time = jst2utcSec(getIntMiddleTimeSec(scan_data));
833 0 : record.interval = obs_header_.IPTIM0;
834 : // std::cout << "TIME=" << record.time << " INTERVAL=" << record.interval
835 : // << std::endl;
836 :
837 0 : Int srctype = (scan_data.SCNTP0 == "ON") ? 0 : 1;
838 0 : record.intent = getIntent(srctype);
839 0 : record.scan = (Int)scan_data.ISCN0;
840 0 : record.subscan = getSubscan(srctype);
841 0 : record.field_id = 0;
842 : // Int ndata_per_ant = obs_header_.NPOL * obs_header_.NSPWIN;
843 0 : record.antenna_id = (Int) getNROArrayBeamId(array_id); //(irow / ndata_per_ant % obs_header_.NBEAM);
844 0 : record.direction_vector(0) = scan_data.SCX0;
845 0 : record.direction_vector(1) = scan_data.SCY0;
846 0 : record.scan_rate = 0.0;
847 0 : record.feed_id = (Int)0;
848 0 : record.spw_id = (Int) getNROArraySpwId(array_id);//(irow % obs_header_.NSPWIN);
849 : // Int ndata_per_scan = obs_header_.NBEAM * obs_header_.NPOL * obs_header_.NSPWIN;
850 0 : record.pol = getNROArrayPol(array_id); //getPolNo(obs_header_.RX0[irow % ndata_per_scan]);
851 0 : record.pol_type = "linear";
852 :
853 : // std::cout << "set data size to " << num_chan_map_[record.spw_id] << " shape "
854 : // << record.data.shape() << std::endl;
855 0 : record.setDataSize(obs_header_.NCH0);
856 0 : auto const spectrum = getSpectrum(irow, scan_data);
857 0 : record.data.resize(spectrum.size());
858 0 : std::copy(spectrum.cbegin(), spectrum.cend(), record.data.begin());
859 0 : size_t flag_len = obs_header_.NCH0;
860 0 : for (size_t i = 0; i < flag_len; ++i) {
861 0 : record.flag(i) = false;
862 : }
863 0 : record.flag_row = false;
864 :
865 : // std::cout << "set tsys size to " << tsys_column_.shape(index)[0]
866 : // << " shape " << record.tsys.shape() << std::endl;
867 0 : record.setTsysSize(1);
868 0 : record.tsys(0) = scan_data.TSYS0;
869 :
870 0 : record.temperature = scan_data.TEMP0 + 273.15; // Celsius to Kelvin
871 0 : record.pressure = scan_data.PATM0;
872 0 : record.rel_humidity = scan_data.PH200;
873 0 : record.wind_speed = scan_data.VWIND0;
874 0 : record.wind_direction = scan_data.DWIND0;
875 :
876 0 : return true;
877 : }
878 :
879 : } //# NAMESPACE CASA - END
|
