Line data Source code
1 : //# Importmiriad: miriad dataset to MeasurementSet conversion
2 : //# Copyright (C) 1997,2000,2001,2002,2013,2015
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify
6 : //# it under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or
8 : //# (at your option) any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful,
11 : //# but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : //# GNU General Public License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this program; if not, write to the Free Software
17 : //# Foundation, Inc., 675 Masve, Cambridge, MA 02139, USA.
18 : //#
19 :
20 : //#Includes
21 : #include <miriad/Filling/Importmiriad.h>
22 :
23 : #include <casacore/casa/Inputs/Input.h>
24 : #include <casacore/casa/OS/File.h>
25 : #include <casacore/casa/Utilities/GenSort.h>
26 :
27 : #include <casacore/casa/Arrays/Cube.h>
28 : #include <casacore/casa/Arrays/Matrix.h>
29 : #include <casacore/casa/Arrays/Vector.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/ArrayUtil.h>
32 : #include <casacore/casa/Arrays/ArrayLogical.h>
33 : #include <casacore/casa/Arrays/MatrixMath.h>
34 :
35 :
36 : #include <casacore/measures/Measures.h>
37 : #include <casacore/measures/Measures/MPosition.h>
38 : #include <casacore/measures/Measures/MeasData.h>
39 : #include <casacore/measures/Measures/Stokes.h>
40 :
41 : #include <casacore/tables/Tables.h>
42 : #include <casacore/tables/Tables/TableInfo.h>
43 :
44 : #include <casacore/ms/MeasurementSets.h>
45 :
46 : #include <casacore/mirlib/maxdimc.h>
47 : #include <casacore/mirlib/miriad.h>
48 : #include <sstream>
49 :
50 : using namespace casa;
51 :
52 : // helper functions
53 :
54 0 : String show_version_info()
55 : {
56 : return "============================================================\n"
57 : "Importmiriad - last few updates:\n"
58 : " Mar 2013 - make it process ATCA/CABB data \n"
59 : " Jul 2015 - deal with multiple zoom setups \n"
60 : " ...\n"
61 0 : "============================================================\n";
62 : }
63 :
64 :
65 : // Convert fits date string of form dd/mm/yy to mjd seconds
66 :
67 0 : Double date2mjd(const String& date)
68 : {
69 : Int day,month,year;
70 :
71 0 : if (date[2] == '/') { // old FITS style (dd/mm/yy)
72 0 : sscanf(date.chars(),"%2d/%2d/%2d",&day,&month,&year);
73 0 : year+=1900;
74 0 : if (year<1950) year+=100; // yuck !
75 : } else { // YEAR-2000 (ISO) convention (ccyy-mm-ddThh:mm:ss.sss)
76 : // cerr<< "Parsing new Year-2000 notation" << endl;
77 0 : sscanf(date.chars(),"%4d-%2d-%2d",&year,&month,&day);
78 : //sscanf(date,"%4d-%2d-%2dT%2d:%2d:%f",&year,&month,&day,&hour,&min,&sec);
79 : }
80 0 : MVTime mjd_date(year,month,(Double)day);
81 0 : return mjd_date.second();
82 0 : }
83 :
84 : // apply CARMA line calibration, the 'linecal' method
85 : // check MIRIADs 'uvcal options=linecal' for the other approach
86 :
87 0 : void linecal(int ndata, float *data, float phi1, float phi2)
88 : {
89 : float x,y,c,s;
90 0 : if (ndata <= 0) return;
91 :
92 0 : c = cos(phi1-phi2);
93 0 : s = sin(phi1-phi2);
94 :
95 0 : for (int i=0; i<ndata*2; i+=2) {
96 0 : x = data[i];
97 0 : y = data[i+1];
98 0 : data[i] = c*x - s*y;
99 0 : data[i+1] = s*x + c*y;
100 : }
101 : }
102 :
103 : // ==============================================================================================
104 1 : Importmiriad::Importmiriad(String& infile, Int debug,
105 1 : Bool Qtsys, Bool Qarrays, Bool Qlinecal)
106 : {
107 1 : infile_p = infile;
108 1 : debug_p = debug;
109 1 : msc_p = 0;
110 1 : nArray_p = 0;
111 1 : nFreqSet_p = 0;
112 1 : nfield = 0; // # mosaiced fields (using offsets?)
113 1 : npoint = 0; // # pointings (using independant RA/DEC?)
114 1 : Qtsys_p = Qtsys;
115 1 : Qarrays_p = Qarrays;
116 1 : Qlinecal_p = Qlinecal;
117 1 : os_p = LogOrigin("Importmiriad");
118 1 : zero_tsys = 0;
119 10001 : for (int i=0; i<MAXFIELD; i++) fcount[i] = 0;
120 65 : for (int i=0; i<MAXANT; i++) phasem1[i] = 0.0;
121 :
122 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::Importmiriad debug_level=" << debug << LogIO::POST;
123 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Opening miriad dataset " << infile_p << LogIO::POST;
124 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << (Qtsys_p ? "tsys weights" : "weights=1") << LogIO::POST;
125 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << (Qarrays_p ? "split arrays" : "single array forced") << LogIO::POST;
126 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << (Qlinecal_p ? "linecal applied" : "linecal not applied") << LogIO::POST;
127 :
128 : if (sizeof(double) != sizeof(Double))
129 : os_p<< LogIO::SEVERE<<"Double != double; importmiriad will probably fail" << LogIO::POST;
130 : if (sizeof(int) != sizeof(Int))
131 : os_p << LogIO::SEVERE<<"int != Int; importmiriad will probably fail" << LogIO::POST;
132 :
133 : // open miriad dataset
134 1 : uvopen_c(&uv_handle_p, infile_p.chars(), "old");
135 :
136 : // preamble data must be UVW (default miriad is UV)
137 1 : uvset_c(uv_handle_p,"preamble","uvw/time/baseline",0,0.0,0.0,0.0);
138 :
139 : // initialize those UV variables that need to be tracked
140 1 : Tracking(-1);
141 1 : }
142 :
143 : // ==============================================================================================
144 1 : Importmiriad::~Importmiriad()
145 : {
146 1 : if (msc_p) { delete msc_p; msc_p=0;}
147 1 : if (Debug(1)) os_p <<LogIO::DEBUG1<< "Importmiriad::~Importmiriad" << LogIO::POST;
148 1 : if (zero_tsys)
149 0 : os_p << "There were " << zero_tsys << " record with no WEIGHT due to zero TSYS" << LogIO::POST;
150 :
151 1 : if (Debug(1))
152 0 : for (int i=0; i<nfield; i++)
153 0 : os_p << LogIO::DEBUG1 << "Field " << i << " = " << fcount[i] << " records" << LogIO::POST;
154 :
155 : // most single MIRIAD files are time ordered, so could check for
156 : // that, and if so, add SORT_ORDER = 'ASCENDING' and COLUMNS = 'TIME'
157 :
158 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "*** Closing " << infile_p << " ***"<<LogIO::POST;
159 : //os_p << "Importmiriad::END" << LogIO::POST;
160 1 : if (Debug(1)) os_p<<LogIO::DEBUG1 << show_version_info()<<LogIO::POST;
161 1 : }
162 :
163 : // ==============================================================================================
164 0 : void Importmiriad::Error(char *msg)
165 : {
166 0 : throw(AipsError(msg));
167 : }
168 :
169 : // ==============================================================================================
170 0 : void Importmiriad::Warning(char *msg)
171 : {
172 0 : os_p << LogIO::WARN<< "### Warning: " << msg << LogIO::POST;
173 0 : }
174 :
175 : // ==============================================================================================
176 913 : Bool Importmiriad::Debug(int level)
177 : {
178 913 : Bool ok=false;
179 913 : if (level <= debug_p) ok=true;
180 913 : return ok;
181 : }
182 :
183 : // ==============================================================================================
184 1 : void Importmiriad::checkInput(Block<Int>& spw, Block<Int>& wide)
185 : {
186 : Int i, nread, nwread, vlen, vupd;
187 : char vtype[10], vdata[256];
188 : Float epoch;
189 :
190 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::checkInput" << LogIO::POST;
191 :
192 : // Let's read one scan and try and derive some basics. If important
193 : // variables not present, bail out (or else scan on)
194 :
195 1 : nvis = 0;
196 : for (;;) { // loop forever until happy or EOF
197 :
198 1 : uvread_c(uv_handle_p, preamble, data, flags, MAXCHAN, &nread);
199 1 : if (nread <= 0) break;
200 1 : nvis++;
201 1 : uvwread_c(uv_handle_p, wdata, wflags, MAXCHAN, &nwread);
202 :
203 1 : if (nvis == 1) {
204 :
205 : // get the initial correllator setup
206 1 : check_window();
207 : // setup the 'keep' array to specify which data we want to keep
208 1 : keep_p.resize(MAXWIN+MAXWIDE); keep_p.set(False);
209 1 : if (spw.nelements()>0) {
210 67 : for (Int i=0; i<MAXWIN+MAXWIDE; i++) keep_p[i]=(spw[0]==-1);
211 2 : for (uInt i=0; i<spw.nelements(); i++) {
212 1 : if (spw[i]>=0 && spw[i]<win[0].nspect) keep_p[spw[i]]=true;
213 : }
214 : }
215 1 : Int n=win[0].nspect;
216 1 : if (wide.nelements()>0) {
217 0 : for (Int i=0; i<win[0].nwide; i++) keep_p[n+i]=(wide[0]==-1);
218 0 : for (uInt i=0; i<wide.nelements(); i++) {
219 0 : if (wide[i]>0 && wide[i]<win[0].nwide) keep_p[n+wide[i]]=true;
220 : }
221 : }
222 : // should store nread + nwread, or handle it as option
223 1 : if (win[freqSet_p].nspect > 0) { // narrow band, with possibly wide band also
224 1 : nchan_p = nread;
225 1 : nwide_p = nwread;
226 : } else { // wide band data only: nread=nwread
227 0 : nchan_p = nread;
228 0 : nwide_p = 0;
229 : }
230 :
231 : // get the initial array configuration
232 :
233 1 : nants_offset_p = 0;
234 1 : uvgetvr_c(uv_handle_p,H_INT, "nants", (char *)&nants_p,1);
235 1 : if (nants_p > MAXANT) {
236 0 : ostringstream msg;
237 0 : msg << "Importmiriad: Too many antennas, "<< nants_p << " > "<< MAXANT<<endl;
238 0 : throw(AipsError(msg.str()));
239 : return;
240 0 : }
241 :
242 1 : uvgetvr_c(uv_handle_p,H_DBLE,"antpos",(char *)antpos,3*nants_p);
243 1 : if (Debug(1)) {
244 0 : os_p << LogIO::DEBUG1 << "Found " << nants_p << " antennas (first scan)" << LogIO::POST;
245 0 : for (int i=0; i<nants_p; i++) {
246 0 : os_p << LogIO::DEBUG1 << antpos[i] << " " <<
247 0 : antpos[i+nants_p] << " " <<
248 0 : antpos[i+nants_p*2] << LogIO::POST;
249 : }
250 : }
251 :
252 : // remember systemp is stored systemp[nants][nwin] in C notation
253 1 : if (win[freqSet_p].nspect > 0) {
254 1 : uvgetvr_c(uv_handle_p,H_REAL,"systemp",(char *)systemp,nants_p*win[freqSet_p].nspect);
255 1 : if (Debug(1)) {
256 0 : os_p << LogIO::DEBUG1 << "Found systemps (first scan)" ;
257 0 : for (Int i=0; i<nants_p; i++) os_p << systemp[i] << " ";
258 0 : os_p << LogIO::POST;
259 : }
260 : } else {
261 0 : uvgetvr_c(uv_handle_p,H_REAL,"wsystemp",(char *)systemp,nants_p);
262 0 : if (Debug(1)) {
263 0 : os_p << LogIO::DEBUG1 << "Found wsystemps (first scan)" ;
264 0 : for (Int i=0; i<nants_p; i++) os_p << systemp[i] << " ";
265 0 : os_p << LogIO::POST;
266 : }
267 : }
268 :
269 1 : if (win[freqSet_p].nspect > 0) {
270 1 : uvgetvr_c(uv_handle_p,H_DBLE,"restfreq",(char *)win[freqSet_p].restfreq,win[freqSet_p].nspect);
271 1 : if (Debug(1)) {
272 0 : os_p << LogIO::DEBUG1 << "Found restfreq (first scan)" ;
273 0 : for (Int i=0; i<win[freqSet_p].nspect; i++) os_p << win[freqSet_p].restfreq[i] << " ";
274 0 : os_p << LogIO::POST;
275 : }
276 : }
277 :
278 : // Note that MIRIAD coordinates are in nanosec, but actual unused
279 : // antennas are filled with -999 values (or sometimes 0!)
280 :
281 1 : uvprobvr_c(uv_handle_p,"project",vtype,&vlen,&vupd);
282 1 : if (vupd) {
283 1 : uvgetvr_c(uv_handle_p,H_BYTE,"project",vdata,32);
284 1 : project_p = vdata;
285 : } else
286 0 : project_p = "unknown";
287 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Project=>" << project_p << "<=" << LogIO::POST;
288 :
289 1 : uvprobvr_c(uv_handle_p,"version",vtype,&vlen,&vupd);
290 1 : if (vupd) {
291 0 : uvgetvr_c(uv_handle_p,H_BYTE,"version",vdata,80);
292 0 : version_p = vdata;
293 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Version=>" << version_p << "<=" << LogIO::POST;
294 : }
295 1 : uvgetvr_c(uv_handle_p,H_BYTE,"source",vdata,16);
296 1 : object_p = vdata;
297 :
298 : // TODO: telescope will now change, so this is not a good idea
299 1 : uvgetvr_c(uv_handle_p,H_BYTE,"telescop",vdata,16);
300 1 : array_p = vdata;
301 :
302 : // array_p = "CARMA";
303 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "First baseline=>" << array_p << "<=" << LogIO::POST;
304 :
305 : // All CARMA (OVRO,BIMA,SZA) & ATCA have this
306 1 : mount_p = 0;
307 : #if 0
308 : if (array_p == "VLA")
309 : mount_p = 1;
310 : uvrdvr_c(uv_handle_p,H_INT,"mount",(char *)&mount_p, (char *)&mount_p, 1);
311 : os_p << "Warning: " << array_p
312 : << " Cannot handle all of this telescope yet" << LogIO::POST;
313 : os_p << "Assumed mount=" << mount_p << LogIO::POST;
314 : #endif
315 :
316 1 : uvprobvr_c(uv_handle_p,"observer",vtype,&vlen,&vupd);
317 1 : if (vupd) {
318 1 : uvgetvr_c(uv_handle_p,H_BYTE,"observer",vdata,16);
319 1 : observer_p = vdata;
320 : } else
321 0 : observer_p = "unknown";
322 :
323 1 : uvgetvr_c(uv_handle_p,H_REAL,"epoch",(char *)&epoch,1);
324 1 : epoch_p = epoch;
325 : // do this globally, we used to do this in the Field table alone
326 1 : epochRef_p=MDirection::J2000;
327 1 : if (nearAbs(epoch_p,1950.0,0.01)) epochRef_p=MDirection::B1950;
328 :
329 1 : uvgetvr_c(uv_handle_p,H_INT,"npol", (char *)&npol_p,1);
330 1 : uvgetvr_c(uv_handle_p,H_INT,"pol",(char *)&pol_p,1);
331 1 : uvgetvr_c(uv_handle_p,H_REAL,"inttime",(char *)&inttime_p,1);
332 1 : uvgetvr_c(uv_handle_p,H_REAL,"jyperk",(char *)&jyperk_p,1);
333 :
334 1 : uvprobvr_c(uv_handle_p,"freq",vtype,&vlen,&vupd);
335 1 : freq_p = 1e9;
336 1 : if (vupd) {
337 0 : uvgetvr_c(uv_handle_p,H_DBLE,"freq",(char *)&freq_p,1);
338 0 : freq_p *= 1e9;
339 : }
340 :
341 1 : uvprobvr_c(uv_handle_p,"ifchain",vtype,&vlen,&vupd);
342 1 : if(!vupd) {
343 0 : for (i=0; i<MAXWIN+MAXWIDE;i++) win[freqSet_p].chain[i]=0;
344 : }
345 : // and initial source position
346 :
347 1 : uvgetvr_c(uv_handle_p,H_DBLE,"ra", (char *)&ra_p, 1);
348 1 : uvgetvr_c(uv_handle_p,H_DBLE,"dec",(char *)&dec_p,1);
349 :
350 : // check if certain calibration tables are present and warn if so,
351 : // since we can't (don't want to) deal with them here; miriad
352 : // programs like uvcat should be used to apply them!
353 :
354 1 : if (hexists_c(uv_handle_p,"gains"))
355 0 : os_p << LogIO::WARN << "gains table present, but cannot apply it" << LogIO::POST;
356 1 : if (hexists_c(uv_handle_p,"bandpass"))
357 0 : os_p << LogIO::WARN << "bandpass table present, but cannot apply it" << LogIO::POST;
358 1 : if (hexists_c(uv_handle_p,"leakage"))
359 0 : os_p << LogIO::WARN << "leakage table present, but cannot apply it" << LogIO::POST;
360 :
361 :
362 1 : if (npol_p > 1) { // read the next npol-1 scans to find the other pols
363 4 : for (i=1; i<npol_p; i++) {
364 3 : uvread_c(uv_handle_p, preamble, data, flags, MAXCHAN, &nread);
365 3 : if (nread <= 0) {
366 0 : break;
367 : }
368 3 : if (Debug(1)) { if (i==1) os_p << LogIO::DEBUG1 << "POL(" << i << ") = " << pol_p[0] << LogIO::POST;}
369 3 : uvgetvr_c(uv_handle_p,H_INT,"pol",(char *)&pol_p[i],1); // FIX
370 3 : if (Debug(1)) os_p << LogIO::DEBUG1 << "POL(" << i+1 << ") = " << pol_p[i] << LogIO::POST;
371 : }
372 : }
373 : // only do one scan
374 1 : break;
375 : }
376 0 : }
377 1 : if (nvis == 0) {
378 0 : throw(AipsError("Importmiriad: Bad first uvread: no narrow or wide band data present"));
379 : return;
380 : }
381 : //os_p << "Importmiriad::checkInput: " << nvis << " records found" << LogIO::POST;
382 : //os_p << "Found " << nvis << " records" << LogIO::POST;
383 1 : uvrewind_c(uv_handle_p);
384 :
385 1 : Int numCorr = npol_p;
386 1 : corrType_p.resize(numCorr);
387 5 : for (i=0; i < numCorr; i++) {
388 : // note: 1-based ref pix
389 4 : corrType_p(i)=pol_p[i];
390 : // convert AIPS-convention Stokes description to CASA enum
391 : // CHECK if these are really the right conversions for CASA
392 4 : if (corrType_p(i)<0) {
393 4 : if (corrType_p(i)==-8) corrType_p(i)=Stokes::YX;
394 4 : if (corrType_p(i)==-7) corrType_p(i)=Stokes::XY;
395 4 : if (corrType_p(i)==-6) corrType_p(i)=Stokes::YY;
396 4 : if (corrType_p(i)==-5) corrType_p(i)=Stokes::XX;
397 4 : if (corrType_p(i)==-4) corrType_p(i)=Stokes::LR;
398 4 : if (corrType_p(i)==-3) corrType_p(i)=Stokes::RL;
399 4 : if (corrType_p(i)==-2) corrType_p(i)=Stokes::LL;
400 4 : if (corrType_p(i)==-1) corrType_p(i)=Stokes::RR;
401 : }
402 : }
403 :
404 1 : Vector<Int> tmp(numCorr); tmp=corrType_p;
405 : // Sort the polarizations to standard order
406 1 : GenSort<Int>::sort(corrType_p);
407 1 : corrIndex_p.resize(numCorr);
408 : // Get the sort indices to rearrange the data to standard order
409 5 : for (i=0;i<numCorr;i++) {
410 20 : for (Int j=0;j<numCorr;j++) {
411 16 : if (corrType_p(j)==tmp(i)) corrIndex_p(i)=j;
412 : }
413 : }
414 :
415 : // Figure out the correlation products from the polarizations
416 1 : corrProduct_p.resize(2,numCorr); corrProduct_p=0;
417 5 : for (i=0; i<numCorr; i++) {
418 4 : Fallible<Int> receptor=Stokes::receptor1(Stokes::type(corrType_p(i)));
419 4 : if (receptor.isValid()) corrProduct_p(0,i)=receptor;
420 4 : receptor=Stokes::receptor2(Stokes::type(corrType_p(i)));
421 4 : if (receptor.isValid()) corrProduct_p(1,i)=receptor;
422 4 : }
423 1 : }
424 :
425 : // ==============================================================================================
426 1 : void Importmiriad::setupMeasurementSet(const String& MSFileName, Bool useTSM)
427 : {
428 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::setupMeasurementSet" << LogIO::POST;
429 :
430 1 : Int nCorr = (array_p=="CARMA" ? 1 : npol_p); // # stokes (1 for CARMA for now)
431 1 : Int nChan = nchan_p; // we are only exporting the narrow channels to the MS
432 :
433 : // Make the MS table
434 1 : TableDesc td = MS::requiredTableDesc();
435 :
436 1 : MS::addColumnToDesc(td, MS::DATA,2);
437 1 : td.removeColumn(MS::columnName(MS::FLAG));
438 1 : MS::addColumnToDesc(td, MS::FLAG,2);
439 :
440 : #if 0
441 : // why does the FITS code do this? We don't need it....
442 : td.removeColumn(MS::columnName(MS::SIGMA));
443 : MS::addColumnToDesc(td, MS::SIGMA, IPosition(1,nCorr),
444 : ColumnDesc::Direct);
445 : #endif
446 :
447 :
448 : // #define OLD_CODE // define this if you want to try the old method again
449 :
450 : #ifdef OLD_CODE
451 : // OLD
452 : if (useTSM) {
453 : td.defineHypercolumn("TiledData",3,
454 : stringToVector(MS::columnName(MS::DATA)+","+
455 : MS::columnName(MS::FLAG)));
456 : }
457 : #else
458 : // NEW
459 1 : if (useTSM) {
460 1 : td.defineHypercolumn("TiledData",3,
461 2 : stringToVector(MS::columnName(MS::DATA)));
462 1 : td.defineHypercolumn("TiledFlag",3,
463 2 : stringToVector(MS::columnName(MS::FLAG)));
464 1 : td.defineHypercolumn("TiledUVW",2,
465 2 : stringToVector(MS::columnName(MS::UVW)));
466 : }
467 : #endif
468 :
469 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Creating MS=" << MSFileName << LogIO::POST;
470 1 : SetupNewTable newtab(MSFileName, td, Table::New);
471 :
472 : // Set the default Storage Manager to be the Incr one
473 1 : IncrementalStMan incrStMan ("ISMData");
474 1 : newtab.bindAll(incrStMan, true);
475 1 : StandardStMan aipsStMan;
476 :
477 :
478 : #ifdef OLD_CODE
479 : // ORIGINAL CODE
480 : newtab.bindColumn(MS::columnName(MS::ANTENNA2), aipsStMan);
481 : if (useTSM) {
482 : // choose a tile size in the channel direction that is <=10
483 : Int tileSize=(nChan+nChan/10)/(nChan/10+1);
484 : // make the tile about 32k big
485 : TiledColumnStMan tiledStMan1("TiledData",
486 : IPosition(3,nCorr,tileSize,
487 : 2000/nCorr/tileSize));
488 : TiledColumnStMan tiledStMan2("TiledWeight",
489 : IPosition(2,tileSize,
490 : 8000/tileSize));
491 : // Bind the DATA and FLAG columns to the tiled stman
492 : newtab.bindColumn(MS::columnName(MS::DATA),tiledStMan1);
493 : newtab.bindColumn(MS::columnName(MS::FLAG),tiledStMan1);
494 : }
495 : // Change some to aipsStMan as they change every row
496 : newtab.bindColumn(MS::columnName(MS::ANTENNA2),aipsStMan);
497 : newtab.bindColumn(MS::columnName(MS::UVW),aipsStMan);
498 : if (!useTSM) {
499 : newtab.bindColumn(MS::columnName(MS::DATA),aipsStMan);
500 : newtab.bindColumn(MS::columnName(MS::FLAG),aipsStMan);
501 : }
502 : MeasurementSet ms(newtab);
503 : #else
504 : // NEW CODE TO ACCOMODATE VARYING SHAPED COLUMNS
505 1 : if (useTSM) {
506 1 : Int tileSize=nChan/10+1;
507 :
508 : TiledShapeStMan tiledStMan1("TiledData",
509 0 : IPosition(3,nCorr,tileSize,
510 2 : 16384/nCorr/tileSize));
511 : TiledShapeStMan tiledStMan1f("TiledFlag",
512 0 : IPosition(3,nCorr,tileSize,
513 2 : 16384/nCorr/tileSize));
514 : //TiledShapeStMan tiledStMan2("TiledWeight",
515 : // IPosition(3,nCorr,tileSize,
516 : // 16384/nCorr/tileSize));
517 : TiledColumnStMan tiledStMan3("TiledUVW",
518 2 : IPosition(2,3,1024));
519 :
520 : // Bind the DATA and FLAG columns to the tiled stman
521 1 : newtab.bindColumn(MS::columnName(MS::DATA),tiledStMan1);
522 1 : newtab.bindColumn(MS::columnName(MS::FLAG),tiledStMan1f);
523 1 : newtab.bindColumn(MS::columnName(MS::UVW),tiledStMan3);
524 1 : } else {
525 0 : newtab.bindColumn(MS::columnName(MS::DATA),aipsStMan);
526 0 : newtab.bindColumn(MS::columnName(MS::FLAG),aipsStMan);
527 0 : newtab.bindColumn(MS::columnName(MS::UVW),aipsStMan);
528 : }
529 1 : TableLock lock(TableLock::AutoLocking);
530 1 : MeasurementSet ms(newtab,lock);
531 : #endif
532 :
533 : // create all subtables
534 : // we make new tables with 0 rows
535 1 : Table::TableOption option=Table::New;
536 :
537 : // Set up the default subtables for the MS
538 1 : ms.createDefaultSubtables(option);
539 :
540 : // Add some optional columns to the required tables
541 : //ms.spectralWindow().addColumn(ArrayColumnDesc<Int>(
542 : // MSSpectralWindow::columnName(MSSpectralWindow::ASSOC_SPW_ID),
543 : // MSSpectralWindow::columnStandardComment(MSSpectralWindow::ASSOC_SPW_ID)));
544 :
545 : //ms.spectralWindow().addColumn(ArrayColumnDesc<String>(
546 : // MSSpectralWindow::columnName(MSSpectralWindow::ASSOC_NATURE),
547 : // MSSpectralWindow::columnStandardComment(MSSpectralWindow::ASSOC_NATURE)));
548 :
549 1 : ms.spectralWindow().addColumn(ScalarColumnDesc<Int>(
550 : MSSpectralWindow::columnName(MSSpectralWindow::DOPPLER_ID),
551 : MSSpectralWindow::columnStandardComment(MSSpectralWindow::DOPPLER_ID)));
552 :
553 : // Now setup some optional columns::
554 :
555 : // the SOURCE table, 2 extra optional columns needed
556 1 : TableDesc sourceDesc = MSSource::requiredTableDesc();
557 1 : MSSource::addColumnToDesc(sourceDesc,MSSourceEnums::REST_FREQUENCY,1);
558 1 : MSSource::addColumnToDesc(sourceDesc,MSSourceEnums::SYSVEL,1);
559 1 : MSSource::addColumnToDesc(sourceDesc,MSSourceEnums::TRANSITION,1);
560 1 : SetupNewTable sourceSetup(ms.sourceTableName(),sourceDesc,option);
561 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SOURCE),
562 2 : Table(sourceSetup));
563 :
564 : // the DOPPLER table, no optional columns needed
565 1 : TableDesc dopplerDesc = MSDoppler::requiredTableDesc();
566 1 : SetupNewTable dopplerSetup(ms.dopplerTableName(),dopplerDesc,option);
567 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::DOPPLER),
568 2 : Table(dopplerSetup));
569 :
570 : // the SYSCAL table, 1 optional column needed
571 1 : TableDesc syscalDesc = MSSysCal::requiredTableDesc();
572 1 : MSSysCal::addColumnToDesc(syscalDesc,MSSysCalEnums::TSYS,1);
573 1 : SetupNewTable syscalSetup(ms.sysCalTableName(),syscalDesc,option);
574 2 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SYSCAL),
575 2 : Table(syscalSetup));
576 :
577 : // update the references to the subtable keywords
578 1 : ms.initRefs();
579 :
580 : { // Set the TableInfo
581 1 : TableInfo& info(ms.tableInfo());
582 1 : info.setType(TableInfo::type(TableInfo::MEASUREMENTSET));
583 1 : info.setSubType(String("MIRIAD"));
584 1 : info.readmeAddLine("Made with Importmiriad");
585 : }
586 :
587 1 : ms_p=ms;
588 1 : msc_p = new MSColumns(ms_p);
589 1 : } // setupMeasurementSet()
590 :
591 : // ==============================================================================================
592 1 : void Importmiriad::fillObsTables()
593 : {
594 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillObsTables" << LogIO::POST;
595 :
596 : char hline[256];
597 : Int heof;
598 :
599 1 : ms_p.observation().addRow();
600 1 : MSObservationColumns msObsCol(ms_p.observation());
601 :
602 1 : msObsCol.telescopeName().put(0,array_p);
603 1 : msObsCol.observer().put(0,observer_p);
604 1 : msObsCol.project().put(0,project_p);
605 1 : if (array_p == "HATCREEK") {
606 0 : Vector<String> blog(1);
607 0 : blog(0) = "See HISTORY for CARMA observing log";
608 0 : msObsCol.log().put(0,blog);
609 0 : }
610 1 : Vector<Double> range(2);
611 1 : MSColumns msCol(ms_p);
612 1 : range(0) = msCol.time()(0)-1;
613 1 : range(1) = msCol.time()(ms_p.nrow()-1)+1;
614 1 : msObsCol.timeRange().put(0,range);
615 :
616 : // should double buffer history, and search for (e.g.)
617 : // GPAVER: Executed on: 96SEP12:15:40:48.0
618 :
619 : // String date("");
620 : // if (date=="") date="01/01/00";
621 : // Double time=date2mjd(date);
622 :
623 1 : MSHistoryColumns msHisCol(ms_p.history());
624 :
625 1 : String history;
626 1 : Int row=-1;
627 1 : hisopen_c(uv_handle_p,"read");
628 : for (;;) {
629 289 : hisread_c(uv_handle_p,hline,256,&heof);
630 289 : if (heof) break;
631 288 : ms_p.history().addRow();
632 288 : row++;
633 288 : msHisCol.observationId().put(row,0);
634 : // msHisCol.time().put(row,time); // fix the "2000/01/01/24:00:00" bug
635 : // nono, better file a report, it appears to be an aips++ problem
636 288 : msHisCol.priority().put(row,"NORMAL");
637 288 : msHisCol.origin().put(row,"Importmiriad::fillObsTables");
638 288 : msHisCol.application().put(row,"importmiriad");
639 288 : msHisCol.message().put(row,hline);
640 : }
641 1 : hisclose_c(uv_handle_p);
642 1 : } // fillObsTables()
643 :
644 : // ==============================================================================================
645 : //
646 : // Loop over the visibility data and fill the main table of the MeasurementSet
647 : // as you find corr/wcorr's
648 : //
649 1 : void Importmiriad::fillMSMainTable()
650 : {
651 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillMSMainTable" << LogIO::POST;
652 :
653 1 : MSColumns& msc(*msc_p); // Get access to the MS columns, new way
654 1 : Int nCorr = (array_p=="CARMA" ? 1 : npol_p); // # stokes (1 for CARMA for now)
655 1 : Int nChan = MAXCHAN; // max # channels to be written
656 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "nCorr = "<<nCorr<<", nChan = "<< nChan <<LogIO::POST;
657 1 : Int nCat = 3; // # initial flagging categories (fixed at 3)
658 1 : Int iscan = 0;
659 1 : Int ifield_old = 0;
660 :
661 1 : Matrix<Complex> vis(nCorr,nChan);
662 1 : Vector<Float> sigma(nCorr);
663 1 : Vector<String> cat(nCat);
664 1 : cat(0)="FLAG_CMD";
665 1 : cat(1)="ORIGINAL";
666 1 : cat(2)="USER";
667 1 : msc.flagCategory().rwKeywordSet().define("CATEGORY",cat);
668 1 : Cube<Bool> flagCat(nCorr,nChan,nCat,false);
669 1 : Matrix<Bool> flag = flagCat.xyPlane(0); // references flagCat's storage
670 1 : Vector<Float> w1(nCorr), w2(nCorr);
671 :
672 1 : uvrewind_c(uv_handle_p);
673 :
674 1 : nAnt_p.resize(1);
675 1 : nAnt_p[0]=0;
676 :
677 1 : receptorAngle_p.resize(1);
678 1 : Int group, row=-1;
679 : Double interval;
680 1 : Bool lastRowFlag = false;
681 :
682 : //os_p << "Found " << win[0].nspect << " spectral window" << (win[0].nspect>1 ? "s":"") << LogIO::POST;
683 :
684 1 : time_p=0;
685 1 : for (group=0; ; group++) { // loop forever until end-of-file
686 : int nread, nwread;
687 211 : uvread_c(uv_handle_p, preamble, data, flags, MAXCHAN, &nread);
688 : // os_p << "UVREAD: " << data[0] << " " << data[1] << LogIO::POST;
689 :
690 211 : Float baseline = preamble[4];
691 211 : Int ant1 = Int(baseline)/256; // baseline = 256*A1 + A2
692 211 : Int ant2 = Int(baseline) - ant1*256; // mostly A1 <= A2
693 :
694 :
695 211 : if (Debug(9)) os_p << LogIO::DEBUG2 << "UVREAD: " << nread << LogIO::POST;
696 211 : if (nread <= 0) break; // done with reading miriad data
697 210 : if (win[freqSet_p].nspect > 0)
698 210 : uvwread_c(uv_handle_p, wdata, wflags, MAXCHAN, &nwread);
699 : else
700 0 : nwread=0;
701 :
702 : // get time in MJD seconds ; input was in JD
703 210 : Double time = (preamble[3] - 2400000.5) * C::day;
704 210 : if (time>time_p) {
705 14 : if (time_p==0) timeFirst_p=time;
706 : // new time slot (assuming time sorted data)
707 : // update tsys data - TODO
708 : }
709 210 : time_p = time;
710 :
711 : // for MIRIAD, this would always cause a single array dataset,
712 : // but we need to count the antpos occurences to find out
713 : // which array configuration we're in.
714 :
715 210 : if (uvupdate_c(uv_handle_p)) { // aha, something important changed
716 210 : if (Debug(4)) {
717 0 : os_p << LogIO::DEBUG2 << "Record " << group+1 << " uvupdate" << LogIO::POST;
718 : }
719 210 : Tracking(group);
720 : } else {
721 0 : if (Debug(5)) os_p << LogIO::DEBUG2 << "Record " << group << LogIO::POST;
722 : }
723 :
724 : // now that phasem1 has been loaded, apply linelength, if needed
725 210 : if (Qlinecal_p) {
726 0 : linecal(nread,data,phasem1[ant1-1],phasem1[ant2-1]);
727 : // linecal(nwread,wdata,phasem1[ant1-1],phasem1[ant2-1]);
728 : }
729 :
730 210 : nAnt_p[nArray_p-1] = max(nAnt_p[nArray_p-1],ant1); // for MIRIAD, and also
731 210 : nAnt_p[nArray_p-1] = max(nAnt_p[nArray_p-1],ant2);
732 210 : ant1--; ant2--; // make them 0-based for CASA
733 :
734 210 : ant1 += nants_offset_p; // correct for different array offsets
735 210 : ant2 += nants_offset_p;
736 :
737 :
738 : // should ant1 and ant2 be offset with (nArray_p-1)*nant_p ???
739 : // in case there are multiple arrays???
740 : // TODO: code should just assuming single array
741 :
742 210 : Vector<Double> uvw(3);
743 210 : uvw(0) = -preamble[0] * 1e-9; // convert to seconds
744 210 : uvw(1) = -preamble[1] * 1e-9; // MIRIAD uses nanosec
745 210 : uvw(2) = -preamble[2] * 1e-9; // note - sign (CASA vs. MIRIAD convention)
746 210 : uvw *= C::c; // Finally convert to meters for CASA
747 :
748 210 : if (group==0 && Debug(1)) {
749 0 : os_p << LogIO::DEBUG1 << "### First record: " << LogIO::POST;
750 0 : os_p << LogIO::DEBUG1 << "### Preamble: " << preamble[0] << " " <<
751 : preamble[1] << " " <<
752 0 : preamble[2] << " nanosec.(MIRIAD convention)" << LogIO::POST;
753 0 : os_p << LogIO::DEBUG1 << "### uvw: " << uvw(0) << " " <<
754 0 : uvw(1) << " " <<
755 0 : uvw(2) << " meter. (CASA convention)" << LogIO::POST;
756 : }
757 :
758 :
759 : // first construct the data (vis & flag) in a single long array
760 : // containing all spectral windows
761 : // In the (optional) loop over all spectral windows, subsets of
762 : // these arrays will be written out
763 :
764 1050 : for (Int i=0; i<nCorr; i++) {
765 840 : Int count = 0; // index into data[] array
766 840 : if (i>0) uvread_c(uv_handle_p, preamble, data, flags, MAXCHAN, &nread);
767 : //if (group==0) {
768 : // os_p << "pol="<< pol<<", nread="<<nread<<LogIO::POST;
769 : // os_p << "data(500)="<< data[1000] <<", "<< data[1001] <<LogIO::POST;
770 : //}
771 58800840 : for (Int chan=0; chan<nChan; chan++) {
772 :
773 : // miriad uses bl=ant1-ant2, FITS/AIPS/CASA use bl=ant2-ant1
774 : // apart from using -(UVW)'s, the visib need to be conjugated as well
775 58800000 : Bool visFlag = (flags[count/2] == 0) ? false : true;
776 58800000 : Float visReal = +data[count]; count++;
777 58800000 : Float visImag = -data[count]; count++;
778 58800000 : Float wt = 1.0;
779 58800000 : if (!visFlag) wt = -wt;
780 :
781 : // check flags array !! need separate counter (count/2)
782 58800000 : Int pol = corrIndex_p(i);
783 58800000 : flag(pol,chan) = (wt<=0);
784 58800000 : vis(pol,chan) = Complex(visReal,visImag);
785 : } // chan
786 : //if (group==0) os_p << "vis = "<<vis(pol,500)<<LogIO::POST;
787 : } // pol
788 :
789 210 : Int ispw=-1;
790 420 : for (Int sno=0; sno < win[freqSet_p].nspect; sno++) {
791 210 : if (not keep_p(sno)) continue;
792 : // IFs go to separate rows in the MS, pol's do not!
793 210 : ms_p.addRow();
794 210 : row++; ispw++;
795 :
796 : // first fill in values for all the unused columns
797 210 : if (row==0) {
798 1 : ifield_old = ifield;
799 1 : msc.feed1().put(row,0);
800 1 : msc.feed2().put(row,0);
801 1 : msc.flagRow().put(row,false);
802 1 : lastRowFlag = false;
803 1 : msc.scanNumber().put(row,iscan);
804 1 : msc.processorId().put(row,-1);
805 1 : msc.observationId().put(row,0);
806 1 : msc.stateId().put(row,-1);
807 : }
808 :
809 210 : Matrix<Complex> tvis(nCorr,win[freqSet_p].nschan[sno]);
810 210 : Cube<Bool> tflagCat(nCorr,win[freqSet_p].nschan[sno],nCat,false);
811 210 : Matrix<Bool> tflag = tflagCat.xyPlane(0); // references flagCat's storage
812 :
813 210 : Int woffset = win[freqSet_p].ischan[sno]-1;
814 210 : Int wsize = win[freqSet_p].nschan[sno];
815 1050 : for (Int j=0; j<nCorr; j++) {
816 1722000 : for (Int i=0; i< wsize; i++) {
817 1721160 : tvis(j,i) = vis(j,i+woffset);
818 1721160 : tflag(j,i) = flag(j,i+woffset);
819 : }
820 : }
821 :
822 : //if (group==0) os_p<<"tvis="<<tvis(0,500)<<", "<<tvis(1,500)<<LogIO::POST;
823 210 : msc.data().put(row,tvis);
824 210 : msc.flag().put(row,tflag);
825 210 : msc.flagCategory().put(row,tflagCat);
826 :
827 210 : Bool rowFlag = allEQ(flag,true);
828 210 : if (rowFlag != lastRowFlag) {
829 0 : msc.flagRow().put(row,rowFlag);
830 0 : lastRowFlag = rowFlag;
831 : }
832 :
833 210 : msc.antenna1().put(row,ant1);
834 210 : msc.antenna2().put(row,ant2);
835 210 : msc.time().put(row,time); // CARMA did begin of scan.., now middle (2009)
836 210 : msc.timeCentroid().put(row,time); // do we really need this ? flagging/blanking ?
837 :
838 210 : interval = inttime_p;
839 210 : msc.exposure().put(row,interval);
840 210 : msc.interval().put(row,interval);
841 210 : Float chnbw = win[freqSet_p].sdf[sno]*1e9;
842 210 : Float factor = interval * abs(chnbw)/jyperk_p/jyperk_p;
843 : // sigma=sqrt(Tx1*Tx2)/sqrt(chnbw*intTime)*JyPerK;
844 210 : if (Qtsys_p) {
845 0 : w2 = 1.0;
846 0 : if( systemp[ant1] == 0 || systemp[ant2] == 0) {
847 0 : zero_tsys++;
848 0 : w1 = 0.0;
849 : } else {
850 0 : w1 = factor/(systemp[ant1]*systemp[ant2]); // see uvio::uvinfo_variance()
851 0 : w2 = sqrt(systemp[ant1]*systemp[ant2])/sqrt(factor);
852 : }
853 : // os_p << w1 << " " << w2 << " " << jyperk_p << " "<< chnbw<< " "<< interval<< LogIO::POST;
854 0 : msc.weight().put(row,w1);
855 0 : msc.sigma().put(row,w2);
856 : } else {
857 210 : w1=factor/50/50; // Use nominal 50K systemp to keep values similar
858 210 : w2=50/sqrt(factor);
859 210 : msc.weight().put(row,w1);
860 210 : msc.sigma().put(row,w2);
861 : }
862 210 : msc.uvw().put(row,uvw);
863 210 : msc.arrayId().put(row,nArray_p-1);
864 : // calc index into table
865 210 : msc.dataDescId().put(row,ddid_p+ispw);
866 210 : msc.fieldId().put(row,ifield);
867 :
868 :
869 210 : if (ifield_old != ifield)
870 0 : iscan++;
871 210 : ifield_old = ifield;
872 210 : msc.scanNumber().put(row,iscan);
873 :
874 210 : } // sno
875 210 : fcount[ifield]++;
876 210 : } // for(grou) : loop over all visibilities
877 1 : show();
878 1 : if (ms_p.nrow()==0) {
879 0 : os_p<<LogIO::SEVERE<<"No data selected, table is empty!" <<LogIO::POST;
880 : }
881 : else {
882 2 : os_p << "Processed " << group << " visibilities from " << infile_p
883 1 : << (Qlinecal_p ? " (applying linecal)." : " (raw)." )
884 1 : << LogIO::POST;
885 1 : os_p << "Found " << npoint << " pointings with "
886 : << nfield << " unique source/fields, "
887 : << source_p.nelements() << " sources and "
888 1 : << nArray_p << " array"<< (nArray_p>1 ? "s":"")<<"."
889 1 : << LogIO::POST;
890 : }
891 1 : if (Debug(1))
892 0 : os_p << LogIO::DEBUG1 << "nAnt_p contains: " << nAnt_p.nelements() << LogIO::POST;
893 :
894 :
895 1 : } // fillMSMainTable()
896 :
897 1 : void Importmiriad::fillAntennaTable()
898 : {
899 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillAntennaTable" << LogIO::POST;
900 1 : Int nAnt=nants_p;
901 :
902 : #if 0
903 : Int array = nArray_p;
904 : // we don't have 'array' yet, and nAnt_p isnt' big enough....
905 : if (nAnt_p[array]>MAXANT)
906 : throw(AipsError("Too many antennas -- should never occur"));
907 : if (nAnt_p[array]>nants_p)
908 : throw(AipsError("Not all antennas found in antenna table:"));
909 :
910 :
911 : receptorAngle_p[array].resize(2*nAnt);
912 : #endif
913 :
914 : // === Here's a recipe to get them in ITRF format in casapy ===
915 : // But there appears to be some GEODETIC vs. GEOCENTRIC issue
916 : // WGS84
917 : // b=me.measure(me.observatory('carma'),'itrf')
918 : // b2=me.measure(me.observatory('carma'),'wgs84')
919 : // lon=me.getvalue(b)['m0']['value']
920 : // lat=me.getvalue(b2)['m1']['value']
921 : // r=me.getvalue(b)['m2']['value']
922 : // x=r*cos(lat)*cos(lon)
923 : // y=r*cos(lat)*sin(lon)
924 : // z=r*sin(lat)
925 : // print 'arrayXYZ_p(0) =',x,';'
926 : // print 'arrayXYZ_p(1) =',y,';'
927 : // print 'arrayXYZ_p(2) =',z,';'
928 1 : arrayXYZ_p.resize(3);
929 1 : if (array_p == "HATCREEK" || array_p == "BIMA") { // Array center:
930 0 : arrayXYZ_p(0) = -2523862.04;
931 0 : arrayXYZ_p(1) = -4123592.80;
932 0 : arrayXYZ_p(2) = 4147750.37;
933 1 : } else if (array_p == "ATCA") {
934 1 : arrayXYZ_p(0) = -4750915.837;
935 1 : arrayXYZ_p(1) = 2792906.182;
936 1 : arrayXYZ_p(2) = -3200483.747;
937 0 : } else if (array_p == "OVRO" || array_p == "CARMA") {
938 0 : arrayXYZ_p(0) = -2397389.65197;
939 0 : arrayXYZ_p(1) = -4482068.56252;
940 0 : arrayXYZ_p(2) = 3843528.41479;
941 : } else {
942 0 : os_p << LogIO::WARN<< "unknown array position for "<<array_p<<LogIO::POST;
943 0 : arrayXYZ_p = 0.0;
944 : }
945 1 : if(Debug(3)) os_p << LogIO::DEBUG2 << "number of antennas ="<<nAnt<<LogIO::POST;
946 1 : if(Debug(3)) os_p << LogIO::DEBUG2 << "array ref pos:"<<arrayXYZ_p<<LogIO::POST;
947 :
948 1 : String timsys = "TAI"; // assume, for now ....
949 :
950 : // store the time keywords ; again, miriad doesn't have this (yet)
951 : // check w/ uvfitsfiller again
952 :
953 : //save value to set time reference frame elsewhere
954 1 : timsys_p=timsys;
955 :
956 : // Antenna diamater:
957 : // Should check the 'antdiam' UV variable, but it doesn't appear to
958 : // exist in our CARMA datasets.
959 : // So, fill in some likely values
960 1 : Float diameter=25; //# most common size (:-)
961 1 : if (array_p=="ATCA") diameter=22; //# only at 'low' freq !!
962 1 : if (array_p=="HATCREEK") diameter=6;
963 1 : if (array_p=="BIMA") diameter=6;
964 1 : if (array_p=="CARMA") diameter=8;
965 1 : if (array_p=="OVRO") diameter=10;
966 :
967 1 : if (nAnt == 15 && array_p=="OVRO") {
968 0 : os_p << "CARMA array (6 OVRO, 9 BIMA) assumed" << LogIO::POST;
969 0 : array_p = "CARMA";
970 1 : } else if (nAnt == 23 && array_p=="OVRO") {
971 0 : os_p << "CARMA array (6 OVRO, 9 BIMA, 8 SZA) assumed" << LogIO::POST;
972 0 : array_p = "CARMA";
973 1 : } else if (array_p=="CARMA") {
974 0 : os_p << "Hurray, CARMA data; version " << version_p << " with " << nAnt
975 0 : << " antennas" << LogIO::POST;
976 1 : } else if (array_p=="ATCA") {
977 1 : os_p <<"Found ATCA data with " << nAnt << " antennas" << LogIO::POST;
978 : } else
979 0 : os_p << "Ant configuration not supported yet" << LogIO::POST;
980 :
981 1 : MSAntennaColumns& ant(msc_p->antenna());
982 1 : Vector<Double> antXYZ(3);
983 :
984 : // add antenna info to table
985 1 : if (nArray_p == 0) { // check if needed
986 1 : ant.setPositionRef(MPosition::ITRF);
987 : //ant.setPositionRef(MPosition::WGS84);
988 : }
989 1 : Int row=ms_p.antenna().nrow()-1;
990 :
991 1 : if (Debug(2)) os_p << LogIO::DEBUG2 << "Importmiriad::fillAntennaTable row=" << row+1
992 0 : << " array " << nArray_p+1 << LogIO::POST;
993 :
994 7 : for (Int i=0; i<nAnt; i++) {
995 :
996 6 : ms_p.antenna().addRow();
997 6 : row++;
998 6 : if (array_p=="OVRO" || array_p=="BIMA" || array_p=="HATCREEK" || array_p=="CARMA") {
999 0 : if (i<6)
1000 0 : ant.dishDiameter().put(row,10.4); // OVRO
1001 0 : else if (i<15)
1002 0 : ant.dishDiameter().put(row,6.1); // BIMA or HATCREEK
1003 : else
1004 0 : ant.dishDiameter().put(row,3.5); // SZA
1005 : } else {
1006 6 : ant.dishDiameter().put(row,diameter); // others
1007 : }
1008 6 : antXYZ(0) = antpos[i]; //# these are now in nano-sec
1009 6 : antXYZ(1) = antpos[i+nAnt];
1010 6 : antXYZ(2) = antpos[i+nAnt*2];
1011 6 : antXYZ *= 1e-9 * C::c; //# and now in meters
1012 6 : if (Debug(2)) os_p << LogIO::DEBUG2 << "Ant " << i+1 << ":" << antXYZ << " (m)." << LogIO::POST;
1013 :
1014 6 : String mount; // really should consult
1015 6 : switch (mount_p) { // the "mount" uv-variable
1016 6 : case 0: mount="ALT-AZ"; break;
1017 0 : case 1: mount="EQUATORIAL"; break;
1018 0 : case 2: mount="X-Y"; break;
1019 0 : case 3: mount="ORBITING"; break;
1020 0 : case 4: mount="BIZARRE"; break;
1021 : // case 5: mount="SPACE-HALCA"; break;
1022 0 : default: mount="UNKNOWN"; break;
1023 : }
1024 6 : ant.mount().put(row,mount);
1025 6 : ant.flagRow().put(row,false);
1026 12 : String antName = "C";
1027 6 : if (array_p=="ATCA") antName="CA0";
1028 6 : antName += String::toString(i+1);
1029 6 : ant.name().put(row,antName);
1030 6 : ant.station().put(row,"ANT" + String::toString(i+1)); // unknown PADs, so for now ANT#
1031 6 : ant.type().put(row,"GROUND-BASED");
1032 :
1033 12 : Vector<Double> offsets(3);
1034 6 : offsets=0.0;
1035 : // store absolute positions, with all offsets 0
1036 :
1037 : #if 1
1038 : // from MirFiller; but why we're rotating this? To reverse miriad's rotation
1039 : // of y-axis to local East
1040 6 : Double arrayLong = atan2(arrayXYZ_p(1),arrayXYZ_p(0));
1041 12 : Matrix<Double> posRot = Rot3D(2,arrayLong);
1042 6 : antXYZ = product(posRot,antXYZ);
1043 : #endif
1044 :
1045 : #if 1
1046 : // This doesn't work because miriad calculated the relative positions with
1047 : // respect to the first antenna with non zero coordinates,
1048 : // not the array reference position. This makes it impossible to invert exactly
1049 6 : ant.position().put(row,arrayXYZ_p+antXYZ);
1050 : #else
1051 : //test
1052 : ant.position().put(row,arrayXYZ_p);
1053 : #endif
1054 6 : ant.offset().put(row,offsets);
1055 :
1056 : // store the angle for use in the feed table
1057 : // receptorAngle_p[array](2*i+0)=polangleA(i)*C::degree;
1058 : // receptorAngle_p[array](2*i+1)=polangleB(i)*C::degree;
1059 6 : }
1060 : // ant.position().rwKeywordSet().define("MEASURE_REFERENCE","ITRF");
1061 :
1062 1 : nArray_p++;
1063 1 : nAnt_p.resize(nArray_p);
1064 1 : nAnt_p[nArray_p-1] = 0;
1065 1 : if (Debug(3) && nArray_p > 1)
1066 0 : os_p << LogIO::DEBUG2 << nAnt_p[nArray_p-2] << LogIO::POST;
1067 :
1068 1 : if (nArray_p > 1) return;
1069 :
1070 : // now do some things which only need to happen the first time around
1071 :
1072 : // store these items in non-standard keywords for now
1073 : //
1074 1 : String arrnam = array_p;
1075 1 : ant.name().rwKeywordSet().define("ARRAY_NAME",arrnam);
1076 1 : ant.position().rwKeywordSet().define("ARRAY_POSITION",arrayXYZ_p);
1077 :
1078 :
1079 : // fill the array table entry
1080 : // this assumes there is one AN table for each (sub)array index encountered.
1081 :
1082 : //PJT ms_p.array().addRow();
1083 : // array is now gone, there is an array_id in the main MS table for
1084 : // id purposes. We store the ARRAY_POSITION as a non-standard keyword
1085 : // with the POSITION collumn in the ANTENNA table (see above)
1086 : #if 0
1087 : MSArrayColumns arr(ms_p.array());
1088 : arr.name().put(array,array_p);
1089 : arr.position().put(array,arrayXYZ);
1090 : arr.position().rwKeywordSet().define("MEASURE_REFERENCE","ITRF");
1091 : #endif
1092 1 : } // fillAntennaTable
1093 :
1094 : // ==============================================================================================
1095 15 : void Importmiriad::fillSyscalTable()
1096 : {
1097 : //if (Debug(1)) os_p << "Importmiriad::fillSyscalTable" << LogIO::POST;
1098 :
1099 15 : MSSysCalColumns& msSys(msc_p->sysCal());
1100 15 : Vector<Float> Systemp(1); // should we set both receptors same?
1101 15 : Int row = ms_p.sysCal().nrow();
1102 :
1103 : // if (Debug(1))
1104 : // for (Int i=0; i<nants_p; i++)
1105 : // os_p << "SYSTEMP: " << i << ": " << systemp[i] << LogIO::POST;
1106 :
1107 30 : for (Int j=0; j<win[freqSet_p].nspect; j++) {
1108 105 : for (Int i=0; i<nants_p; i++) {
1109 90 : ms_p.sysCal().addRow();
1110 90 : msSys.antennaId().put(row,i); // i, or i+nants_offset_p ????
1111 90 : msSys.feedId().put(row,0);
1112 90 : msSys.spectralWindowId().put(row,j); // all of them for now .....
1113 90 : msSys.time().put(row,time_p);
1114 90 : msSys.interval().put(row,-1.0);
1115 :
1116 90 : Systemp(0) = systemp[i+j*nants_p];
1117 90 : msSys.tsys().put(row,Systemp);
1118 90 : row++;
1119 : }
1120 : }
1121 :
1122 :
1123 :
1124 : // this may actually be a nasty problem for MIRIAD datasets that are not
1125 : // timesorted. A temporary table needs to be written with all records,
1126 : // which then needs to be sorted and 'recomputed'
1127 15 : }
1128 :
1129 : // ==============================================================================================
1130 1 : void Importmiriad::fillSpectralWindowTable(String vel)
1131 : {
1132 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillSpectralWindowTable" << LogIO::POST;
1133 :
1134 1 : MSSpWindowColumns& msSpW(msc_p->spectralWindow());
1135 1 : MSDataDescColumns& msDD(msc_p->dataDescription());
1136 1 : MSPolarizationColumns& msPol(msc_p->polarization());
1137 1 : MSDopplerColumns& msDop(msc_p->doppler());
1138 :
1139 1 : Int nCorr = ( array_p=="CARMA" ? 1 : npol_p); // CARMA wants 1 polarization
1140 : Int i, side;
1141 1 : Double BW = 0.0;
1142 :
1143 1 : MDirection::Types dirtype = epochRef_p; // MDirection::B1950 or MDirection::J2000;
1144 2 : MEpoch ep(Quantity(time_p, "s"), MEpoch::UTC);
1145 : // ERROR:: type specifier omitted for parameter in older AIPS++, now works in CASA
1146 1 : MPosition obspos(MVPosition(arrayXYZ_p), MPosition::ITRF);
1147 : //MPosition obspos(MVPosition(arrayXYZ_p), MPosition::WGS84);
1148 2 : MDirection dir(Quantity(ra_p, "rad"), Quantity(dec_p, "rad"), dirtype);
1149 1 : MeasFrame frame(ep, obspos, dir);
1150 :
1151 1 : String velsys = vel;
1152 : // Keep previous default, for ATCA leave at TOPO (multi-source data)
1153 1 : if (array_p!="ATCA" && velsys=="") velsys="LSRK";
1154 :
1155 1 : Bool convert=true;
1156 : MFrequency::Types freqsys_p;
1157 1 : if (velsys=="LSRK") {
1158 0 : freqsys_p = MFrequency::LSRK; // LSRD vs. LSRK
1159 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << "USE_LSRK" << LogIO::POST;
1160 1 : } else if (velsys=="LSRD"){
1161 0 : freqsys_p = MFrequency::LSRD; // LSRD vs. LSRK
1162 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << "USE_LSRD" << LogIO::POST;
1163 : } else {
1164 1 : freqsys_p = MFrequency::TOPO; // use TOPO if unspecified
1165 1 : convert=false;
1166 : }
1167 :
1168 : MFrequency::Convert tolsr(MFrequency::TOPO,
1169 1 : MFrequency::Ref(freqsys_p, frame)); // LSRD vs. LSRK
1170 : // fill out the polarization info (only 1 entry allowed for now)
1171 1 : ms_p.polarization().addRow();
1172 1 : msPol.numCorr().put(0,nCorr);
1173 1 : msPol.corrType().put(0,corrType_p);
1174 1 : msPol.corrProduct().put(0,corrProduct_p);
1175 1 : msPol.flagRow().put(0,false);
1176 :
1177 : // fill out doppler table (only 1 entry needed, CARMA data only identify 1 line :-(
1178 1 : if (array_p=="CARMA") {
1179 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad:: now writing Doppler table " << LogIO::POST;
1180 0 : for (i=0; i<win[0].nspect; i++) {
1181 0 : ms_p.doppler().addRow();
1182 0 : msDop.dopplerId().put(i,i);
1183 0 : msDop.sourceId().put(i,-1); // how the heck..... for all i guess
1184 0 : msDop.transitionId().put(i,-1);
1185 0 : msDop.velDefMeas().put(i,MDoppler(Quantity(0),MDoppler::RADIO));
1186 : }
1187 : }
1188 1 : Int ddid=-1;
1189 2 : for (Int k=0; k < nFreqSet_p; k++) {
1190 2 : for (Int i=0; i < win[k].nspect; i++) {
1191 1 : if(not keep_p(i)) continue;
1192 1 : ddid++;
1193 1 : Int n = win[k].nschan[i];
1194 1 : Vector<Double> f(n), w(n),cs(n);
1195 :
1196 1 : ms_p.spectralWindow().addRow();
1197 1 : ms_p.dataDescription().addRow();
1198 :
1199 1 : msDD.spectralWindowId().put(ddid,ddid);
1200 1 : msDD.polarizationId().put(ddid,0);
1201 1 : msDD.flagRow().put(ddid,false);
1202 :
1203 1 : msSpW.numChan().put(ddid,win[k].nschan[i]);
1204 1 : BW = 0.0;
1205 1 : Double fwin = win[k].sfreq[i]*1e9;
1206 1 : if (convert) {
1207 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Fwin: OBS=" << fwin/1e9;
1208 0 : fwin = tolsr(fwin).getValue().getValue();
1209 0 : if (Debug(1)) os_p << LogIO::DEBUG1 << " LSR=" << fwin/1e9 << LogIO::POST;
1210 : }
1211 2050 : for (Int j=0; j < win[k].nschan[i]; j++) {
1212 2049 : f(j) = fwin + j * win[k].sdf[i] * 1e9;
1213 2049 : w(j) = abs(win[k].sdf[i]*1e9);
1214 2049 : cs(j) = win[k].sdf[i]*1e9;
1215 2049 : BW += w(j);
1216 : }
1217 :
1218 1 : msSpW.chanFreq().put(ddid,f);
1219 1 : if (i<win[k].nspect) {
1220 : // I think restfreq should just be in source table,
1221 : // but leave for now if it makes sense
1222 1 : if (win[k].restfreq[i]>0 && freqsys_p!=MFrequency::TOPO) {
1223 0 : msSpW.refFrequency().put(ddid,win[k].restfreq[i]*1e9);
1224 : } else {
1225 1 : msSpW.refFrequency().put(ddid,win[k].sfreq[i]*1e9);
1226 : }
1227 : } else
1228 0 : msSpW.refFrequency().put(ddid,freq_p); // no reference for wide band???
1229 :
1230 1 : msSpW.resolution().put(ddid,w);
1231 1 : msSpW.chanWidth().put(ddid,cs); // MSSelection wants width<0 if spectrum inverted
1232 1 : msSpW.effectiveBW().put(ddid,w);
1233 1 : msSpW.totalBandwidth().put(ddid,BW);
1234 1 : Int ifchain = win[k].chain[i];
1235 1 : msSpW.ifConvChain().put(ddid,ifchain);
1236 : // can also do it implicitly via Measures you give to the freq's
1237 1 : msSpW.measFreqRef().put(ddid,freqsys_p);
1238 1 : if (i<win[k].nspect && array_p=="CARMA")
1239 0 : msSpW.dopplerId().put(ddid,i); // CARMA has only 1 ref freq line
1240 : else
1241 1 : msSpW.dopplerId().put(ddid,-1); // no ref
1242 :
1243 1 : if (win[k].sdf[i] > 0) side = 1;
1244 1 : else if (win[k].sdf[i] < 0) side = -1;
1245 0 : else side = 0;
1246 :
1247 1 : switch (win[k].code[i]) {
1248 1 : case 'N':
1249 1 : msSpW.netSideband().put(ddid,side);
1250 1 : msSpW.freqGroup().put(ddid,1);
1251 1 : msSpW.freqGroupName().put(ddid,"MULTI-CHANNEL-DATA");
1252 1 : break;
1253 0 : case 'W':
1254 0 : msSpW.netSideband().put(ddid,side);
1255 0 : msSpW.freqGroup().put(ddid,3);
1256 0 : msSpW.freqGroupName().put(ddid,"SIDE-BAND-AVERAGE");
1257 0 : break;
1258 0 : case 'S':
1259 0 : msSpW.netSideband().put(ddid,side);
1260 0 : msSpW.freqGroup().put(ddid,2);
1261 0 : msSpW.freqGroupName().put(ddid,"MULTI-CHANNEL-AVG");
1262 0 : break;
1263 0 : default:
1264 0 : throw(AipsError("Bad code for a spectral window"));
1265 : break;
1266 : }
1267 1 : }
1268 : }
1269 1 : }
1270 :
1271 : // ==============================================================================================
1272 1 : void Importmiriad::fillFieldTable()
1273 : {
1274 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillFieldTable" << LogIO::POST;
1275 :
1276 : // set the DIRECTION MEASURE REFERENCE for appropriate columns
1277 : // but note we're not varying them accross rows
1278 : /// MDirection::Types epochRef=MDirection::J2000;
1279 : /// if (nearAbs(epoch_p,1950.0,0.01)) epochRef=MDirection::B1950;
1280 1 : msc_p->setDirectionRef(epochRef_p);
1281 :
1282 1 : MSFieldColumns& msField(msc_p->field());
1283 :
1284 1 : Vector<Double> radec(2), pm(2);
1285 1 : Vector<MDirection> radecMeas(1);
1286 : Int fld;
1287 : Double cosdec;
1288 :
1289 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillFieldTable() adding " << nfield
1290 0 : << " fields" << LogIO::POST;
1291 :
1292 1 : pm = 0; // Proper motion is zero
1293 :
1294 1 : if (nfield == 0) { // if no pointings found, say there is 1
1295 0 : os_p << "Warning: no dra/ddec pointings found, creating 1." << LogIO::POST;
1296 0 : nfield = npoint = 1;
1297 0 : dra[0] = ddec[0] = 0.0;
1298 : }
1299 :
1300 2 : for (fld = 0; fld < nfield; fld++) {
1301 :
1302 1 : ms_p.field().addRow();
1303 :
1304 1 : if (Debug(1))
1305 0 : os_p << LogIO::DEBUG1 << "FLD: " << fld << " " << field[fld] << " " << source_p[field[fld]] << LogIO::POST;
1306 :
1307 1 : msField.sourceId().put(fld,field[fld]);
1308 1 : msField.name().put(fld,source_p[field[fld]]); // this is the source name
1309 1 : msField.code().put(fld,purpose_p[field[fld]]);
1310 1 : msField.numPoly().put(fld,0);
1311 :
1312 1 : cosdec = cos(dec[fld]);
1313 1 : radec(0) = ra[fld] + dra[fld]/cosdec; // RA, in radians
1314 1 : radec(1) = dec[fld] + ddec[fld]; // DEC, in radians
1315 1 : radecMeas(0).set(MVDirection(radec(0), radec(1)), MDirection::Ref(epochRef_p));
1316 :
1317 1 : msField.delayDirMeasCol().put(fld,radecMeas);
1318 1 : msField.phaseDirMeasCol().put(fld,radecMeas);
1319 1 : msField.referenceDirMeasCol().put(fld,radecMeas);
1320 : // put in best guess for time of position, this is not the coord epoch,
1321 : // rather it is meant to cope with moving objects which have a rate.
1322 1 : msField.time().put(fld,timeFirst_p);
1323 : }
1324 1 : }
1325 :
1326 : // ==============================================================================================
1327 1 : void Importmiriad::fillSourceTable()
1328 : {
1329 : // According to MS2 specs we should have a source table with
1330 : // an entry for every source/spectral window combination that occurs.
1331 : // For the moment we ignore spectral window (and TIME) as an index and
1332 : // just use -1. This means there is only a single rest frequency.
1333 :
1334 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillSourceTable" << LogIO::POST;
1335 1 : Int ns = 0;
1336 : Int skip;
1337 :
1338 1 : MSSourceColumns& msSource(msc_p->source());
1339 :
1340 1 : Vector<Double> radec(2);
1341 1 : Vector<Double> restFreq(1),sysvel(1);
1342 1 : sysvel(0)=0;
1343 1 : Int m=win[0].nspect;
1344 1 : restFreq(0)=0;
1345 : // pick first non zero restFreq
1346 2 : for (Int i=0; i<m; i++) {
1347 1 : if (win[0].restfreq[i]>0) {
1348 0 : restFreq(0) = win[0].restfreq[i] * 1e9; // convert from GHz to Hz
1349 0 : break;
1350 : }
1351 : }
1352 :
1353 : //String key("MEASURE_REFERENCE");
1354 : //msSource.restFrequency().rwKeywordSet().define(key,"REST");
1355 :
1356 1 : if (Debug(1)) {
1357 0 : os_p << LogIO::DEBUG1 << "Importmiriad::fillSourceTable() querying " << source_p.nelements()
1358 0 : << " sources" << LogIO::POST;
1359 0 : os_p << LogIO::DEBUG1 << source_p << LogIO::POST;
1360 : }
1361 :
1362 :
1363 : //
1364 2 : for (uInt src=0; src < source_p.nelements(); src++) {
1365 :
1366 1 : skip = 0;
1367 1 : for (uInt i=0; i<src; i++) { // loop over sources to avoid duplicates
1368 0 : if (source_p[i] == source_p[src]) {
1369 0 : skip=1; cerr<<"Found duplicate source name! Fix code!"<<LogIO::POST;
1370 0 : break;
1371 : }
1372 : }
1373 :
1374 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "source : " << source_p[src] << " " << skip << LogIO::POST;
1375 :
1376 1 : if (skip) continue; // if seen before, don't add it again
1377 1 : ms_p.source().addRow();
1378 :
1379 1 : radec(0) = ras_p[src];
1380 1 : radec(1) = decs_p[src];
1381 :
1382 1 : msSource.sourceId().put(ns,src);
1383 1 : msSource.name().put(ns,source_p[src]);
1384 1 : msSource.spectralWindowId().put(src,-1); // set valid for all
1385 1 : msSource.direction().put(ns,radec);
1386 1 : msSource.numLines().put(ns,1);
1387 1 : msSource.restFrequency().put(ns,restFreq);
1388 1 : msSource.time().put(ns,0.0); // valid for all times
1389 1 : msSource.interval().put(ns,0); // valid forever
1390 1 : msSource.sysvel().put(ns,sysvel);
1391 : // TODO?
1392 : // missing position/sysvel/transition in the produced MS/SOURCE sub-table ??
1393 :
1394 : // listobs complains:
1395 : // No systemic velocity information found in SOURCE table.
1396 :
1397 1 : ns++;
1398 : }
1399 :
1400 : // TODO: #sources wrong if you take raw miriad before noise taken out
1401 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillSourceTable() added " << ns << " sources" << LogIO::POST;
1402 1 : }
1403 :
1404 : // ==============================================================================================
1405 1 : void Importmiriad::fillFeedTable()
1406 : {
1407 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fillFeedTable" << LogIO::POST;
1408 :
1409 1 : MSFeedColumns msfc(ms_p.feed());
1410 :
1411 : // find out the POLARIZATION_TYPE
1412 : // In the fits files we handle there can be only a single, uniform type
1413 : // of polarization so the following should work.
1414 1 : MSPolarizationColumns& msPolC(msc_p->polarization());
1415 :
1416 1 : Int numCorr=msPolC.numCorr()(0);
1417 1 : Vector<String> rec_type(2); rec_type="";
1418 1 : if (corrType_p(0)>=Stokes::RR && corrType_p(numCorr-1)<=Stokes::LL) {
1419 0 : rec_type(0)="R"; rec_type(1)="L";
1420 : }
1421 1 : if (corrType_p(0)>=Stokes::XX && corrType_p(numCorr-1)<=Stokes::YY) {
1422 1 : rec_type(0)="X"; rec_type(1)="Y";
1423 : }
1424 :
1425 1 : Matrix<Complex> polResponse(2,2);
1426 1 : polResponse=0.; polResponse(0,0)=polResponse(1,1)=1.;
1427 1 : Matrix<Double> offset(2,2); offset=0.;
1428 1 : Vector<Double> position(3); position=0.;
1429 1 : Vector<Double> ra(2);
1430 1 : ra = 0.0;
1431 1 : if (array_p == "ATCA") {
1432 1 : ra(0)=C::pi/4;
1433 : // 7mm feed is different; assumes all spectral windows are in same band for now...
1434 1 : if (win[freqSet_p].sfreq[0]>30. && win[freqSet_p].sfreq[0]<50.) ra(0)+=C::pi/2;
1435 1 : ra(1)=ra(0)+C::pi/2;
1436 : }
1437 :
1438 : // fill the feed table
1439 : // will only do UP TO the largest antenna referenced in the dataset
1440 1 : Int row=-1;
1441 1 : if (Debug(3)) os_p << LogIO::DEBUG2 << nAnt_p.nelements() << LogIO::POST;
1442 2 : for (Int arr=0; arr< (Int)nAnt_p.nelements(); arr++) {
1443 1 : if (Debug(3)) os_p << LogIO::DEBUG2 << nAnt_p[arr] << LogIO::POST;
1444 7 : for (Int ant=0; ant<nAnt_p[arr]; ant++) {
1445 6 : ms_p.feed().addRow(); row++;
1446 :
1447 6 : msfc.antennaId().put(row,ant);
1448 6 : msfc.beamId().put(row,-1);
1449 6 : msfc.feedId().put(row,0);
1450 6 : msfc.interval().put(row,DBL_MAX);
1451 :
1452 : // msfc.phasedFeedId().put(row,-1); // now optional
1453 6 : msfc.spectralWindowId().put(row,-1);
1454 6 : msfc.time().put(row,0.);
1455 6 : msfc.numReceptors().put(row,2);
1456 6 : msfc.beamOffset().put(row,offset);
1457 6 : msfc.polarizationType().put(row,rec_type);
1458 6 : msfc.polResponse().put(row,polResponse);
1459 6 : msfc.position().put(row,position);
1460 : // fix these when incremental array building is ok.
1461 : // although for CARMA this would never change ....
1462 6 : msfc.receptorAngle().put(row,ra);
1463 : // msfc.receptorAngle().put(row,receptorAngle_p[arr](Slice(2*ant,2)));
1464 : }
1465 : }
1466 1 : }
1467 :
1468 : // ==============================================================================================
1469 1 : void Importmiriad::fixEpochReferences() {
1470 :
1471 1 : if (Debug(1)) os_p << LogIO::DEBUG1 << "Importmiriad::fixEpochReferences" << LogIO::POST;
1472 :
1473 1 : if (timsys_p=="IAT") timsys_p="TAI";
1474 1 : if (timsys_p=="UTC" || timsys_p=="TAI") {
1475 1 : String key("MEASURE_REFERENCE");
1476 1 : MSColumns msc(ms_p);
1477 1 : msc.time().rwKeywordSet().define(key,timsys_p);
1478 1 : msc.feed().time().rwKeywordSet().define(key,timsys_p);
1479 1 : msc.field().time().rwKeywordSet().define(key,timsys_p);
1480 : // Fits obslog time is probably local time instead of TAI or UTC
1481 : //PJT msc.obsLog().time().rwKeywordSet().define(key,timsys_p);
1482 1 : } else {
1483 0 : if (timsys_p!="")
1484 0 : cerr << "Unhandled time reference frame: "<<timsys_p<<LogIO::POST;
1485 : }
1486 1 : }
1487 :
1488 : //
1489 : // track some important uv variables to get notified when they change
1490 : //
1491 : // ==============================================================================================
1492 211 : void Importmiriad::Tracking(int record)
1493 : {
1494 211 : if (Debug(3)) os_p << LogIO::DEBUG2 << "Importmiriad::Tracking" << LogIO::POST;
1495 :
1496 : char vtype[10], vdata[10];
1497 : int vlen, vupd, idat, vupd1, vupd2, vupd3, j, k;
1498 : // float dx, dy;
1499 : // Float rdat;
1500 : // Double ddat;
1501 :
1502 211 : if (record < 0) { // first time around: set variables to track
1503 1 : uvtrack_c(uv_handle_p,"nschan","u"); // narrow lines
1504 1 : uvtrack_c(uv_handle_p,"nspect","u"); // window averages
1505 1 : uvtrack_c(uv_handle_p,"ischan","u");
1506 1 : uvtrack_c(uv_handle_p,"sdf","u");
1507 1 : uvtrack_c(uv_handle_p,"sfreq","u"); // changes a lot (doppler)
1508 :
1509 1 : uvtrack_c(uv_handle_p,"restfreq","u"); // never really changes....
1510 1 : uvtrack_c(uv_handle_p,"freq","u"); // never really changes....
1511 1 : uvtrack_c(uv_handle_p,"ifchain","u"); // optional
1512 :
1513 :
1514 1 : uvtrack_c(uv_handle_p,"nwide","u");
1515 1 : uvtrack_c(uv_handle_p,"wfreq","u");
1516 1 : uvtrack_c(uv_handle_p,"wwidth","u");
1517 :
1518 1 : uvtrack_c(uv_handle_p,"antpos","u"); // array's
1519 1 : uvtrack_c(uv_handle_p,"pol","u"); // pol's
1520 1 : uvtrack_c(uv_handle_p,"dra","u"); // fields
1521 1 : uvtrack_c(uv_handle_p,"ddec","u"); // fields
1522 :
1523 1 : uvtrack_c(uv_handle_p,"ra","u"); // source position
1524 1 : uvtrack_c(uv_handle_p,"dec","u"); // source position
1525 :
1526 1 : uvtrack_c(uv_handle_p,"inttime","u");
1527 :
1528 1 : if (Qlinecal_p)
1529 0 : uvtrack_c(uv_handle_p,"phasem1","u"); // linelength meaurements
1530 :
1531 : // weather:
1532 : // uvtrack_c(uv_handle_p,"airtemp","u");
1533 : // uvtrack_c(uv_handle_p,"dewpoint","u");
1534 : // uvtrack_c(uv_handle_p,"relhumid","u");
1535 : // uvtrack_c(uv_handle_p,"winddir","u");
1536 : // uvtrack_c(uv_handle_p,"windmph","u");
1537 :
1538 1 : return;
1539 : }
1540 :
1541 : // here is all the special tracking code...
1542 :
1543 210 : check_window(); // check if the freq setup has changed
1544 :
1545 210 : uvprobvr_c(uv_handle_p,"pol",vtype,&vlen,&vupd);
1546 210 : if (vupd && npol_p==1) {
1547 0 : uvrdvr_c(uv_handle_p,H_INT,"pol",(char *)&idat, NULL, 1);
1548 0 : if (idat != pol_p[0])
1549 0 : os_p << LogIO::WARN<<"polarization changed to " << *pol_p << LogIO::POST;
1550 0 : pol_p[0] = idat;
1551 : }
1552 :
1553 210 : uvprobvr_c(uv_handle_p,"npol",vtype,&vlen,&vupd);
1554 210 : if (vupd) {
1555 1 : uvrdvr_c(uv_handle_p,H_INT,"npol",(char *)&idat, NULL, 1);
1556 1 : if (idat != npol_p)
1557 0 : throw(AipsError("Cannot handle a changing npol yet"));
1558 : }
1559 :
1560 210 : uvprobvr_c(uv_handle_p,"inttime",vtype,&vlen,&vupd);
1561 210 : if (vupd) {
1562 1 : uvgetvr_c(uv_handle_p,H_REAL,"inttime",(char *)&inttime_p,1);
1563 : }
1564 210 : uvprobvr_c(uv_handle_p,"antpos",vtype,&vlen,&vupd);
1565 210 : if (vupd && record) {
1566 0 : if (Qarrays_p)
1567 0 : nants_offset_p += nants_p; // increment from size of previous array
1568 0 : uvgetvr_c(uv_handle_p,H_INT, "nants", (char *)&nants_p,1);
1569 0 : uvgetvr_c(uv_handle_p,H_DBLE,"antpos",(char *)antpos,3*nants_p);
1570 0 : if (Debug(2)) {
1571 0 : os_p << LogIO::DEBUG2 << "Found " << nants_p << " antennas for array "
1572 0 : << nArray_p << LogIO::POST;
1573 0 : for (int i=0; i<nants_p; i++) {
1574 0 : os_p << antpos[i] << " " <<
1575 0 : antpos[i+nants_p] << " " <<
1576 0 : antpos[i+nants_p*2] << LogIO::POST;
1577 : }
1578 : }
1579 0 : if (Debug(2)) os_p << LogIO::DEBUG2 << "Warning: antpos changed at record " << record << LogIO::POST;
1580 0 : if (Qarrays_p)
1581 0 : fillAntennaTable();
1582 : }
1583 :
1584 210 : if (Qlinecal_p) {
1585 0 : uvprobvr_c(uv_handle_p,"phasem1",vtype,&vlen,&vupd);
1586 0 : if (vupd) {
1587 : // lets assume (hope?) that nants_p didn't change, it better not.
1588 0 : uvgetvr_c(uv_handle_p,H_REAL,"phasem1",(char *)phasem1,nants_p);
1589 : // os_p << "PHASEM1: " << phasem1[0] << " " << phasem1[1] << " ...\n";
1590 : }
1591 : }
1592 :
1593 210 : if (win[freqSet_p].nspect > 0) {
1594 210 : uvprobvr_c(uv_handle_p,"systemp",vtype,&vlen,&vupd);
1595 210 : if (vupd) {
1596 14 : uvgetvr_c(uv_handle_p,H_REAL,"systemp",(char *)systemp,nants_p*win[freqSet_p].nspect);
1597 14 : if (Debug(3)) {
1598 0 : os_p << LogIO::DEBUG2 << "Found systemps (new scan)" ;
1599 0 : for (Int i=0; i<nants_p; i++) os_p << systemp[i] << " ";
1600 0 : os_p << LogIO::POST;
1601 : }
1602 14 : fillSyscalTable();
1603 : }
1604 : } else {
1605 0 : uvprobvr_c(uv_handle_p,"wsystemp",vtype,&vlen,&vupd);
1606 0 : if (vupd) {
1607 0 : uvgetvr_c(uv_handle_p,H_REAL,"wsystemp",(char *)systemp,nants_p);
1608 0 : if (Debug(3)) {
1609 0 : os_p << LogIO::DEBUG2 << "Found wsystemps (new scan)" ;
1610 0 : for (Int i=0; i<nants_p; i++) os_p << systemp[i] << " ";
1611 0 : os_p << LogIO::POST;
1612 : }
1613 : }
1614 : }
1615 :
1616 : // Go after a new pointing (where {source,ra,dec} was changed)
1617 : // SOURCE and DRA/DDEC are mixed together they define a row in the FIELD table
1618 : // so we need to build a field index here as well
1619 :
1620 210 : uvprobvr_c(uv_handle_p,"source",vtype,&vlen,&vupd);
1621 210 : if (vupd) {
1622 1 : uvgetvr_c(uv_handle_p,H_BYTE,"source",vdata,16);
1623 1 : object_p = vdata;
1624 :
1625 1 : j=-1;
1626 1 : for (uInt i=0; i<source_p.nelements(); i++) { // find first matching source name
1627 0 : if (source_p[i] == object_p) {
1628 0 : j = i ;
1629 0 : os_p << "Found old source: " << object_p << LogIO::POST;
1630 :
1631 0 : break;
1632 : }
1633 : }
1634 1 : if (j==-1) {
1635 1 : os_p << "Found new source: " << object_p << LogIO::POST;
1636 1 : source_p.resize(source_p.nelements()+1, true); // need to copy the old values
1637 1 : source_p[source_p.nelements()-1] = object_p;
1638 :
1639 1 : ras_p.resize(ras_p.nelements()+1, true);
1640 1 : decs_p.resize(decs_p.nelements()+1, true);
1641 1 : ras_p[ras_p.nelements()-1] = 0.0; // if no source at (0,0) offset
1642 1 : decs_p[decs_p.nelements()-1] = 0.0; // these would never be initialized
1643 : }
1644 1 : uvprobvr_c(uv_handle_p,"purpose",vtype,&vlen,&vupd3);
1645 1 : purpose_p.resize(purpose_p.nelements()+1, true); // need to copy the old values
1646 1 : purpose_p[purpose_p.nelements()-1] = " ";
1647 1 : if (vupd3) {
1648 0 : uvgetvr_c(uv_handle_p,H_BYTE,"purpose",vdata,16);
1649 0 : purpose_p[purpose_p.nelements()-1] = vdata;
1650 : }
1651 : }
1652 :
1653 210 : uvprobvr_c(uv_handle_p,"dra", vtype,&vlen,&vupd1);
1654 210 : uvprobvr_c(uv_handle_p,"ddec",vtype,&vlen,&vupd2);
1655 :
1656 210 : uvgetvr_c(uv_handle_p,H_DBLE,"ra", (char *)&ra_p, 1);
1657 210 : uvgetvr_c(uv_handle_p,H_DBLE,"dec",(char *)&dec_p,1);
1658 :
1659 210 : if (vupd || vupd1 || vupd2) { // if either source or offset changed, find FIELD_ID
1660 1 : npoint++;
1661 1 : if (vupd1) uvgetvr_c(uv_handle_p,H_REAL,"dra", (char *)&dra_p, 1);
1662 1 : if (vupd2) uvgetvr_c(uv_handle_p,H_REAL,"ddec",(char *)&ddec_p, 1);
1663 :
1664 1 : j=-1;
1665 1 : for (uInt i=0; i<source_p.nelements(); i++) { // find first matching source name
1666 1 : if (source_p[i] == object_p) {
1667 1 : j = i ;
1668 1 : break;
1669 : }
1670 : }
1671 : // j should always be >= 0 now, and is the source index
1672 1 : k=-1;
1673 1 : for (Int i=0; i<nfield; i++) { // check if we had this pointing/source before
1674 0 : if (dra[i] == dra_p && ddec[i] == ddec_p && field[i] == j) {
1675 0 : k = i;
1676 0 : break;
1677 : }
1678 : }
1679 : // k could be -1, when a new field/source is found
1680 : // else it is >=0 and the index into the field array
1681 :
1682 1 : if (Debug(1)) {
1683 0 : os_p << LogIO::DEBUG1 << "POINTING: " << npoint
1684 0 : << " source: " << object_p << " [" << j << "," << k << "] "
1685 0 : << " dra/ddec: " << dra_p << " " << ddec_p << LogIO::POST;
1686 : }
1687 :
1688 1 : if (k<0) { // we have a new source/field combination
1689 1 : ifield = nfield;
1690 1 : nfield++;
1691 1 : if (Debug(2)) os_p << LogIO::DEBUG2 << "Adding new field " << ifield
1692 0 : << " for " << object_p << " " << source_p[j]
1693 : << " at "
1694 0 : << dra_p *206264.8062 << " "
1695 0 : << ddec_p*206264.8062 << " arcsec." << LogIO::POST;
1696 1 : if (Debug(1)) show();
1697 :
1698 1 : if (nfield >= MAXFIELD) {
1699 0 : os_p << "Cannot handle more than " << MAXFIELD << " fields." << LogIO::POST;
1700 0 : exit(1);
1701 : }
1702 1 : ra[ifield] = ra_p;
1703 1 : dec[ifield] = dec_p;
1704 1 : dra[ifield] = dra_p;
1705 1 : ddec[ifield] = ddec_p;
1706 1 : field[ifield] = j;
1707 1 : if (dra_p == 0.0 && ddec_p==0.0) { // store ra/dec for SOURCE table as well
1708 1 : ras_p[j] = ra_p;
1709 1 : decs_p[j] = dec_p;
1710 : }
1711 : } else {
1712 0 : ifield = k;
1713 : }
1714 :
1715 1 : if (Debug(3)) os_p << LogIO::DEBUG2 << "Warning: pointing " << j
1716 : << " (dra/ddec) changed at record " << record << " : "
1717 0 : << dra_p *206264.8062 << " "
1718 0 : << ddec_p*206264.8062 << LogIO::POST;
1719 : }
1720 : } // Tracking()
1721 :
1722 :
1723 :
1724 : // ==============================================================================================
1725 : // Check for changes of the frequency setup and store them all
1726 211 : void Importmiriad::check_window()
1727 : {
1728 : int vlen,vupd;
1729 : char vtype[10];
1730 211 : Int next = nFreqSet_p;
1731 211 : if (Debug(2)) os_p << LogIO::DEBUG2 << "Importmiriad::check_window" << LogIO::POST;
1732 211 : Int idx, nchan=0, nspect=0, nwide=0;
1733 :
1734 211 : uvrdvr_c(uv_handle_p,H_INT,"nchan",(char *)&nchan, (char *)&nchan, 1);
1735 211 : uvrdvr_c(uv_handle_p,H_INT,"nspect",(char *)&nspect,(char *)&nspect, 1);
1736 211 : win[next].nspect = nspect;
1737 211 : uvrdvr_c(uv_handle_p,H_INT,"nwide",(char *)&nwide, (char *)&nwide, 1);
1738 211 : win[next].nwide = nwide;
1739 :
1740 211 : if (nspect > 0 && nspect <= MAXWIN) {
1741 :
1742 211 : uvgetvr_c(uv_handle_p,H_INT,"ischan",(char *)win[next].ischan, nspect);
1743 211 : uvgetvr_c(uv_handle_p,H_INT,"nschan",(char *)win[next].nschan, nspect);
1744 211 : uvgetvr_c(uv_handle_p,H_DBLE,"restfreq",(char *)win[next].restfreq, nspect);
1745 211 : uvgetvr_c(uv_handle_p,H_DBLE,"sdf",(char *)win[next].sdf, nspect);
1746 211 : uvgetvr_c(uv_handle_p,H_DBLE,"sfreq",(char *)win[next].sfreq, nspect);
1747 211 : uvprobvr_c(uv_handle_p,"ifchain", vtype,&vlen,&vupd);
1748 211 : if (vtype[0]=='i')
1749 211 : uvgetvr_c(uv_handle_p,H_INT,"ifchain",(char *)win[next].chain, nspect);
1750 : }
1751 211 : if (nwide>0) {
1752 0 : uvgetvr_c(uv_handle_p,H_REAL,"wfreq",(char *)win[next].wfreq, nwide);
1753 0 : uvgetvr_c(uv_handle_p,H_REAL,"wwidth",(char *)win[next].wwidth, nwide);
1754 : }
1755 :
1756 211 : if (nwide >0 && nspect != nwide) {
1757 : // don't know how to handle this
1758 : // we could assume the smaller one is the one we should deal with
1759 : // but there's no way to check how select=win() was used...
1760 : // also, if you've used uvcat options=nowide, nwide=0 and nspect non-zero.
1761 : // we really don't care about the wide bands anymore
1762 0 : if (nwide < nspect)
1763 0 : throw(AipsError("nspect != nwide"));
1764 : else {
1765 0 : nwide = nspect;
1766 : }
1767 : }
1768 :
1769 422 : for (Int i=0; i<nspect; i++) {
1770 211 : win[next].code[i] = 'N';
1771 : }
1772 :
1773 211 : idx = (nspect > 0 ? nspect : 0); // idx points into the combined win.xxx[] elements
1774 211 : for (Int i=0; i<nwide; i++) {
1775 0 : Int side = (win[next].sdf[i] < 0 ? -1 : 1);
1776 0 : win[next].code[idx] = 'S';
1777 0 : win[next].ischan[idx] = nchan + i + 1;
1778 0 : win[next].nschan[idx] = 1;
1779 0 : win[next].sfreq[idx] = win[next].wfreq[i];
1780 0 : win[next].sdf[idx] = side * win[next].wwidth[i];
1781 0 : win[next].restfreq[idx] = -1.0; // no meaning
1782 0 : idx++;
1783 : }
1784 :
1785 211 : if (nspect>0) {
1786 : // Got all the window details, now check if we already have this one
1787 211 : Bool found=false;
1788 211 : for (Int i=0; i<nFreqSet_p; i++) {
1789 : // compare win[i] and win[next]
1790 210 : found=compareWindows(win[i],win[next]);
1791 210 : if (found) {
1792 210 : freqSet_p=i;
1793 210 : break;
1794 : }
1795 : }
1796 :
1797 211 : if (not found && Debug(1)) {
1798 0 : os_p << LogIO::DEBUG1 << "Layout of spectral windows (check_window): nspect=" << nspect
1799 0 : << " nwide=" << nwide << LogIO::POST;
1800 0 : os_p << LogIO::DEBUG1 << "(N=narrow W=wide, S=spectral window averages)" << LogIO::POST;
1801 :
1802 0 : for (Int i=0; i<nspect+nwide; i++)
1803 0 : os_p << LogIO::DEBUG1 << win[next].code[i] << ": " << i+1 << " " << keep_p(i) << " "
1804 : << win[next].nschan[i] << " " << win[next].ischan[i] << " "
1805 : << win[next].sfreq[i] << " " << win[next].sdf[i] << " " << win[next].restfreq[i]
1806 0 : << LogIO::POST;
1807 : }
1808 :
1809 :
1810 211 : if (not found) {
1811 1 : os_p << "New frequency setting with "<<nspect<<" spectral windows"<<LogIO::POST;
1812 1 : nFreqSet_p=nFreqSet_p+1;
1813 1 : if (nFreqSet_p>=MAXFSET) throw(AipsError("Too many frequency settings"));
1814 1 : freqSet_p=next;
1815 : }
1816 : // Calculate datadesc_id offset
1817 211 : ddid_p=0;
1818 211 : for (Int i=0; i<freqSet_p; i++){
1819 0 : for (Int j=0; j<win[i].nspect+win[i].nwide; j++) {
1820 0 : if (keep_p(j)) ddid_p+=1;
1821 : }
1822 : }
1823 : }
1824 211 : }
1825 :
1826 : // ==============================================================================================
1827 210 : Bool Importmiriad::compareWindows(WINDOW& win1,WINDOW& win2)
1828 : {
1829 : // Check if two freq/corr windows are the same (within tolerance)
1830 210 : if (win1.nspect!= win2.nspect || win1.nwide!=win2.nwide) return false;
1831 420 : for (Int i=0; i<win1.nspect; i++){
1832 210 : if (win1.nschan[i]!=win2.nschan[i]) return false;
1833 210 : Double w = abs(win1.sdf[i]);
1834 210 : if (abs(win1.sdf[i]-win2.sdf[i])>0.01*w) return false;
1835 210 : if (abs(win1.sfreq[i]-win2.sfreq[i])>0.5*w) return false;
1836 210 : if (abs(win1.restfreq[i]-win2.restfreq[i])>0.5*w) return false;
1837 210 : if (win1.chain[i]!=win2.chain[i]) return false;
1838 : }
1839 : // could check wides, but since they are not written..
1840 210 : return true;
1841 : }
1842 :
1843 : // ==============================================================================================
1844 1 : void Importmiriad::show()
1845 : {
1846 : #if 0
1847 : os_p << "Importmiriad::show()" << LogIO::POST;
1848 : for (int i=0; i<source_p.nelements(); i++) os_p << "SOURCE_P_1: " << source_p[i] << LogIO::POST;
1849 : #endif
1850 1 : }
1851 :
1852 : // ==============================================================================================
1853 1 : void Importmiriad::close()
1854 : {
1855 : // close the file
1856 1 : uvclose_c(uv_handle_p);
1857 1 : }
1858 :
|
