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

update from livedata CVS

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.