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