Changeset 1452 for trunk/external


Ignore:
Timestamp:
11/19/08 20:41:16 (16 years ago)
Author:
Malte Marquarding
Message:

update from livedata CVS

Location:
trunk/external/atnf
Files:
5 added
2 deleted
22 edited

Legend:

Unmodified
Added
Removed
  • trunk/external/atnf/PKSIO/FITSreader.cc

    r1325 r1452  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: FITSreader.cc,v 19.2 2006/05/19 02:14:50 mcalabre Exp $
     29//# $Id: FITSreader.cc,v 19.3 2008-09-12 02:37:37 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The FITSreader class is an abstract base class for the Parkes Multibeam
     
    5353        const int getSpectra,
    5454        const int getXPol,
    55         const int getFeedPos)
     55        const int coordSys)
    5656{
    5757  int maxNChan = 0;
     
    8383  cGetSpectra = getSpectra && cHaveSpectra;
    8484  cGetXPol    = getXPol    && cGetXPol;
    85   cGetFeedPos = getFeedPos;
     85  cCoordSys   = coordSys;
    8686
    8787  return maxNChan;
  • trunk/external/atnf/PKSIO/FITSreader.h

    r1399 r1452  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: FITSreader.h,v 19.6 2007/11/12 03:37:56 cal103 Exp $
     29//# $Id: FITSreader.h,v 19.9 2008-11-17 06:28:04 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The FITSreader class is an abstract base class for the Parkes Multibeam
     
    3838#define ATNF_FITSREADER_H
    3939
    40 #include <atnf/PKSIO/PKSMBrecord.h>
     40#include <atnf/PKSIO/MBrecord.h>
     41#include <atnf/PKSIO/PKSmsg.h>
     42
     43using namespace std;
    4144
    4245// <summary>
     
    4447// </summary>
    4548
    46 class FITSreader
     49class FITSreader : public PKSmsg
    4750{
    4851  public:
     
    9194    // Set data selection criteria.  Channel numbering is 1-relative, zero or
    9295    // negative channel numbers are taken to be offsets from the last channel.
     96    // Coordinate systems are
     97    //   0: equatorial (RA,Dec),
     98    //   1: vertical (Az,El),
     99    //   2: feed-plane.
    93100    int select(
    94101        const int startChan[],
     
    96103        const int refChan[],
    97104        const int getSpectra = 1,
    98         const int getXPol = 0,
    99         const int getFeedPos = 0);
     105        const int getXPol  = 0,
     106        const int coordSys = 0);
    100107
    101108    // Find the range in time and position of the data selected.
     
    109116    // Read the next data record.
    110117    virtual int read(
    111         PKSMBrecord &record) = 0;
     118        MBrecord &record) = 0;
    112119
    113120    // Close the RPFITS file.
     
    115122
    116123  protected:
    117     int    *cBeams, *cEndChan, cGetFeedPos, cGetSpectra, cGetXPol, cHaveBase,
    118            cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF, *cNPol,
    119            *cRefChan, *cStartChan;
     124    int  *cBeams, *cEndChan, cCoordSys, cGetSpectra, cGetXPol, cHaveBase,
     125         cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF, *cNPol,
     126         *cRefChan, *cStartChan;
     127
     128    // For use in constructing messages.
     129    char cMsg[256];
    120130};
    121131
  • trunk/external/atnf/PKSIO/MBFITSreader.cc

    r1427 r1452  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBFITSreader.cc,v 19.38 2008-06-26 02:24:22 cal103 Exp $
     29//# $Id: MBFITSreader.cc,v 19.54 2008-11-17 06:51:55 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
     
    3535//#---------------------------------------------------------------------------
    3636
     37#include <atnf/pks/pks_maths.h>
    3738#include <atnf/PKSIO/MBFITSreader.h>
    38 #include <atnf/PKSIO/PKSMBrecord.h>
    39 
    40 #include <RPFITS.h>
     39#include <atnf/PKSIO/MBrecord.h>
    4140
    4241#include <casa/math.h>
     
    4746#include <unistd.h>
    4847
     48#include <RPFITS.h>
     49
    4950using namespace std;
    5051
     
    5253const double PI = 3.141592653589793238462643;
    5354const double TWOPI = 2.0 * PI;
     55const double HALFPI = PI / 2.0;
    5456const double R2D = 180.0 / PI;
    5557
     
    9395
    9496  cMBopen = 0;
     97
     98  // Tell RPFITSIN not to report errors directly.
     99  iostat_.errlun = -1;
     100
     101  // By default, messages are written to stderr.
     102  initMsg();
    95103}
    96104
     
    121129        int  &extraSysCal)
    122130{
     131  // Clear the message stack.
     132  clearMsg();
     133
    123134  if (cMBopen) {
    124135    close();
     
    129140  // Open the RPFITS file.
    130141  int jstat = -3;
    131   rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
    132             &cBin, &cIFno, &cSrcNo);
    133 
    134   if (jstat) {
    135     fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname);
     142  if (rpfitsin(jstat)) {
     143    sprintf(cMsg, "ERROR: failed to open MBFITS file\n       %s", rpname);
     144    logMsg(cMsg);
    136145    return 1;
    137146  }
     
    149158  // Read the first header.
    150159  jstat = -1;
    151   rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
    152             &cBin, &cIFno, &cSrcNo);
    153 
    154   if (jstat) {
    155     fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname);
     160  if (rpfitsin(jstat)) {
     161    sprintf(cMsg, "ERROR: failed to read MBFITS header in file\n"
     162                  "       %s", rpname);
     163    logMsg(cMsg);
    156164    close();
    157165    return 1;
     
    163171  // Non-ATNF data may not store the position in (u,v,w).
    164172  if (strncmp(names_.sta, "tid", 3) == 0) {
    165     fprintf(stderr, "WARNING, found Tidbinbilla data");
     173    sprintf(cMsg, "WARNING: Found Tidbinbilla data");
    166174    cSUpos = 1;
    167175  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
    168     fprintf(stderr, "WARNING, found Hobart data");
     176    sprintf(cMsg, "WARNING: Found Hobart data");
    169177    cSUpos = 1;
    170178  } else if (strncmp(names_.sta, "CED", 3) == 0) {
    171     fprintf(stderr, "WARNING, found Ceduna data");
     179    sprintf(cMsg, "WARNING: Found Ceduna data");
    172180    cSUpos = 1;
    173181  } else {
     
    176184
    177185  if (cSUpos) {
    178     fprintf(stderr, ", using telescope position from SU table.\n");
     186    strcat(cMsg, ", using telescope position\n         from SU table.");
     187    logMsg(cMsg);
    179188    cInterp = 0;
    180189  }
     190
     191  // Mean scan rate (for timestamp repairs).
     192  cNRate = 0;
     193  cAvRate[0] = 0.0;
     194  cAvRate[1] = 0.0;
     195  cCode5 = 0;
    181196
    182197
     
    190205
    191206  if (cNBeam <= 0) {
    192     fprintf(stderr, "ERROR, couldn't determine number of beams.\n");
     207    logMsg("ERROR, couldn't determine number of beams.");
    193208    close();
    194209    return 1;
     
    203218  // ...beams present in the data.
    204219  for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
    205     cBeams[anten_.ant_num[iBeam] - 1] = 1;
     220    // Guard against dubious beam numbers, e.g. zeroes in
     221    // 1999-09-29_1632_024848p14_071b.hpf and the four scans following.
     222    // Note that the actual beam number is decoded from the 'baseline' random
     223    // parameter for each spectrum and is only used for beam selection.
     224    int beamNo = anten_.ant_num[iBeam];
     225    if (beamNo != iBeam+1) {
     226      char sta[8];
     227      strncpy(sta, names_.sta+(8*iBeam), 8);
     228      char *cp = sta + 7;
     229      while (*cp == ' ') *(cp--) = '\0';
     230
     231      sprintf(cMsg,
     232        "WARNING: RPFITSIN returned beam number %2d for AN table\n"
     233        "         entry %2d with name '%.8s'", beamNo, iBeam+1, sta);
     234
     235      char text[8];
     236      sprintf(text, "MB%2.2d", iBeam+1);
     237      cp = cMsg + strlen(cMsg);
     238      if (strncmp(sta, text, 8) == 0) {
     239        beamNo = iBeam + 1;
     240        sprintf(cp, "; using beam number %2d.", beamNo);
     241      } else {
     242        sprintf(cp, ".");
     243      }
     244
     245      logMsg(cMsg);
     246    }
     247
     248    if (0 < beamNo && beamNo <= cNBeam) {
     249      cBeams[beamNo-1] = 1;
     250    }
    206251  }
    207252
     
    292337  // Read the first syscal record.
    293338  if (rpget(1, cEOS)) {
    294     fprintf(stderr, "ERROR, failed to read first syscal record.\n");
     339    logMsg("ERROR, failed to read first syscal record.");
    295340    close();
    296341    return 1;
     
    328373{
    329374  if (!cMBopen) {
    330     fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
     375    logMsg("ERROR, an MBFITS file has not been opened.");
    331376    return 1;
    332377  }
     
    348393    antPos[1] =  2816759.046;
    349394    antPos[2] = -3454035.950;
     395
    350396  } else if (strncmp(names_.sta, "HOH", 3) == 0) {
    351397    // Parkes HOH receiver.
     
    354400    antPos[1] =  2816759.046;
    355401    antPos[2] = -3454035.950;
     402
    356403  } else if (strncmp(names_.sta, "CA0", 3) == 0) {
    357404    // An ATCA antenna, use the array centre position.
     
    360407    antPos[1] =  2792906.182;
    361408    antPos[2] = -3200483.747;
     409
     410    // ATCA-104.  Updated position at epoch 2007/06/24 from Chris Phillips.
     411    // antPos[0] = -4751640.182; // ± 0.008
     412    // antPos[1] =  2791700.322; // ± 0.006
     413    // antPos[2] = -3200490.668; // ± 0.007
     414    //
    362415  } else if (strncmp(names_.sta, "MOP", 3) == 0) {
    363     // Mopra.
     416    // Mopra.  Updated position at epoch 2007/06/24 from Chris Phillips.
    364417    sprintf(telescope, "%-16.16s", "ATMOPRA");
    365     antPos[0] = -4682768.630;
    366     antPos[1] =  2802619.060;
    367     antPos[2] = -3291759.900;
     418    antPos[0] = -4682769.444; // ± 0.009
     419    antPos[1] =  2802618.963; // ± 0.006
     420    antPos[2] = -3291758.864; // ± 0.008
     421
    368422  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
    369423    // Hobart.
     
    372426    antPos[1] =  2522347.567;
    373427    antPos[2] = -4311562.569;
     428
    374429  } else if (strncmp(names_.sta, "CED", 3) == 0) {
    375     // Ceduna.
     430    // Ceduna.  Updated position at epoch 2007/06/24 from Chris Phillips.
    376431    sprintf(telescope, "%-16.16s", "CEDUNA");
    377     antPos[0] = -3749943.657;
    378     antPos[1] =  3909017.709;
    379     antPos[2] = -3367518.309;
     432    antPos[0] = -3753443.168; // ± 0.017
     433    antPos[1] =  3912709.794; // ± 0.017
     434    antPos[2] = -3348067.060; // ± 0.016
     435
    380436  } else if (strncmp(names_.sta, "tid", 3) == 0) {
    381437    // DSS.
     
    448504//--------------------------------------------------------- MBFITSreader::read
    449505
    450 // Read the next data record.
     506// Read the next data record (if you're feeling lucky).
    451507
    452508int MBFITSreader::read(
    453         PKSMBrecord &MBrec)
     509        MBrecord &MBrec)
    454510{
    455511  int beamNo = -1;
    456   int haveData, status;
    457   PKSMBrecord *iMBuff = 0x0;
     512  int haveData, pCode = 0, status;
     513  double raRate = 0.0, decRate = 0.0, paRate = 0.0;
     514  MBrecord *iMBuff = 0x0;
    458515
    459516  if (!cMBopen) {
    460     fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
     517    logMsg("ERROR, an MBFITS file has not been opened.");
    461518    return 1;
    462519  }
    463520
    464   // Positions recorded in the input records do not coincide with the midpoint
    465   // of the integration and hence the input must be buffered so that true
    466   // positions may be interpolated.
     521  // Positions recorded in the input records usually do not coincide with the
     522  // midpoint of the integration and hence the input must be buffered so that
     523  // true positions may be interpolated.
    467524  //
    468525  // On the first call nBeamSel buffers of length nBin, are allocated and
     
    494551
    495552      // Read the next record.
     553      pCode = 0;
    496554      if ((status = rpget(0, cEOS)) == -1) {
    497555        // EOF.
     
    502560
    503561#ifdef PKSIO_DEBUG
    504         printf("\nEnd-of-file detected, flushing last scan.\n");
     562        fprintf(stderr, "\nEnd-of-file detected, flushing last cycle.\n");
    505563#endif
    506564
     
    586644
    587645          if (cNBin > 1 && cNBeamSel > 1) {
    588             fprintf(stderr, "ERROR, cannot handle binning mode for multiple "
    589                             "beams.\n");
     646            logMsg("ERROR, cannot handle binning mode for multiple beams.");
    590647            close();
    591648            return 1;
    592649          }
    593650
    594           // Allocate buffer data storage; the PKSMBrecord constructor zeroes
     651          // Allocate buffer data storage; the MBrecord constructor zeroes
    595652          // class members such as cycleNo that are tested in the first pass
    596653          // below.
    597654          int nBuff = cNBeamSel * cNBin;
    598           cBuffer = new PKSMBrecord[nBuff];
     655          cBuffer = new MBrecord[nBuff];
    599656
    600657          // Allocate memory for spectral arrays.
     
    602659            cBuffer[ibuff].setNIFs(simulIF);
    603660            cBuffer[ibuff].allocate(0, maxChan, maxXpol);
     661
     662            // Signal that this IF in this buffer has been flushed.
     663            for (int iIF = 0; iIF < simulIF; iIF++) {
     664              cBuffer[ibuff].IFno[iIF] = 0;
     665            }
    604666          }
    605667
     
    609671          cScanNo  = 1;
    610672          cCycleNo = 0;
    611           cPrevUTC = 0.0;
     673          cPrevUTC = -1.0;
    612674        }
    613675
     
    616678          cScanNo++;
    617679          cCycleNo = 0;
    618           cPrevUTC = 0.0;
     680          cPrevUTC = -1.0;
    619681        }
    620682
    621683        // Check for change-of-day.
    622         if (cUTC < cPrevUTC - 85800.0) {
     684        double cod = 0.0;
     685        if ((cUTC + 86400.0) < (cPrevUTC + 600.0)) {
     686          // cUTC should continue to increase past 86400 during a single scan.
     687          // However, if the RPFITS file contains multiple scans that straddle
     688          // midnight then cUTC can jump backwards from the end of one scan to
     689          // the start of the next.
     690#ifdef PKSIO_DEBUG
     691          fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f",
     692            cPrevUTC, cUTC);
     693#endif
     694          // Can't change the recorded value of cUTC directly (without also
     695          // changing dateobs) so change-of-day must be recorded separately as
     696          // an offset to be applied when comparing integration timestamps.
     697          cod = 86400.0;
     698
     699        } else if (cUTC < cPrevUTC - 1.0) {
     700          sprintf(cMsg,
     701            "WARNING: Cycle %d:%03d-%03d, UTC went backwards from\n"
     702            "         %.1f to %.1f!  Incrementing day number,\n"
     703            "         positions may be unreliable.", cScanNo, cCycleNo,
     704            cCycleNo+1, cPrevUTC, cUTC);
     705          logMsg(cMsg);
    623706          cUTC += 86400.0;
    624707        }
     
    630713
    631714        // New integration cycle?
    632         if (cUTC > cPrevUTC) {
     715        if ((cUTC+cod) > cPrevUTC) {
    633716          cCycleNo++;
    634717          cPrevUTC = cUTC + 0.0001;
     
    637720        // Apply beam selection.
    638721        beamNo = int(cBaseline / 256.0);
     722        if (beamNo == 1) {
     723          // Store the position of beam 1 for grid convergence corrections.
     724          cRA0  = cU;
     725          cDec0 = cV;
     726        }
    639727        iBeamSel = cBeamSel[beamNo-1];
    640728        if (iBeamSel < 0) continue;
     
    648736
    649737        sprintf(cDateObs, "%-10.10s", names_.datobs);
     738        cDateObs[10] = '\0';
    650739
    651740        // Compute buffer number.
     
    660749
    661750        // Begin flush cycle?
    662         if (cEOS || (iMBuff->nIF && cUTC > iMBuff->utc + 0.0001)) {
     751        if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
    663752          cFlushing = 1;
    664753          cFlushBin = 0;
     
    667756
    668757#ifdef PKSIO_DEBUG
    669         printf("\n In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
    670         if (cEOS) printf("Start of new scan, flushing previous scan.\n");
     758        char rel = '=';
     759        double dt = utcDiff(cUTC, cW);
     760        if (dt < 0.0) {
     761          rel = '<';
     762        } else if (dt > 0.0) {
     763          rel = '>';
     764        }
     765
     766        fprintf(stderr, "\n In:%4d%4d%3d%3d  %.3f %c %.3f (%+.3fs) - "
     767          "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
     768          cFlushing ? "" : "not ");
     769        if (cEOS) {
     770          fprintf(stderr, "Start of new scan, flushing previous scan.\n");
     771        }
    671772#endif
    672773      }
     
    733834      // as the 'u' and 'v' visibility coordinates, must be interpolated to
    734835      // the integration time which RPFITSIN returns as 'cUTC', this usually
    735       // being a second or two later.
     836      // being a second or two later.  The interpolation method used here is
     837      // based on the scan rate.
    736838      //
    737839      // "This" RA, Dec, and UTC refers to the position currently stored in
    738       // the buffer marked for output (iMBuff).  This position will be
    739       // interpolated to the midpoint of that integration using the position
    740       // recorded in the "next" integration which is currently sitting in the
    741       // RPFITS commons.  The interpolation method used here is based on the
    742       // scan rate.  At the end of a scan, or if the next position has not
    743       // been updated, the most recent determination of the scan rate will be
    744       // used for extrapolation.
     840      // the buffer marked for output (iMBuff).  This position is interpolated
     841      // to the midpoint of that integration using either
     842      //   a) the rate currently sitting in iMBuff, which was computed from
     843      //      the previous integration, otherwise
     844      //   b) from the position recorded in the "next" integration which is
     845      //      currently sitting in the RPFITS commons,
     846      // so that the position timestamps straddle the midpoint of the
     847      // integration and is thereby interpolated rather than extrapolated.
    745848      //
    746       // The rate "age" is the offset from "this" integration (in iMBuff) of
    747       // the earliest integration in the pair used to compute the rate.  A
    748       // rate "age" of 0 thus refers to the normal situation where the rate
    749       // is determined from "this" integration and the "next" one.  An age
    750       // of 1 cycle means that it is determined from "this" integration and
    751       // the one preceding it, which should be equally reliable.  An age
    752       // of 2 cycles means that the rate is determined from the previous
    753       // integration and the one before that, so the extrapolation spans one
    754       // integration cycle.  Thus it has a "staleness" of 1.
    755 
     849      // At the end of a scan, or if the next position has not been updated
     850      // or its timestamp does not advance sufficiently, the most recent
     851      // determination of the scan rate will be used for extrapolation which
     852      // is quantified by the "rate age" measured in seconds beyond the
     853      // interval defined by the position timestamps.
     854
     855      // At this point, iMBuff contains cU, cV, cW, parAngle and focusRot
     856      // stored from the previous call to rpget() for this beam (i.e. "this"),
     857      // and also raRate, decRate and paRate computed from that integration
     858      // and the previous one.
    756859      double thisRA  = iMBuff->ra;
    757860      double thisDec = iMBuff->dec;
    758861      double thisUTC = cPosUTC[iBeamSel];
     862      double thisPA  = iMBuff->parAngle + iMBuff->focusRot;
     863
     864#ifdef PKSIO_DEBUG
     865      fprintf(stderr, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
     866        iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
     867#endif
    759868
    760869      if (cEOF || cEOS) {
    761         iMBuff->rateAge++;
    762         iMBuff->rateson = 0;
    763        
     870        // Use rates from the last cycle.
     871        raRate  = iMBuff->raRate;
     872        decRate = iMBuff->decRate;
     873        paRate  = iMBuff->paRate;
     874
    764875      } else {
    765         // Note that the time recorded as the 'w' visibility coordinate
    766         // cycles through 86400 back to 0 at midnight, whereas that in 'cUTC'
    767         // continues to increase past 86400.
    768 
    769         double nextRA  = cU;
    770         double nextDec = cV;
    771         double nextUTC = cW;
    772 
    773         if (nextUTC < thisUTC) {
    774           // Must have cycled through midnight.
    775           nextUTC += 86400.0;
    776         }
    777 
    778         // Guard against RA cycling through 24h in either direction.
    779         if (fabs(nextRA - thisRA) > PI) {
    780           if (nextRA < thisRA) {
    781             nextRA += TWOPI;
     876        if (cW == thisUTC) {
     877          // The control system at Mopra typically does not update the
     878          // positions between successive integration cycles at the end of a
     879          // scan (nor are they flagged).  In this case we use the previously
     880          // computed rates, even if from the previous scan since these are
     881          // likely to be a better guess than anything else.
     882          raRate  = iMBuff->raRate;
     883          decRate = iMBuff->decRate;
     884          paRate  = iMBuff->paRate;
     885
     886          if (cU == thisRA && cV == thisDec) {
     887            // Position and timestamp unchanged.
     888            pCode = 1;
     889
     890          } else if (fabs(cU-thisRA) < 0.0001 && fabs(cV-thisDec) < 0.0001) {
     891            // Allow small rounding errors (seen infrequently).
     892            pCode = 1;
     893
    782894          } else {
    783             nextRA -= TWOPI;
     895            // (cU,cV) are probably rubbish (not yet seen in practice).
     896            pCode = 2;
     897            cU = thisRA;
     898            cV = thisDec;
    784899          }
    785         }
    786900
    787901#ifdef PKSIO_DEBUG
    788         printf("Previous ra, dec, UTC:  %8.4f  %8.4f  %7.1f\n", thisRA*R2D,
    789           thisDec*R2D, thisUTC);
    790         printf("Current  ra, dec, UTC:  %8.4f  %8.4f  %7.1f\n", nextRA*R2D,
    791           nextDec*R2D, nextUTC);
     902          fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
     903            "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
    792904#endif
    793905
    794         // The control system at Mopra typically does not update the
    795         // positions between successive integration cycles at the end of a
    796         // scan (nor are they flagged).  In this case we use the previously
    797         // computed rates, even if from the previous scan since these are
    798         // likely to be a better guess than anything else.
    799 
    800         double dUTC = nextUTC - thisUTC;
    801 
    802         // Scan rate for this beam.
    803         if (dUTC > 0.0) {
    804           iMBuff->raRate  = (nextRA  - thisRA)  / dUTC;
    805           iMBuff->decRate = (nextDec - thisDec) / dUTC;
    806           iMBuff->rateAge = 0;
    807           iMBuff->rateson = 0;
    808 
    809           if (cInterp == 2) {
    810             // Use the same interpolation scheme as the original pksmbfits
    811             // client.  This incorrectly assumed that (nextUTC - thisUTC) is
    812             // equal to the integration time and interpolated by computing a
    813             // weighted sum of the positions before and after the required
    814             // time.
    815 
    816             double utc = iMBuff->utc;
    817             if (utc - thisUTC > 100.0) {
    818               // Must have cycled through midnight.
    819               utc -= 86400.0;
     906        } else {
     907          double nextRA  = cU;
     908          double nextDec = cV;
     909
     910          // Check and, if necessary, repair the position timestamp,
     911          // remembering that pCode refers to the NEXT cycle.
     912          pCode = fixw(cDateObs, cCycleNo, beamNo, cAvRate, thisRA, thisDec,
     913                       thisUTC, nextRA, nextDec, cW);
     914          if (pCode > 0) pCode += 3;
     915          double nextUTC = cW;
     916
     917#ifdef PKSIO_DEBUG
     918          fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
     919            "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
     920            utcDiff(nextUTC, thisUTC));
     921#endif
     922
     923          // Compute the scan rate for this beam.
     924          double dUTC = utcDiff(nextUTC, thisUTC);
     925          if ((0.0 < dUTC) && (dUTC < 600.0)) {
     926            scanRate(cRA0, cDec0, thisRA, thisDec, nextRA, nextDec, dUTC,
     927                     raRate, decRate);
     928
     929            // Update the mean scan rate.
     930            cAvRate[0] = (cAvRate[0]*cNRate +  raRate) / (cNRate + 1);
     931            cAvRate[1] = (cAvRate[1]*cNRate + decRate) / (cNRate + 1);
     932            cNRate++;
     933
     934            // Rate of change of position angle.
     935            if (sc_.sc_ant <= anten_.nant) {
     936              paRate = 0.0;
     937            } else {
     938              int iOff = sc_.sc_q * (sc_.sc_ant - 1) - 1;
     939              double nextPA = sc_.sc_cal[iOff + 4] + sc_.sc_cal[iOff + 7];
     940              double paDiff = nextPA - thisPA;
     941              if (paDiff > PI) {
     942                paDiff -= TWOPI;
     943              } else if (paDiff < -PI) {
     944                paDiff += TWOPI;
     945              }
     946              paRate = paDiff / dUTC;
    820947            }
    821948
    822             double tw1 = 1.0 - (utc - thisUTC) / iMBuff->exposure;
    823             double tw2 = 1.0 - (nextUTC - utc) / iMBuff->exposure;
    824             double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);
    825 
    826             iMBuff->raRate  *= gamma;
    827             iMBuff->decRate *= gamma;
    828           }
    829 
    830         } else {
    831           iMBuff->rateAge++;
    832 
    833           // Staleness codes.
    834           if (dUTC < 0.0) {
    835             iMBuff->rateson = 3;
     949            if (cInterp == 2) {
     950              // Use the same interpolation scheme as the original pksmbfits
     951              // client.  This incorrectly assumed that (nextUTC - thisUTC) is
     952              // equal to the integration time and interpolated by computing a
     953              // weighted sum of the positions before and after the required
     954              // time.
     955
     956              double utc = iMBuff->utc;
     957              double tw1 = 1.0 - utcDiff(utc, thisUTC) / iMBuff->exposure;
     958              double tw2 = 1.0 - utcDiff(nextUTC, utc) / iMBuff->exposure;
     959              double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);
     960
     961              // Guard against RA cycling through 24h in either direction.
     962              if (fabs(nextRA - thisRA) > PI) {
     963                if (nextRA < thisRA) {
     964                  nextRA += TWOPI;
     965                } else {
     966                  nextRA -= TWOPI;
     967                }
     968              }
     969
     970              raRate  = gamma * (nextRA  - thisRA)  / dUTC;
     971              decRate = gamma * (nextDec - thisDec) / dUTC;
     972            }
     973
    836974          } else {
    837             if (nextRA != thisRA || nextDec != thisDec) {
    838               iMBuff->rateson = 2;
     975            if (cCycleNo == 2 && fabs(utcDiff(cUTC,cW)) < 600.0) {
     976              // thisUTC (i.e. cW for the first cycle) is rubbish, and
     977              // probably the position as well (extremely rare in practice,
     978              // e.g. 97-12-19_1029_235708-18_586e.hpf which actually has the
     979              // t/1000 scaling bug in the first cycle).
     980              iMBuff->pCode = 3;
     981              thisRA  = cU;
     982              thisDec = cV;
     983              thisUTC = cW;
     984              raRate  = 0.0;
     985              decRate = 0.0;
     986              paRate  = 0.0;
     987
    839988            } else {
    840               iMBuff->rateson = 1;
     989              // cW is rubbish and probably (cU,cV), and possibly the
     990              // parallactic angle and everything else as well (rarely seen
     991              // in practice, e.g. 97-12-09_0743_235707-58_327c.hpf and
     992              // 97-09-01_0034_123717-42_242b.hpf, the latter with bad
     993              // parallactic angle).
     994              pCode = 3;
     995              cU = thisRA;
     996              cV = thisDec;
     997              cW = thisUTC;
     998              raRate  = iMBuff->raRate;
     999              decRate = iMBuff->decRate;
     1000              paRate  = iMBuff->paRate;
    8411001            }
    8421002          }
     
    8441004      }
    8451005
     1006
     1007      // Choose the closest rate determination.
     1008      if (cCycleNo == 1) {
     1009        // Scan containing a single integration.
     1010        iMBuff->raRate  = 0.0;
     1011        iMBuff->decRate = 0.0;
     1012        iMBuff->paRate  = 0.0;
     1013
     1014      } else {
     1015        double dUTC = iMBuff->utc - cPosUTC[iBeamSel];
     1016
     1017        if (dUTC >= 0.0) {
     1018          // In HIPASS/ZOA, the position timestamp, which should always occur
     1019          // on the whole second, normally precedes an integration midpoint
     1020          // falling on the half-second.  Consequently, positive ages are
     1021          // always half-integral.
     1022          dUTC = utcDiff(iMBuff->utc, cW);
     1023          if (dUTC > 0.0) {
     1024            iMBuff->rateAge = dUTC;
     1025          } else {
     1026            iMBuff->rateAge = 0.0f;
     1027          }
     1028
     1029          iMBuff->raRate  =  raRate;
     1030          iMBuff->decRate = decRate;
     1031          iMBuff->paRate  =  paRate;
     1032
     1033        } else {
     1034          // In HIPASS/ZOA, negative ages occur when the integration midpoint,
     1035          // occurring on the whole second, precedes the position timestamp.
     1036          // Thus negative ages are always an integral number of seconds.
     1037          // They have only been seen to occur sporadically in the period
     1038          // 1999/05/31 to 1999/11/01, e.g. 1999-07-26_1821_005410-74_007c.hpf
     1039          //
     1040          // In recent (2008/10/07) Mopra data, small negative ages (~10ms,
     1041          // occasionally up to ~300ms) seem to be the norm, with both the
     1042          // position timestamp and integration midpoint falling close to but
     1043          // not on the integral second.
     1044          if (cCycleNo == 2) {
     1045            // We have to start with something!
     1046            iMBuff->rateAge = dUTC;
     1047
     1048          } else {
     1049            // Although we did not record the relevant position timestamp
     1050            // explicitly, it can easily be deduced.
     1051            double w = iMBuff->utc - utcDiff(cUTC, iMBuff->utc) -
     1052                       iMBuff->rateAge;
     1053            dUTC = utcDiff(iMBuff->utc, w);
     1054
     1055            if (dUTC > 0.0) {
     1056              iMBuff->rateAge = 0.0f;
     1057            } else {
     1058              iMBuff->rateAge = dUTC;
     1059            }
     1060          }
     1061
     1062          iMBuff->raRate  =  raRate;
     1063          iMBuff->decRate = decRate;
     1064          iMBuff->paRate  =  paRate;
     1065        }
     1066      }
     1067
    8461068#ifdef PKSIO_DEBUG
    847       printf("Doing position interpolation for beam %d.\n", iMBuff->beamNo);
    848       printf("RA and Dec rates and age: %7.4f  %7.4f  %d\n",
    849         iMBuff->raRate*R2D, iMBuff->decRate*R2D, iMBuff->rateAge);
     1069      double avRate = sqrt(cAvRate[0]*cAvRate[0] + cAvRate[1]*cAvRate[1]);
     1070      fprintf(stderr, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
     1071        "pCode %d\n", raRate*R2D, decRate*R2D, avRate*R2D, paRate*R2D, pCode);
    8501072#endif
     1073
    8511074
    8521075      // Compute the position of this beam for all bins.
     
    8561079        cBuffer[jbuff].raRate  = iMBuff->raRate;
    8571080        cBuffer[jbuff].decRate = iMBuff->decRate;
    858 
    859         double dutc = cBuffer[jbuff].utc - thisUTC;
    860         if (dutc > 100.0) {
     1081        cBuffer[jbuff].paRate  = iMBuff->paRate;
     1082
     1083        double dUTC = utcDiff(cBuffer[jbuff].utc, thisUTC);
     1084        if (dUTC > 100.0) {
    8611085          // Must have cycled through midnight.
    862           dutc -= 86400.0;
    863         }
    864 
    865         cBuffer[jbuff].ra  = thisRA  + cBuffer[jbuff].raRate  * dutc;
    866         cBuffer[jbuff].dec = thisDec + cBuffer[jbuff].decRate * dutc;
    867         if (cBuffer[jbuff].ra < 0.0) {
    868           cBuffer[jbuff].ra += TWOPI;
    869         } else if (cBuffer[jbuff].ra > TWOPI) {
    870           cBuffer[jbuff].ra -= TWOPI;
    871         }
     1086          dUTC -= 86400.0;
     1087        }
     1088
     1089        applyRate(cRA0, cDec0, thisRA, thisDec,
     1090          cBuffer[jbuff].raRate, cBuffer[jbuff].decRate, dUTC,
     1091          cBuffer[jbuff].ra, cBuffer[jbuff].dec);
     1092
     1093#ifdef PKSIO_DEBUG
     1094        fprintf(stderr, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
     1095          "age: %d %.1fs)\n", iMBuff->cycleNo, cBuffer[jbuff].ra*R2D,
     1096          cBuffer[jbuff].dec*R2D, cBuffer[jbuff].utc, iMBuff->pCode,
     1097          iMBuff->rateAge);
     1098#endif
    8721099      }
    8731100    }
     
    8801107
    8811108#ifdef PKSIO_DEBUG
    882       printf("Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo, MBrec.beamNo,
    883         MBrec.IFno[0]);
     1109      fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
     1110        MBrec.beamNo, MBrec.IFno[0]);
    8841111#endif
    8851112
     
    9231150        // Sanity check on the number of IFs in the new scan.
    9241151        if (if_.n_if != cNIF) {
    925           fprintf(stderr, "WARNING, scan %d has %d IFs instead of %d, "
    926             "continuing.\n", cScanNo, if_.n_if, cNIF);
     1152          sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
     1153            "continuing.", cScanNo, if_.n_if, cNIF);
     1154          logMsg(cMsg);
    9271155        }
    9281156      }
     
    9351163      }
    9361164
    937       iMBuff->scanNo  = cScanNo;
    938       iMBuff->cycleNo = cCycleNo;
    939 
    940       // Times.
    941       strncpy(iMBuff->datobs, cDateObs, 10);
    942       iMBuff->utc = cUTC;
    943       iMBuff->exposure = param_.intbase;
    944 
    945       // Source identification.
    946       sprintf(iMBuff->srcName, "%-16.16s",
    947               names_.su_name + (cSrcNo-1)*16);
    948       iMBuff->srcRA  = doubles_.su_ra[cSrcNo-1];
    949       iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
    950 
    951       // Rest frequency of the line of interest.
    952       iMBuff->restFreq = doubles_.rfreq;
    953       if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
    954         // Fix the HI rest frequency recorded for Parkes multibeam data.
    955         double reffreq  = doubles_.freq;
    956         double restfreq = doubles_.rfreq;
    957         if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
    958              fabs(reffreq - 1420.40575e6) < 100.0) {
    959           iMBuff->restFreq = 1420.40575e6;
    960         }
    961       }
    962 
    963       // Observation type.
    964       int j;
    965       for (j = 0; j < 15; j++) {
    966         iMBuff->obsType[j] = names_.card[11+j];
    967         if (iMBuff->obsType[j] == '\'') break;
    968       }
    969       iMBuff->obsType[j] = '\0';
    970 
    971       // Beam-dependent parameters.
    972       iMBuff->beamNo = beamNo;
    973 
    974       // Beam position at the specified time.
    975       if (cSUpos) {
    976         // Non-ATNF data that does not store the position in (u,v,w).
    977         iMBuff->ra  = doubles_.su_ra[cSrcNo-1];
    978         iMBuff->dec = doubles_.su_dec[cSrcNo-1];
    979       } else {
    980         iMBuff->ra  = cU;
    981         iMBuff->dec = cV;
    982       }
    983       cPosUTC[iBeamSel] = cW;
     1165#ifdef PKSIO_DEBUG
     1166      fprintf(stderr, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
     1167#endif
     1168
     1169      // Store IF-independent parameters only for the first IF of a new cycle,
     1170      // particularly because this is the only one for which the scan rates
     1171      // are computed above.
     1172      int firstIF = (iMBuff->nIF == 0);
     1173      if (firstIF) {
     1174        iMBuff->scanNo  = cScanNo;
     1175        iMBuff->cycleNo = cCycleNo;
     1176
     1177        // Times.
     1178        strcpy(iMBuff->datobs, cDateObs);
     1179        iMBuff->utc = cUTC;
     1180        iMBuff->exposure = param_.intbase;
     1181
     1182        // Source identification.
     1183        sprintf(iMBuff->srcName, "%-16.16s",
     1184                names_.su_name + (cSrcNo-1)*16);
     1185        iMBuff->srcName[16] = '\0';
     1186        iMBuff->srcRA  = doubles_.su_ra[cSrcNo-1];
     1187        iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
     1188
     1189        // Rest frequency of the line of interest.
     1190        iMBuff->restFreq = doubles_.rfreq;
     1191        if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
     1192          // Fix the HI rest frequency recorded for Parkes multibeam data.
     1193          double reffreq  = doubles_.freq;
     1194          double restfreq = doubles_.rfreq;
     1195          if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
     1196               fabs(reffreq - 1420.405752e6) < 100.0) {
     1197            iMBuff->restFreq = 1420.405752e6;
     1198          }
     1199        }
     1200
     1201        // Observation type.
     1202        int j;
     1203        for (j = 0; j < 15; j++) {
     1204          iMBuff->obsType[j] = names_.card[11+j];
     1205          if (iMBuff->obsType[j] == '\'') break;
     1206        }
     1207        iMBuff->obsType[j] = '\0';
     1208
     1209        // Beam-dependent parameters.
     1210        iMBuff->beamNo = beamNo;
     1211
     1212        // Beam position at the specified time.
     1213        if (cSUpos) {
     1214          // Non-ATNF data that does not store the position in (u,v,w).
     1215          iMBuff->ra  = doubles_.su_ra[cSrcNo-1];
     1216          iMBuff->dec = doubles_.su_dec[cSrcNo-1];
     1217        } else {
     1218          iMBuff->ra  = cU;
     1219          iMBuff->dec = cV;
     1220        }
     1221        cPosUTC[iBeamSel] = cW;
     1222        iMBuff->pCode = pCode;
     1223
     1224        // Store rates for next time.
     1225        iMBuff->raRate  =  raRate;
     1226        iMBuff->decRate = decRate;
     1227        iMBuff->paRate  =  paRate;
     1228      }
    9841229
    9851230      // IF-dependent parameters.
     
    9921237
    9931238      iIFSel = cIFSel[iIF];
    994       iMBuff->nIF++;
    995       iMBuff->IFno[iIFSel]  = cIFno;
     1239      if (iMBuff->IFno[iIFSel] == 0) {
     1240        iMBuff->nIF++;
     1241        iMBuff->IFno[iIFSel] = cIFno;
     1242      } else {
     1243        // Integration cycle written to the output file twice (the only known
     1244        // example is 1999-05-22_1914_000-031805_03v.hpf).
     1245        sprintf(cMsg, "WARNING: Integration cycle %d:%d, beam %2d, \n"
     1246                      "         IF %d was duplicated.", cScanNo, cCycleNo-1,
     1247                      beamNo, cIFno);
     1248        logMsg(cMsg);
     1249      }
    9961250      iMBuff->nChan[iIFSel] = nChan;
    9971251      iMBuff->nPol[iIFSel]  = cNPol[iIF];
     
    10931347
    10941348
    1095       // Parallactic angle.
    1096       iMBuff->parAngle = sc_.sc_cal[scq*iBeam + 11];
    1097 
    10981349      // Calibration factor applied to the data by the correlator.
    10991350      if (scq > 14) {
     
    11061357      }
    11071358
    1108       if (sc_.sc_ant <= anten_.nant) {
    1109         // No extra syscal information present.
    1110         iMBuff->extraSysCal = 0;
    1111         iMBuff->azimuth   = 0.0f;
    1112         iMBuff->elevation = 0.0f;
    1113         iMBuff->parAngle  = 0.0f;
    1114         iMBuff->focusAxi  = 0.0f;
    1115         iMBuff->focusTan  = 0.0f;
    1116         iMBuff->focusRot  = 0.0f;
    1117         iMBuff->temp      = 0.0f;
    1118         iMBuff->pressure  = 0.0f;
    1119         iMBuff->humidity  = 0.0f;
    1120         iMBuff->windSpeed = 0.0f;
    1121         iMBuff->windAz    = 0.0f;
    1122         strcpy(iMBuff->tcalTime, "                ");
    1123         iMBuff->refBeam = 0;
    1124 
    1125       } else {
    1126         // Additional information for Parkes Multibeam data.
    1127         int iOff = scq*(sc_.sc_ant - 1) - 1;
    1128         iMBuff->extraSysCal = 1;
    1129         iMBuff->azimuth   = sc_.sc_cal[iOff + 2];
    1130         iMBuff->elevation = sc_.sc_cal[iOff + 3];
    1131         iMBuff->parAngle  = sc_.sc_cal[iOff + 4];
    1132         iMBuff->focusAxi  = sc_.sc_cal[iOff + 5] * 1e-3;
    1133         iMBuff->focusTan  = sc_.sc_cal[iOff + 6] * 1e-3;
    1134         iMBuff->focusRot  = sc_.sc_cal[iOff + 7];
    1135         iMBuff->temp      = sc_.sc_cal[iOff + 8];
    1136         iMBuff->pressure  = sc_.sc_cal[iOff + 9];
    1137         iMBuff->humidity  = sc_.sc_cal[iOff + 10];
    1138         iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
    1139         iMBuff->windAz    = sc_.sc_cal[iOff + 12];
    1140 
    1141         char *tcalTime = iMBuff->tcalTime;
    1142         sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
     1359      if (firstIF) {
     1360        if (sc_.sc_ant <= anten_.nant) {
     1361          // No extra syscal information present.
     1362          iMBuff->extraSysCal = 0;
     1363          iMBuff->azimuth   = 0.0f;
     1364          iMBuff->elevation = 0.0f;
     1365          iMBuff->parAngle  = 0.0f;
     1366          iMBuff->focusAxi  = 0.0f;
     1367          iMBuff->focusTan  = 0.0f;
     1368          iMBuff->focusRot  = 0.0f;
     1369          iMBuff->temp      = 0.0f;
     1370          iMBuff->pressure  = 0.0f;
     1371          iMBuff->humidity  = 0.0f;
     1372          iMBuff->windSpeed = 0.0f;
     1373          iMBuff->windAz    = 0.0f;
     1374          strcpy(iMBuff->tcalTime, "                ");
     1375          iMBuff->refBeam = 0;
     1376
     1377        } else {
     1378          // Additional information for Parkes Multibeam data.
     1379          int iOff = scq*(sc_.sc_ant - 1) - 1;
     1380          iMBuff->extraSysCal = 1;
     1381
     1382          iMBuff->azimuth   = sc_.sc_cal[iOff + 2];
     1383          iMBuff->elevation = sc_.sc_cal[iOff + 3];
     1384          iMBuff->parAngle  = sc_.sc_cal[iOff + 4];
     1385
     1386          iMBuff->focusAxi  = sc_.sc_cal[iOff + 5] * 1e-3;
     1387          iMBuff->focusTan  = sc_.sc_cal[iOff + 6] * 1e-3;
     1388          iMBuff->focusRot  = sc_.sc_cal[iOff + 7];
     1389
     1390          iMBuff->temp      = sc_.sc_cal[iOff + 8];
     1391          iMBuff->pressure  = sc_.sc_cal[iOff + 9];
     1392          iMBuff->humidity  = sc_.sc_cal[iOff + 10];
     1393          iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
     1394          iMBuff->windAz    = sc_.sc_cal[iOff + 12];
     1395
     1396          char *tcalTime = iMBuff->tcalTime;
     1397          sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
     1398          tcalTime[16] = '\0';
    11431399
    11441400#ifndef AIPS_LITTLE_ENDIAN
    1145         // Do byte swapping on the ASCII date string.
    1146         for (int j = 0; j < 16; j += 4) {
    1147           char ctmp;
    1148           ctmp = tcalTime[j];
    1149           tcalTime[j]   = tcalTime[j+3];
    1150           tcalTime[j+3] = ctmp;
    1151           ctmp = tcalTime[j+1];
    1152           tcalTime[j+1] = tcalTime[j+2];
    1153           tcalTime[j+2] = ctmp;
    1154         }
     1401          // Do byte swapping on the ASCII date string.
     1402          for (int j = 0; j < 16; j += 4) {
     1403            char ctmp;
     1404            ctmp = tcalTime[j];
     1405            tcalTime[j]   = tcalTime[j+3];
     1406            tcalTime[j+3] = ctmp;
     1407            ctmp = tcalTime[j+1];
     1408            tcalTime[j+1] = tcalTime[j+2];
     1409            tcalTime[j+2] = ctmp;
     1410          }
    11551411#endif
    11561412
    1157         // Reference beam number.
    1158         float refbeam = sc_.sc_cal[iOff + 17];
    1159         if (refbeam > 0.0f || refbeam < 100.0f) {
    1160           iMBuff->refBeam = int(refbeam);
    1161         } else {
    1162           iMBuff->refBeam = 0;
     1413          // Reference beam number.
     1414          float refbeam = sc_.sc_cal[iOff + 17];
     1415          if (refbeam > 0.0f || refbeam < 100.0f) {
     1416            iMBuff->refBeam = int(refbeam);
     1417          } else {
     1418            iMBuff->refBeam = 0;
     1419          }
    11631420        }
    11641421      }
     
    11851442  while (numErr < 10) {
    11861443    int lastjstat = jstat;
    1187     rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
    1188               &cBin, &cIFno, &cSrcNo);
    1189 
    1190     switch(jstat) {
     1444
     1445    switch(rpfitsin(jstat)) {
    11911446    case -1:
    11921447      // Read failed; retry.
    11931448      numErr++;
    1194       fprintf(stderr, "RPFITS read failed - retrying.\n");
     1449      logMsg("WARNING: RPFITS read failed - retrying.");
    11951450      jstat = 0;
    11961451      break;
     
    12481503    default:
    12491504      // Shouldn't reach here.
    1250       fprintf(stderr, "Unrecognized RPFITSIN return code: %d (retrying)\n",
    1251               jstat);
     1505      sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
     1506                    "(retrying).", jstat);
     1507      logMsg(cMsg);
    12521508      jstat = 0;
    12531509      break;
     
    12551511  }
    12561512
    1257   fprintf(stderr, "ERROR, RPFITS read failed too many times.\n");
     1513  logMsg("ERROR: RPFITS read failed too many times.");
    12581514  return 2;
     1515}
     1516
     1517//----------------------------------------------------- MBFITSreader::rpfitsin
     1518
     1519// Wrapper around RPFITSIN that reports errors.  Returned RPFITSIN subroutine
     1520// arguments are captured as MBFITSreader member variables.
     1521
     1522int MBFITSreader::rpfitsin(int &jstat)
     1523
     1524{
     1525  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
     1526            &cBin, &cIFno, &cSrcNo);
     1527
     1528  // Handle messages from RPFITSIN.
     1529  if (names_.errmsg[0] != ' ') {
     1530    int i;
     1531    for (i = 80; i > 0; i--) {
     1532      if (names_.errmsg[i-1] != ' ') break;
     1533    }
     1534
     1535    sprintf(cMsg, "WARNING: Cycle %d:%03d, RPFITSIN reported -\n"
     1536                  "         %.*s", cScanNo, cCycleNo, i, names_.errmsg);
     1537    logMsg(cMsg);
     1538  }
     1539
     1540  return jstat;
     1541}
     1542
     1543//------------------------------------------------------- MBFITSreader::fixPos
     1544
     1545// Check and, if necessary, repair a position timestamp.
     1546//
     1547// Problems with the position timestamp manifest themselves via the scan rate:
     1548//
     1549//   1) Zero scan rate pairs, 1997/02/28 to 1998/01/07
     1550//
     1551//      These occur because the position timestamp for the first integration
     1552//      of the pair is erroneous; the value recorded is t/1000, where t is the
     1553//      true value.
     1554//        Earliest known: 97-02-28_1725_132653-42_258a.hpf
     1555//          Latest known: 98-01-02_1923_095644-50_165c.hpf
     1556//        (time range chosen to encompass observing runs).
     1557//
     1558//   2) Slow-fast scan rate pairs (0.013 - 0.020 deg/s),
     1559//        1997/03/28 to 1998/01/07.
     1560//
     1561//      The UTC position timestamp is 1.0s later than it should be (never
     1562//      earlier), almost certainly arising from an error in the telescope
     1563//      control system.
     1564//        Earliest known: 97-03-28_0150_010420-74_008d.hpf
     1565//          Latest known: 98-01-04_1502_065150-02_177c.hpf
     1566//        (time range chosen to encompass observing runs).
     1567//
     1568//   3) Slow-fast scan rate pairs (0.015 - 0.018 deg/s),
     1569//        1999/05/20 to 2001/07/12 (HIPASS and ZOA),
     1570//        2001/09/02 to 2001/12/04 (HIPASS and ZOA),
     1571//        2002/03/28 to 2002/05/13 (ZOA only),
     1572//        2003/04/26 to 2003/06/09 (ZOA only).
     1573//        Earliest known: 1999-05-20_1818_175720-50_297e.hpf
     1574//          Latest known: 2001-12-04_1814_065531p14_173e.hpf (HIPASS)
     1575//                        2003-06-09_1924_352-085940_-6c.hpf (ZOA)
     1576//
     1577//      Caused by the Linux signalling NaN problem.  IEEE "signalling" NaNs
     1578//      are silently transformed to "quiet" NaNs during assignment by setting
     1579//      bit 22.  This affected RPFITS because of its use of VAX-format
     1580//      floating-point numbers which, with their permuted bytes, may sometimes
     1581//      appear as signalling NaNs.
     1582//
     1583//      The problem arose when the linux correlator came online and was
     1584//      fixed with a workaround to the RPFITS library (repeated episodes
     1585//      are probably due to use of an older version of the library).  It
     1586//      should not have affected the data significantly because of the
     1587//      low relative error, which ranges from 0.0000038 to 0.0000076, but
     1588//      it is important for the computation of scan rates which requires
     1589//      taking the difference of two large UTC timestamps, one or other
     1590//      of which will have 0.5s added to it.
     1591//
     1592// The return value identifies which, if any, of these problems was repaired.
     1593
     1594int MBFITSreader::fixw(
     1595  const char *datobs,
     1596  int    cycleNo,
     1597  int    beamNo,
     1598  double avRate[2],
     1599  double thisRA,
     1600  double thisDec,
     1601  double thisUTC,
     1602  double nextRA,
     1603  double nextDec,
     1604  float &nextUTC)
     1605{
     1606  if (strcmp(datobs, "2003-06-09") > 0) {
     1607    return 0;
     1608
     1609  } else if (strcmp(datobs, "1998-01-07") <= 0) {
     1610    if (nextUTC < thisUTC && (nextUTC + 86400.0) > (thisUTC + 600.0)) {
     1611      // Possible scaling problem.
     1612      double diff = nextUTC*1000.0 - thisUTC;
     1613      if (0.0 < diff && diff < 600.0) {
     1614        nextUTC *= 1000.0;
     1615        return 1;
     1616      } else {
     1617        // Irreparable.
     1618        return -1;
     1619      }
     1620    }
     1621
     1622    if (cycleNo > 2) {
     1623      if (beamNo == 1) {
     1624        // This test is only reliable for beam 1.
     1625        double dUTC = nextUTC - thisUTC;
     1626        if (dUTC < 0.0) dUTC += 86400.0;
     1627
     1628        // Guard against RA cycling through 24h in either direction.
     1629        if (fabs(nextRA - thisRA) > PI) {
     1630          if (nextRA < thisRA) {
     1631            nextRA += TWOPI;
     1632          } else {
     1633            nextRA -= TWOPI;
     1634          }
     1635        }
     1636
     1637        double  dRA = (nextRA  - thisRA) * cos(nextDec);
     1638        double dDec =  nextDec - thisDec;
     1639        double  arc = sqrt(dRA*dRA + dDec*dDec);
     1640
     1641        double averate = sqrt(avRate[0]*avRate[0] + avRate[1]*avRate[1]);
     1642        double diff1 = fabs(averate - arc/(dUTC-1.0));
     1643        double diff2 = fabs(averate - arc/dUTC);
     1644        if ((diff1 < diff2) && (diff1 < 0.05*averate)) {
     1645          nextUTC -= 1.0;
     1646          cCode5 = cycleNo;
     1647          return 2;
     1648        } else {
     1649          cCode5 = 0;
     1650        }
     1651
     1652      } else {
     1653        if (cycleNo == cCode5) {
     1654          nextUTC -= 1.0;
     1655          return 2;
     1656        }
     1657      }
     1658    }
     1659
     1660  } else if ((strcmp(datobs, "1999-05-20") >= 0 &&
     1661              strcmp(datobs, "2001-07-12") <= 0) ||
     1662             (strcmp(datobs, "2001-09-02") >= 0 &&
     1663              strcmp(datobs, "2001-12-04") <= 0) ||
     1664             (strcmp(datobs, "2002-03-28") >= 0 &&
     1665              strcmp(datobs, "2002-05-13") <= 0) ||
     1666             (strcmp(datobs, "2003-04-26") >= 0 &&
     1667              strcmp(datobs, "2003-06-09") <= 0)) {
     1668    // Signalling NaN problem, e.g. 1999-07-26_1839_011106-74_009c.hpf.
     1669    // Position timestamps should always be an integral number of seconds.
     1670    double resid = nextUTC - int(nextUTC);
     1671    if (resid == 0.5) {
     1672      nextUTC -= 0.5;
     1673      return 3;
     1674    }
     1675  }
     1676
     1677  return 0;
    12591678}
    12601679
     
    12921711  }
    12931712}
     1713
     1714//-------------------------------------------------------------------- utcDiff
     1715
     1716// Subtract two UTCs (s) allowing for any plausible number of cycles through
     1717// 86400s, returning a result in the range [-43200, +43200]s.
     1718
     1719double MBFITSreader::utcDiff(double utc1, double utc2)
     1720{
     1721  double diff = utc1 - utc2;
     1722
     1723  if (diff > 43200.0) {
     1724    diff -= 86400.0;
     1725    while (diff > 43200.0) diff -= 86400.0;
     1726  } else if (diff < -43200.0) {
     1727    diff += 86400.0;
     1728    while (diff < -43200.0) diff += 86400.0;
     1729  }
     1730
     1731  return diff;
     1732}
     1733
     1734//------------------------------------------------------- scanRate & applyRate
     1735
     1736// Compute and apply the scan rate corrected for grid convergence.  (ra0,dec0)
     1737// are the coordinates of the central beam, assumed to be the tracking centre.
     1738// The rate computed in RA will be a rate of change of angular distance in the
     1739// direction of increasing RA at the position of the central beam.  Similarly
     1740// for declination.  Angles in radian, time in s.
     1741
     1742void MBFITSreader::scanRate(
     1743  double ra0,
     1744  double dec0,
     1745  double ra1,
     1746  double dec1,
     1747  double ra2,
     1748  double dec2,
     1749  double dt,
     1750  double &raRate,
     1751  double &decRate)
     1752{
     1753  // Transform to a system where the central beam lies on the equator at 12h.
     1754  eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
     1755  eulerx(ra2, dec2, ra0+HALFPI, -dec0, -HALFPI, ra2, dec2);
     1756
     1757  raRate  = (ra2  - ra1)  / dt;
     1758  decRate = (dec2 - dec1) / dt;
     1759}
     1760
     1761
     1762void MBFITSreader::applyRate(
     1763  double ra0,
     1764  double dec0,
     1765  double ra1,
     1766  double dec1,
     1767  double raRate,
     1768  double decRate,
     1769  double dt,
     1770  double &ra2,
     1771  double &dec2)
     1772{
     1773  // Transform to a system where the central beam lies on the equator at 12h.
     1774  eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
     1775
     1776  ra2  = ra1  + (raRate  * dt);
     1777  dec2 = dec1 + (decRate * dt);
     1778
     1779  // Transform back.
     1780  eulerx(ra2, dec2, -HALFPI, dec0, ra0+HALFPI, ra2, dec2);
     1781}
     1782
     1783//--------------------------------------------------------------------- eulerx
     1784
     1785void MBFITSreader::eulerx(
     1786  double lng0,
     1787  double lat0,
     1788  double phi0,
     1789  double theta,
     1790  double phi,
     1791  double &lng1,
     1792  double &lat1)
     1793
     1794// Applies the Euler angle based transformation of spherical coordinates.
     1795//
     1796//     phi0  Longitude of the ascending node in the old system, radians.  The
     1797//           ascending node is the point of intersection of the equators of
     1798//           the two systems such that the equator of the new system crosses
     1799//           from south to north as viewed in the old system.
     1800//
     1801//    theta  Angle between the poles of the two systems, radians.  THETA is
     1802//           positive for a positive rotation about the ascending node.
     1803//
     1804//      phi  Longitude of the ascending node in the new system, radians.
     1805
     1806{
     1807  // Compute intermediaries.
     1808  double lng0p  = lng0 - phi0;
     1809  double slng0p = sin(lng0p);
     1810  double clng0p = cos(lng0p);
     1811  double slat0  = sin(lat0);
     1812  double clat0  = cos(lat0);
     1813  double ctheta = cos(theta);
     1814  double stheta = sin(theta);
     1815
     1816  double x = clat0*clng0p;
     1817  double y = clat0*slng0p*ctheta + slat0*stheta;
     1818
     1819  // Longitude in the new system.
     1820  if (x != 0.0 || y != 0.0) {
     1821    lng1 = phi + atan2(y, x);
     1822  } else {
     1823    // Longitude at the poles in the new system is consistent with that
     1824    // specified in the old system.
     1825    lng1 = phi + lng0p;
     1826  }
     1827  lng1 = fmod(lng1, TWOPI);
     1828  if (lng1 < 0.0) lng1 += TWOPI;
     1829
     1830  lat1 = asin(slat0*ctheta - clat0*stheta*slng0p);
     1831}
  • trunk/external/atnf/PKSIO/MBFITSreader.h

    r1427 r1452  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBFITSreader.h,v 19.15 2008-06-26 02:14:36 cal103 Exp $
     29//# $Id: MBFITSreader.h,v 19.21 2008-11-17 06:33:10 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
     
    3939
    4040#include <atnf/PKSIO/FITSreader.h>
    41 #include <atnf/PKSIO/PKSMBrecord.h>
     41#include <atnf/PKSIO/MBrecord.h>
     42
     43using namespace std;
    4244
    4345// <summary>
     
    102104
    103105    // Read the next data record.
    104     virtual int read(PKSMBrecord &record);
     106    virtual int read(MBrecord &record);
    105107
    106108    // Close the RPFITS file.
     
    112114    float cUTC, cU, cV, cW, *cVis, *cWgt;
    113115
    114     char   cDateObs[10];
     116    char   cDateObs[12];
    115117    int    *cBeamSel, *cChanOff, cFirst, *cIFSel, cInterp, cIntTime, cMBopen,
    116118           cMopra, cNBeamSel, cNBin, cRetry, cSUpos, *cXpolOff;
     
    119121    int    cEOF, cEOS, cFlushBin, cFlushIF, cFlushing;
    120122    double *cPosUTC;
    121     PKSMBrecord *cBuffer;
     123    MBrecord *cBuffer;
    122124
     125    // Scan and cycle number bookkeeping.
    123126    int    cCycleNo, cScanNo;
    124127    double cPrevUTC;
     
    126129    // Read the next data record from the RPFITS file.
    127130    int rpget(int syscalonly, int &EOS);
     131    int rpfitsin(int &jstat);
     132
     133    // Check and, if necessary, repair a position timestamp.
     134    int    cCode5, cNRate;
     135    double cAvRate[2];
     136    int fixw(const char *datobs, int cycleNo, int beamNo, double avRate[2],
     137             double thisRA, double thisDec, double thisUTC,
     138             double nextRA, double nextDec, float &nextUTC);
     139
     140    // Subtract two UTCs (s).
     141    double utcDiff(double utc1, double utc2);
     142
     143    // Compute and apply the scan rate corrected for grid convergence.
     144    double cRA0, cDec0;
     145    void  scanRate(double ra0, double dec0,
     146                   double ra1, double dec1,
     147                   double ra2, double dec2, double dt,
     148                   double &raRate, double &decRate);
     149    void applyRate(double ra0, double dec0,
     150                   double ra1, double dec1,
     151                   double raRate, double decRate, double dt,
     152                   double &ra2, double &dec2);
     153    void eulerx(double lng0, double lat0, double phi0, double theta,
     154                double phi, double &lng1, double &lat1);
    128155};
    129156
  • trunk/external/atnf/PKSIO/PKSFITSreader.cc

    r1427 r1452  
    2525//#                        Charlottesville, VA 22903-2475 USA
    2626//#
    27 //# $Id: PKSFITSreader.cc,v 19.14 2008-06-26 02:07:21 cal103 Exp $
     27//# $Id: PKSFITSreader.cc,v 19.21 2008-11-17 06:54:25 cal103 Exp $
    2828//#---------------------------------------------------------------------------
    2929//# Original: 2000/08/02, Mark Calabretta, ATNF
    3030//#---------------------------------------------------------------------------
    3131
     32#include <atnf/PKSIO/PKSmsg.h>
    3233#include <atnf/PKSIO/MBFITSreader.h>
    3334#include <atnf/PKSIO/SDFITSreader.h>
    3435#include <atnf/PKSIO/PKSFITSreader.h>
    35 #include <atnf/PKSIO/PKSMBrecord.h>
    36 
     36#include <atnf/PKSIO/PKSrecord.h>
     37
     38#include <casa/stdio.h>
    3739#include <casa/Arrays/Array.h>
    3840#include <casa/BasicMath/Math.h>
    3941#include <casa/Quanta/MVTime.h>
    40 
    41 #include <casa/stdio.h>
    4242
    4343//----------------------------------------------- PKSFITSreader::PKSFITSreader
     
    5050        const Bool   interpolate)
    5151{
    52   cFITSMBrec.setNIFs(1);
     52  cMBrec.setNIFs(1);
    5353
    5454  if (fitsType == "SDFITS") {
     
    5757    cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
    5858  }
     59
     60  // By default, messages are written to stderr.
     61  initMsg();
    5962}
    6063
     
    6770  close();
    6871  delete cReader;
     72}
     73
     74//------------------------------------------------------ PKSFITSreader::setMsg
     75
     76// Set message disposition.  If fd is non-zero messages will be written
     77// to that file descriptor, else stored for retrieval by getMsg().
     78
     79Int PKSFITSreader::setMsg(FILE *fd)
     80{
     81  PKSmsg::setMsg(fd);
     82  cReader->setMsg(fd);
     83
     84  return 0;
    6985}
    7086
     
    8399        Bool   &haveSpectra)
    84100{
     101  clearMsg();
     102
    85103  int    extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
    86104         nIF, *nPol_, status;
    87   if ((status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF,
    88                               cIFs, nChan_, nPol_, haveXPol_, haveBase_,
    89                               haveSpectra_, extraSysCal))) {
     105  status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF, cIFs,
     106                         nChan_, nPol_, haveXPol_, haveBase_, haveSpectra_,
     107                         extraSysCal);
     108  logMsg(cReader->getMsg());
     109  cReader->clearMsg();
     110  if (status) {
    90111    return status;
    91112  }
     
    146167  char   bunit_[32], datobs[32], dopplerFrame_[32], observer_[32],
    147168         obsType_[32], project_[32], radecsys[32], telescope[32];
     169  int    status;
    148170  float  equinox_;
    149171  double antPos[3], utc;
    150172
    151   if (cReader->getHeader(observer_, project_, telescope, antPos, obsType_,
    152                          bunit_, equinox_, radecsys, dopplerFrame_,
    153                          datobs, utc, refFreq, bandwidth)) {
     173  status = cReader->getHeader(observer_, project_, telescope, antPos,
     174                              obsType_, bunit_, equinox_, radecsys,
     175                              dopplerFrame_, datobs, utc, refFreq, bandwidth);
     176  logMsg(cReader->getMsg());
     177  cReader->clearMsg();
     178  if (status) {
    154179    return 1;
    155180  }
     
    187212  double *startfreq, *endfreq;
    188213
    189   Int status;
    190   if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
     214  Int status = cReader->getFreqInfo(nIF, startfreq, endfreq);
     215
     216  logMsg(cReader->getMsg());
     217  cReader->clearMsg();
     218  if (!status) {
    191219    startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
    192220    endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
     
    208236        const Bool getSpectra,
    209237        const Bool getXPol,
    210         const Bool getFeedPos)
     238        const Int  coordSys)
    211239{
    212240  // Apply beam selection.
     
    275303  cGetSpectra = getSpectra;
    276304  cGetXPol    = getXPol;
    277   cGetFeedPos = getFeedPos;
     305  cCoordSys   = coordSys;
    278306
    279307  uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
    280                                   cGetFeedPos);
     308                                  cCoordSys);
     309  logMsg(cReader->getMsg());
     310  cReader->clearMsg();
    281311
    282312  delete [] end;
     
    301331  double* posns;
    302332
    303   Int status;
    304   if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
     333  Int status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns);
     334  logMsg(cReader->getMsg());
     335  cReader->clearMsg();
     336
     337  if (!status) {
    305338    timeSpan.resize(2);
    306339
     
    321354// Read the next data record.
    322355
    323 Int PKSFITSreader::read(MBrecord &MBrec)
    324 {
    325   Int status;
    326 
    327   if ((status = cReader->read(cFITSMBrec))) {
     356Int PKSFITSreader::read(PKSrecord &pksrec)
     357{
     358  Int status = cReader->read(cMBrec);
     359  logMsg(cReader->getMsg());
     360  cReader->clearMsg();
     361
     362  if (status) {
    328363    if (status != -1) {
    329364      status = 1;
     
    334369
    335370
    336   uInt nChan = cFITSMBrec.nChan[0];
    337   uInt nPol  = cFITSMBrec.nPol[0];
    338 
    339   MBrec.scanNo  = cFITSMBrec.scanNo;
    340   MBrec.cycleNo = cFITSMBrec.cycleNo;
     371  uInt nChan = cMBrec.nChan[0];
     372  uInt nPol  = cMBrec.nPol[0];
     373
     374  pksrec.scanNo  = cMBrec.scanNo;
     375  pksrec.cycleNo = cMBrec.cycleNo;
    341376
    342377  // Extract MJD.
    343378  Int day, month, year;
    344   sscanf(cFITSMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
    345   MBrec.mjd = MVTime(year, month, Double(day)).day() + cFITSMBrec.utc/86400.0;
    346 
    347   MBrec.interval  = cFITSMBrec.exposure;
    348 
    349   MBrec.fieldName = trim(cFITSMBrec.srcName);
    350   MBrec.srcName   = MBrec.fieldName;
    351 
    352   MBrec.srcDir.resize(2);
    353   MBrec.srcDir(0) = cFITSMBrec.srcRA;
    354   MBrec.srcDir(1) = cFITSMBrec.srcDec;
    355 
    356   MBrec.srcPM.resize(2);
    357   MBrec.srcPM(0)  = 0.0;
    358   MBrec.srcPM(1)  = 0.0;
    359   MBrec.srcVel    = 0.0;
    360   MBrec.obsType   = trim(cFITSMBrec.obsType);
    361 
    362   MBrec.IFno = cFITSMBrec.IFno[0];
    363   Double chanWidth = fabs(cFITSMBrec.fqDelt[0]);
    364   MBrec.refFreq   = cFITSMBrec.fqRefVal[0];
    365   MBrec.bandwidth = chanWidth * nChan;
    366   MBrec.freqInc   = cFITSMBrec.fqDelt[0];
    367   MBrec.restFreq  = cFITSMBrec.restFreq;
    368 
    369   MBrec.tcal.resize(nPol);
     379  sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
     380  pksrec.mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
     381
     382  pksrec.interval  = cMBrec.exposure;
     383
     384  pksrec.fieldName = trim(cMBrec.srcName);
     385  pksrec.srcName   = pksrec.fieldName;
     386
     387  pksrec.srcDir.resize(2);
     388  pksrec.srcDir(0) = cMBrec.srcRA;
     389  pksrec.srcDir(1) = cMBrec.srcDec;
     390
     391  pksrec.srcPM.resize(2);
     392  pksrec.srcPM(0)  = 0.0;
     393  pksrec.srcPM(1)  = 0.0;
     394  pksrec.srcVel    = 0.0;
     395  pksrec.obsType   = trim(cMBrec.obsType);
     396
     397  pksrec.IFno = cMBrec.IFno[0];
     398  Double chanWidth = fabs(cMBrec.fqDelt[0]);
     399  pksrec.refFreq   = cMBrec.fqRefVal[0];
     400  pksrec.bandwidth = chanWidth * nChan;
     401  pksrec.freqInc   = cMBrec.fqDelt[0];
     402  pksrec.restFreq  = cMBrec.restFreq;
     403
     404  pksrec.tcal.resize(nPol);
    370405  for (uInt ipol = 0; ipol < nPol; ipol++) {
    371     MBrec.tcal(ipol) = cFITSMBrec.tcal[0][ipol];
    372   }
    373   MBrec.tcalTime  = trim(cFITSMBrec.tcalTime);
    374   MBrec.azimuth   = cFITSMBrec.azimuth;
    375   MBrec.elevation = cFITSMBrec.elevation;
    376   MBrec.parAngle  = cFITSMBrec.parAngle;
    377   MBrec.focusAxi  = cFITSMBrec.focusAxi;
    378   MBrec.focusTan  = cFITSMBrec.focusTan;
    379   MBrec.focusRot  = cFITSMBrec.focusRot;
    380 
    381   MBrec.temperature = cFITSMBrec.temp;
    382   MBrec.pressure    = cFITSMBrec.pressure;
    383   MBrec.humidity    = cFITSMBrec.humidity;
    384   MBrec.windSpeed   = cFITSMBrec.windSpeed;
    385   MBrec.windAz      = cFITSMBrec.windAz;
    386 
    387   MBrec.refBeam = cFITSMBrec.refBeam;
    388   MBrec.beamNo  = cFITSMBrec.beamNo;
    389 
    390   MBrec.direction.resize(2);
    391   MBrec.direction(0) = cFITSMBrec.ra;
    392   MBrec.direction(1) = cFITSMBrec.dec;
    393 
    394   MBrec.scanRate.resize(2);
    395   MBrec.scanRate(0)  = cFITSMBrec.raRate;
    396   MBrec.scanRate(1)  = cFITSMBrec.decRate;
    397   MBrec.rateAge      = cFITSMBrec.rateAge;
    398   MBrec.rateson      = cFITSMBrec.rateson;
    399 
    400   MBrec.tsys.resize(nPol);
    401   MBrec.sigma.resize(nPol);
    402   MBrec.calFctr.resize(nPol);
     406    pksrec.tcal(ipol) = cMBrec.tcal[0][ipol];
     407  }
     408  pksrec.tcalTime  = trim(cMBrec.tcalTime);
     409  pksrec.azimuth   = cMBrec.azimuth;
     410  pksrec.elevation = cMBrec.elevation;
     411  pksrec.parAngle  = cMBrec.parAngle;
     412
     413  pksrec.focusAxi  = cMBrec.focusAxi;
     414  pksrec.focusTan  = cMBrec.focusTan;
     415  pksrec.focusRot  = cMBrec.focusRot;
     416
     417  pksrec.temperature = cMBrec.temp;
     418  pksrec.pressure    = cMBrec.pressure;
     419  pksrec.humidity    = cMBrec.humidity;
     420  pksrec.windSpeed   = cMBrec.windSpeed;
     421  pksrec.windAz      = cMBrec.windAz;
     422
     423  pksrec.refBeam = cMBrec.refBeam;
     424  pksrec.beamNo  = cMBrec.beamNo;
     425
     426  pksrec.direction.resize(2);
     427  pksrec.direction(0) = cMBrec.ra;
     428  pksrec.direction(1) = cMBrec.dec;
     429  pksrec.pCode        = cMBrec.pCode;
     430  pksrec.rateAge      = cMBrec.rateAge;
     431  pksrec.scanRate.resize(2);
     432  pksrec.scanRate(0)  = cMBrec.raRate;
     433  pksrec.scanRate(1)  = cMBrec.decRate;
     434  pksrec.paRate       = cMBrec.paRate;
     435
     436  pksrec.tsys.resize(nPol);
     437  pksrec.sigma.resize(nPol);
     438  pksrec.calFctr.resize(nPol);
    403439  for (uInt ipol = 0; ipol < nPol; ipol++) {
    404     MBrec.tsys(ipol)  = cFITSMBrec.tsys[0][ipol];
    405     MBrec.sigma(ipol) = (MBrec.tsys(ipol) / 0.81) /
    406                          sqrt(MBrec.interval * chanWidth);
    407     MBrec.calFctr(ipol) = cFITSMBrec.calfctr[0][ipol];
    408   }
    409 
    410   if (cFITSMBrec.haveBase) {
    411     MBrec.baseLin.resize(2,nPol);
    412     MBrec.baseSub.resize(9,nPol);
     440    pksrec.tsys(ipol)  = cMBrec.tsys[0][ipol];
     441    pksrec.sigma(ipol) = (pksrec.tsys(ipol) / 0.81) /
     442                            sqrt(pksrec.interval * chanWidth);
     443    pksrec.calFctr(ipol) = cMBrec.calfctr[0][ipol];
     444  }
     445
     446  if (cMBrec.haveBase) {
     447    pksrec.baseLin.resize(2,nPol);
     448    pksrec.baseSub.resize(9,nPol);
    413449
    414450    for (uInt ipol = 0; ipol < nPol; ipol++) {
    415       MBrec.baseLin(0,ipol) = cFITSMBrec.baseLin[0][ipol][0];
    416       MBrec.baseLin(1,ipol) = cFITSMBrec.baseLin[0][ipol][1];
     451      pksrec.baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
     452      pksrec.baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
    417453
    418454      for (uInt j = 0; j < 9; j++) {
    419         MBrec.baseSub(j,ipol) = cFITSMBrec.baseSub[0][ipol][j];
     455        pksrec.baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
    420456      }
    421457    }
    422458
    423459  } else {
    424     MBrec.baseLin.resize(0,0);
    425     MBrec.baseSub.resize(0,0);
    426   }
    427 
    428   if (cGetSpectra && cFITSMBrec.haveSpectra) {
    429     MBrec.spectra.resize(nChan,nPol);
    430     MBrec.spectra.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.spectra[0],
     460    pksrec.baseLin.resize(0,0);
     461    pksrec.baseSub.resize(0,0);
     462  }
     463
     464  if (cGetSpectra && cMBrec.haveSpectra) {
     465    pksrec.spectra.resize(nChan,nPol);
     466    pksrec.spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0],
    431467      SHARE);
    432468
    433     MBrec.flagged.resize(nChan,nPol);
    434     MBrec.flagged.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.flagged[0],
     469    pksrec.flagged.resize(nChan,nPol);
     470    pksrec.flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0],
    435471      SHARE);
    436472
    437473  } else {
    438     MBrec.spectra.resize(0,0);
    439     MBrec.flagged.resize(0,0);
     474    pksrec.spectra.resize(0,0);
     475    pksrec.flagged.resize(0,0);
    440476  }
    441477
    442478  if (cGetXPol) {
    443     MBrec.xCalFctr = Complex(cFITSMBrec.xcalfctr[0][0],
    444                              cFITSMBrec.xcalfctr[0][1]);
    445     MBrec.xPol.resize(nChan);
    446     MBrec.xPol.takeStorage(IPosition(1,nChan), (Complex *)cFITSMBrec.xpol[0],
     479    pksrec.xCalFctr = Complex(cMBrec.xcalfctr[0][0],
     480                             cMBrec.xcalfctr[0][1]);
     481    pksrec.xPol.resize(nChan);
     482    pksrec.xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0],
    447483      SHARE);
    448484  }
     
    451487}
    452488
    453 //-------------------------------------------------------- PKSFITSreader::read
    454 
    455 // Read the next data record, just the basics.
    456 
    457 Int PKSFITSreader::read(
    458         Int           &IFno,
    459         Vector<Float> &tsys,
    460         Vector<Float> &calFctr,
    461         Matrix<Float> &baseLin,
    462         Matrix<Float> &baseSub,
    463         Matrix<Float> &spectra,
    464         Matrix<uChar> &flagged)
    465 {
    466   Int status;
    467 
    468   if ((status = cReader->read(cFITSMBrec))) {
    469     if (status != -1) {
    470       status = 1;
    471     }
    472 
    473     return status;
    474   }
    475 
    476   IFno = cFITSMBrec.IFno[0];
    477 
    478   uInt nChan = cFITSMBrec.nChan[0];
    479   uInt nPol  = cFITSMBrec.nPol[0];
    480 
    481   tsys.resize(nPol);
    482   calFctr.resize(nPol);
    483   for (uInt ipol = 0; ipol < nPol; ipol++) {
    484     tsys(ipol) = cFITSMBrec.tsys[0][ipol];
    485     calFctr(ipol) = cFITSMBrec.calfctr[0][ipol];
    486   }
    487 
    488   if (cFITSMBrec.haveBase) {
    489     baseLin.resize(2,nPol);
    490     baseSub.resize(9,nPol);
    491 
    492     for (uInt ipol = 0; ipol < nPol; ipol++) {
    493       baseLin(0,ipol) = cFITSMBrec.baseLin[0][ipol][0];
    494       baseLin(1,ipol) = cFITSMBrec.baseLin[0][ipol][1];
    495 
    496       for (uInt j = 0; j < 9; j++) {
    497         baseSub(j,ipol) = cFITSMBrec.baseSub[0][ipol][j];
    498       }
    499     }
    500 
    501   } else {
    502     baseLin.resize(0,0);
    503     baseSub.resize(0,0);
    504   }
    505 
    506   if (cGetSpectra && cFITSMBrec.haveSpectra) {
    507     spectra.resize(nChan,nPol);
    508     spectra.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.spectra[0],
    509       SHARE);
    510 
    511     flagged.resize(nChan,nPol);
    512     flagged.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.flagged[0],
    513       SHARE);
    514 
    515   } else {
    516     spectra.resize(0,0);
    517     flagged.resize(0,0);
    518   }
    519 
    520   return 0;
    521 }
    522 
    523489//------------------------------------------------------- PKSFITSreader::close
    524490
     
    528494{
    529495  cReader->close();
     496  logMsg(cReader->getMsg());
     497  cReader->clearMsg();
    530498}
    531499
  • trunk/external/atnf/PKSIO/PKSFITSreader.h

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSFITSreader.h,v 19.13 2008-06-26 02:02:43 cal103 Exp $
     28//# $Id: PKSFITSreader.h,v 19.17 2008-11-17 06:38:05 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# This class is basically a wrapper class for reading data from either an
     
    3939
    4040#include <atnf/PKSIO/FITSreader.h>
     41#include <atnf/PKSIO/PKSrecord.h>
    4142#include <atnf/PKSIO/PKSreader.h>
    4243
    4344#include <casa/aips.h>
     45#include <casa/stdio.h>
    4446#include <casa/Arrays/Vector.h>
    4547#include <casa/Arrays/Matrix.h>
     
    4749#include <casa/BasicSL/String.h>
    4850
     51#include <casa/namespace.h>
     52
    4953// <summary>
    5054// Class to read Parkes Multibeam data from a FITS file.
    5155// </summary>
    5256
    53 #include <casa/namespace.h>
    5457class PKSFITSreader : public PKSreader
    5558{
     
    6366    // Destructor.
    6467    virtual ~PKSFITSreader();
     68
     69    // Set message disposition.
     70    virtual Int setMsg(
     71        FILE *fd = 0x0);
    6572
    6673    // Open the FITS file for reading.
     
    104111        const Bool getSpectra = True,
    105112        const Bool getXPol    = False,
    106         const Bool getFeedPos = False);
     113        const Int  coordSys   = 0);
    107114
    108115    // Find the range of the data selected in time and position.
     
    114121
    115122    // Read the next data record.
    116     virtual Int read(MBrecord &mbrec);
    117 
    118     // Read the next data record, just the basics.
    119     virtual Int read(
    120         Int           &IFno,
    121         Vector<Float> &tsys,
    122         Vector<Float> &calFctr,
    123         Matrix<Float> &baseLin,
    124         Matrix<Float> &baseSub,
    125         Matrix<Float> &spectra,
    126         Matrix<uChar> &flagged);
     123    virtual Int read(PKSrecord &pksrec);
    127124
    128125    // Close the FITS file.
     
    132129    Int    *cBeams, *cIFs;
    133130    uInt   cNBeam, cNIF;
    134     PKSMBrecord cFITSMBrec;
     131    MBrecord cMBrec;
    135132    FITSreader  *cReader;
    136133
  • trunk/external/atnf/PKSIO/PKSMS2reader.cc

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2reader.cc,v 19.13 2008-06-26 02:08:02 cal103 Exp $
     28//# $Id: PKSMS2reader.cc,v 19.21 2008-11-17 06:55:18 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/03, Mark Calabretta, ATNF
    3131//#---------------------------------------------------------------------------
    3232
    33 
    34 // AIPS++ includes.
     33#include <atnf/pks/pks_maths.h>
     34#include <atnf/PKSIO/PKSmsg.h>
     35#include <atnf/PKSIO/PKSrecord.h>
     36#include <atnf/PKSIO/PKSMS2reader.h>
     37
    3538#include <casa/stdio.h>
    3639#include <casa/Arrays/ArrayMath.h>
     
    3942#include <tables/Tables.h>
    4043
    41 // Parkes includes.
    42 #include <atnf/pks/pks_maths.h>
    43 #include <atnf/PKSIO/PKSMS2reader.h>
    44 
    45 
    4644//------------------------------------------------- PKSMS2reader::PKSMS2reader
    4745
     
    5149{
    5250  cMSopen = False;
     51
     52  // By default, messages are written to stderr.
     53  initMsg();
    5354}
    5455
     
    353354        const Bool getSpectra,
    354355        const Bool getXPol,
    355         const Bool getFeedPos)
     356        const Int  coordSys)
    356357{
    357358  if (!cMSopen) {
     
    432433  cGetXPol = cGetXPol && getXPol;
    433434
    434   // Get feed positions?  (Not available.)
    435   cGetFeedPos = False;
     435  // Coordinate system?  (Only equatorial available.)
     436  cCoordSys = 0;
    436437
    437438  return maxNChan;
     
    486487// Read the next data record.
    487488
    488 Int PKSMS2reader::read(MBrecord &MBrec)
     489Int PKSMS2reader::read(PKSrecord &pksrec)
    489490{
    490491  if (!cMSopen) {
     
    514515
    515516  // Renumerate scan no. Here still is 1-based
    516   MBrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
    517 
    518   if (MBrec.scanNo != cScanNo) {
     517  pksrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
     518
     519  if (pksrec.scanNo != cScanNo) {
    519520    // Start of new scan.
    520     cScanNo  = MBrec.scanNo;
     521    cScanNo  = pksrec.scanNo;
    521522    cCycleNo = 1;
    522523    cTime    = cTimeCol(cIdx);
     
    524525
    525526  Double time = cTimeCol(cIdx);
    526   MBrec.mjd      = time/86400.0;
    527   MBrec.interval = cIntervalCol(cIdx);
     527  pksrec.mjd      = time/86400.0;
     528  pksrec.interval = cIntervalCol(cIdx);
    528529
    529530  // Reconstruct the integration cycle number; due to small latencies the
    530531  // integration time is usually slightly less than the time between cycles,
    531532  // resetting cTime will prevent the difference from accumulating.
    532   cCycleNo += nint((time - cTime)/MBrec.interval);
    533   MBrec.cycleNo = cCycleNo;
    534   cTime   = time;
     533  cCycleNo += nint((time - cTime)/pksrec.interval);
     534  pksrec.cycleNo = cCycleNo;
     535  cTime = time;
    535536
    536537  Int fieldId = cFieldIdCol(cIdx);
    537   MBrec.fieldName = cFieldNameCol(fieldId);
     538  pksrec.fieldName = cFieldNameCol(fieldId);
    538539
    539540  Int srcId = cSrcIdCol(fieldId);
    540   MBrec.srcName = cSrcNameCol(srcId);
    541   MBrec.srcDir  = cSrcDirCol(srcId);
    542   MBrec.srcPM   = cSrcPMCol(srcId);
     541  pksrec.srcName = cSrcNameCol(srcId);
     542  pksrec.srcDir  = cSrcDirCol(srcId);
     543  pksrec.srcPM   = cSrcPMCol(srcId);
    543544
    544545  // Systemic velocity.
    545546  if (!cHaveSrcVel) {
    546     MBrec.srcVel = 0.0f;
     547    pksrec.srcVel = 0.0f;
    547548  } else {
    548     MBrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
     549    pksrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
    549550  }
    550551
    551552  // Observation type.
    552553  Int stateId = cStateIdCol(cIdx);
    553   MBrec.obsType = cObsModeCol(stateId);
    554 
    555   MBrec.IFno = iIF + 1;
     554  pksrec.obsType = cObsModeCol(stateId);
     555
     556  pksrec.IFno = iIF + 1;
    556557  Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
    557558
     
    560561  if (nChan == 1) {
    561562    cout << "The input is continuum data. "<< endl;
    562     MBrec.freqInc  = chanFreq(0);
    563     MBrec.refFreq  = chanFreq(0);
    564     MBrec.restFreq = 0.0f;
     563    pksrec.freqInc  = chanFreq(0);
     564    pksrec.refFreq  = chanFreq(0);
     565    pksrec.restFreq = 0.0f;
    565566  } else {
    566567    if (cStartChan(iIF) <= cEndChan(iIF)) {
    567       MBrec.freqInc = chanFreq(1) - chanFreq(0);
     568      pksrec.freqInc = chanFreq(1) - chanFreq(0);
    568569    } else {
    569       MBrec.freqInc = chanFreq(0) - chanFreq(1);
    570     }
    571 
    572     MBrec.refFreq  = chanFreq(cRefChan(iIF)-1);
    573     MBrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
    574   }
    575   MBrec.bandwidth = abs(MBrec.freqInc * nChan);
    576 
    577   MBrec.tcal.resize(cNPol(iIF));
    578   MBrec.tcal      = 0.0f;
    579   MBrec.tcalTime  = "";
    580   MBrec.azimuth   = 0.0f;
    581   MBrec.elevation = 0.0f;
    582   MBrec.parAngle  = 0.0f;
    583   MBrec.focusAxi  = 0.0f;
    584   MBrec.focusTan  = 0.0f;
    585   MBrec.focusRot  = 0.0f;
     570      pksrec.freqInc = chanFreq(0) - chanFreq(1);
     571    }
     572
     573    pksrec.refFreq  = chanFreq(cRefChan(iIF)-1);
     574    pksrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
     575  }
     576  pksrec.bandwidth = abs(pksrec.freqInc * nChan);
     577
     578  pksrec.tcal.resize(cNPol(iIF));
     579  pksrec.tcal      = 0.0f;
     580  pksrec.tcalTime  = "";
     581  pksrec.azimuth   = 0.0f;
     582  pksrec.elevation = 0.0f;
     583  pksrec.parAngle  = 0.0f;
     584
     585  pksrec.focusAxi  = 0.0f;
     586  pksrec.focusTan  = 0.0f;
     587  pksrec.focusRot  = 0.0f;
    586588
    587589  // Find the appropriate entry in the WEATHER subtable.
     
    596598  if (weatherIdx < 0) {
    597599    // No appropriate WEATHER entry.
    598     MBrec.temperature = 0.0f;
    599     MBrec.pressure    = 0.0f;
    600     MBrec.humidity    = 0.0f;
     600    pksrec.temperature = 0.0f;
     601    pksrec.pressure    = 0.0f;
     602    pksrec.humidity    = 0.0f;
    601603  } else {
    602     MBrec.temperature = cTemperatureCol(weatherIdx);
    603     MBrec.pressure    = cPressureCol(weatherIdx);
    604     MBrec.humidity    = cHumidityCol(weatherIdx);
    605   }
    606 
    607   MBrec.windSpeed = 0.0f;
    608   MBrec.windAz    = 0.0f;
    609 
    610   MBrec.refBeam = 0;
    611   MBrec.beamNo  = ibeam + 1;
     604    pksrec.temperature = cTemperatureCol(weatherIdx);
     605    pksrec.pressure    = cPressureCol(weatherIdx);
     606    pksrec.humidity    = cHumidityCol(weatherIdx);
     607  }
     608
     609  pksrec.windSpeed = 0.0f;
     610  pksrec.windAz    = 0.0f;
     611
     612  pksrec.refBeam = 0;
     613  pksrec.beamNo  = ibeam + 1;
    612614
    613615  Matrix<Double> pointingDir = cPointingCol(fieldId);
    614   MBrec.direction = pointingDir.column(0);
     616  pksrec.direction = pointingDir.column(0);
     617  pksrec.pCode = 0;
     618  pksrec.rateAge = 0.0f;
    615619  uInt ncols = pointingDir.ncolumn();
    616620  if (ncols == 1) {
    617     MBrec.scanRate = 0.0f;
     621    pksrec.scanRate = 0.0f;
    618622  } else {
    619     MBrec.scanRate  = pointingDir.column(1);
    620   }
    621   MBrec.rateAge = 0;
    622   MBrec.rateson = 0;
     623    pksrec.scanRate(0) = pointingDir.column(1)(0);
     624    pksrec.scanRate(1) = pointingDir.column(1)(1);
     625  }
     626  pksrec.paRate = 0.0f;
    623627
    624628  // Get Tsys assuming that entries in the SYSCAL table match the main table.
     
    630634  }
    631635  if (cHaveTsys) {
    632     cTsysCol.get(cIdx, MBrec.tsys, True);
     636    cTsysCol.get(cIdx, pksrec.tsys, True);
    633637  } else {
    634638    Int numReceptor;
    635639    cNumReceptorCol.get(0, numReceptor);
    636     MBrec.tsys.resize(numReceptor);
    637     MBrec.tsys = 1.0f;
    638   }
    639   cSigmaCol.get(cIdx, MBrec.sigma, True);
     640    pksrec.tsys.resize(numReceptor);
     641    pksrec.tsys = 1.0f;
     642  }
     643  cSigmaCol.get(cIdx, pksrec.sigma, True);
    640644
    641645  // Calibration factors (if available).
    642   MBrec.calFctr.resize(cNPol(iIF));
     646  pksrec.calFctr.resize(cNPol(iIF));
    643647  if (cHaveCalFctr) {
    644     cCalFctrCol.get(cIdx, MBrec.calFctr);
     648    cCalFctrCol.get(cIdx, pksrec.calFctr);
    645649  } else {
    646     MBrec.calFctr = 0.0f;
     650    pksrec.calFctr = 0.0f;
    647651  }
    648652
    649653  // Baseline parameters (if available).
    650654  if (cHaveBaseLin) {
    651     MBrec.baseLin.resize(2,cNPol(iIF));
    652     cBaseLinCol.get(cIdx, MBrec.baseLin);
    653 
    654     MBrec.baseSub.resize(9,cNPol(iIF));
    655     cBaseSubCol.get(cIdx, MBrec.baseSub);
     655    pksrec.baseLin.resize(2,cNPol(iIF));
     656    cBaseLinCol.get(cIdx, pksrec.baseLin);
     657
     658    pksrec.baseSub.resize(9,cNPol(iIF));
     659    cBaseSubCol.get(cIdx, pksrec.baseSub);
    656660
    657661  } else {
    658     MBrec.baseLin.resize(0,0);
    659     MBrec.baseSub.resize(0,0);
     662    pksrec.baseLin.resize(0,0);
     663    pksrec.baseSub.resize(0,0);
    660664  }
    661665
     
    670674    // Transpose spectra.
    671675    Int nPol = tmpData.nrow();
    672     MBrec.spectra.resize(nChan, nPol);
    673     MBrec.flagged.resize(nChan, nPol);
     676    pksrec.spectra.resize(nChan, nPol);
     677    pksrec.flagged.resize(nChan, nPol);
    674678    if (cEndChan(iIF) >= cStartChan(iIF)) {
    675679      // Simple transposition.
    676680      for (Int ipol = 0; ipol < nPol; ipol++) {
    677681        for (Int ichan = 0; ichan < nChan; ichan++) {
    678           MBrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
    679           MBrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
     682          pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
     683          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
    680684        }
    681685      }
     
    686690      for (Int ipol = 0; ipol < nPol; ipol++) {
    687691        for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
    688           MBrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
    689           MBrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
     692          pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
     693          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
    690694        }
    691695      }
     
    696700  if (cGetXPol) {
    697701    if (cHaveXCalFctr) {
    698       cXCalFctrCol.get(cIdx, MBrec.xCalFctr);
     702      cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
    699703    } else {
    700       MBrec.xCalFctr = Complex(0.0f, 0.0f);
    701     }
    702 
    703     cDataCol.get(cIdx, MBrec.xPol, True);
     704      pksrec.xCalFctr = Complex(0.0f, 0.0f);
     705    }
     706
     707    cDataCol.get(cIdx, pksrec.xPol, True);
    704708
    705709    if (cEndChan(iIF) < cStartChan(iIF)) {
     
    707711      Int jchan = nChan - 1;
    708712      for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
    709         ctmp = MBrec.xPol(ichan);
    710         MBrec.xPol(ichan) = MBrec.xPol(jchan);
    711         MBrec.xPol(jchan) = ctmp;
    712       }
    713     }
    714   }
    715 
    716   cIdx++;
    717 
    718   return 0;
    719 }
    720 
    721 //--------------------------------------------------------- PKSMS2reader::read
    722 
    723 // Read the next data record, just the basics.
    724 
    725 Int PKSMS2reader::read(
    726         Int           &IFno,
    727         Vector<Float> &tsys,
    728         Vector<Float> &calFctr,
    729         Matrix<Float> &baseLin,
    730         Matrix<Float> &baseSub,
    731         Matrix<Float> &spectra,
    732         Matrix<uChar> &flagged)
    733 {
    734   if (!cMSopen) {
    735     return 1;
    736   }
    737 
    738   // Check for EOF.
    739   if (cIdx >= cNRow) {
    740     return -1;
    741   }
    742 
    743   // Find the next selected beam and IF.
    744   Int ibeam;
    745   Int iIF;
    746   while (True) {
    747     ibeam = cBeamNoCol(cIdx);
    748     iIF   = cDataDescIdCol(cIdx);
    749     if (cBeams(ibeam) && cIFs(iIF)) {
    750       break;
    751     }
    752 
    753     // Check for EOF.
    754     if (++cIdx >= cNRow) {
    755       return -1;
    756     }
    757   }
    758 
    759   IFno = iIF + 1;
    760 
    761   // Get Tsys assuming that entries in the SYSCAL table match the main table.
    762   cTsysCol.get(cIdx, tsys, True);
    763 
    764   // Calibration factors (if available).
    765   if (cHaveCalFctr) {
    766     cCalFctrCol.get(cIdx, calFctr, True);
    767   } else {
    768     calFctr.resize(cNPol(iIF));
    769     calFctr = 0.0f;
    770   }
    771 
    772   // Baseline parameters (if available).
    773   if (cHaveBaseLin) {
    774     baseLin.resize(2,cNPol(iIF));
    775     cBaseLinCol.get(cIdx, baseLin);
    776 
    777     baseSub.resize(9,cNPol(iIF));
    778     cBaseSubCol.get(cIdx, baseSub);
    779 
    780   } else {
    781     baseLin.resize(0,0);
    782     baseSub.resize(0,0);
    783   }
    784 
    785   if (cGetSpectra) {
    786     // Get spectral data.
    787     Matrix<Float> tmpData;
    788     Matrix<Bool>  tmpFlag;
    789     cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
    790     cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
    791 
    792     // Transpose spectra.
    793     Int nChan = tmpData.ncolumn();
    794     Int nPol  = tmpData.nrow();
    795     spectra.resize(nChan, nPol);
    796     flagged.resize(nChan, nPol);
    797     if (cEndChan(iIF) >= cStartChan(iIF)) {
    798       // Simple transposition.
    799       for (Int ipol = 0; ipol < nPol; ipol++) {
    800         for (Int ichan = 0; ichan < nChan; ichan++) {
    801           spectra(ichan,ipol) = tmpData(ipol,ichan);
    802           flagged(ichan,ipol) = tmpFlag(ipol,ichan);
    803         }
    804       }
    805 
    806     } else {
    807       // Transpose with inversion.
    808       Int jchan = nChan - 1;
    809       for (Int ipol = 0; ipol < nPol; ipol++) {
    810         for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
    811           spectra(ichan,ipol) = tmpData(ipol,jchan);
    812           flagged(ichan,ipol) = tmpFlag(ipol,jchan);
    813         }
     713        ctmp = pksrec.xPol(ichan);
     714        pksrec.xPol(ichan) = pksrec.xPol(jchan);
     715        pksrec.xPol(jchan) = ctmp;
    814716      }
    815717    }
  • trunk/external/atnf/PKSIO/PKSMS2reader.h

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2reader.h,v 19.14 2008-06-26 02:03:51 cal103 Exp $
     28//# $Id: PKSMS2reader.h,v 19.17 2008-11-17 06:38:53 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/03, Mark Calabretta, ATNF
     
    3535
    3636#include <atnf/PKSIO/PKSreader.h>
     37#include <atnf/PKSIO/PKSrecord.h>
    3738
    3839#include <casa/aips.h>
     
    4647#include <tables/Tables/ScalarColumn.h>
    4748
     49#include <casa/namespace.h>
     50
    4851// <summary>
    4952// Class to read Parkes Multibeam data from a v2 MS.
    5053// </summary>
    51 
    52 #include <casa/namespace.h>
    5354
    5455class PKSMS2reader : public PKSreader
     
    101102        const Bool getSpectra = True,
    102103        const Bool getXPol    = False,
    103         const Bool getFeedPos = False);
     104        const Int  coordSys   = 0);
    104105
    105106    // Find the range of the data selected in time and position.
     
    111112
    112113    // Read the next data record.
    113     virtual Int read(MBrecord &mbrec);
    114 
    115     // Read the next data record, just the basics.
    116     virtual Int read(
    117         Int           &IFno,
    118         Vector<Float> &tsys,
    119         Vector<Float> &calFctr,
    120         Matrix<Float> &baseLin,
    121         Matrix<Float> &baseSub,
    122         Matrix<Float> &spectra,
    123         Matrix<uChar> &flagged);
     114    virtual Int read(PKSrecord &pksrec);
    124115
    125116    // Close the MS.
  • trunk/external/atnf/PKSIO/PKSMS2writer.cc

    r1399 r1452  
    22//# PKSMS2writer.cc: Class to write Parkes multibeam data to a measurementset.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2007
     4//# Copyright (C) 2000-2008
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2writer.cc,v 19.12 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: PKSMS2writer.cc,v 19.14 2008-11-17 06:56:13 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
     31#include <atnf/PKSIO/PKSrecord.h>
    3132#include <atnf/PKSIO/PKSMS2writer.h>
    3233
     
    5455PKSMS2writer::PKSMS2writer()
    5556{
     57  cPKSMS = 0x0;
     58
     59  // By default, messages are written to stderr.
     60  initMsg();
    5661}
    5762
     
    8489        const Bool   haveBase)
    8590{
     91  if (cPKSMS) {
     92    logMsg("ERROR: Output MS already open, close it first.");
     93    return 1;
     94  }
     95
     96  // Clear the message stack.
     97  clearMsg();
     98
    8699  // Open a MS table.
    87100  TableDesc pksDesc = MS::requiredTableDesc();
     
    369382
    370383Int PKSMS2writer::write(
    371         const Int             scanNo,
    372         const Int             cycleNo,
    373         const Double          mjd,
    374         const Double          interval,
    375         const String          fieldName,
    376         const String          srcName,
    377         const Vector<Double>  srcDir,
    378         const Vector<Double>  srcPM,
    379         const Double          srcVel,
    380         const String          obsMode,
    381         const Int             IFno,
    382         const Double          refFreq,
    383         const Double          bandwidth,
    384         const Double          freqInc,
    385         const Double          restFreq,
    386         const Vector<Float>   tcal,
    387         const String          tcalTime,
    388         const Float           azimuth,
    389         const Float           elevation,
    390         const Float           parAngle,
    391         const Float           focusAxi,
    392         const Float           focusTan,
    393         const Float           focusRot,
    394         const Float           temperature,
    395         const Float           pressure,
    396         const Float           humidity,
    397         const Float           windSpeed,
    398         const Float           windAz,
    399         const Int             refBeam,
    400         const Int             beamNo,
    401         const Vector<Double>  direction,
    402         const Vector<Double>  scanRate,
    403         const Vector<Float>   tsys,
    404         const Vector<Float>   sigma,
    405         const Vector<Float>   calFctr,
    406         const Matrix<Float>   baseLin,
    407         const Matrix<Float>   baseSub,
    408         const Matrix<Float>   &spectra,
    409         const Matrix<uChar>   &flagged,
    410         const Complex         xCalFctr,
    411         const Vector<Complex> &xPol)
     384        const PKSrecord &pksrec)
    412385{
    413386  // Extend the time range in the OBSERVATION subtable.
    414387  Vector<Double> timerange(2);
    415388  cObservationCols->timeRange().get(0, timerange);
    416   Double time = mjd*86400.0;
     389  Double time = pksrec.mjd*86400.0;
    417390  if (timerange(0) == 0.0) {
    418391    timerange(0) = time;
     
    421394  cObservationCols->timeRange().put(0, timerange);
    422395
    423   Int iIF = IFno - 1;
     396  Int iIF = pksrec.IFno - 1;
    424397  Int nChan = cNChan(iIF);
    425398  Int nPol  = cNPol(iIF);
     
    427400  // IFno is the 1-relative row number in the DATA_DESCRIPTION,
    428401  // SPECTRAL_WINDOW, and POLARIZATION subtables.
    429   if (Int(cDataDescription.nrow()) < IFno) {
     402  if (Int(cDataDescription.nrow()) < pksrec.IFno) {
    430403    // Add a new entry to each subtable.
    431     addDataDescriptionEntry(IFno);
    432     addSpectralWindowEntry(IFno, nChan, refFreq, bandwidth, freqInc);
    433     addPolarizationEntry(IFno, nPol);
     404    addDataDescriptionEntry(pksrec.IFno);
     405    addSpectralWindowEntry(pksrec.IFno, nChan, pksrec.refFreq,
     406      pksrec.bandwidth, pksrec.freqInc);
     407    addPolarizationEntry(pksrec.IFno, nPol);
    434408  }
    435409
    436410  // Find or add the source to the SOURCE subtable.
    437   Int srcId = addSourceEntry(srcName, srcDir, srcPM, restFreq, srcVel);
    438 
    439   // Find or add the obsMode to the STATE subtable.
    440   Int stateId = addStateEntry(obsMode);
     411  Int srcId = addSourceEntry(pksrec.srcName, pksrec.srcDir, pksrec.srcPM,
     412    pksrec.restFreq, pksrec.srcVel);
     413
     414  // Find or add the obsType to the STATE subtable.
     415  Int stateId = addStateEntry(pksrec.obsType);
    441416
    442417  // FIELD subtable.
    443   Int fieldId = addFieldEntry(fieldName, time, direction, scanRate, srcId);
     418  Vector<Double> scanRate(2);
     419  scanRate(0) = pksrec.scanRate(0);
     420  scanRate(1) = pksrec.scanRate(1);
     421  Int fieldId = addFieldEntry(pksrec.fieldName, time, pksrec.direction,
     422    scanRate, srcId);
    444423
    445424  // POINTING subtable.
    446   addPointingEntry(time, interval, fieldName, direction, scanRate);
     425  addPointingEntry(time, pksrec.interval, pksrec.fieldName, pksrec.direction,
     426    scanRate);
    447427
    448428  // SYSCAL subtable.
    449   addSysCalEntry(beamNo, iIF, time, interval, tcal, tsys);
     429  addSysCalEntry(pksrec.beamNo, iIF, time, pksrec.interval, pksrec.tcal,
     430    pksrec.tsys);
    450431
    451432  // Handle weather information.
     
    453434  Int nWeather = wTime.nrow();
    454435  if (nWeather == 0 || time > wTime(nWeather-1)) {
    455     addWeatherEntry(time, interval, pressure, humidity, temperature);
     436    addWeatherEntry(time, pksrec.interval, pksrec.pressure, pksrec.humidity,
     437      pksrec.temperature);
    456438  }
    457439
     
    465447  cMSCols->antenna1().put(irow, 0);
    466448  cMSCols->antenna2().put(irow, 0);
    467   cMSCols->feed1().put(irow, beamNo-1);
    468   cMSCols->feed2().put(irow, beamNo-1);
     449  cMSCols->feed1().put(irow, pksrec.beamNo-1);
     450  cMSCols->feed2().put(irow, pksrec.beamNo-1);
    469451  cMSCols->dataDescId().put(irow, iIF);
    470452  cMSCols->processorId().put(irow, 0);
     
    472454
    473455  // Non-key attributes.
    474   cMSCols->interval().put(irow, interval);
    475   cMSCols->exposure().put(irow, interval);
     456  cMSCols->interval().put(irow, pksrec.interval);
     457  cMSCols->exposure().put(irow, pksrec.interval);
    476458  cMSCols->timeCentroid().put(irow, time);
    477   cMSCols->scanNumber().put(irow, scanNo);
     459  cMSCols->scanNumber().put(irow, pksrec.scanNo);
    478460  cMSCols->arrayId().put(irow, 0);
    479461  cMSCols->observationId().put(irow, 0);
     
    485467  // Baseline fit parameters.
    486468  if (cHaveBase) {
    487     cBaseLinCol->put(irow, baseLin);
    488 
    489     if (baseSub.nrow() == 9) {
    490       cBaseSubCol->put(irow, baseSub);
     469    cBaseLinCol->put(irow, pksrec.baseLin);
     470
     471    if (pksrec.baseSub.nrow() == 9) {
     472      cBaseSubCol->put(irow, pksrec.baseSub);
    491473
    492474    } else {
    493475      Matrix<Float> tmp(9, 2, 0.0f);
    494476      for (Int ipol = 0; ipol < nPol; ipol++) {
    495         for (uInt j = 0; j < baseSub.nrow(); j++) {
    496           tmp(j,ipol) = baseSub(j,ipol);
     477        for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) {
     478          tmp(j,ipol) = pksrec.baseSub(j,ipol);
    497479        }
    498480      }
     
    506488  for (Int ipol = 0; ipol < nPol; ipol++) {
    507489    for (Int ichan = 0; ichan < nChan; ichan++) {
    508       tmpData(ipol,ichan) = spectra(ichan,ipol);
    509       tmpFlag(ipol,ichan) = flagged(ichan,ipol);
     490      tmpData(ipol,ichan) = pksrec.spectra(ichan,ipol);
     491      tmpFlag(ipol,ichan) = pksrec.flagged(ichan,ipol);
    510492    }
    511493  }
    512494
    513   cCalFctrCol->put(irow, calFctr);
     495  cCalFctrCol->put(irow, pksrec.calFctr);
    514496  cMSCols->floatData().put(irow, tmpData);
    515497  cMSCols->flag().put(irow, tmpFlag);
     
    517499  // Cross-polarization spectra.
    518500  if (cHaveXPol(iIF)) {
    519     cXCalFctrCol->put(irow, xCalFctr);
    520     cMSCols->data().put(irow, xPol);
    521   }
    522 
    523   cMSCols->sigma().put(irow, sigma);
     501    cXCalFctrCol->put(irow, pksrec.xCalFctr);
     502    cMSCols->data().put(irow, pksrec.xPol);
     503  }
     504
     505  cMSCols->sigma().put(irow, pksrec.sigma);
    524506
    525507  Vector<Float> weight(1, 1.0f);
     
    586568  cSysCal          = MSSysCal();
    587569  cWeather         = MSWeather();
     570
    588571  // Release the main table.
    589   delete cPKSMS; cPKSMS=0;
     572  delete cPKSMS;
     573  cPKSMS = 0x0;
    590574}
    591575
     
    984968
    985969Int PKSMS2writer::addStateEntry(
    986         const String obsMode)
     970        const String obsType)
    987971{
    988972  // Look for an entry in the STATE subtable.
    989973  for (uInt n = 0; n < cStateCols->nrow(); n++) {
    990     if (cStateCols->obsMode()(n) == obsMode) {
     974    if (cStateCols->obsMode()(n) == obsType) {
    991975      return n;
    992976    }
     
    998982
    999983  // Data.
    1000   if (obsMode.contains("RF")) {
     984  if (obsType.contains("RF")) {
    1001985    cStateCols->sig().put(n, False);
    1002986    cStateCols->ref().put(n, True);
    1003   } else if (!obsMode.contains("PA")) {
     987  } else if (!obsType.contains("PA")) {
    1004988    // Signal and reference are both false for "paddle" data.
    1005989    cStateCols->sig().put(n, True);
     
    1010994  cStateCols->cal().put(n, 0.0);
    1011995  cStateCols->subScan().put(n, 0);
    1012   cStateCols->obsMode().put(n, obsMode);
     996  cStateCols->obsMode().put(n, obsType);
    1013997
    1014998  // Flags.
  • trunk/external/atnf/PKSIO/PKSMS2writer.h

    r1399 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2writer.h,v 19.12 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: PKSMS2writer.h,v 19.13 2008-11-17 06:40:25 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
     
    3232#define ATNF_PKSMS2WRITER_H
    3333
     34#include <atnf/PKSIO/PKSrecord.h>
    3435#include <atnf/PKSIO/PKSwriter.h>
    3536
     
    3839#include <casa/Arrays/Vector.h>
    3940#include <casa/BasicSL/Complex.h>
     41#include <casa/BasicSL/String.h>
    4042#include <ms/MeasurementSets/MeasurementSet.h>
    4143#include <ms/MeasurementSets/MSColumns.h>
    42 #include <casa/BasicSL/String.h>
     44
     45#include <casa/namespace.h>
    4346
    4447// <summary>
     
    4649// </summary>
    4750
    48 #include <casa/namespace.h>
    4951class PKSMS2writer : public PKSwriter
    5052{
     
    7476    // Write the next data record.
    7577    virtual Int write(
    76         const Int             scanNo,
    77         const Int             cycleNo,
    78         const Double          mjd,
    79         const Double          interval,
    80         const String          fieldName,
    81         const String          srcName,
    82         const Vector<Double>  srcDir,
    83         const Vector<Double>  srcPM,
    84         const Double          srcVel,
    85         const String          obsMode,
    86         const Int             IFno,
    87         const Double          refFreq,
    88         const Double          bandwidth,
    89         const Double          freqInc,
    90         const Double          restFreq,
    91         const Vector<Float>   tcal,
    92         const String          tcalTime,
    93         const Float           azimuth,
    94         const Float           elevation,
    95         const Float           parAngle,
    96         const Float           focusAxi,
    97         const Float           focusTan,
    98         const Float           focusRot,
    99         const Float           temperature,
    100         const Float           pressure,
    101         const Float           humidity,
    102         const Float           windSpeed,
    103         const Float           windAz,
    104         const Int             refBeam,
    105         const Int             beamNo,
    106         const Vector<Double>  direction,
    107         const Vector<Double>  scanRate,
    108         const Vector<Float>   tsys,
    109         const Vector<Float>   sigma,
    110         const Vector<Float>   calFctr,
    111         const Matrix<Float>   baseLin,
    112         const Matrix<Float>   baseSub,
    113         const Matrix<Float>   &spectra,
    114         const Matrix<uChar>   &flagged,
    115         const Complex         xCalFctr,
    116         const Vector<Complex> &xPol);
     78        const PKSrecord &pksrec);
    11779
    11880    // Close the MS, flushing all associated Tables.
  • trunk/external/atnf/PKSIO/PKSSDwriter.cc

    r1399 r1452  
    22//# PKSSDwriter.cc: Class to write Parkes multibeam data to an SDFITS file.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2007
     4//# Copyright (C) 2000-2008
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSSDwriter.cc,v 19.13 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: PKSSDwriter.cc,v 19.15 2008-11-17 06:56:50 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
    31 #include <atnf/PKSIO/PKSMBrecord.h>
     31#include <atnf/PKSIO/MBrecord.h>
    3232#include <atnf/PKSIO/PKSSDwriter.h>
    3333
     34#include <casa/stdio.h>
    3435#include <casa/Quanta/MVTime.h>
    3536
    36 
    3737//--------------------------------------------------- PKSSDwriter::PKSSDwriter
    3838
     
    4141PKSSDwriter::PKSSDwriter()
    4242{
     43  // By default, messages are written to stderr.
     44  initMsg();
    4345}
    4446
     
    5052{
    5153  close();
     54}
     55
     56//-------------------------------------------------------- PKSSDwriter::setMsg
     57
     58// Set message disposition.  If fd is non-zero messages will be written
     59// to that file descriptor, else stored for retrieval by getMsg().
     60
     61Int PKSSDwriter::setMsg(FILE *fd)
     62{
     63  PKSmsg::setMsg(fd);
     64  cSDwriter.setMsg(fd);
     65
     66  return 0;
    5267}
    5368
     
    7186        const Bool   haveBase)
    7287{
     88  // Clear the message stack.
     89  clearMsg();
     90
    7391  double antPos[3];
    7492  antPos[0] = antPosition(0);
     
    107125        (int *)cNPol.getStorage(deleteIt),
    108126        (int *)cHaveXPol.getStorage(deleteIt), (int)cHaveBase, 1);
     127  logMsg(cSDwriter.getMsg());
     128  cSDwriter.clearMsg();
    109129  if (status) {
    110     cSDwriter.reportError();
    111130    cSDwriter.deleteFile();
    112131    close();
     
    121140
    122141Int PKSSDwriter::write(
    123         const Int             scanNo,
    124         const Int             cycleNo,
    125         const Double          mjd,
    126         const Double          interval,
    127         const String          fieldName,
    128         const String          srcName,
    129         const Vector<Double>  srcDir,
    130         const Vector<Double>  srcPM,
    131         const Double          srcVel,
    132         const String          obsMode,
    133         const Int             IFno,
    134         const Double          refFreq,
    135         const Double          bandwidth,
    136         const Double          freqInc,
    137         const Double          restFreq,
    138         const Vector<Float>   tcal,
    139         const String          tcalTime,
    140         const Float           azimuth,
    141         const Float           elevation,
    142         const Float           parAngle,
    143         const Float           focusAxi,
    144         const Float           focusTan,
    145         const Float           focusRot,
    146         const Float           temperature,
    147         const Float           pressure,
    148         const Float           humidity,
    149         const Float           windSpeed,
    150         const Float           windAz,
    151         const Int             refBeam,
    152         const Int             beamNo,
    153         const Vector<Double>  direction,
    154         const Vector<Double>  scanRate,
    155         const Vector<Float>   tsys,
    156         const Vector<Float>   sigma,
    157         const Vector<Float>   calFctr,
    158         const Matrix<Float>   baseLin,
    159         const Matrix<Float>   baseSub,
    160         const Matrix<Float>   &spectra,
    161         const Matrix<uChar>   &flagged,
    162         const Complex         xCalFctr,
    163         const Vector<Complex> &xPol)
     142        const PKSrecord &pksrec)
    164143{
    165144  // Do basic checks.
     145  Int IFno = pksrec.IFno;
    166146  uInt iIF = IFno - 1;
    167147  if (IFno < 1 || Int(cNIF) < IFno) {
     
    172152  }
    173153
    174   uInt nChan = spectra.nrow();
     154  uInt nChan = pksrec.spectra.nrow();
    175155  if (nChan != cNChan(iIF)) {
    176156    cerr << "PKSDwriter::write: "
     
    181161  }
    182162
    183   uInt nPol = spectra.ncolumn();
     163  uInt nPol = pksrec.spectra.ncolumn();
    184164  if (nPol != cNPol(iIF)) {
    185165    cerr << "PKSDwriter::write: "
     
    191171
    192172  // Extract calendar information from mjd.
    193   MVTime time(mjd);
     173  MVTime time(pksrec.mjd);
    194174  Int year  = time.year();
    195175  Int month = time.month();
    196176  Int day   = time.monthday();
    197177
    198   // Transfer data to a single-IF PKSMBrecord.
    199   PKSMBrecord mbrec(1);
     178  // Transfer data to a single-IF MBrecord.
     179  MBrecord mbrec(1);
    200180
    201181  // Start with basic beam- and IF-independent bookkeeping information.
    202   mbrec.scanNo  = scanNo;
    203   mbrec.cycleNo = cycleNo;
     182  mbrec.scanNo  = pksrec.scanNo;
     183  mbrec.cycleNo = pksrec.cycleNo;
    204184
    205185  sprintf(mbrec.datobs, "%4.4d-%2.2d-%2.2d", year, month, day);
    206   mbrec.utc      = fmod(mjd, 1.0) * 86400.0;
    207   mbrec.exposure = float(interval);
    208 
    209   strncpy(mbrec.srcName, (char *)srcName.chars(), 17);
    210   mbrec.srcRA    = srcDir(0);
    211   mbrec.srcDec   = srcDir(1);
    212 
    213   mbrec.restFreq = restFreq;
    214 
    215   strncpy(mbrec.obsType, (char *)obsMode.chars(), 16);
     186  mbrec.utc      = fmod(pksrec.mjd, 1.0) * 86400.0;
     187  mbrec.exposure = float(pksrec.interval);
     188
     189  strncpy(mbrec.srcName, (char *)pksrec.srcName.chars(), 17);
     190  mbrec.srcRA    = pksrec.srcDir(0);
     191  mbrec.srcDec   = pksrec.srcDir(1);
     192
     193  mbrec.restFreq = pksrec.restFreq;
     194
     195  strncpy(mbrec.obsType, (char *)pksrec.obsType.chars(), 16);
    216196
    217197  // Now beam-dependent parameters.
    218   mbrec.beamNo   = beamNo;
    219   mbrec.ra       = direction(0);
    220   mbrec.dec      = direction(1);
    221   mbrec.raRate   = scanRate(0);
    222   mbrec.decRate  = scanRate(1);
     198  mbrec.beamNo   = pksrec.beamNo;
     199  mbrec.ra       = pksrec.direction(0);
     200  mbrec.dec      = pksrec.direction(1);
     201  mbrec.raRate   = pksrec.scanRate(0);
     202  mbrec.decRate  = pksrec.scanRate(1);
    223203
    224204  // Now IF-dependent parameters.
     
    229209
    230210  mbrec.fqRefPix[0] = (nChan/2) + 1;
    231   mbrec.fqRefVal[0] = refFreq;
    232   mbrec.fqDelt[0]   = freqInc;
     211  mbrec.fqRefVal[0] = pksrec.refFreq;
     212  mbrec.fqDelt[0]   = pksrec.freqInc;
    233213
    234214  // Now the data itself.
    235   for (uInt i = 0; i < tsys.nelements(); i++) {
    236     mbrec.tsys[0][i] = tsys(i);
     215  for (uInt i = 0; i < pksrec.tsys.nelements(); i++) {
     216    mbrec.tsys[0][i] = pksrec.tsys(i);
    237217  }
    238218
    239219  for (uInt ipol = 0; ipol < nPol; ipol++) {
    240     mbrec.calfctr[0][ipol] = calFctr(ipol);
     220    mbrec.calfctr[0][ipol] = pksrec.calFctr(ipol);
    241221  }
    242222
    243223  if (cHaveXPol(iIF)) {
    244     mbrec.xcalfctr[0][0] = xCalFctr.real();
    245     mbrec.xcalfctr[0][1] = xCalFctr.imag();
     224    mbrec.xcalfctr[0][0] = pksrec.xCalFctr.real();
     225    mbrec.xcalfctr[0][1] = pksrec.xCalFctr.imag();
    246226  } else {
    247227    mbrec.xcalfctr[0][0] = 0.0f;
     
    253233
    254234    for (uInt ipol = 0; ipol < nPol; ipol++) {
    255       mbrec.baseLin[0][ipol][0] = baseLin(0,ipol);
    256       mbrec.baseLin[0][ipol][1] = baseLin(1,ipol);
    257 
    258       for (uInt j = 0; j < baseSub.nrow(); j++) {
    259         mbrec.baseSub[0][ipol][j] = baseSub(j,ipol);
     235      mbrec.baseLin[0][ipol][0] = pksrec.baseLin(0,ipol);
     236      mbrec.baseLin[0][ipol][1] = pksrec.baseLin(1,ipol);
     237
     238      for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) {
     239        mbrec.baseSub[0][ipol][j] = pksrec.baseSub(j,ipol);
    260240      }
    261       for (uInt j = baseSub.nrow(); j < 9; j++) {
     241      for (uInt j = pksrec.baseSub.nrow(); j < 9; j++) {
    262242        mbrec.baseSub[0][ipol][j] = 0.0f;
    263243      }
     
    269249
    270250  Bool delSpectra = False;
    271   const Float *specstor = spectra.getStorage(delSpectra);
     251  const Float *specstor = pksrec.spectra.getStorage(delSpectra);
    272252  mbrec.spectra[0] = (float *)specstor;
    273253
    274254  Bool delFlagged = False;
    275   const uChar *flagstor = flagged.getStorage(delFlagged);
     255  const uChar *flagstor = pksrec.flagged.getStorage(delFlagged);
    276256  mbrec.flagged[0] = (unsigned char *)flagstor;
    277257
     
    279259  const Complex *xpolstor;
    280260  if (cHaveXPol(iIF)) {
    281     xpolstor = xPol.getStorage(delXPol);
     261    xpolstor = pksrec.xPol.getStorage(delXPol);
    282262  } else {
    283263    xpolstor = 0;
     
    287267  // Finish off with system calibration parameters.
    288268  mbrec.extraSysCal = 1;
    289   mbrec.refBeam     = refBeam;
    290   for (uInt i = 0; i < tcal.nelements(); i++) {
    291     mbrec.tcal[0][i] = tcal(i);
    292   }
    293   strncpy(mbrec.tcalTime, (char *)tcalTime.chars(), 16);
    294   mbrec.azimuth   = azimuth;
    295   mbrec.elevation = elevation;
    296   mbrec.parAngle  = parAngle;
    297   mbrec.focusAxi  = focusAxi;
    298   mbrec.focusTan  = focusTan;
    299   mbrec.focusRot  = focusRot;
    300   mbrec.temp      = temperature;
    301   mbrec.pressure  = pressure;
    302   mbrec.humidity  = humidity;
    303   mbrec.windSpeed = windSpeed;
    304   mbrec.windAz    = windAz;
     269  mbrec.refBeam     = pksrec.refBeam;
     270  for (uInt i = 0; i < pksrec.tcal.nelements(); i++) {
     271    mbrec.tcal[0][i] = pksrec.tcal(i);
     272  }
     273  strncpy(mbrec.tcalTime, (char *)pksrec.tcalTime.chars(), 16);
     274  mbrec.azimuth   = pksrec.azimuth;
     275  mbrec.elevation = pksrec.elevation;
     276  mbrec.parAngle  = pksrec.parAngle;
     277  mbrec.focusAxi  = pksrec.focusAxi;
     278  mbrec.focusTan  = pksrec.focusTan;
     279  mbrec.focusRot  = pksrec.focusRot;
     280  mbrec.temp      = pksrec.temperature;
     281  mbrec.pressure  = pksrec.pressure;
     282  mbrec.humidity  = pksrec.humidity;
     283  mbrec.windSpeed = pksrec.windSpeed;
     284  mbrec.windAz    = pksrec.windAz;
    305285
    306286  Int status = cSDwriter.write(mbrec);
     287  logMsg(cSDwriter.getMsg());
     288  cSDwriter.clearMsg();
    307289  if (status) {
    308     cSDwriter.reportError();
    309290    status = 1;
    310291  }
    311292
    312   spectra.freeStorage(specstor, delSpectra);
    313   flagged.freeStorage(flagstor, delFlagged);
    314   xPol.freeStorage(xpolstor, delXPol);
     293  pksrec.spectra.freeStorage(specstor, delSpectra);
     294  pksrec.flagged.freeStorage(flagstor, delFlagged);
     295  pksrec.xPol.freeStorage(xpolstor, delXPol);
    315296
    316297  return status;
     
    338319{
    339320  cSDwriter.close();
    340 }
     321  logMsg(cSDwriter.getMsg());
     322  cSDwriter.clearMsg();
     323}
  • trunk/external/atnf/PKSIO/PKSSDwriter.h

    r1399 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSSDwriter.h,v 19.14 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: PKSSDwriter.h,v 19.16 2008-11-17 06:41:48 cal103 Exp $
    2929//# Original: 2000/07/21, Mark Calabretta, ATNF
    3030//#---------------------------------------------------------------------------
     
    3434
    3535#include <atnf/PKSIO/PKSwriter.h>
     36#include <atnf/PKSIO/PKSrecord.h>
    3637#include <atnf/PKSIO/SDFITSwriter.h>
    3738
    3839#include <casa/aips.h>
     40#include <casa/stdio.h>
    3941#include <casa/Arrays/Vector.h>
    4042#include <casa/Arrays/Matrix.h>
     
    4244#include <casa/BasicSL/String.h>
    4345
     46#include <casa/namespace.h>
     47
    4448// <summary>
    4549// Class to write Parkes multibeam data to an SDFITS file.
    4650// </summary>
    4751
    48 #include <casa/namespace.h>
    4952class PKSSDwriter : public PKSwriter
    5053{
     
    5558    // Destructor.
    5659    virtual ~PKSSDwriter();
     60
     61    // Set message disposition.
     62    virtual Int setMsg(
     63        FILE *fd = 0x0);
    5764
    5865    // Create the SDFITS file and write static data.
     
    7481    // Write the next data record.
    7582    virtual Int write(
    76         const Int             scanNo,
    77         const Int             cycleNo,
    78         const Double          mjd,
    79         const Double          interval,
    80         const String          fieldName,
    81         const String          srcName,
    82         const Vector<Double>  srcDir,
    83         const Vector<Double>  srcPM,
    84         const Double          srcVel,
    85         const String          obsMode,
    86         const Int             IFno,
    87         const Double          refFreq,
    88         const Double          bandwidth,
    89         const Double          freqInc,
    90         const Double          restFreq,
    91         const Vector<Float>   tcal,
    92         const String          tcalTime,
    93         const Float           azimuth,
    94         const Float           elevation,
    95         const Float           parAngle,
    96         const Float           focusAxi,
    97         const Float           focusTan,
    98         const Float           focusRot,
    99         const Float           temperature,
    100         const Float           pressure,
    101         const Float           humidity,
    102         const Float           windSpeed,
    103         const Float           windAz,
    104         const Int             refBeam,
    105         const Int             beamNo,
    106         const Vector<Double>  direction,
    107         const Vector<Double>  scanRate,
    108         const Vector<Float>   tsys,
    109         const Vector<Float>   sigma,
    110         const Vector<Float>   calFctr,
    111         const Matrix<Float>   baselin,
    112         const Matrix<Float>   basesub,
    113         const Matrix<Float>   &spectra,
    114         const Matrix<uChar>   &flagged,
    115         const Complex         xCalFctr,
    116         const Vector<Complex> &xPol);
     83        const PKSrecord &pksrec);
    11784
    11885    // Write a history record.
  • trunk/external/atnf/PKSIO/PKSreader.cc

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSreader.cc,v 19.6 2008-06-26 01:54:08 cal103 Exp $
     28//# $Id: PKSreader.cc,v 19.12 2008-11-17 06:58:07 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/23, Mark Calabretta, ATNF
     
    4141#include <casa/OS/File.h>
    4242
    43 
    4443//--------------------------------------------------------------- getPKSreader
    4544
     
    5049        const Int retry,
    5150        const Int interpolate,
    52         String &format,
    53         Vector<Bool> &beams,
    54         Vector<Bool> &IFs,
    55         Vector<uInt> &nChan,
    56         Vector<uInt> &nPol,
    57         Vector<Bool> &haveXPol,
    58         Bool   &haveBase,
    59         Bool   &haveSpectra)
     51        String &format)
    6052{
    6153  // Check accessibility of the input.
     
    6355  if (!inFile.exists()) {
    6456    format = "DATASET NOT FOUND";
    65     return 0;
     57    return 0x0;
    6658  }
    6759
    6860  if (!inFile.isReadable()) {
    6961    format = "DATASET UNREADABLE";
    70     return 0;
     62    return 0x0;
    7163  }
    7264
    7365  // Determine the type of input.
    74   PKSreader *reader = 0;
     66  PKSreader *reader = 0x0;
    7567  if (inFile.isRegular()) {
    7668    // Is it MBFITS or SDFITS?
     
    112104  }
    113105
     106  return reader;
     107}
     108
     109//--------------------------------------------------------------- getPKSreader
     110
     111// Search a list of directories for a Parkes Multibeam dataset and return an
     112// appropriate PKSreader for it.
     113
     114PKSreader* getPKSreader(
     115        const String name,
     116        const Vector<String> directories,
     117        const Int retry,
     118        const Int interpolate,
     119        Int    &iDir,
     120        String &format)
     121{
     122  PKSreader *reader = 0x0;
     123
     124  iDir = -1;
     125  Int nDir = directories.nelements();
     126  for (Int i = 0; i < nDir; i++) {
     127    String inName = directories(i) + "/" + name;
     128    reader = getPKSreader(inName, retry, interpolate, format);
     129    if (reader) {
     130      iDir = i;
     131      break;
     132    }
     133  }
     134
     135  return reader;
     136}
     137
     138//--------------------------------------------------------------- getPKSreader
     139
     140// Open an appropriate PKSreader for a Parkes Multibeam dataset.
     141
     142PKSreader* getPKSreader(
     143        const String name,
     144        const Int retry,
     145        const Int interpolate,
     146        String &format,
     147        Vector<Bool> &beams,
     148        Vector<Bool> &IFs,
     149        Vector<uInt> &nChan,
     150        Vector<uInt> &nPol,
     151        Vector<Bool> &haveXPol,
     152        Bool   &haveBase,
     153        Bool   &haveSpectra)
     154{
     155  PKSreader *reader = getPKSreader(name, retry, interpolate, format);
    114156
    115157  // Try to open it.
     
    119161      format += " OPEN ERROR";
    120162      delete reader;
    121     } else {
    122       return reader;
    123     }
    124   }
    125 
    126   return 0;
    127 }
    128 
    129 
    130 //--------------------------------------------------------------- getPKSreader
    131 
    132 // Search a list of directories for a Parkes Multibeam dataset and return an
     163      reader = 0x0;
     164    }
     165  }
     166
     167  return reader;
     168}
     169
     170//--------------------------------------------------------------- getPKSreader
     171
     172// Search a list of directories for a Parkes Multibeam dataset and open an
    133173// appropriate PKSreader for it.
    134174
     
    148188        Bool   &haveSpectra)
    149189{
    150   Int nDir = directories.nelements();
    151   for (iDir = 0; iDir < nDir; iDir++) {
    152     String inName = directories(iDir) + "/" + name;
    153     PKSreader *reader = getPKSreader(inName, retry, interpolate, format,
    154                                      beams, IFs, nChan, nPol, haveXPol,
    155                                      haveBase, haveSpectra);
    156     if (reader != 0) {
    157       return reader;
    158     }
    159   }
    160 
    161   iDir = -1;
    162   return 0;
    163 }
    164 
    165 
    166 //-------------------------------------------------------- PKSFITSreader::read
    167 
    168 // Read the next data record.
    169 
    170 Int PKSreader::read(
    171         Int             &scanNo,
    172         Int             &cycleNo,
    173         Double          &mjd,
    174         Double          &interval,
    175         String          &fieldName,
    176         String          &srcName,
    177         Vector<Double>  &srcDir,
    178         Vector<Double>  &srcPM,
    179         Double          &srcVel,
    180         String          &obsType,
    181         Int             &IFno,
    182         Double          &refFreq,
    183         Double          &bandwidth,
    184         Double          &freqInc,
    185         Double          &restFreq,
    186         Vector<Float>   &tcal,
    187         String          &tcalTime,
    188         Float           &azimuth,
    189         Float           &elevation,
    190         Float           &parAngle,
    191         Float           &focusAxi,
    192         Float           &focusTan,
    193         Float           &focusRot,
    194         Float           &temperature,
    195         Float           &pressure,
    196         Float           &humidity,
    197         Float           &windSpeed,
    198         Float           &windAz,
    199         Int             &refBeam,
    200         Int             &beamNo,
    201         Vector<Double>  &direction,
    202         Vector<Double>  &scanRate,
    203         Vector<Float>   &tsys,
    204         Vector<Float>   &sigma,
    205         Vector<Float>   &calFctr,
    206         Matrix<Float>   &baseLin,
    207         Matrix<Float>   &baseSub,
    208         Matrix<Float>   &spectra,
    209         Matrix<uChar>   &flagged,
    210         Complex         &xCalFctr,
    211         Vector<Complex> &xPol)
    212 {
    213   Int status;
    214   MBrecord MBrec;
    215 
    216   if ((status = read(MBrec))) {
    217     if (status != -1) {
    218       status = 1;
    219     }
    220 
    221     return status;
    222   }
    223 
    224   scanNo      = MBrec.scanNo;
    225   cycleNo     = MBrec.cycleNo;
    226   mjd         = MBrec.mjd;
    227   interval    = MBrec.interval;
    228   fieldName   = MBrec.fieldName;
    229   srcName     = MBrec.srcName;
    230   srcDir      = MBrec.srcDir;
    231   srcPM       = MBrec.srcPM;
    232   srcVel      = MBrec.srcVel;
    233   obsType     = MBrec.obsType;
    234   IFno        = MBrec.IFno;
    235   refFreq     = MBrec.refFreq;
    236   bandwidth   = MBrec.bandwidth;
    237   freqInc     = MBrec.freqInc;
    238   restFreq    = MBrec.restFreq;
    239   tcal        = MBrec.tcal;
    240   tcalTime    = MBrec.tcalTime;
    241   azimuth     = MBrec.azimuth;
    242   elevation   = MBrec.elevation;
    243   parAngle    = MBrec.parAngle;
    244   focusAxi    = MBrec.focusAxi;
    245   focusTan    = MBrec.focusTan;
    246   focusRot    = MBrec.focusRot;
    247   temperature = MBrec.temperature;
    248   pressure    = MBrec.pressure;
    249   humidity    = MBrec.humidity;
    250   windSpeed   = MBrec.windSpeed;
    251   windAz      = MBrec.windAz;
    252   refBeam     = MBrec.refBeam;
    253   beamNo      = MBrec.beamNo;
    254   direction   = MBrec.direction;
    255   scanRate    = MBrec.scanRate;
    256   tsys        = MBrec.tsys;
    257   sigma       = MBrec.sigma;
    258   calFctr     = MBrec.calFctr;
    259   baseLin     = MBrec.baseLin;
    260   baseSub     = MBrec.baseSub;
    261   spectra     = MBrec.spectra;
    262   flagged     = MBrec.flagged;
    263   xCalFctr    = MBrec.xCalFctr;
    264   xPol        = MBrec.xPol;
    265 
    266   return 0;
    267 }
     190  PKSreader *reader = getPKSreader(name, directories, retry, interpolate,
     191                                   iDir, format);
     192
     193  // Try to open it.
     194  if (reader) {
     195    if (reader->open(name, beams, IFs, nChan, nPol, haveXPol, haveBase,
     196                     haveSpectra)) {
     197      format += " OPEN ERROR";
     198      delete reader;
     199      reader = 0x0;
     200    }
     201  }
     202
     203  return reader;
     204}
  • trunk/external/atnf/PKSIO/PKSreader.h

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSreader.h,v 19.13 2008-06-26 01:50:24 cal103 Exp $
     28//# $Id: PKSreader.h,v 19.22 2008-11-17 06:44:34 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/02, Mark Calabretta, ATNF
     
    3434#define ATNF_PKSREADER_H
    3535
     36#include <atnf/PKSIO/PKSmsg.h>
     37#include <atnf/PKSIO/PKSrecord.h>
     38
    3639#include <casa/aips.h>
    3740#include <casa/Arrays/Matrix.h>
    3841#include <casa/Arrays/Vector.h>
    39 #include <casa/BasicSL/Complex.h>
    4042#include <casa/BasicSL/String.h>
    4143
     
    4547// Class to read Parkes multibeam data.
    4648// </summary>
     49
     50// Return an appropriate PKSreader for a Parkes Multibeam dataset.
     51class PKSreader* getPKSreader(
     52        const String name,
     53        const Int retry,
     54        const Int interpolate,
     55        String &format);
     56
     57// As above, but search a list of directories for it.
     58class PKSreader* getPKSreader(
     59        const String name,
     60        const Vector<String> directories,
     61        const Int retry,
     62        const Int interpolate,
     63        Int    &iDir,
     64        String &format);
    4765
    4866// Open an appropriate PKSreader for a Parkes Multibeam dataset.
     
    7694        Bool   &haveSpectra);
    7795
    78 class MBrecord;
    79 
    80 class PKSreader
     96class PKSreader : public PKSmsg
    8197{
    8298  public:
     
    116132    // Set data selection criteria.  Channel numbering is 1-relative, zero or
    117133    // negative channel numbers are taken to be offsets from the last channel.
     134    // Coordinate system selection (only supported for SDFITS input):
     135    //   0: equatorial (RA,Dec),
     136    //   1: vertical (Az,El),
     137    //   2: feed-plane.
    118138    virtual uInt select(
    119139        const Vector<Bool> beamSel,
     
    124144        const Bool getSpectra = True,
    125145        const Bool getXPol    = False,
    126         const Bool getFeedPos = False) = 0;
     146        const Int  coordSys   = 0) = 0;
    127147
    128148    // Find the range of the data selected in time and position.
     
    133153        Matrix<Double> &positions) = 0;
    134154
    135     // Read the next data record (MBrecord is defined below).
    136     virtual Int read(MBrecord &mbrec) = 0;
    137 
    138     // Read the next data record (for backwards compatibility, do not use).
    139     virtual Int read(
    140         Int             &scanNo,
    141         Int             &cycleNo,
    142         Double          &mjd,
    143         Double          &interval,
    144         String          &fieldName,
    145         String          &srcName,
    146         Vector<Double>  &srcDir,
    147         Vector<Double>  &srcPM,
    148         Double          &srcVel,
    149         String          &obsType,
    150         Int             &IFno,
    151         Double          &refFreq,
    152         Double          &bandwidth,
    153         Double          &freqInc,
    154         Double          &restFreq,
    155         Vector<Float>   &tcal,
    156         String          &tcalTime,
    157         Float           &azimuth,
    158         Float           &elevation,
    159         Float           &parAngle,
    160         Float           &focusAxi,
    161         Float           &focusTan,
    162         Float           &focusRot,
    163         Float           &temperature,
    164         Float           &pressure,
    165         Float           &humidity,
    166         Float           &windSpeed,
    167         Float           &windAz,
    168         Int             &refBeam,
    169         Int             &beamNo,
    170         Vector<Double>  &direction,
    171         Vector<Double>  &scanRate,
    172         Vector<Float>   &tsys,
    173         Vector<Float>   &sigma,
    174         Vector<Float>   &calFctr,
    175         Matrix<Float>   &baseLin,
    176         Matrix<Float>   &baseSub,
    177         Matrix<Float>   &spectra,
    178         Matrix<uChar>   &flagged,
    179         Complex         &xCalFctr,
    180         Vector<Complex> &xPol);
    181 
    182     // Read the next data record, just the basics.
    183     virtual Int read(
    184         Int           &IFno,
    185         Vector<Float> &tsys,
    186         Vector<Float> &calFctr,
    187         Matrix<Float> &baseLin,
    188         Matrix<Float> &baseSub,
    189         Matrix<Float> &spectra,
    190         Matrix<uChar> &flagged) = 0;
     155    // Read the next data record.
     156    virtual Int read(PKSrecord &pksrec) = 0;
    191157
    192158    // Close the input file.
     
    194160
    195161  protected:
    196     Bool   cGetFeedPos, cGetSpectra, cGetXPol;
     162    Bool  cGetSpectra, cGetXPol;
     163    Int   cCoordSys;
    197164
    198165    Vector<uInt> cNChan, cNPol;
     
    200167};
    201168
    202 
    203 // Essentially just a struct used as a function argument.
    204 class MBrecord
    205 {
    206   public:
    207     Int             scanNo;
    208     Int             cycleNo;
    209     Double          mjd;
    210     Double          interval;
    211     String          fieldName;
    212     String          srcName;
    213     Vector<Double>  srcDir;
    214     Vector<Double>  srcPM;
    215     Double          srcVel;
    216     String          obsType;
    217     Int             IFno;
    218     Double          refFreq;
    219     Double          bandwidth;
    220     Double          freqInc;
    221     Double          restFreq;
    222     Vector<Float>   tcal;
    223     String          tcalTime;
    224     Float           azimuth;
    225     Float           elevation;
    226     Float           parAngle;
    227     Float           focusAxi;
    228     Float           focusTan;
    229     Float           focusRot;
    230     Float           temperature;
    231     Float           pressure;
    232     Float           humidity;
    233     Float           windSpeed;
    234     Float           windAz;
    235     Int             refBeam;
    236     Int             beamNo;
    237     Vector<Double>  direction;
    238     Vector<Double>  scanRate;
    239     Int             rateAge;
    240     Int             rateson;
    241     Vector<Float>   tsys;
    242     Vector<Float>   sigma;
    243     Vector<Float>   calFctr;
    244     Matrix<Float>   baseLin;
    245     Matrix<Float>   baseSub;
    246     Matrix<Float>   spectra;
    247     Matrix<uChar>   flagged;
    248     Complex         xCalFctr;
    249     Vector<Complex> xPol;
    250 };
    251 
    252169#endif
  • trunk/external/atnf/PKSIO/PKSwriter.h

    r1399 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSwriter.h,v 19.14 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: PKSwriter.h,v 19.16 2008-11-17 06:46:36 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
    3131#ifndef ATNF_PKSWRITER_H
    3232#define ATNF_PKSWRITER_H
     33
     34#include <atnf/PKSIO/PKSmsg.h>
     35#include <atnf/PKSIO/PKSrecord.h>
    3336
    3437#include <casa/aips.h>
     
    3841#include <casa/BasicSL/String.h>
    3942
     43#include <casa/namespace.h>
     44
    4045// <summary>
    4146// Class to write out Parkes multibeam data.
    4247// </summary>
    4348
    44 #include <casa/namespace.h>
    45 class PKSwriter
     49class PKSwriter : public PKSmsg
    4650{
    4751  public:
     
    6771    // Write the next data record.
    6872    virtual Int write (
    69         const Int             scanNo,
    70         const Int             cycleNo,
    71         const Double          mjd,
    72         const Double          interval,
    73         const String          fieldName,
    74         const String          srcName,
    75         const Vector<Double>  srcDir,
    76         const Vector<Double>  srcPM,
    77         const Double          srcVel,
    78         const String          obsMode,
    79         const Int             IFno,
    80         const Double          refFreq,
    81         const Double          bandwidth,
    82         const Double          freqInc,
    83         const Double          restFreq,
    84         const Vector<Float>   tcal,
    85         const String          tcalTime,
    86         const Float           azimuth,
    87         const Float           elevation,
    88         const Float           parAngle,
    89         const Float           focusAxi,
    90         const Float           focusTan,
    91         const Float           focusRot,
    92         const Float           temperature,
    93         const Float           pressure,
    94         const Float           humidity,
    95         const Float           windSpeed,
    96         const Float           windAz,
    97         const Int             refBeam,
    98         const Int             beamNo,
    99         const Vector<Double>  direction,
    100         const Vector<Double>  scanRate,
    101         const Vector<Float>   tsys,
    102         const Vector<Float>   sigma,
    103         const Vector<Float>   calFctr,
    104         const Matrix<Float>   baseLin,
    105         const Matrix<Float>   baseSub,
    106         const Matrix<Float>   &spectra,
    107         const Matrix<uChar>   &flagged,
    108         const Complex         xCalFctr,
    109         const Vector<Complex> &xPol) = 0;
     73        const PKSrecord &pksrec) = 0;
    11074
    11175    // Write a history record.
  • trunk/external/atnf/PKSIO/SDFITSreader.cc

    r1427 r1452  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: SDFITSreader.cc,v 19.24 2008-06-26 02:13:11 cal103 Exp $
     28//# $Id: SDFITSreader.cc,v 19.33 2008-11-17 06:58:34 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# The SDFITSreader class reads single dish FITS files such as those written
     
    3434//#---------------------------------------------------------------------------
    3535
     36#include <atnf/pks/pks_maths.h>
     37#include <atnf/PKSIO/PKSmsg.h>
     38#include <atnf/PKSIO/MBrecord.h>
     39#include <atnf/PKSIO/SDFITSreader.h>
     40
     41#include <casa/math.h>
     42#include <casa/stdio.h>
     43
    3644#include <algorithm>
    3745#include <strings.h>
    38 
    39 // AIPS++ includes.
    40 #include <casa/iostream.h>
    41 #include <casa/math.h>
    42 #include <casa/stdio.h>
    43 
    44 // ATNF includes.
    45 #include <atnf/pks/pks_maths.h>
    46 #include <atnf/PKSIO/PKSMBrecord.h>
    47 #include <atnf/PKSIO/SDFITSreader.h>
    48 
    4946
    5047class FITSparm
     
    8683  cEndChan   = 0x0;
    8784  cRefChan   = 0x0;
     85
     86  // By default, messages are written to stderr.
     87  initMsg();
    8888}
    8989
     
    114114        int    &extraSysCal)
    115115{
     116  // Clear the message stack.
     117  clearMsg();
     118
    116119  if (cSDptr) {
    117120    close();
     
    121124  cStatus = 0;
    122125  if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
    123     cerr << "Failed to open SDFITS file: " << sdName << endl;
    124     reportError();
     126    sprintf(cMsg, "ERROR: Failed to open SDFITS file\n       %s", sdName);
     127    logMsg(cMsg);
    125128    return 1;
    126129  }
     
    139142        cALFA_CIMA = 1;
    140143
     144        // Check for later versions of CIMAFITS.
     145        float version;
     146        readParm("VERSION", TFLOAT, &version);
     147        if (version >= 2.0f) cALFA_CIMA = int(version);
     148
    141149      } else {
    142         cerr << "Failed to locate SDFITS binary table." << endl;
    143         reportError();
     150        logMsg("ERROR: Failed to locate SDFITS binary table.");
    144151        close();
    145152        return 1;
     
    177184    cNAxis = 5;
    178185    if (readDim(DATA, 1, &cNAxis, cNAxes)) {
    179       reportError();
     186      logMsg();
    180187      close();
    181188      return 1;
     
    196203    if (cNAxis < 4) {
    197204      // Need at least four axes (for now).
    198       cerr << "DATA array contains fewer than four axes." << endl;
     205      logMsg("ERROR: DATA array contains fewer than four axes.");
    199206      close();
    200207      return 1;
    201208    } else if (cNAxis > 5) {
    202209      // We support up to five axes.
    203       cerr << "DATA array contains more than five axes." << endl;
     210      logMsg("ERROR: DATA array contains more than five axes.");
    204211      close();
    205212      return 1;
     
    212219    findData(DATAXED, "DATAXED", TSTRING);
    213220    if (cData[DATAXED].colnum < 0) {
    214       cerr << "DATA array column absent from binary table." << endl;
     221      logMsg("ERROR: DATA array column absent from binary table.");
    215222      close();
    216223      return 1;
     
    242249
    243250  if (cStatus) {
    244     reportError();
     251    logMsg();
    245252    close();
    246253    return 1;
     
    295302  for (int iaxis = 0; iaxis < 4; iaxis++) {
    296303    if (cReqax[iaxis] < 0) {
    297       cerr << "Could not find required DATA array axes." << endl;
     304      logMsg("ERROR: Could not find required DATA array axes.");
    298305      close();
    299306      return 1;
     
    345352
    346353  if (cStatus) {
    347     reportError();
     354    logMsg();
    348355    close();
    349356    return 1;
     
    356363    cALFAscan = 0;
    357364    cScanNo = 0;
    358     if (cALFA_BD) {
     365    if (cALFA_CIMA) {
     366      findData(SCAN,  "SCAN_ID", TINT);
     367      if (cALFA_CIMA > 1) {
     368        findData(CYCLE, "RECNUM", TINT);
     369      } else {
     370        findData(CYCLE, "SUBSCAN", TINT);
     371      }
     372    } else if (cALFA_BD) {
    359373      findData(SCAN,  "SCAN_NUMBER", TINT);
    360374      findData(CYCLE, "PATTERN_NUMBER", TINT);
    361     } else if (cALFA_CIMA) {
    362       findData(SCAN,  "SCAN_ID", TINT);
    363       findData(CYCLE, "SUBSCAN", TINT);
    364375    }
    365376  } else {
     
    372383  // Beam number, 1-relative by default.
    373384  cBeam_1rel = 1;
    374   if (cData[BEAM].colnum < 0) {
     385  if (cALFA) {
     386    // ALFA INPUT_ID, 0-relative (overrides BEAM column if present).
     387    findData(BEAM, "INPUT_ID", TSHORT);
     388    cBeam_1rel = 0;
     389
     390  } else if (cData[BEAM].colnum < 0) {
    375391    if (beamCRVAL) {
    376392      // There is a BEAM axis.
    377393      findData(BEAM, beamCRVAL, TDOUBLE);
    378 
    379394    } else {
    380       if (cALFA) {
    381         // ALFA data, 0-relative.
    382         findData(BEAM, "INPUT_ID", TSHORT);
    383       } else {
    384         // ms2sdfits output, 0-relative "feed" number.
    385         findData(BEAM, "MAIN_FEED1", TSHORT);
    386       }
    387 
     395      // ms2sdfits output, 0-relative "feed" number.
     396      findData(BEAM, "MAIN_FEED1", TSHORT);
    388397      cBeam_1rel = 0;
    389398    }
     
    394403  if (cALFA && cData[IF].colnum < 0) {
    395404    // ALFA data, 0-relative.
    396     findData(IF, "IFVAL", TSHORT);
     405    if (cALFA_CIMA > 1) {
     406      findData(IF, "IFN", TSHORT);
     407    } else {
     408      findData(IF, "IFVAL", TSHORT);
     409    }
    397410    cIF_1rel = 0;
    398411  }
     
    477490  fits_get_num_rows(cSDptr, &cNRow, &cStatus);
    478491  if (!cNRow) {
    479     cerr << "Table contains no entries." << endl;
     492    logMsg("ERROR: Table contains no entries.");
    480493    close();
    481494    return 1;
     
    490503    if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
    491504                      &beamNul, beamCol, &anynul, &cStatus)) {
    492       reportError();
    493505      delete [] beamCol;
     506      logMsg();
    494507      close();
    495508      return 1;
     
    505518      // Check validity.
    506519      if (beamCol[irow] < cBeam_1rel) {
    507         cerr << "SDFITS file contains invalid beam number." << endl;
    508520        delete [] beamCol;
     521        logMsg("ERROR: SDFITS file contains invalid beam number.");
    509522        close();
    510523        return 1;
     
    546559    if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
    547560                      &IFNul, IFCol, &anynul, &cStatus)) {
    548       reportError();
    549561      delete [] IFCol;
     562      logMsg();
    550563      close();
    551564      return 1;
     
    561574      // Check validity.
    562575      if (IFCol[irow] < cIF_1rel) {
    563         cerr << "SDFITS file contains invalid IF number." << endl;
    564576        delete [] IFCol;
     577        logMsg("ERROR: SDFITS file contains invalid IF number.");
    565578        close();
    566579        return 1;
     
    594607            // Variable dimension array.
    595608            if (readDim(DATA, irow+1, &cNAxis, cNAxes)) {
    596               reportError();
     609              logMsg();
    597610              close();
    598611              return 1;
     
    622635
    623636          if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
    624             reportError();
     637            logMsg();
    625638            close();
    626639            return 1;
     
    656669  }
    657670
    658   if (cALFA) {
    659     // ALFA labels each polarization as a separate IF.
     671  if (cALFA && cALFA_CIMA < 2) {
     672    // Older ALFA data labels each polarization as a separate IF.
    660673    cNPol[0] = cNIF;
    661674    cNIF = 1;
     
    776789
    777790  // Brightness unit.
    778   strcpy(bunit, cData[DATA].units);
     791  if (cData[DATAXED].colnum >= 0) {
     792    strcpy(bunit, "Jy");
     793  } else {
     794    strcpy(bunit, cData[DATA].units);
     795  }
     796
    779797  if (strcmp(bunit, "JY") == 0) {
    780798    bunit[1] = 'y';
     
    836854
    837855  if (cStatus) {
    838     reportError();
     856    logMsg();
    839857    return 1;
    840858  }
     
    849867
    850868  if (cStatus) {
    851     reportError();
     869    logMsg();
    852870    return 1;
    853871  }
     
    902920    if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
    903921                      &IFNul, IFCol, &anynul, &cStatus)) {
    904       reportError();
    905922      delete [] IFCol;
     923      logMsg();
    906924      close();
    907925      return 1;
     
    966984        double* &positions)
    967985{
    968   int anynul;
    969 
    970986  // Has the file been opened?
    971987  if (!cSDptr) {
     
    981997  }
    982998
     999  int anynul;
    9831000  if (cData[BEAM].colnum > 0) {
    9841001    short *beamCol = new short[cNRow];
     
    9861003    if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
    9871004                      &beamNul, beamCol, &anynul, &cStatus)) {
    988       reportError();
    9891005      delete [] beamCol;
    9901006      delete [] sel;
     1007      logMsg();
    9911008      return 1;
    9921009    }
     
    10061023    if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
    10071024                      &IFNul, IFCol, &anynul, &cStatus)) {
    1008       reportError();
    10091025      delete [] IFCol;
    10101026      delete [] sel;
     1027      logMsg();
    10111028      return 1;
    10121029    }
     
    10661083
    10671084  // Retrieve positions for selected data.
    1068   double *ra  = new double[cNRow];
    1069   double *dec = new double[cNRow];
    1070   fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum,  1, 1, nRow, 0, ra,
    1071                 &anynul, &cStatus);
    1072   fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
    1073                 &anynul, &cStatus);
    1074 
    1075   if (cALFA_BD) {
    1076     for (int irow = 0; irow < nRow; irow++) {
    1077       // Convert hours to degrees.
    1078       ra[irow] *= 15.0;
    1079     }
    1080   }
    1081 
    10821085  int isel = 0;
    10831086  positions = new double[2*nSel];
    10841087
    1085   // Parameters needed to compute feed-plane coordinates.
    1086   double *srcRA, *srcDec;
    1087   float  *par, *rot;
    1088   if (cGetFeedPos) {
    1089     srcRA  = new double[cNRow];
    1090     srcDec = new double[cNRow];
    1091     par    = new float[cNRow];
    1092     rot    = new float[cNRow];
    1093     fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum,   1, 1, nRow, 0,
    1094                   srcRA,  &anynul, &cStatus);
    1095     fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum,  1, 1, nRow, 0,
    1096                   srcDec, &anynul, &cStatus);
    1097     fits_read_col(cSDptr, TFLOAT,  cData[PARANGLE].colnum, 1, 1, nRow, 0,
    1098                   par,    &anynul, &cStatus);
    1099     fits_read_col(cSDptr, TFLOAT,  cData[FOCUSROT].colnum, 1, 1, nRow, 0,
    1100                   rot,    &anynul, &cStatus);
    1101 
    1102     for (int irow = 0; irow < nRow; irow++) {
    1103       if (sel[irow]) {
    1104         // Convert to feed-plane coordinates.
    1105         Double dist, pa;
    1106         distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R, srcDec[irow]*D2R,
    1107                dist, pa);
    1108 
    1109         Double spin = (par[irow] + rot[irow])*D2R - pa + PI;
    1110         if (spin > 2.0*PI) spin -= 2.0*PI;
    1111         Double squint = PI/2.0 - dist;
    1112 
    1113         positions[isel++] = spin;
    1114         positions[isel++] = squint;
    1115       }
    1116     }
    1117 
    1118     delete [] srcRA;
    1119     delete [] srcDec;
    1120     delete [] par;
    1121     delete [] rot;
     1088  if (cCoordSys == 1) {
     1089    // Vertical (Az,El).
     1090    if (cData[AZIMUTH].colnum  < 1 ||
     1091        cData[ELEVATIO].colnum < 1) {
     1092      logMsg("WARNING: Azimuth/elevation information absent.");
     1093      cStatus = -1;
     1094
     1095    } else {
     1096      float *az = new float[cNRow];
     1097      float *el = new float[cNRow];
     1098      fits_read_col(cSDptr, TFLOAT, cData[AZIMUTH].colnum,  1, 1, nRow, 0, az,
     1099                    &anynul, &cStatus);
     1100      fits_read_col(cSDptr, TFLOAT, cData[ELEVATIO].colnum, 1, 1, nRow, 0, el,
     1101                    &anynul, &cStatus);
     1102
     1103      if (!cStatus) {
     1104        for (int irow = 0; irow < nRow; irow++) {
     1105          if (sel[irow]) {
     1106            positions[isel++] = az[irow] * D2R;
     1107            positions[isel++] = el[irow] * D2R;
     1108          }
     1109        }
     1110      }
     1111
     1112      delete [] az;
     1113      delete [] el;
     1114    }
    11221115
    11231116  } else {
    1124     for (int irow = 0; irow < nRow; irow++) {
    1125       if (sel[irow]) {
    1126         positions[isel++] =  ra[irow] * D2R;
    1127         positions[isel++] = dec[irow] * D2R;
    1128       }
    1129     }
    1130   }
    1131 
     1117    double *ra  = new double[cNRow];
     1118    double *dec = new double[cNRow];
     1119    fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum,  1, 1, nRow, 0, ra,
     1120                  &anynul, &cStatus);
     1121    fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
     1122                  &anynul, &cStatus);
     1123    if (cStatus) {
     1124      delete [] ra;
     1125      delete [] dec;
     1126      goto cleanup;
     1127    }
     1128
     1129    if (cALFA_BD) {
     1130      for (int irow = 0; irow < nRow; irow++) {
     1131        // Convert hours to degrees.
     1132        ra[irow] *= 15.0;
     1133      }
     1134    }
     1135
     1136    if (cCoordSys == 0) {
     1137      // Equatorial (RA,Dec).
     1138      for (int irow = 0; irow < nRow; irow++) {
     1139        if (sel[irow]) {
     1140          positions[isel++] =  ra[irow] * D2R;
     1141          positions[isel++] = dec[irow] * D2R;
     1142        }
     1143      }
     1144
     1145    } else if (cCoordSys == 2) {
     1146      // Feed-plane.
     1147      if (cData[OBJ_RA].colnum   < 0 ||
     1148          cData[OBJ_DEC].colnum  < 0 ||
     1149          cData[PARANGLE].colnum < 1 ||
     1150          cData[FOCUSROT].colnum < 1) {
     1151        logMsg("WARNING: Insufficient information to compute feed-plane\n"
     1152               "         coordinates.");
     1153        cStatus = -1;
     1154
     1155      } else {
     1156        double *srcRA  = new double[cNRow];
     1157        double *srcDec = new double[cNRow];
     1158        float  *par = new float[cNRow];
     1159        float  *rot = new float[cNRow];
     1160
     1161        if (cData[OBJ_RA].colnum == 0) {
     1162          // Header keyword.
     1163          readData(OBJ_RA, 0, srcRA);
     1164          for (int irow = 1; irow < nRow; irow++) {
     1165            srcRA[irow] = *srcRA;
     1166          }
     1167        } else {
     1168          // Table column.
     1169          fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum,   1, 1, nRow,
     1170                        0, srcRA,  &anynul, &cStatus);
     1171        }
     1172
     1173        if (cData[OBJ_DEC].colnum == 0) {
     1174          // Header keyword.
     1175          readData(OBJ_DEC, 0, srcDec);
     1176          for (int irow = 1; irow < nRow; irow++) {
     1177            srcDec[irow] = *srcDec;
     1178          }
     1179        } else {
     1180          // Table column.
     1181          fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum,  1, 1, nRow,
     1182                        0, srcDec, &anynul, &cStatus);
     1183        }
     1184
     1185        fits_read_col(cSDptr, TFLOAT,  cData[PARANGLE].colnum, 1, 1, nRow, 0,
     1186                      par,    &anynul, &cStatus);
     1187        fits_read_col(cSDptr, TFLOAT,  cData[FOCUSROT].colnum, 1, 1, nRow, 0,
     1188                      rot,    &anynul, &cStatus);
     1189
     1190        if (!cStatus) {
     1191          for (int irow = 0; irow < nRow; irow++) {
     1192            if (sel[irow]) {
     1193              // Convert to feed-plane coordinates.
     1194              Double dist, pa;
     1195              distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R,
     1196                     srcDec[irow]*D2R, dist, pa);
     1197
     1198              Double spin = (par[irow] + rot[irow])*D2R - pa + PI;
     1199              if (spin > 2.0*PI) spin -= 2.0*PI;
     1200              Double squint = PI/2.0 - dist;
     1201
     1202              positions[isel++] = spin;
     1203              positions[isel++] = squint;
     1204            }
     1205          }
     1206        }
     1207
     1208        delete [] srcRA;
     1209        delete [] srcDec;
     1210        delete [] par;
     1211        delete [] rot;
     1212      }
     1213    }
     1214
     1215    delete [] ra;
     1216    delete [] dec;
     1217  }
     1218
     1219cleanup:
    11321220  delete [] sel;
    1133   delete [] ra;
    1134   delete [] dec;
    1135 
    1136   return cStatus;
     1221
     1222  if (cStatus) {
     1223    nSel = 0;
     1224    delete [] positions;
     1225    logMsg();
     1226    cStatus = 0;
     1227    return 1;
     1228  }
     1229
     1230  return 0;
    11371231}
    11381232
     
    11431237
    11441238int SDFITSreader::read(
    1145         PKSMBrecord &mbrec)
     1239        MBrecord &mbrec)
    11461240{
    11471241  // Has the file been opened?
     
    11751269          readData(OBSMODE, cRow, chars);
    11761270          if (strcmp(chars, "CAL") == 0) {
    1177             // iIF is really the polarization in ALFA data.
    1178             alfaCal(iBeam, iIF);
    1179             continue;
     1271            if (cALFA_CIMA > 1) {
     1272              for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
     1273                alfaCal(iBeam, iIF, iPol);
     1274              }
     1275              continue;
     1276            } else {
     1277              // iIF is really the polarization in older ALFA data.
     1278              alfaCal(iBeam, 0, iIF);
     1279              continue;
     1280            }
    11801281          }
    11811282        }
     
    12891390    mbrec.decRate = scanrate[1] * D2R;
    12901391  }
    1291   mbrec.rateAge = 0;
    1292   mbrec.rateson = 0;
     1392  mbrec.paRate = 0.0f;
    12931393
    12941394  // IF-dependent parameters.
     
    13281428
    13291429  if (cStatus) {
    1330     reportError();
     1430    logMsg();
    13311431    return 1;
    13321432  }
     
    13651465
    13661466  if (cStatus) {
    1367     reportError();
     1467    logMsg();
    13681468    return 1;
    13691469  }
     
    13921492      trc[cReqax[1]] = ipol+1;
    13931493
    1394       if (cALFA) {
     1494      if (cALFA && cALFA_CIMA < 2) {
    13951495        // ALFA data: polarizations are stored in successive rows.
    13961496        blc[cReqax[1]] = 1;
     
    14111511
    14121512        if ((status = readDim(DATA, cRow, &naxis, cNAxes))) {
    1413           reportError();
     1513          logMsg();
    14141514
    14151515        } else if ((status = (naxis != cNAxis))) {
    1416           cerr << "DATA array dimensions changed." << endl;
     1516          logMsg("ERROR: DATA array dimensions changed.");
    14171517        }
    14181518
     
    14281528          blc, trc, inc, 0, mbrec.spectra[0] + ipol*nChan, &anynul,
    14291529          &cStatus)) {
    1430         reportError();
     1530        logMsg();
    14311531        delete [] blc;
    14321532        delete [] trc;
     
    14741574            cNAxes, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,
    14751575            &cStatus)) {
    1476           reportError();
     1576          logMsg();
    14771577          delete [] blc;
    14781578          delete [] trc;
     
    15251625    if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
    15261626        blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
    1527       reportError();
     1627      logMsg();
    15281628      delete [] blc;
    15291629      delete [] trc;
     
    15561656
    15571657  if (cStatus) {
    1558     reportError();
     1658    logMsg();
    15591659    return 1;
    15601660  }
     
    15641664  readData(TCAL,     cRow, &mbrec.tcal[0]);
    15651665  readData(TCALTIME, cRow,  mbrec.tcalTime);
     1666
    15661667  readData(AZIMUTH,  cRow, &mbrec.azimuth);
    15671668  readData(ELEVATIO, cRow, &mbrec.elevation);
    15681669  readData(PARANGLE, cRow, &mbrec.parAngle);
     1670
    15691671  readData(FOCUSAXI, cRow, &mbrec.focusAxi);
    15701672  readData(FOCUSTAN, cRow, &mbrec.focusTan);
    15711673  readData(FOCUSROT, cRow, &mbrec.focusRot);
     1674
    15721675  readData(TAMBIENT, cRow, &mbrec.temp);
    15731676  readData(PRESSURE, cRow, &mbrec.pressure);
     
    15881691
    15891692  if (cStatus) {
    1590     reportError();
     1693    logMsg();
    15911694    return 1;
    15921695  }
    15931696
    15941697  return 0;
    1595 }
    1596 
    1597 
    1598 //------------------------------------------------------ SDFITSreader::alfaCal
    1599 
    1600 // Process ALFA calibration data.
    1601 
    1602 int SDFITSreader::alfaCal(
    1603         short iBeam,
    1604         short iPol)
    1605 {
    1606   int  calOn;
    1607   char chars[32];
    1608   if (cALFA_BD) {
    1609     readData("OBS_NAME", TSTRING, cRow, chars);
    1610   } else {
    1611     readData("SCANTYPE", TSTRING, cRow, chars);
    1612   }
    1613 
    1614   if (strcmp(chars, "ON") == 0) {
    1615     calOn = 1;
    1616   } else if (strcmp(chars, "OFF") == 0) {
    1617     calOn = 0;
    1618   } else {
    1619     return 1;
    1620   }
    1621 
    1622   // Read cal data.
    1623   long *blc = new long[cNAxis+1];
    1624   long *trc = new long[cNAxis+1];
    1625   long *inc = new long[cNAxis+1];
    1626   for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
    1627     blc[iaxis] = 1;
    1628     trc[iaxis] = 1;
    1629     inc[iaxis] = 1;
    1630   }
    1631 
    1632   // User channel selection.
    1633   int startChan = cStartChan[0];
    1634   int endChan   = cEndChan[0];
    1635 
    1636   blc[cNAxis] = cRow;
    1637   trc[cNAxis] = cRow;
    1638   blc[cReqax[0]] = std::min(startChan, endChan);
    1639   trc[cReqax[0]] = std::max(startChan, endChan);
    1640   blc[cReqax[1]] = 1;
    1641   trc[cReqax[1]] = 1;
    1642 
    1643   float spectrum[endChan];
    1644   int anynul;
    1645   if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
    1646       blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
    1647     reportError();
    1648     delete [] blc;
    1649     delete [] trc;
    1650     delete [] inc;
    1651     return 1;
    1652   }
    1653 
    1654   // Average the spectrum.
    1655   float mean = 1e9f;
    1656   for (int k = 0; k < 2; k++) {
    1657     float discrim = 2.0f * mean;
    1658 
    1659     int nChan = 0;
    1660     float sum = 0.0f;
    1661 
    1662     float *chanN = spectrum + abs(endChan - startChan) + 1;
    1663     for (float *chan = spectrum; chan < chanN; chan++) {
    1664       // Simple discriminant that eliminates strong radar interference.
    1665       if (*chan < discrim) {
    1666         nChan++;
    1667         sum += *chan;
    1668       }
    1669     }
    1670 
    1671     mean = sum / nChan;
    1672   }
    1673 
    1674   if (calOn) {
    1675     cALFAcalOn[iBeam][iPol]  += mean;
    1676   } else {
    1677     cALFAcalOff[iBeam][iPol] += mean;
    1678   }
    1679 
    1680   if (cALFAcalOn[iBeam][iPol] != 0.0f &&
    1681       cALFAcalOff[iBeam][iPol] != 0.0f) {
    1682     // Tcal should come from the TCAL table, it varies weakly with beam,
    1683     // polarization, and frequency.  However, TCAL is not written properly.
    1684     float Tcal = 12.0f;
    1685     cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
    1686                                     cALFAcalOff[iBeam][iPol]);
    1687 
    1688     // Scale from K to Jy; the gain also varies weakly with beam,
    1689     // polarization, frequency, and zenith angle.
    1690     float fluxCal = 10.0f;
    1691     cALFAcal[iBeam][iPol] /= fluxCal;
    1692 
    1693     cALFAcalOn[iBeam][iPol]  = 0.0f;
    1694     cALFAcalOff[iBeam][iPol] = 0.0f;
    1695   }
    1696 
    1697   return 0;
    1698 }
    1699 
    1700 
    1701 //-------------------------------------------------- SDFITSreader::reportError
    1702 
    1703 // Print the error message corresponding to the input status value and all the
    1704 // messages on the CFITSIO error stack to stderr.
    1705 
    1706 void SDFITSreader::reportError()
    1707 {
    1708   fits_report_error(stderr, cStatus);
    17091698}
    17101699
     
    17251714    if (cEndChan)   delete [] cEndChan;
    17261715    if (cRefChan)   delete [] cRefChan;
     1716  }
     1717}
     1718
     1719//------------------------------------------------------- SDFITSreader::logMsg
     1720
     1721// Log a message.  If the current CFITSIO status value is non-zero, also log
     1722// the corresponding error message and the CFITSIO message stack.
     1723
     1724void SDFITSreader::logMsg(const char *msg)
     1725{
     1726  FITSreader::logMsg(msg);
     1727
     1728  if (cStatus > 0) {
     1729    fits_get_errstatus(cStatus, cMsg);
     1730    FITSreader::logMsg(cMsg);
     1731
     1732    while (fits_read_errmsg(cMsg)) {
     1733      FITSreader::logMsg(cMsg);
     1734    }
    17271735  }
    17281736}
     
    20112019  }
    20122020}
     2021
     2022//------------------------------------------------------ SDFITSreader::alfaCal
     2023
     2024// Process ALFA calibration data.
     2025
     2026int SDFITSreader::alfaCal(
     2027        short iBeam,
     2028        short iIF,
     2029        short iPol)
     2030{
     2031  int  calOn;
     2032  char chars[32];
     2033  if (cALFA_BD) {
     2034    readData("OBS_NAME", TSTRING, cRow, chars);
     2035  } else {
     2036    readData("SCANTYPE", TSTRING, cRow, chars);
     2037  }
     2038
     2039  if (strcmp(chars, "ON") == 0) {
     2040    calOn = 1;
     2041  } else if (strcmp(chars, "OFF") == 0) {
     2042    calOn = 0;
     2043  } else {
     2044    return 1;
     2045  }
     2046
     2047  // Read cal data.
     2048  long *blc = new long[cNAxis+1];
     2049  long *trc = new long[cNAxis+1];
     2050  long *inc = new long[cNAxis+1];
     2051  for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
     2052    blc[iaxis] = 1;
     2053    trc[iaxis] = 1;
     2054    inc[iaxis] = 1;
     2055  }
     2056
     2057  // User channel selection.
     2058  int startChan = cStartChan[iIF];
     2059  int endChan   = cEndChan[iIF];
     2060
     2061  blc[cNAxis] = cRow;
     2062  trc[cNAxis] = cRow;
     2063  blc[cReqax[0]] = std::min(startChan, endChan);
     2064  trc[cReqax[0]] = std::max(startChan, endChan);
     2065  if (cALFA_CIMA > 1) {
     2066    // CIMAFITS 2.x has a legitimate STOKES axis...
     2067    blc[cReqax[1]] = iPol+1;
     2068    trc[cReqax[1]] = iPol+1;
     2069  } else {
     2070    // ...older ALFA data does not.
     2071    blc[cReqax[1]] = 1;
     2072    trc[cReqax[1]] = 1;
     2073  }
     2074
     2075  float spectrum[endChan];
     2076  int anynul;
     2077  if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
     2078      blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
     2079    logMsg();
     2080    delete [] blc;
     2081    delete [] trc;
     2082    delete [] inc;
     2083    return 1;
     2084  }
     2085
     2086  // Average the spectrum.
     2087  float mean = 1e9f;
     2088  for (int k = 0; k < 2; k++) {
     2089    float discrim = 2.0f * mean;
     2090
     2091    int nChan = 0;
     2092    float sum = 0.0f;
     2093
     2094    float *chanN = spectrum + abs(endChan - startChan) + 1;
     2095    for (float *chan = spectrum; chan < chanN; chan++) {
     2096      // Simple discriminant that eliminates strong radar interference.
     2097      if (*chan < discrim) {
     2098        nChan++;
     2099        sum += *chan;
     2100      }
     2101    }
     2102
     2103    mean = sum / nChan;
     2104  }
     2105
     2106  if (calOn) {
     2107    cALFAcalOn[iBeam][iPol]  += mean;
     2108  } else {
     2109    cALFAcalOff[iBeam][iPol] += mean;
     2110  }
     2111
     2112  if (cALFAcalOn[iBeam][iPol] != 0.0f &&
     2113      cALFAcalOff[iBeam][iPol] != 0.0f) {
     2114    // Tcal should come from the TCAL table, it varies weakly with beam,
     2115    // polarization, and frequency.  However, TCAL is not written properly.
     2116    float Tcal = 12.0f;
     2117    cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
     2118                                    cALFAcalOff[iBeam][iPol]);
     2119
     2120    // Scale from K to Jy; the gain also varies weakly with beam,
     2121    // polarization, frequency, and zenith angle.
     2122    float fluxCal = 10.0f;
     2123    cALFAcal[iBeam][iPol] /= fluxCal;
     2124
     2125    cALFAcalOn[iBeam][iPol]  = 0.0f;
     2126    cALFAcalOff[iBeam][iPol] = 0.0f;
     2127  }
     2128
     2129  return 0;
     2130}
  • trunk/external/atnf/PKSIO/SDFITSreader.h

    r1399 r1452  
    22//# SDFITSreader.h: ATNF CFITSIO interface class for SDFITS input.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2007
     4//# Copyright (C) 2000-2008
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: SDFITSreader.h,v 19.13 2007/11/12 03:37:56 cal103 Exp $
     28//# $Id: SDFITSreader.h,v 19.16 2008-11-17 06:47:05 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# The SDFITSreader class reads single dish FITS files such as those written
     
    3838
    3939#include <atnf/PKSIO/FITSreader.h>
    40 #include <atnf/PKSIO/PKSMBrecord.h>
     40#include <atnf/PKSIO/MBrecord.h>
    4141
    4242#include <fitsio.h>
     43
     44using namespace std;
    4345
    4446// <summary>
     
    100102
    101103    // Read the next data record.
    102     virtual int read(PKSMBrecord &record);
    103 
    104     // Print out CFITSIO error messages.
    105     void reportError(void);
     104    virtual int read(MBrecord &record);
    106105
    107106    // Close the SDFITS file.
     
    125124          HUMIDITY, WINDSPEE, WINDDIRE, NDATA};
    126125
     126    // Message handling.
     127    virtual void logMsg(const char *msg = 0x0);
     128
    127129    void findData(int iData, char *name, int type);
    128130    int   readDim(int iData, long iRow, int *naxis, long naxes[]);
     
    132134    void  findCol(char *name, int *colnum);
    133135
    134     // These are for ALFA data, "BDFITS" or "CIMAFITS".
     136    // These are for ALFA data: "BDFITS" or "CIMAFITS".
    135137    int   cALFA, cALFA_BD, cALFA_CIMA, cALFAscan, cScanNo;
    136138    float cALFAcal[8][2], cALFAcalOn[8][2], cALFAcalOff[8][2];
    137     int   alfaCal(short iBeam, short iIF);
     139    int   alfaCal(short iBeam, short iIF, short iPol);
    138140
    139141    // These are for GBT data.
  • trunk/external/atnf/PKSIO/SDFITSwriter.cc

    r1399 r1452  
    22//# SDFITSwriter.cc: ATNF CFITSIO interface class for SDFITS output.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2007
     4//# Copyright (C) 2000-2008
    55//# Mark Calabretta, ATNF
    66//#
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: SDFITSwriter.cc,v 19.12 2007/11/12 03:37:56 cal103 Exp $
     29//# $Id: SDFITSwriter.cc,v 19.15 2008-11-17 06:58:58 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# Original: 2000/07/24, Mark Calabretta, ATNF
    3232//#---------------------------------------------------------------------------
    3333
     34#include <atnf/PKSIO/MBrecord.h>
     35#include <atnf/PKSIO/PKSmsg.h>
     36#include <atnf/PKSIO/SDFITSwriter.h>
     37
     38#include <casa/iostream.h>
     39
    3440#include <algorithm>
    3541#include <math.h>
    36 
    37 // AIPS++ includes.
    38 #include <casa/iostream.h>
    39 
    40 // ATNF includes.
    41 #include <atnf/PKSIO/PKSMBrecord.h>
    42 #include <atnf/PKSIO/SDFITSwriter.h>
    4342
    4443using namespace std;
     
    5554{
    5655  // Default constructor.
    57   cSDptr = 0;
     56  cSDptr = 0x0;
     57
     58  // By default, messages are written to stderr.
     59  initMsg();
    5860}
    5961
     
    8688        int    extraSysCal)
    8789{
     90  if (cSDptr) {
     91    logMsg("ERROR: Output file already open, close it first.");
     92    return 1;
     93  }
     94
     95  // Clear the message stack.
     96  clearMsg();
     97
    8898  // Prepend an '!' to the output name to force it to be overwritten.
    8999  char sdname[80];
     
    94104  cStatus = 0;
    95105  if (fits_create_file(&cSDptr, sdname, &cStatus)) {
     106    sprintf(cMsg, "ERROR: Failed to create SDFITS file\n       %s", sdName);
     107    logMsg(cMsg);
    96108    return cStatus;
    97109  }
     
    141153  // Write required primary header keywords.
    142154  if (fits_write_imghdr(cSDptr, 8, 0, 0, &cStatus)) {
     155    logMsg("ERROR: Failed to write required primary header keywords.");
    143156    return cStatus;
    144157  }
     
    160173  char version[7];
    161174  char date[11];
    162   sscanf("$Revision: 19.12 $", "%*s%s", version);
    163   sscanf("$Date: 2007/11/12 03:37:56 $", "%*s%s", date);
     175  sscanf("$Revision: 19.15 $", "%*s%s", version);
     176  sscanf("$Date: 2008-11-17 06:58:58 $", "%*s%s", date);
    164177  sprintf(text, "SDFITSwriter (v%s, %s)", version, date);
    165178  fits_write_key_str(cSDptr, "ORIGIN", text, "output class", &cStatus);
     
    171184  fits_write_comment(cSDptr, text, &cStatus);
    172185
     186  if (cStatus) {
     187    logMsg("ERROR: Failed in writing primary header.");
     188    return cStatus;
     189  }
     190
     191
    173192  // Create an SDFITS extension.
    174193  long nrow = 0;
     
    176195  if (fits_create_tbl(cSDptr, BINARY_TBL, nrow, ncol, NULL, NULL, NULL,
    177196      "SINGLE DISH", &cStatus)) {
     197    logMsg("ERROR: Failed to create a binary table extension.");
    178198    return 1;
    179199  }
     
    524544  }
    525545
     546  if (cStatus) {
     547    logMsg("ERROR: Failed in writing binary table header.");
     548  }
     549
    526550  return cStatus;
    527551}
     
    531555// Write a record to the SDFITS file.
    532556
    533 int SDFITSwriter::write(PKSMBrecord &mbrec)
     557int SDFITSwriter::write(MBrecord &mbrec)
    534558{
    535559  char *cptr;
     
    740764  }
    741765
     766  if (cStatus) {
     767    logMsg("ERROR: Failed in writing binary table entry.");
     768  }
     769
    742770  return cStatus;
    743771}
     
    751779
    752780{
    753   fits_write_history(cSDptr, text, &cStatus);
     781  if (!cSDptr) {
     782    return 1;
     783  }
     784
     785  if (fits_write_history(cSDptr, text, &cStatus)) {
     786    logMsg("ERROR: Failed in writing HISTORY records.");
     787  }
    754788
    755789  return cStatus;
    756 }
    757 
    758 //-------------------------------------------------- SDFITSwriter::reportError
    759 
    760 // Print the error message corresponding to the input status value and all the
    761 // messages on the CFITSIO error stack to stderr.
    762 
    763 void SDFITSwriter::reportError()
    764 {
    765   fits_report_error(stderr, cStatus);
    766790}
    767791
     
    774798  if (cSDptr) {
    775799    cStatus = 0;
    776     fits_close_file(cSDptr, &cStatus);
     800    if (fits_close_file(cSDptr, &cStatus)) {
     801      logMsg("ERROR: Failed to close file.");
     802    }
     803
    777804    cSDptr = 0;
    778805  }
     
    787814  if (cSDptr) {
    788815    cStatus = 0;
    789     fits_delete_file(cSDptr, &cStatus);
     816    if (fits_delete_file(cSDptr, &cStatus)) {
     817      logMsg("ERROR: Failed to close and delete file.");
     818    }
     819
    790820    cSDptr = 0;
    791821  }
    792822}
     823
     824//------------------------------------------------------- SDFITSwriter::logMsg
     825
     826// Log a message.  If the current CFITSIO status value is non-zero, also log
     827// the corresponding error message and dump the CFITSIO message stack.
     828
     829void SDFITSwriter::logMsg(const char *msg)
     830{
     831  PKSmsg::logMsg(msg);
     832
     833  if (cStatus) {
     834    fits_get_errstatus(cStatus, cMsg);
     835    PKSmsg::logMsg(cMsg);
     836
     837    while (fits_read_errmsg(cMsg)) {
     838      PKSmsg::logMsg(cMsg);
     839    }
     840  }
     841}
  • trunk/external/atnf/PKSIO/SDFITSwriter.h

    r1399 r1452  
    22//# SDFITSwriter.h: ATNF CFITSIO interface class for SDFITS output.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2007
     4//# Copyright (C) 2000-2008
    55//# Mark Calabretta, ATNF
    66//#
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: SDFITSwriter.h,v 19.7 2007/11/12 03:37:56 cal103 Exp $
     29//# $Id: SDFITSwriter.h,v 19.9 2008-11-17 06:48:32 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# Original: 2000/07/24, Mark Calabretta, ATNF
     
    3535#define ATNF_SDFITSWRITER_H
    3636
    37 #include <atnf/PKSIO/PKSMBrecord.h>
     37#include <atnf/PKSIO/MBrecord.h>
     38#include <atnf/PKSIO/PKSmsg.h>
    3839
    3940#include <fitsio.h>
     41
     42using namespace std;
    4043
    4144// <summary>
     
    4346// </summary>
    4447
    45 class SDFITSwriter
     48class SDFITSwriter : public PKSmsg
    4649{
    4750  public:
     
    5053
    5154    // Destructor.
    52     ~SDFITSwriter();
     55    virtual ~SDFITSwriter();
    5356
    5457    // Create a new SDFITSwriter and store static data.
     
    7174
    7275    // Store time-variable data.
    73     int write(PKSMBrecord &record);
     76    int write(MBrecord &record);
    7477
    7578    // Write a history record.
    7679    int history(char* text);
    77 
    78     // Print out CFITSIO error messages.
    79     void reportError();
    8080
    8181    // Close the SDFITS file.
     
    9090         *cNChan, cNIF, *cNPol, cStatus;
    9191    long cRow;
     92
     93    // Message handling.
     94    char cMsg[256];
     95    virtual void logMsg(const char *msg = 0x0);
    9296};
    9397
  • trunk/external/atnf/PKSIO/makefile

    r1325 r1452  
    1 # $Id: makefile,v 19.0 2003/07/16 03:34:05 aips2adm Exp $
     1# $Id: makefile,v 19.0 2003-07-16 03:34:05 aips2adm Exp $
    22
    33XLIBLIST := CFITSIO RPFITS
  • trunk/external/atnf/pks/makefile

    r1325 r1452  
    1 # $Id: makefile,v 19.0 2003/07/16 03:33:47 aips2adm Exp $
     1# $Id: makefile,v 19.0 2003-07-16 03:33:47 aips2adm Exp $
    22
    33# Use the generic AIPS++ class implementation makefile.
  • trunk/external/atnf/pks/pks_maths.cc

    r1427 r1452  
    22//# pks_maths.cc: Mathematical functions for Parkes single-dish data reduction
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 1994-2006
     4//# Copyright (C) 1994-2008
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2727//#
    2828//# Original: Mark Calabretta
    29 //# $Id: pks_maths.cc,v 1.5 2006-05-19 00:12:35 mcalabre Exp $
     29//# $Id: pks_maths.cc,v 1.6 2008-10-17 02:36:16 cal103 Exp $
    3030//----------------------------------------------------------------------------
    3131
     
    295295//----------------------------------------------------------------------- azel
    296296
    297 // Convert (ra,dec) to (az,el), from
    298 // http://aa.usno.navy.mil/faq/docs/Alt_Az.html.  Position as a Cartesian
    299 // triplet in m, UT1 in MJD form, and all angles in radian.
     297// Convert (ra,dec) to (az,el).  Position as a Cartesian triplet in m, UT1 in
     298// MJD form, and all angles in radian.
    300299
    301300void azel(const Vector<Double> position, Double ut1, Double ra, Double dec,
    302301          Double &az, Double &el)
    303302{
    304   // Get gocentric longitude and latitude (rad).
     303  // Get geocentric longitude and latitude (rad).
    305304  Double x = position(0);
    306305  Double y = position(1);
     
    318317
    319318  // Azimuth and elevation (rad).
    320   az = atan2(cos(dec)*sin(ha), cos(dec)*sin(lat)*cos(ha) - sin(dec)*cos(lat));
     319  az = atan2(-cos(dec)*sin(ha),
     320            sin(dec)*cos(lat) - cos(dec)*sin(lat)*cos(ha));
     321  el = asin(sin(dec)*sin(lat) + cos(dec)*cos(lat)*cos(ha));
     322
    321323  if (az < 0.0) az += C::_2pi;
    322   el = asin(cos(dec)*cos(lat)*cos(ha) + sin(dec)*sin(lat));
    323324}
    324325
Note: See TracChangeset for help on using the changeset viewer.