Changeset 1453 for branches/alma/external/atnf/PKSIO/MBFITSreader.cc

Timestamp:

11/27/08 12:47:15 (16 years ago)

Author:

TakTsutsumi

Message:

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree

File:

• branches/alma/external/atnf/PKSIO/MBFITSreader.cc (modified) (34 diffs)

branches/alma/external/atnf/PKSIO/MBFITSreader.cc

@@ -2 +2 @@
 //# MBFITSreader.cc: ATNF single-dish RPFITS reader.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2006
 //# Mark Calabretta, ATNF
 //#
@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: MBFITSreader.cc,v 19.34 2007/07/02 06:12:18 cal103 Exp $
+//# $Id$
 //#---------------------------------------------------------------------------
 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes
@@ -81 +81 @@
   cRefChan   = 0x0;
 
-  cVis = 0x0;
-  cWgt = 0x0;
+  cVis = new float[2*4*8163];
 
   cBeamSel   = 0x0;
@@ -92 +91 @@
 
   cMBopen = 0;
+  jstat = -3;
 }
 
@@ -127 +127 @@
 
   // Open the RPFITS file.
-  int jstat = -3;
-  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-            &cBin, &cIFno, &cSrcNo);
+  rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
+            &if_no, &sourceno);
 
   if (jstat) {
-    fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname);
+    fprintf(stderr, "Failed to open MBFITS file: %s\n", rpname);
     return 1;
   }
@@ -148 +147 @@
   // Read the first header.
   jstat = -1;
-  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-            &cBin, &cIFno, &cSrcNo);
+  rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
+            &if_no, &sourceno);
 
   if (jstat) {
-    fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname);
+    fprintf(stderr, "Failed to read MBFITS header: %s\n", rpname);
     close();
     return 1;
@@ -160 +159 @@
   cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0;
 
-  // Non-ATNF data may not store the position in (u,v,w).
-  if (strncmp(names_.sta, "tid", 3) == 0) {
-    fprintf(stderr, "WARNING, found Tidbinbilla data");
-    cSUpos = 1;
-  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
-    fprintf(stderr, "WARNING, found Hobart data");
-    cSUpos = 1;
-  } else if (strncmp(names_.sta, "CED", 3) == 0) {
-    fprintf(stderr, "WARNING, found Ceduna data");
-    cSUpos = 1;
-  } else {
-    cSUpos = 0;
-  }
-
-  if (cSUpos) {
-    fprintf(stderr, ", using telescope position from SU table.\n");
+  // Tidbinbilla data has some more.
+  cTid = strncmp(names_.sta, "tid", 3) == 0;
+  if (cTid) {
+    // Telescope position is stored in the source table.
     cInterp = 0;
   }
@@ -189 +176 @@
 
   if (cNBeam <= 0) {
-    fprintf(stderr, "ERROR, couldn't determine number of beams.\n");
+    fprintf(stderr, "Couldn't determine number of beams.\n");
     close();
     return 1;
@@ -250 +237 @@
   }
 
-  // Allocate memory for RPFITSIN subroutine arguments.
-  if (cVis) delete [] cVis;
-  if (cWgt) delete [] cWgt;
-  cVis = new float[2*maxProd];
-  cWgt = new float[maxProd];
+  // Is the vis array declared by RPFITS.h large enough?
+  if (8*8193 < maxProd) {
+    // Need to allocate more memory for RPFITSIN.
+    cVis = new float[2*maxProd];
+  }
 
   nChan    = cNChan;
@@ -291 +278 @@
   // Read the first syscal record.
   if (rpget(1, cEOS)) {
-    fprintf(stderr, "ERROR, failed to read first syscal record.\n");
+    fprintf(stderr, "Error: Failed to read first syscal record.\n");
     close();
     return 1;
@@ -326 +313 @@
 {
   if (!cMBopen) {
-    fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
+    fprintf(stderr, "An MBFITS file has not been opened.\n");
     return 1;
   }
@@ -399 +386 @@
   // Time at start of observation.
   sprintf(datobs, "%-10.10s", names_.datobs);
-  utc = cUTC;
+  utc = ut;
 
   // Spectral parameters.
@@ -448 +435 @@
 
   if (!cMBopen) {
-    fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
+    fprintf(stderr, "An MBFITS file has not been opened.\n");
     return 1;
   }
@@ -576 +563 @@
 
           if (cNBin > 1 && cNBeamSel > 1) {
-            fprintf(stderr, "ERROR, cannot handle binning mode for multiple "
+            fprintf(stderr, "Cannot handle binning mode for multiple "
                             "beams.\n");
             close();
@@ -582 +569 @@
           }
 
-          // Allocate buffer data storage; the PKSMBrecord constructor zeroes
-          // class members such as cycleNo that are tested in the first pass
-          // below.
+          // Allocate buffer data storage.
           int nBuff = cNBeamSel * cNBin;
           cBuffer = new PKSMBrecord[nBuff];
@@ -599 +584 @@
           cScanNo  = 1;
           cCycleNo = 0;
-          cPrevUTC = 0.0;
+          cUTC = 0.0;
           cStaleness = new int[cNBeamSel];
           for (int iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
@@ -610 +595 @@
           cScanNo++;
           cCycleNo = 0;
-          cPrevUTC = 0.0;
+          cUTC = 0.0;
         }
 
         // Apply beam selection.
-        beamNo = int(cBaseline / 256.0);
+        beamNo = int(baseline / 256.0);
         iBeamSel = cBeamSel[beamNo-1];
         if (iBeamSel < 0) continue;
 
         // Sanity check (mainly for MOPS).
-        if (cIFno > cNIF) continue;
+        if (if_no > cNIF) continue;
 
         // Apply IF selection.
-        iIFSel = cIFSel[cIFno - 1];
+        iIFSel = cIFSel[if_no - 1];
         if (iIFSel < 0) continue;
 
         sprintf(cDateObs, "%-10.10s", names_.datobs);
 
-        // Check for change-of-day.
-        if (cUTC < cPrevUTC - 85800.0) {
-          cUTC += 86400.0;
+        // Change-of-day; note that the ut variable from RPFITS.h is global
+        // and will be preserved between calls to this function.
+        if (ut < cUTC - 85800.0) {
+          ut += 86400.0;
+        }
+
+        // New integration cycle?
+        if (ut > cUTC) {
+          cCycleNo++;
+          cUTC = ut + 0.0001;
         }
 
         if (cNBin > 1) {
           // Binning mode: correct the time.
-          cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
-        }
-
-        // New integration cycle?
-        if (cUTC > cPrevUTC) {
-          cCycleNo++;
-          cPrevUTC = cUTC + 0.0001;
+          ut += param_.intbase * (bin - (cNBin + 1)/2.0);
         }
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
 
         if (cCycleNo < iMBuff->cycleNo) {
@@ -654 +640 @@
 
         // Begin flush cycle?
-        if (cEOS || (iMBuff->nIF && cUTC > iMBuff->utc + 0.0001)) {
+        if (cEOS || (iMBuff->nIF && ut > iMBuff->utc + 0.0001)) {
           cFlushing = 1;
           cFlushBin = 0;
@@ -661 +647 @@
 
 #ifdef PKSIO_DEBUG
-        printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
+        printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, if_no);
         if (cEOS) printf("Start of new scan, flushing previous scan.\n");
 #endif
@@ -710 +696 @@
 
         // The last record read must have been the first of a new cycle.
-        beamNo = int(cBaseline / 256.0);
+        beamNo = int(baseline / 256.0);
         iBeamSel = cBeamSel[beamNo-1];
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
       }
     }
@@ -735 +721 @@
         // returned as the 'u' and 'v' visibility coordinates, must be
         // interpolated to the integration time which RPFITSIN returns as
-        // 'cUTC', this usually being a second or two later.
+        // 'ut', this usually being a second or two later.
         //
         // Note that the time recorded as the 'w' visibility coordinate
-        // cycles through 86400 back to 0 at midnight, whereas that in 'cUTC'
+        // cycles through 86400 back to 0 at midnight, whereas that in 'ut'
         // continues to increase past 86400.
 
-        double thisRA  = cU;
-        double thisDec = cV;
-        double thisUTC = cW;
+        double thisRA  = u;
+        double thisDec = v;
+        double thisUTC = w;
 
         if (thisUTC < prevUTC) {
@@ -865 +851 @@
         // Signal that all IFs in this buffer location have been flushed.
         iMBuff->nIF = 0;
-
-        // Stop cEOS being set when the next integration is read.
-        iMBuff->cycleNo = 0;
-
       } else {
         // Carry on flushing the other IFs.
@@ -877 +859 @@
       if (cFlushBin == cNBin - 1) {
         if (cEOS || cEOF) {
+          // Stop cEOS being set when the next integration is read.
+          iMBuff->cycleNo = 0;
+
           // Carry on flushing other buffers.
           cFlushIF = 0;
@@ -884 +869 @@
         cFlushing = 0;
 
-        beamNo = int(cBaseline / 256.0);
+        beamNo = int(baseline / 256.0);
         iBeamSel = cBeamSel[beamNo-1];
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
       }
     }
@@ -903 +888 @@
       }
 
-      // Sanity check on incomplete integrations within a scan.
-      if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {
-        // Force the incomplete integration to be flushed before proceeding.
-        cFlushing = 1;
-        continue;
-      }
-
       iMBuff->scanNo  = cScanNo;
       iMBuff->cycleNo = cCycleNo;
@@ -915 +893 @@
       // Times.
       strncpy(iMBuff->datobs, cDateObs, 10);
-      iMBuff->utc = cUTC;
+      iMBuff->utc = ut;
       iMBuff->exposure = param_.intbase;
 
       // Source identification.
       sprintf(iMBuff->srcName, "%-16.16s",
-              names_.su_name + (cSrcNo-1)*16);
-      iMBuff->srcRA  = doubles_.su_ra[cSrcNo-1];
-      iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
+              names_.su_name + (sourceno-1)*16);
+      iMBuff->srcRA  = doubles_.su_ra[sourceno-1];
+      iMBuff->srcDec = doubles_.su_dec[sourceno-1];
 
       // Rest frequency of the line of interest.
@@ -948 +926 @@
 
       // Beam position at the specified time.
-      if (cSUpos) {
-        // Non-ATNF data that does not store the position in (u,v,w).
-        iMBuff->ra  = doubles_.su_ra[cSrcNo-1];
-        iMBuff->dec = doubles_.su_dec[cSrcNo-1];
+      if (cTid) {
+        // Tidbinbilla data.
+        iMBuff->ra  = doubles_.su_ra[sourceno-1];
+        iMBuff->dec = doubles_.su_dec[sourceno-1];
       } else {
-        iMBuff->ra  = cU;
-        iMBuff->dec = cV;
-      }
-      cPosUTC[iBeamSel] = cW;
+        iMBuff->ra  = u;
+        iMBuff->dec = v;
+      }
+      cPosUTC[iBeamSel] = w;
 
       // IF-dependent parameters.
-      int iIF = cIFno - 1;
+      int iIF = if_no - 1;
       int startChan = cStartChan[iIF];
       int endChan   = cEndChan[iIF];
@@ -968 +946 @@
       iIFSel = cIFSel[iIF];
       iMBuff->nIF++;
-      iMBuff->IFno[iIFSel]  = cIFno;
+      iMBuff->IFno[iIFSel]  = if_no;
       iMBuff->nChan[iIFSel] = nChan;
       iMBuff->nPol[iIFSel]  = cNPol[iIF];
@@ -1027 +1005 @@
 
       // The baseline flag may be set independently.
-      if (rpflag == 0) rpflag = cFlag;
+      if (rpflag == 0) rpflag = flag;
 
       // Copy and scale data.
@@ -1157 +1135 @@
   int numErr = 0;
 
-  int jstat = 0;
+  jstat = 0;
   while (numErr < 10) {
     int lastjstat = jstat;
-    rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-              &cBin, &cIFno, &cSrcNo);
+    rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
+              &if_no, &sourceno);
 
     switch(jstat) {
@@ -1174 +1152 @@
       // Successful read.
       if (lastjstat == 0) {
-        if (cBaseline == -1) {
+        if (baseline == -1) {
           // Syscal data.
           if (syscalonly) {
@@ -1230 +1208 @@
   }
 
-  fprintf(stderr, "ERROR, RPFITS read failed too many times.\n");
+  fprintf(stderr, "RPFITS read failed too many times.\n");
   return 2;
 }
@@ -1241 +1219 @@
 {
   if (cMBopen) {
-    int jstat = 1;
-    rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-              &cBin, &cIFno, &cSrcNo);
+    jstat = 1;
+    rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
+              &if_no, &sourceno);
 
     if (cBeams)     delete [] cBeams;
@@ -1254 +1232 @@
     if (cRefChan)   delete [] cRefChan;
 
-    if (cVis) delete [] cVis;
-    if (cWgt) delete [] cWgt;
+    if (cVis)       delete [] cVis;
 
     if (cBeamSel)   delete [] cBeamSel;