//#---------------------------------------------------------------------------
//# SDFITSreader.cc: ATNF interface class for SDFITS input using CFITSIO.
//#---------------------------------------------------------------------------
//# livedata - processing pipeline for single-dish, multibeam spectral data.
//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
//#
//# This file is part of livedata.
//#
//# livedata is free software: you can redistribute it and/or modify it under
//# the terms of the GNU General Public License as published by the Free
//# Software Foundation, either version 3 of the License, or (at your option)
//# any later version.
//#
//# livedata is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
//# more details.
//#
//# You should have received a copy of the GNU General Public License along
//# with livedata. If not, see .
//#
//# Correspondence concerning livedata may be directed to:
//# Internet email: mcalabre@atnf.csiro.au
//# Postal address: Dr. Mark Calabretta
//# Australia Telescope National Facility, CSIRO
//# PO Box 76
//# Epping NSW 1710
//# AUSTRALIA
//#
//# http://www.atnf.csiro.au/computing/software/livedata.html
//# $Id: SDFITSreader.cc,v 19.45 2009-09-30 07:23:48 cal103 Exp $
//#---------------------------------------------------------------------------
//# The SDFITSreader class reads single dish FITS files such as those written
//# by SDFITSwriter containing Parkes Multibeam data.
//#
//# Original: 2000/08/09, Mark Calabretta, ATNF
//#---------------------------------------------------------------------------
#include
#include
#include
#include
#include
#include
#include
#include
class FITSparm
{
public:
char *name; // Keyword or column name.
int type; // Expected keyvalue or column data type.
int colnum; // Column number; 0 for keyword; -1 absent.
int coltype; // Column data type, as found.
long nelem; // Column data repeat count; < 0 for vardim.
int tdimcol; // TDIM column number; 0 for keyword; -1 absent.
char units[32]; // Units from TUNITn keyword.
};
// Numerical constants.
const double PI = 3.141592653589793238462643;
// Factor to convert radians to degrees.
const double D2R = PI / 180.0;
//---------------------------------------------------- SDFITSreader::(statics)
int SDFITSreader::sInit = 1;
int SDFITSreader::sReset = 0;
int (*SDFITSreader::sALFAcalNon)[2] = (int (*)[2])(new float[16]);
int (*SDFITSreader::sALFAcalNoff)[2] = (int (*)[2])(new float[16]);
float (*SDFITSreader::sALFAcalOn)[2] = (float (*)[2])(new float[16]);
float (*SDFITSreader::sALFAcalOff)[2] = (float (*)[2])(new float[16]);
float (*SDFITSreader::sALFAcal)[2] = (float (*)[2])(new float[16]);
//------------------------------------------------- SDFITSreader::SDFITSreader
SDFITSreader::SDFITSreader()
{
// Default constructor.
cSDptr = 0x0;
// Allocate space for data descriptors.
cData = new FITSparm[NDATA];
for (int iData = 0; iData < NDATA; iData++) {
cData[iData].colnum = -1;
}
// Initialize pointers.
cBeams = 0x0;
cIFs = 0x0;
cStartChan = 0x0;
cEndChan = 0x0;
cRefChan = 0x0;
// By default, messages are written to stderr.
initMsg();
}
//------------------------------------------------ SDFITSreader::~SDFITSreader
SDFITSreader::~SDFITSreader()
{
close();
delete [] cData;
}
//--------------------------------------------------------- SDFITSreader::open
// Open an SDFITS file for reading.
int SDFITSreader::open(
char* sdName,
int &nBeam,
int* &beams,
int &nIF,
int* &IFs,
int* &nChan,
int* &nPol,
int* &haveXPol,
int &haveBase,
int &haveSpectra,
int &extraSysCal)
{
// Clear the message stack.
clearMsg();
if (cSDptr) {
close();
}
// Open the SDFITS file.
cStatus = 0;
if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
sprintf(cMsg, "ERROR: Failed to open SDFITS file\n %s", sdName);
logMsg(cMsg);
return 1;
}
// Move to the SDFITS extension.
cALFA = cALFA_BD = cALFA_CIMA = 0;
if (fits_movnam_hdu(cSDptr, BINARY_TBL, "SINGLE DISH", 0, &cStatus)) {
// No SDFITS table, look for BDFITS or CIMAFITS.
cStatus = 0;
if (fits_movnam_hdu(cSDptr, BINARY_TBL, "BDFITS", 0, &cStatus) == 0) {
cALFA_BD = 1;
} else {
cStatus = 0;
if (fits_movnam_hdu(cSDptr, BINARY_TBL, "CIMAFITS", 0, &cStatus) == 0) {
cALFA_CIMA = 1;
// Check for later versions of CIMAFITS.
float version;
readParm("VERSION", TFLOAT, &version);
if (version >= 2.0f) cALFA_CIMA = int(version);
} else {
logMsg("ERROR: Failed to locate SDFITS binary table.");
close();
return 1;
}
}
// Arecibo ALFA data of some kind.
cALFA = 1;
if (sInit) {
for (int iBeam = 0; iBeam < 8; iBeam++) {
for (int iPol = 0; iPol < 2; iPol++) {
sALFAcalOn[iBeam][iPol] = 0.0f;
sALFAcalOff[iBeam][iPol] = 0.0f;
// Nominal factor to calibrate spectra in Jy.
sALFAcal[iBeam][iPol] = 3.0f;
}
}
sInit = 0;
}
}
// GBT data.
char telescope[32];
readParm("TELESCOP", TSTRING, telescope); // Core.
cGBT = strncmp(telescope, "GBT", 3) == 0 ||
strncmp(telescope, "NRAO_GBT", 8) == 0;
// Check that the DATA array column is present.
findData(DATA, "DATA", TFLOAT);
haveSpectra = cHaveSpectra = cData[DATA].colnum > 0;
cNAxisTime = 0;
if (cHaveSpectra) {
// Find the number of data axes (must be the same for each IF).
cNAxes = 5;
if (readDim(DATA, 1, &cNAxes, cNAxis)) {
logMsg();
close();
return 1;
}
if (cALFA_BD) {
// ALFA BDFITS: variable length arrays don't actually vary and there is
// no TDIM (or MAXISn) card; use the LAGS_IN value.
cNAxes = 5;
readParm("LAGS_IN", TLONG, cNAxis);
cNAxis[1] = 1;
cNAxis[2] = 1;
cNAxis[3] = 1;
cNAxis[4] = 1;
cData[DATA].nelem = cNAxis[0];
}
if (cNAxes < 4) {
// Need at least four axes (for now).
logMsg("ERROR: DATA array contains fewer than four axes.");
close();
return 1;
} else if (cNAxes > 5) {
// We support up to five axes.
logMsg("ERROR: DATA array contains more than five axes.");
close();
return 1;
}
findData(FLAGGED, "FLAGGED", TBYTE);
} else {
// DATA column not present, check for a DATAXED keyword.
findData(DATAXED, "DATAXED", TSTRING);
if (cData[DATAXED].colnum < 0) {
logMsg("ERROR: DATA array column absent from binary table.");
close();
return 1;
}
// Determine the number of axes and their length.
char dataxed[32];
readParm("DATAXED", TSTRING, dataxed);
for (int iaxis = 0; iaxis < 5; iaxis++) cNAxis[iaxis] = 0;
sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1, cNAxis+2,
cNAxis+3, cNAxis+4);
for (int iaxis = 4; iaxis > -1; iaxis--) {
if (cNAxis[iaxis] == 0) cNAxes = iaxis;
}
}
char *CTYPE[5] = {"CTYPE1", "CTYPE2", "CTYPE3", "CTYPE4", "CTYPE5"};
char *CRPIX[5] = {"CRPIX1", "CRPIX2", "CRPIX3", "CRPIX4", "CRPIX5"};
char *CRVAL[5] = {"CRVAL1", "CRVAL2", "CRVAL3", "CRVAL4", "CRVAL5"};
char *CDELT[5] = {"CDELT1", "CDELT2", "CDELT3", "CDELT4", "CDELT5"};
// Find required DATA array axes.
char ctype[5][72];
for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
strcpy(ctype[iaxis], "");
readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]); // Core.
}
if (cStatus) {
logMsg();
close();
return 1;
}
char *fqCRVAL = 0;
char *fqCDELT = 0;
char *fqCRPIX = 0;
char *raCRVAL = 0;
char *decCRVAL = 0;
char *timeCRVAL = 0;
char *timeCDELT = 0;
char *timeCRPIX = 0;
char *beamCRVAL = 0;
cFreqAxis = -1;
cStokesAxis = -1;
cRaAxis = -1;
cDecAxis = -1;
cTimeAxis = -1;
cBeamAxis = -1;
for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
if (strncmp(ctype[iaxis], "FREQ", 4) == 0) {
cFreqAxis = iaxis;
fqCRVAL = CRVAL[iaxis];
fqCDELT = CDELT[iaxis];
fqCRPIX = CRPIX[iaxis];
} else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
cStokesAxis = iaxis;
} else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
cRaAxis = iaxis;
raCRVAL = CRVAL[iaxis];
} else if (strncmp(ctype[iaxis], "DEC", 3) == 0) {
cDecAxis = iaxis;
decCRVAL = CRVAL[iaxis];
} else if (strcmp(ctype[iaxis], "TIME") == 0) {
// TIME (UTC seconds since midnight); axis type, if present, takes
// precedence over keyword.
cTimeAxis = iaxis;
timeCRVAL = CRVAL[iaxis];
// Check for non-degeneracy.
if ((cNAxisTime = cNAxis[iaxis]) > 1) {
timeCDELT = CDELT[iaxis];
timeCRPIX = CRPIX[iaxis];
sprintf(cMsg, "DATA array contains a TIME axis of length %ld.",
cNAxisTime);
logMsg(cMsg);
}
} else if (strcmp(ctype[iaxis], "BEAM") == 0) {
// BEAM can be a keyword or axis type.
cBeamAxis = iaxis;
beamCRVAL = CRVAL[iaxis];
}
}
if (cALFA_BD) {
// Fixed in ALFA CIMAFITS.
cRaAxis = 2;
raCRVAL = "CRVAL2A";
cDecAxis = 3;
decCRVAL = "CRVAL3A";
}
// Check that required axes are present.
if (cFreqAxis < 0 || cStokesAxis < 0 || cRaAxis < 0 || cDecAxis < 0) {
logMsg("ERROR: Could not find required DATA array axes.");
close();
return 1;
}
// Set up machinery for data retrieval.
findData(SCAN, "SCAN", TINT); // Shared.
findData(CYCLE, "CYCLE", TINT); // Additional.
findData(DATE_OBS, "DATE-OBS", TSTRING); // Core.
if (cTimeAxis >= 0) {
// The DATA array has a TIME axis.
if (cNAxisTime > 1) {
// Non-degenerate.
findData(TimeRefVal, timeCRVAL, TDOUBLE); // Time reference value.
findData(TimeDelt, timeCDELT, TDOUBLE); // Time increment.
findData(TimeRefPix, timeCRPIX, TFLOAT); // Time reference pixel.
} else {
// Degenerate, treat its like a simple TIME keyword.
findData(TIME, timeCRVAL, TDOUBLE);
}
} else {
findData(TIME, "TIME", TDOUBLE); // Core.
}
findData(EXPOSURE, "EXPOSURE", TFLOAT); // Core.
findData(OBJECT, "OBJECT", TSTRING); // Core.
findData(OBJ_RA, "OBJ-RA", TDOUBLE); // Additional.
findData(OBJ_DEC, "OBJ-DEC", TDOUBLE); // Additional.
findData(RESTFRQ, "RESTFRQ", TDOUBLE); // Additional.
findData(OBSMODE, "OBSMODE", TSTRING); // Shared.
findData(BEAM, "BEAM", TSHORT); // Additional.
findData(IF, "IF", TSHORT); // Additional.
findData(FqRefVal, fqCRVAL, TDOUBLE); // Frequency reference value.
findData(FqDelt, fqCDELT, TDOUBLE); // Frequency increment.
findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel.
findData(RA, raCRVAL, TDOUBLE); // Right ascension.
findData(DEC, decCRVAL, TDOUBLE); // Declination.
findData(SCANRATE, "SCANRATE", TFLOAT); // Additional.
findData(TSYS, "TSYS", TFLOAT); // Core.
findData(CALFCTR, "CALFCTR", TFLOAT); // Additional.
findData(XCALFCTR, "XCALFCTR", TFLOAT); // Additional.
findData(BASELIN, "BASELIN", TFLOAT); // Additional.
findData(BASESUB, "BASESUB", TFLOAT); // Additional.
findData(XPOLDATA, "XPOLDATA", TFLOAT); // Additional.
findData(REFBEAM, "REFBEAM", TSHORT); // Additional.
findData(TCAL, "TCAL", TFLOAT); // Shared.
findData(TCALTIME, "TCALTIME", TSTRING); // Additional.
findData(AZIMUTH, "AZIMUTH", TFLOAT); // Shared.
findData(ELEVATIO, "ELEVATIO", TFLOAT); // Shared.
findData(PARANGLE, "PARANGLE", TFLOAT); // Additional.
findData(FOCUSAXI, "FOCUSAXI", TFLOAT); // Additional.
findData(FOCUSTAN, "FOCUSTAN", TFLOAT); // Additional.
findData(FOCUSROT, "FOCUSROT", TFLOAT); // Additional.
findData(TAMBIENT, "TAMBIENT", TFLOAT); // Shared.
findData(PRESSURE, "PRESSURE", TFLOAT); // Shared.
findData(HUMIDITY, "HUMIDITY", TFLOAT); // Shared.
findData(WINDSPEE, "WINDSPEE", TFLOAT); // Shared.
findData(WINDDIRE, "WINDDIRE", TFLOAT); // Shared.
if (cStatus) {
logMsg();
close();
return 1;
}
// Check for alternative column names.
if (cALFA) {
// ALFA data.
cALFAscan = 0;
cScanNo = 0;
if (cALFA_CIMA) {
findData(SCAN, "SCAN_ID", TINT);
if (cALFA_CIMA > 1) {
// Note that RECNUM increases by cNAxisTime per row.
findData(CYCLE, "RECNUM", TINT);
} else {
findData(CYCLE, "SUBSCAN", TINT);
}
} else if (cALFA_BD) {
findData(SCAN, "SCAN_NUMBER", TINT);
findData(CYCLE, "PATTERN_NUMBER", TINT);
}
} else {
readData(SCAN, 1, &cFirstScanNo);
}
cCycleNo = 0;
cLastUTC = 0.0;
// Beam number, 1-relative by default.
cBeam_1rel = 1;
if (cALFA) {
// ALFA INPUT_ID, 0-relative (overrides BEAM column if present).
findData(BEAM, "INPUT_ID", TSHORT);
cBeam_1rel = 0;
} else if (cData[BEAM].colnum < 0) {
if (beamCRVAL) {
// There is a BEAM axis.
findData(BEAM, beamCRVAL, TDOUBLE);
} else {
// ms2sdfits output, 0-relative "feed" number.
findData(BEAM, "MAIN_FEED1", TSHORT);
cBeam_1rel = 0;
}
}
// IF number, 1-relative by default.
cIF_1rel = 1;
if (cALFA && cData[IF].colnum < 0) {
// ALFA data, 0-relative.
if (cALFA_CIMA > 1) {
findData(IF, "IFN", TSHORT);
} else {
findData(IF, "IFVAL", TSHORT);
}
cIF_1rel = 0;
}
// ms2sdfits writes a scalar "TSYS" column that averages the polarizations.
int colnum;
findCol("SYSCAL_TSYS", &colnum);
if (colnum > 0) {
// This contains the vector Tsys.
findData(TSYS, "SYSCAL_TSYS", TFLOAT);
}
// XPOLDATA?
if (cData[SCANRATE].colnum < 0) {
findData(SCANRATE, "FIELD_POINTING_DIR_RATE", TFLOAT);
}
if (cData[RESTFRQ].colnum < 0) {
findData(RESTFRQ, "RESTFREQ", TDOUBLE);
if (cData[RESTFRQ].colnum < 0) {
findData(RESTFRQ, "SPECTRAL_WINDOW_REST_FREQUENCY", TDOUBLE);
}
}
if (cData[OBJ_RA].colnum < 0) {
findData(OBJ_RA, "SOURCE_DIRECTION", TDOUBLE);
}
if (cData[OBJ_DEC].colnum < 0) {
findData(OBJ_DEC, "SOURCE_DIRECTION", TDOUBLE);
}
// REFBEAM?
if (cData[TCAL].colnum < 0) {
findData(TCAL, "SYSCAL_TCAL", TFLOAT);
} else if (cALFA_BD) {
// ALFA BDFITS has a different TCAL with 64 elements - kill it!
findData(TCAL, "NO NO NO", TFLOAT);
}
if (cALFA_BD) {
// ALFA BDFITS.
findData(AZIMUTH, "CRVAL2B", TFLOAT);
findData(ELEVATIO, "CRVAL3B", TFLOAT);
}
if (cALFA) {
// ALFA data.
findData(PARANGLE, "PARA_ANG", TFLOAT);
}
if (cData[TAMBIENT].colnum < 0) {
findData(TAMBIENT, "WEATHER_TEMPERATURE", TFLOAT);
}
if (cData[PRESSURE].colnum < 0) {
findData(PRESSURE, "WEATHER_PRESSURE", TFLOAT);
}
if (cData[HUMIDITY].colnum < 0) {
findData(HUMIDITY, "WEATHER_REL_HUMIDITY", TFLOAT);
}
if (cData[WINDSPEE].colnum < 0) {
findData(WINDSPEE, "WEATHER_WIND_SPEED", TFLOAT);
}
if (cData[WINDDIRE].colnum < 0) {
findData(WINDDIRE, "WEATHER_WIND_DIRECTION", TFLOAT);
}
// Find the number of rows.
fits_get_num_rows(cSDptr, &cNRow, &cStatus);
if (!cNRow) {
logMsg("ERROR: Table contains no entries.");
close();
return 1;
}
// Determine which beams are present in the data.
if (cData[BEAM].colnum > 0) {
short *beamCol = new short[cNRow];
short beamNul = 1;
int anynul;
if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
&beamNul, beamCol, &anynul, &cStatus)) {
delete [] beamCol;
logMsg();
close();
return 1;
}
// Find the maximum beam number.
cNBeam = cBeam_1rel - 1;
for (int irow = 0; irow < cNRow; irow++) {
if (beamCol[irow] > cNBeam) {
cNBeam = beamCol[irow];
}
// Check validity.
if (beamCol[irow] < cBeam_1rel) {
delete [] beamCol;
logMsg("ERROR: SDFITS file contains invalid beam number.");
close();
return 1;
}
}
if (!cBeam_1rel) cNBeam++;
// Find all beams present in the data.
cBeams = new int[cNBeam];
for (int ibeam = 0; ibeam < cNBeam; ibeam++) {
cBeams[ibeam] = 0;
}
for (int irow = 0; irow < cNRow; irow++) {
cBeams[beamCol[irow] - cBeam_1rel] = 1;
}
delete [] beamCol;
} else {
// No BEAM column.
cNBeam = 1;
cBeams = new int[1];
cBeams[0] = 1;
}
// Passing back the address of the array allows PKSFITSreader::select() to
// modify its elements directly.
nBeam = cNBeam;
beams = cBeams;
// Determine which IFs are present in the data.
if (cData[IF].colnum > 0) {
short *IFCol = new short[cNRow];
short IFNul = 1;
int anynul;
if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
&IFNul, IFCol, &anynul, &cStatus)) {
delete [] IFCol;
logMsg();
close();
return 1;
}
// Find the maximum IF number.
cNIF = cIF_1rel - 1;
for (int irow = 0; irow < cNRow; irow++) {
if (IFCol[irow] > cNIF) {
cNIF = IFCol[irow];
}
// Check validity.
if (IFCol[irow] < cIF_1rel) {
delete [] IFCol;
logMsg("ERROR: SDFITS file contains invalid IF number.");
close();
return 1;
}
}
if (!cIF_1rel) cNIF++;
// Find all IFs present in the data.
cIFs = new int[cNIF];
cNChan = new int[cNIF];
cNPol = new int[cNIF];
cHaveXPol = new int[cNIF];
cGetXPol = 0;
for (int iIF = 0; iIF < cNIF; iIF++) {
cIFs[iIF] = 0;
cNChan[iIF] = 0;
cNPol[iIF] = 0;
cHaveXPol[iIF] = 0;
}
for (int irow = 0; irow < cNRow; irow++) {
int iIF = IFCol[irow] - cIF_1rel;
if (cIFs[iIF] == 0) {
cIFs[iIF] = 1;
// Find the axis lengths.
if (cHaveSpectra) {
if (cData[DATA].nelem < 0) {
// Variable dimension array.
if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
logMsg();
close();
return 1;
}
}
} else {
if (cData[DATAXED].colnum > 0) {
char dataxed[32];
readParm("DATAXED", TSTRING, dataxed);
sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1,
cNAxis+2, cNAxis+3, cNAxis+4);
}
}
// Number of channels and polarizations.
cNChan[iIF] = cNAxis[cFreqAxis];
cNPol[iIF] = cNAxis[cStokesAxis];
cHaveXPol[iIF] = 0;
// Is cross-polarization data present?
if (cData[XPOLDATA].colnum > 0) {
// Check that it conforms.
int nAxis;
long nAxes[2];
if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
logMsg();
close();
return 1;
}
// Default is to get it if we have it.
if (nAxis == 2 &&
nAxes[0] == 2 &&
nAxes[1] == cNChan[iIF]) {
cGetXPol = cHaveXPol[iIF] = 1;
}
}
}
}
delete [] IFCol;
} else {
// No IF column.
cNIF = 1;
cIFs = new int[1];
cIFs[0] = 1;
cNChan = new int[1];
cNPol = new int[1];
cHaveXPol = new int[1];
cGetXPol = 0;
// Number of channels and polarizations.
cNChan[0] = cNAxis[cFreqAxis];
cNPol[0] = cNAxis[cStokesAxis];
cHaveXPol[0] = 0;
}
if (cALFA && cALFA_CIMA < 2) {
// Older ALFA data labels each polarization as a separate IF.
cNPol[0] = cNIF;
cNIF = 1;
}
// Passing back the address of the array allows PKSFITSreader::select() to
// modify its elements directly.
nIF = cNIF;
IFs = cIFs;
nChan = cNChan;
nPol = cNPol;
haveXPol = cHaveXPol;
// Default channel range selection.
cStartChan = new int[cNIF];
cEndChan = new int[cNIF];
cRefChan = new int[cNIF];
for (int iIF = 0; iIF < cNIF; iIF++) {
cStartChan[iIF] = 1;
cEndChan[iIF] = cNChan[iIF];
cRefChan[iIF] = cNChan[iIF]/2 + 1;
}
// Default is to get it if we have it.
cGetSpectra = cHaveSpectra;
// Are baseline parameters present?
cHaveBase = 0;
if (cData[BASELIN].colnum) {
// Check that it conforms.
int nAxis, status = 0;
long nAxes[2];
if (fits_read_tdim(cSDptr, cData[BASELIN].colnum, 2, &nAxis, nAxes,
&status) == 0) {
cHaveBase = (nAxis == 2);
}
}
haveBase = cHaveBase;
// Is extra system calibration data available?
cExtraSysCal = 0;
for (int iparm = REFBEAM; iparm < NDATA; iparm++) {
if (cData[iparm].colnum >= 0) {
cExtraSysCal = 1;
break;
}
}
extraSysCal = cExtraSysCal;
// Extras for ALFA data.
cALFAacc = 0.0f;
if (cALFA_CIMA > 1) {
// FFTs per second when the Mock correlator operates in RFI blanking mode.
readData("PHFFTACC", TFLOAT, 0, &cALFAacc);
}
cRow = 0;
cTimeIdx = cNAxisTime;
return 0;
}
//---------------------------------------------------- SDFITSreader::getHeader
// Get parameters describing the data.
int SDFITSreader::getHeader(
char observer[32],
char project[32],
char telescope[32],
double antPos[3],
char obsMode[32],
char bunit[32],
float &equinox,
char radecsys[32],
char dopplerFrame[32],
char datobs[32],
double &utc,
double &refFreq,
double &bandwidth)
{
// Has the file been opened?
if (!cSDptr) {
return 1;
}
// Read parameter values.
readParm("OBSERVER", TSTRING, observer); // Shared.
readParm("PROJID", TSTRING, project); // Shared.
readParm("TELESCOP", TSTRING, telescope); // Core.
antPos[0] = 0.0;
antPos[1] = 0.0;
antPos[2] = 0.0;
if (readParm("ANTENNA_POSITION", TDOUBLE, antPos)) {
readParm("OBSGEO-X", TDOUBLE, antPos); // Additional.
readParm("OBSGEO-Y", TDOUBLE, antPos + 1); // Additional.
readParm("OBSGEO-Z", TDOUBLE, antPos + 2); // Additional.
}
if (antPos[0] == 0.0) {
if (strncmp(telescope, "ATPKS", 5) == 0) {
// Parkes coordinates.
antPos[0] = -4554232.087;
antPos[1] = 2816759.046;
antPos[2] = -3454035.950;
} else if (strncmp(telescope, "ATMOPRA", 7) == 0) {
// Mopra coordinates.
antPos[0] = -4682768.630;
antPos[1] = 2802619.060;
antPos[2] = -3291759.900;
} else if (strncmp(telescope, "ARECIBO", 7) == 0) {
// Arecibo coordinates.
antPos[0] = 2390486.900;
antPos[1] = -5564731.440;
antPos[2] = 1994720.450;
}
}
readData(OBSMODE, 1, obsMode); // Shared.
// Brightness unit.
if (cData[DATAXED].colnum >= 0) {
strcpy(bunit, "Jy");
} else {
strcpy(bunit, cData[DATA].units);
}
if (strcmp(bunit, "JY") == 0) {
bunit[1] = 'y';
} else if (strcmp(bunit, "JY/BEAM") == 0) {
strcpy(bunit, "Jy/beam");
}
readParm("EQUINOX", TFLOAT, &equinox); // Shared.
if (cStatus == 405) {
// EQUINOX was written as string value in early versions.
cStatus = 0;
char strtmp[32];
readParm("EQUINOX", TSTRING, strtmp);
sscanf(strtmp, "%f", &equinox);
}
if (readParm("RADESYS", TSTRING, radecsys)) { // Additional.
if (readParm("RADECSYS", TSTRING, radecsys)) { // Additional.
strcpy(radecsys, "");
}
}
if (readParm("SPECSYS", TSTRING, dopplerFrame)) { // Additional.
// Fallback value.
strcpy(dopplerFrame, "TOPOCENT");
// Look for VELFRAME, written by earlier versions of Livedata.
if (readParm("VELFRAME", TSTRING, dopplerFrame)) { // Additional.
// No, try digging it out of the CTYPE card (AIPS convention).
char keyw[9], ctype[9];
sprintf(keyw, "CTYPE%ld", cFreqAxis+1);
readParm(keyw, TSTRING, ctype);
if (strncmp(ctype, "FREQ-", 5) == 0) {
strcpy(dopplerFrame, ctype+5);
if (strcmp(dopplerFrame, "LSR") == 0) {
// LSR unqualified usually means LSR (kinematic).
strcpy(dopplerFrame, "LSRK");
} else if (strcmp(dopplerFrame, "HEL") == 0) {
// Almost certainly barycentric.
strcpy(dopplerFrame, "BARYCENT");
}
} else {
strcpy(dopplerFrame, "");
}
}
// Translate to FITS standard names.
if (strncmp(dopplerFrame, "TOP", 3) == 0) {
strcpy(dopplerFrame, "TOPOCENT");
} else if (strncmp(dopplerFrame, "GEO", 3) == 0) {
strcpy(dopplerFrame, "GEOCENTR");
} else if (strncmp(dopplerFrame, "HEL", 3) == 0) {
strcpy(dopplerFrame, "HELIOCEN");
} else if (strncmp(dopplerFrame, "BARY", 4) == 0) {
strcpy(dopplerFrame, "BARYCENT");
}
}
if (cStatus) {
logMsg();
return 1;
}
// Get parameters from first row of table.
readTime(1, 1, datobs, utc);
readData(FqRefVal, 1, &refFreq);
readParm("BANDWID", TDOUBLE, &bandwidth); // Core.
if (cStatus) {
logMsg();
return 1;
}
return 0;
}
//-------------------------------------------------- SDFITSreader::getFreqInfo
// Get frequency parameters for each IF.
int SDFITSreader::getFreqInfo(
int &nIF,
double* &startFreq,
double* &endFreq)
{
float fqRefPix;
double fqDelt, fqRefVal;
nIF = cNIF;
startFreq = new double[nIF];
endFreq = new double[nIF];
if (cData[IF].colnum > 0) {
short *IFCol = new short[cNRow];
short IFNul = 1;
int anynul;
if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
&IFNul, IFCol, &anynul, &cStatus)) {
delete [] IFCol;
logMsg();
close();
return 1;
}
for (int iIF = 0; iIF < nIF; iIF++) {
if (cIFs[iIF]) {
// Find the first occurrence of this IF in the table.
int IFno = iIF + cIF_1rel;
for (int irow = 0; irow < cNRow;) {
if (IFCol[irow++] == IFno) {
readData(FqRefPix, irow, &fqRefPix);
readData(FqRefVal, irow, &fqRefVal);
readData(FqDelt, irow, &fqDelt);
if (cALFA_BD) {
unsigned char invert;
readData("UPPERSB", TBYTE, irow, &invert);
if (invert) {
fqDelt = -fqDelt;
}
}
startFreq[iIF] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
endFreq[iIF] = fqRefVal + (cNChan[iIF] - fqRefPix) * fqDelt;
break;
}
}
} else {
startFreq[iIF] = 0.0;
endFreq[iIF] = 0.0;
}
}
delete [] IFCol;
} else {
// No IF column, read the first table entry.
readData(FqRefPix, 1, &fqRefPix);
readData(FqRefVal, 1, &fqRefVal);
readData(FqDelt, 1, &fqDelt);
startFreq[0] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
endFreq[0] = fqRefVal + (cNChan[0] - fqRefPix) * fqDelt;
}
return cStatus;
}
//---------------------------------------------------- SDFITSreader::findRange
// Find the range of the data in time and position.
int SDFITSreader::findRange(
int &nRow,
int &nSel,
char dateSpan[2][32],
double utcSpan[2],
double* &positions)
{
// Has the file been opened?
if (!cSDptr) {
return 1;
}
nRow = cNRow;
// Find the number of rows selected.
short *sel = new short[cNRow];
for (int irow = 0; irow < cNRow; irow++) {
sel[irow] = 1;
}
int anynul;
if (cData[BEAM].colnum > 0) {
short *beamCol = new short[cNRow];
short beamNul = 1;
if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
&beamNul, beamCol, &anynul, &cStatus)) {
delete [] beamCol;
delete [] sel;
logMsg();
return 1;
}
for (int irow = 0; irow < cNRow; irow++) {
if (!cBeams[beamCol[irow]-cBeam_1rel]) {
sel[irow] = 0;
}
}
delete [] beamCol;
}
if (cData[IF].colnum > 0) {
short *IFCol = new short[cNRow];
short IFNul = 1;
if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
&IFNul, IFCol, &anynul, &cStatus)) {
delete [] IFCol;
delete [] sel;
logMsg();
return 1;
}
for (int irow = 0; irow < cNRow; irow++) {
if (!cIFs[IFCol[irow]-cIF_1rel]) {
sel[irow] = 0;
}
}
delete [] IFCol;
}
nSel = 0;
for (int irow = 0; irow < cNRow; irow++) {
nSel += sel[irow];
}
// Find the time range assuming the data is in chronological order.
readTime(1, 1, dateSpan[0], utcSpan[0]);
readTime(cNRow, cNAxisTime, dateSpan[1], utcSpan[1]);
// Retrieve positions for selected data.
int isel = 0;
positions = new double[2*nSel];
if (cCoordSys == 1) {
// Horizontal (Az,El).
if (cData[AZIMUTH].colnum < 0 ||
cData[ELEVATIO].colnum < 0) {
logMsg("WARNING: Azimuth/elevation information absent.");
cStatus = -1;
} else {
float *az = new float[cNRow];
float *el = new float[cNRow];
readCol(AZIMUTH, az);
readCol(ELEVATIO, el);
if (!cStatus) {
for (int irow = 0; irow < cNRow; irow++) {
if (sel[irow]) {
positions[isel++] = az[irow] * D2R;
positions[isel++] = el[irow] * D2R;
}
}
}
delete [] az;
delete [] el;
}
} else if (cCoordSys == 3) {
// ZPA-EL.
if (cData[BEAM].colnum < 0 ||
cData[FOCUSROT].colnum < 0 ||
cData[ELEVATIO].colnum < 0) {
logMsg("WARNING: ZPA/elevation information absent.");
cStatus = -1;
} else {
short *beam = new short[cNRow];
float *rot = new float[cNRow];
float *el = new float[cNRow];
readCol(BEAM, beam);
readCol(FOCUSROT, rot);
readCol(ELEVATIO, el);
if (!cStatus) {
for (int irow = 0; irow < cNRow; irow++) {
if (sel[irow]) {
Int beamNo = beam[irow];
Double zpa = rot[irow];
if (beamNo > 1) {
// Beam geometry for the Parkes multibeam.
if (beamNo < 8) {
zpa += -60.0 + 60.0*(beamNo-2);
} else {
zpa += -90.0 + 60.0*(beamNo-8);
}
if (zpa < -180.0) {
zpa += 360.0;
} else if (zpa > 180.0) {
zpa -= 360.0;
}
}
positions[isel++] = zpa * D2R;
positions[isel++] = el[irow] * D2R;
}
}
}
delete [] beam;
delete [] rot;
delete [] el;
}
} else {
double *ra = new double[cNRow];
double *dec = new double[cNRow];
readCol(RA, ra);
readCol(DEC, dec);
if (cStatus) {
delete [] ra;
delete [] dec;
goto cleanup;
}
if (cALFA_BD) {
for (int irow = 0; irow < cNRow; irow++) {
// Convert hours to degrees.
ra[irow] *= 15.0;
}
}
if (cCoordSys == 0) {
// Equatorial (RA,Dec).
for (int irow = 0; irow < cNRow; irow++) {
if (sel[irow]) {
positions[isel++] = ra[irow] * D2R;
positions[isel++] = dec[irow] * D2R;
}
}
} else if (cCoordSys == 2) {
// Feed-plane.
if (cData[OBJ_RA].colnum < 0 ||
cData[OBJ_DEC].colnum < 0 ||
cData[PARANGLE].colnum < 0 ||
cData[FOCUSROT].colnum < 0) {
logMsg("WARNING: Insufficient information to compute feed-plane\n"
" coordinates.");
cStatus = -1;
} else {
double *srcRA = new double[cNRow];
double *srcDec = new double[cNRow];
float *par = new float[cNRow];
float *rot = new float[cNRow];
readCol(OBJ_RA, srcRA);
readCol(OBJ_DEC, srcDec);
readCol(PARANGLE, par);
readCol(FOCUSROT, rot);
if (!cStatus) {
for (int irow = 0; irow < cNRow; irow++) {
if (sel[irow]) {
// Convert to feed-plane coordinates.
Double dist, pa;
distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R,
srcDec[irow]*D2R, dist, pa);
Double spin = (par[irow] + rot[irow])*D2R - pa;
if (spin > 2.0*PI) spin -= 2.0*PI;
Double squint = PI/2.0 - dist;
positions[isel++] = spin;
positions[isel++] = squint;
}
}
}
delete [] srcRA;
delete [] srcDec;
delete [] par;
delete [] rot;
}
}
delete [] ra;
delete [] dec;
}
cleanup:
delete [] sel;
if (cStatus) {
nSel = 0;
delete [] positions;
logMsg();
cStatus = 0;
return 1;
}
return 0;
}
//--------------------------------------------------------- SDFITSreader::read
// Read the next data record.
int SDFITSreader::read(
MBrecord &mbrec)
{
// Has the file been opened?
if (!cSDptr) {
return 1;
}
// Find the next selected beam and IF.
short iBeam = 0, iIF = 0;
while (1) {
if (++cTimeIdx > cNAxisTime) {
if (++cRow > cNRow) break;
cTimeIdx = 1;
}
if (cData[BEAM].colnum > 0) {
readData(BEAM, cRow, &iBeam);
// Convert to 0-relative.
if (cBeam_1rel) iBeam--;
}
if (cBeams[iBeam]) {
if (cData[IF].colnum > 0) {
readData(IF, cRow, &iIF);
// Convert to 0-relative.
if (cIF_1rel) iIF--;
}
if (cIFs[iIF]) {
if (cALFA) {
// ALFA data, check for calibration data.
char chars[32];
readData(OBSMODE, cRow, chars);
if (strcmp(chars, "DROP") == 0) {
// Completely flagged integration.
continue;
} else if (strcmp(chars, "CAL") == 0) {
sReset = 1;
if (cALFA_CIMA > 1) {
for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
alfaCal(iBeam, iIF, iPol);
}
continue;
} else {
// iIF is really the polarization in older ALFA data.
alfaCal(iBeam, 0, iIF);
continue;
}
} else {
// Reset for the next CAL record.
if (sReset) {
for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
sALFAcalNon[iBeam][iPol] = 0;
sALFAcalNoff[iBeam][iPol] = 0;
sALFAcalOn[iBeam][iPol] = 0.0f;
sALFAcalOff[iBeam][iPol] = 0.0f;
}
sReset = 0;
sprintf(cMsg, "ALFA cal factors for beam %d: %.3e, %.3e",
iBeam+1, sALFAcal[iBeam][0], sALFAcal[iBeam][1]);
logMsg(cMsg);
}
}
}
break;
}
}
}
// EOF?
if (cRow > cNRow) {
return -1;
}
if (cALFA) {
int scanNo;
readData(SCAN, cRow, &scanNo);
if (scanNo != cALFAscan) {
cScanNo++;
cALFAscan = scanNo;
}
mbrec.scanNo = cScanNo;
} else {
readData(SCAN, cRow, &mbrec.scanNo);
// Ensure that scan number is 1-relative.
mbrec.scanNo -= (cFirstScanNo - 1);
}
// Times.
char datobs[32];
readTime(cRow, cTimeIdx, datobs, mbrec.utc);
strcpy(mbrec.datobs, datobs);
if (cData[CYCLE].colnum > 0) {
readData(CYCLE, cRow, &mbrec.cycleNo);
mbrec.cycleNo += cTimeIdx - 1;
if (cALFA_BD) mbrec.cycleNo++;
} else {
// Cycle number not recorded, must do our own bookkeeping.
if (mbrec.utc != cLastUTC) {
mbrec.cycleNo = ++cCycleNo;
cLastUTC = mbrec.utc;
}
}
readData(EXPOSURE, cRow, &mbrec.exposure);
// Source identification.
readData(OBJECT, cRow, mbrec.srcName);
readData(OBJ_RA, cRow, &mbrec.srcRA);
if (strcmp(cData[OBJ_RA].name, "OBJ-RA") == 0) {
mbrec.srcRA *= D2R;
}
if (strcmp(cData[OBJ_DEC].name, "OBJ-DEC") == 0) {
readData(OBJ_DEC, cRow, &mbrec.srcDec);
mbrec.srcDec *= D2R;
}
// Line rest frequency (Hz).
readData(RESTFRQ, cRow, &mbrec.restFreq);
if (mbrec.restFreq == 0.0 && cALFA_BD) {
mbrec.restFreq = 1420.40575e6;
}
// Observation mode.
readData(OBSMODE, cRow, mbrec.obsType);
// Beam-dependent parameters.
mbrec.beamNo = iBeam + 1;
readData(RA, cRow, &mbrec.ra);
readData(DEC, cRow, &mbrec.dec);
mbrec.ra *= D2R;
mbrec.dec *= D2R;
if (cALFA_BD) mbrec.ra *= 15.0;
float scanrate[2];
readData(SCANRATE, cRow, &scanrate);
if (strcmp(cData[SCANRATE].name, "SCANRATE") == 0) {
mbrec.raRate = scanrate[0] * D2R;
mbrec.decRate = scanrate[1] * D2R;
}
mbrec.paRate = 0.0f;
// IF-dependent parameters.
int startChan = cStartChan[iIF];
int endChan = cEndChan[iIF];
int refChan = cRefChan[iIF];
// Allocate data storage.
int nChan = abs(endChan - startChan) + 1;
int nPol = cNPol[iIF];
if (cGetSpectra || cGetXPol) {
int nxpol = cGetXPol ? 2*nChan : 0;
mbrec.allocate(0, nChan*nPol, nxpol);
}
mbrec.nIF = 1;
mbrec.IFno[0] = iIF + 1;
mbrec.nChan[0] = nChan;
mbrec.nPol[0] = nPol;
readData(FqRefPix, cRow, mbrec.fqRefPix);
readData(FqRefVal, cRow, mbrec.fqRefVal);
readData(FqDelt, cRow, mbrec.fqDelt);
if (cALFA_BD) {
unsigned char invert;
int anynul, colnum;
findCol("UPPERSB", &colnum);
fits_read_col(cSDptr, TBYTE, colnum, cRow, 1, 1, 0, &invert, &anynul,
&cStatus);
if (invert) {
mbrec.fqDelt[0] = -mbrec.fqDelt[0];
}
}
if (cStatus) {
logMsg();
return 1;
}
// Adjust for channel selection.
if (mbrec.fqRefPix[0] != refChan) {
mbrec.fqRefVal[0] += (refChan - mbrec.fqRefPix[0]) * mbrec.fqDelt[0];
mbrec.fqRefPix[0] = refChan;
}
if (endChan < startChan) {
mbrec.fqDelt[0] = -mbrec.fqDelt[0];
}
// The data may only have a scalar Tsys value.
mbrec.tsys[0][0] = 0.0f;
mbrec.tsys[0][1] = 0.0f;
if (cData[TSYS].nelem >= nPol) {
readData(TSYS, cRow, mbrec.tsys[0]);
}
for (int j = 0; j < 2; j++) {
mbrec.calfctr[0][j] = 0.0f;
}
if (cData[CALFCTR].colnum > 0) {
readData(CALFCTR, cRow, mbrec.calfctr);
}
if (cHaveBase) {
mbrec.haveBase = 1;
readData(BASELIN, cRow, mbrec.baseLin);
readData(BASESUB, cRow, mbrec.baseSub);
} else {
mbrec.haveBase = 0;
}
if (cStatus) {
logMsg();
return 1;
}
// Read data, sectioning and transposing it in the process.
long *blc = new long[cNAxes+1];
long *trc = new long[cNAxes+1];
long *inc = new long[cNAxes+1];
for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
blc[iaxis] = 1;
trc[iaxis] = 1;
inc[iaxis] = 1;
}
blc[cFreqAxis] = std::min(startChan, endChan);
trc[cFreqAxis] = std::max(startChan, endChan);
if (cTimeAxis >= 0) {
blc[cTimeAxis] = cTimeIdx;
trc[cTimeAxis] = cTimeIdx;
}
blc[cNAxes] = cRow;
trc[cNAxes] = cRow;
mbrec.haveSpectra = cGetSpectra;
if (cGetSpectra) {
int anynul;
for (int iPol = 0; iPol < nPol; iPol++) {
blc[cStokesAxis] = iPol+1;
trc[cStokesAxis] = iPol+1;
if (cALFA && cALFA_CIMA < 2) {
// ALFA data: polarizations are stored in successive rows.
blc[cStokesAxis] = 1;
trc[cStokesAxis] = 1;
if (iPol) {
if (++cRow > cNRow) {
return -1;
}
blc[cNAxes] = cRow;
trc[cNAxes] = cRow;
}
} else if (cData[DATA].nelem < 0) {
// Variable dimension array; get axis lengths.
int naxes = 5, status;
if ((status = readDim(DATA, cRow, &naxes, cNAxis))) {
logMsg();
} else if ((status = (naxes != cNAxes))) {
logMsg("ERROR: DATA array dimensions changed.");
}
if (status) {
delete [] blc;
delete [] trc;
delete [] inc;
return 1;
}
}
if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
&cStatus)) {
logMsg();
delete [] blc;
delete [] trc;
delete [] inc;
return 1;
}
if (endChan < startChan) {
// Reverse the spectrum.
float *iptr = mbrec.spectra[0] + iPol*nChan;
float *jptr = iptr + nChan - 1;
float *mid = iptr + nChan/2;
while (iptr < mid) {
float tmp = *iptr;
*(iptr++) = *jptr;
*(jptr--) = tmp;
}
}
if (cALFA) {
// ALFA data, rescale the spectrum.
float el, zd;
readData(ELEVATIO, cRow, &el);
zd = 90.0f - el;
float factor = sALFAcal[iBeam][iPol] / alfaGain(zd);
if (cALFA_CIMA > 1) {
// Rescale according to the number of unblanked accumulations.
int colnum, naccum;
findCol("STAT", &colnum);
fits_read_col(cSDptr, TINT, colnum, cRow, 10*(cTimeIdx-1)+2, 1, 0,
&naccum, &anynul, &cStatus);
factor *= cALFAacc / naccum;
}
float *chan = mbrec.spectra[0] + iPol*nChan;
float *chanN = chan + nChan;
while (chan < chanN) {
// Approximate conversion to Jy.
*(chan++) *= factor;
}
}
if (mbrec.tsys[0][iPol] == 0.0) {
// Compute Tsys as the average across the spectrum.
float *chan = mbrec.spectra[0] + iPol*nChan;
float *chanN = chan + nChan;
float *tsys = mbrec.tsys[0] + iPol;
while (chan < chanN) {
*tsys += *(chan++);
}
*tsys /= nChan;
}
// Read data flags.
if (cData[FLAGGED].colnum > 0) {
if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxes,
cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
&cStatus)) {
logMsg();
delete [] blc;
delete [] trc;
delete [] inc;
return 1;
}
if (endChan < startChan) {
// Reverse the flag vector.
unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
unsigned char *jptr = iptr + nChan - 1;
for (int ichan = 0; ichan < nChan/2; ichan++) {
unsigned char tmp = *iptr;
*(iptr++) = *jptr;
*(jptr--) = tmp;
}
}
} else {
// All channels are unflagged by default.
unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
for (int ichan = 0; ichan < nChan; ichan++) {
*(iptr++) = 0;
}
}
}
}
// Read cross-polarization data.
if (cGetXPol) {
int anynul;
for (int j = 0; j < 2; j++) {
mbrec.xcalfctr[0][j] = 0.0f;
}
if (cData[XCALFCTR].colnum > 0) {
readData(XCALFCTR, cRow, mbrec.xcalfctr);
}
blc[0] = 1;
trc[0] = 2;
blc[1] = std::min(startChan, endChan);
trc[1] = std::max(startChan, endChan);
blc[2] = cRow;
trc[2] = cRow;
int nAxis = 2;
long nAxes[] = {2, nChan};
if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
logMsg();
delete [] blc;
delete [] trc;
delete [] inc;
return 1;
}
if (endChan < startChan) {
// Invert the cross-polarization spectrum.
float *iptr = mbrec.xpol[0];
float *jptr = iptr + nChan - 2;
for (int ichan = 0; ichan < nChan/2; ichan++) {
float tmp = *iptr;
*iptr = *jptr;
*jptr = tmp;
tmp = *(iptr+1);
*(iptr+1) = *(jptr+1);
*(jptr+1) = tmp;
iptr += 2;
jptr -= 2;
}
}
}
delete [] blc;
delete [] trc;
delete [] inc;
if (cStatus) {
logMsg();
return 1;
}
mbrec.extraSysCal = cExtraSysCal;
readData(REFBEAM, cRow, &mbrec.refBeam);
readData(TCAL, cRow, &mbrec.tcal[0]);
readData(TCALTIME, cRow, mbrec.tcalTime);
readData(AZIMUTH, cRow, &mbrec.azimuth);
readData(ELEVATIO, cRow, &mbrec.elevation);
readData(PARANGLE, cRow, &mbrec.parAngle);
readData(FOCUSAXI, cRow, &mbrec.focusAxi);
readData(FOCUSTAN, cRow, &mbrec.focusTan);
readData(FOCUSROT, cRow, &mbrec.focusRot);
readData(TAMBIENT, cRow, &mbrec.temp);
readData(PRESSURE, cRow, &mbrec.pressure);
readData(HUMIDITY, cRow, &mbrec.humidity);
readData(WINDSPEE, cRow, &mbrec.windSpeed);
readData(WINDDIRE, cRow, &mbrec.windAz);
if (cALFA_BD) {
// ALFA BDFITS stores zenith angle rather than elevation.
mbrec.elevation = 90.0 - mbrec.elevation;
}
mbrec.azimuth *= D2R;
mbrec.elevation *= D2R;
mbrec.parAngle *= D2R;
mbrec.focusRot *= D2R;
mbrec.windAz *= D2R;
if (cStatus) {
logMsg();
return 1;
}
return 0;
}
//-------------------------------------------------------- SDFITSreader::close
// Close the SDFITS file.
void SDFITSreader::close()
{
if (cSDptr) {
int status = 0;
fits_close_file(cSDptr, &status);
cSDptr = 0x0;
if (cBeams) delete [] cBeams;
if (cIFs) delete [] cIFs;
if (cStartChan) delete [] cStartChan;
if (cEndChan) delete [] cEndChan;
if (cRefChan) delete [] cRefChan;
}
}
//------------------------------------------------------- SDFITSreader::logMsg
// Log a message. If the current CFITSIO status value is non-zero, also log
// the corresponding error message and the CFITSIO message stack.
void SDFITSreader::logMsg(const char *msg)
{
FITSreader::logMsg(msg);
if (cStatus > 0) {
fits_get_errstatus(cStatus, cMsg);
FITSreader::logMsg(cMsg);
while (fits_read_errmsg(cMsg)) {
FITSreader::logMsg(cMsg);
}
}
}
//----------------------------------------------------- SDFITSreader::findData
// Locate a data item in the SDFITS file.
void SDFITSreader::findData(
int iData,
char *name,
int type)
{
cData[iData].name = name;
cData[iData].type = type;
int colnum;
findCol(name, &colnum);
cData[iData].colnum = colnum;
// Determine the number of data elements.
if (colnum > 0) {
int coltype;
long nelem, width;
fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
fits_get_bcolparms(cSDptr, colnum, 0x0, cData[iData].units, 0x0, 0x0, 0x0,
0x0, 0x0, 0x0, &cStatus);
// Look for a TDIMnnn keyword or column.
char tdim[8];
sprintf(tdim, "TDIM%d", colnum);
findCol(tdim, &cData[iData].tdimcol);
if (coltype < 0) {
// CFITSIO returns coltype < 0 for variable length arrays.
cData[iData].coltype = -coltype;
cData[iData].nelem = -nelem;
} else {
cData[iData].coltype = coltype;
// Is there a TDIMnnn column?
if (cData[iData].tdimcol > 0) {
// Yes, dimensions of the fixed-length array could still vary.
cData[iData].nelem = -nelem;
} else {
cData[iData].nelem = nelem;
}
}
} else if (colnum == 0) {
// Keyword.
cData[iData].coltype = 0;
cData[iData].nelem = 1;
cData[iData].tdimcol = -1;
}
}
//------------------------------------------------------ SDFITSreader::findCol
// Locate a parameter in the SDFITS file.
void SDFITSreader::findCol(
char *name,
int *colnum)
{
*colnum = 0;
int status = 0;
fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
if (status) {
// Not a real column - maybe it's virtual.
char card[81];
status = 0;
fits_read_card(cSDptr, name, card, &status);
if (status) {
// Not virtual either.
*colnum = -1;
}
// Clear error messages.
fits_clear_errmsg();
}
}
//------------------------------------------------------ SDFITSreader::readDim
// Determine the dimensions of an array in the SDFITS file.
int SDFITSreader::readDim(
int iData,
long iRow,
int *naxes,
long naxis[])
{
int colnum = cData[iData].colnum;
if (colnum <= 0) {
return 1;
}
int maxdim = *naxes;
if (cData[iData].tdimcol < 0) {
// No TDIMnnn column for this array.
if (cData[iData].nelem < 0) {
// Variable length array; read the array descriptor.
*naxes = 1;
long dummy;
if (fits_read_descript(cSDptr, colnum, iRow, naxis, &dummy, &cStatus)) {
return 1;
}
} else {
// Read the repeat count from TFORMnnn.
if (fits_read_tdim(cSDptr, colnum, maxdim, naxes, naxis, &cStatus)) {
return 1;
}
}
} else {
// Read the TDIMnnn value from the header or table.
char tdim[8], tdimval[64];
sprintf(tdim, "TDIM%d", colnum);
readData(tdim, TSTRING, iRow, tdimval);
// fits_decode_tdim() checks that the TDIMnnn value is within the length
// of the array in the specified column number but unfortunately doesn't
// recognize variable-length arrays. Hence we must decode it here.
char *tp = tdimval;
if (*tp != '(') return 1;
tp++;
*naxes = 0;
for (size_t j = 1; j < strlen(tdimval); j++) {
if (tdimval[j] == ',' || tdimval[j] == ')') {
sscanf(tp, "%ld", naxis + (*naxes)++);
if (tdimval[j] == ')') break;
tp = tdimval + j + 1;
}
}
}
return 0;
}
//----------------------------------------------------- SDFITSreader::readParm
// Read a parameter value from the SDFITS file.
int SDFITSreader::readParm(
char *name,
int type,
void *value)
{
return readData(name, type, 1, value);
}
//----------------------------------------------------- SDFITSreader::readData
// Read a data value from the SDFITS file.
int SDFITSreader::readData(
char *name,
int type,
long iRow,
void *value)
{
int colnum;
findCol(name, &colnum);
if (colnum > 0 && iRow > 0) {
// Read the first value from the specified row of the table.
int coltype;
long nelem, width;
fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
int anynul;
if (type == TSTRING) {
if (nelem) {
fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
&cStatus);
} else {
strcpy((char *)value, "");
}
} else {
if (nelem) {
fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, value, &anynul,
&cStatus);
} else {
if (type == TSHORT) {
*((short *)value) = 0;
} else if (type == TINT) {
*((int *)value) = 0;
} else if (type == TFLOAT) {
*((float *)value) = 0.0f;
} else if (type == TDOUBLE) {
*((double *)value) = 0.0;
}
}
}
} else if (colnum == 0) {
// Read keyword value.
fits_read_key(cSDptr, type, name, value, 0, &cStatus);
} else {
// Not present.
if (type == TSTRING) {
strcpy((char *)value, "");
} else if (type == TSHORT) {
*((short *)value) = 0;
} else if (type == TINT) {
*((int *)value) = 0;
} else if (type == TFLOAT) {
*((float *)value) = 0.0f;
} else if (type == TDOUBLE) {
*((double *)value) = 0.0;
}
}
return colnum < 0;
}
//----------------------------------------------------- SDFITSreader::readData
// Read data from the SDFITS file.
int SDFITSreader::readData(
int iData,
long iRow,
void *value)
{
int type = cData[iData].type;
int colnum = cData[iData].colnum;
if (colnum > 0 && iRow > 0) {
// Read the required number of values from the specified row of the table.
long nelem = cData[iData].nelem;
int anynul;
if (type == TSTRING) {
if (nelem) {
fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
&cStatus);
} else {
strcpy((char *)value, "");
}
} else {
if (nelem) {
fits_read_col(cSDptr, type, colnum, iRow, 1, abs(nelem), 0, value,
&anynul, &cStatus);
} else {
if (type == TSHORT) {
*((short *)value) = 0;
} else if (type == TINT) {
*((int *)value) = 0;
} else if (type == TFLOAT) {
*((float *)value) = 0.0f;
} else if (type == TDOUBLE) {
*((double *)value) = 0.0;
}
}
}
} else if (colnum == 0) {
// Read keyword value.
char *name = cData[iData].name;
fits_read_key(cSDptr, type, name, value, 0, &cStatus);
} else {
// Not present.
if (type == TSTRING) {
strcpy((char *)value, "");
} else if (type == TSHORT) {
*((short *)value) = 0;
} else if (type == TINT) {
*((int *)value) = 0;
} else if (type == TFLOAT) {
*((float *)value) = 0.0f;
} else if (type == TDOUBLE) {
*((double *)value) = 0.0;
}
}
return colnum < 0;
}
//------------------------------------------------------ SDFITSreader::readCol
// Read a scalar column from the SDFITS file.
int SDFITSreader::readCol(
int iData,
void *value)
{
int type = cData[iData].type;
if (cData[iData].colnum > 0) {
// Table column.
int anynul;
fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRow, 0,
value, &anynul, &cStatus);
} else {
// Header keyword.
readData(iData, 0, value);
for (int irow = 1; irow < cNRow; irow++) {
if (type == TSHORT) {
((short *)value)[irow] = *((short *)value);
} else if (type == TINT) {
((int *)value)[irow] = *((int *)value);
} else if (type == TFLOAT) {
((float *)value)[irow] = *((float *)value);
} else if (type == TDOUBLE) {
((double *)value)[irow] = *((double *)value);
}
}
}
return cData[iData].colnum < 0;
}
//----------------------------------------------------- SDFITSreader::readTime
// Read the time from the SDFITS file.
int SDFITSreader::readTime(
long iRow,
int iPix,
char *datobs,
double &utc)
{
readData(DATE_OBS, iRow, datobs);
if (cData[TIME].colnum >= 0) {
readData(TIME, iRow, &utc);
} else if (cNAxisTime > 1) {
double timeDelt, timeRefPix, timeRefVal;
readData(TimeRefVal, iRow, &timeRefVal);
readData(TimeDelt, iRow, &timeDelt);
readData(TimeRefPix, iRow, &timeRefPix);
utc = timeRefVal + (iPix - timeRefPix) * timeDelt;
}
if (cALFA_BD) utc *= 3600.0;
// Check DATE-OBS format.
if (datobs[2] == '/') {
// Translate an old-format DATE-OBS.
datobs[9] = datobs[1];
datobs[8] = datobs[0];
datobs[2] = datobs[6];
datobs[5] = datobs[3];
datobs[3] = datobs[7];
datobs[6] = datobs[4];
datobs[7] = '-';
datobs[4] = '-';
datobs[1] = '9';
datobs[0] = '1';
} else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
// Dig UTC out of a new-format DATE-OBS.
int hh, mm;
float ss;
sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
utc = (hh*60 + mm)*60 + ss;
}
datobs[10] = '\0';
return 0;
}
//------------------------------------------------------ SDFITSreader::alfaCal
// Process ALFA calibration data.
int SDFITSreader::alfaCal(
short iBeam,
short iIF,
short iPol)
{
int calOn;
char chars[32];
if (cALFA_BD) {
readData("OBS_NAME", TSTRING, cRow, chars);
} else {
readData("SCANTYPE", TSTRING, cRow, chars);
}
if (strcmp(chars, "ON") == 0) {
calOn = 1;
} else if (strcmp(chars, "OFF") == 0) {
calOn = 0;
} else {
return 1;
}
// Read cal data.
long *blc = new long[cNAxes+1];
long *trc = new long[cNAxes+1];
long *inc = new long[cNAxes+1];
for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
blc[iaxis] = 1;
trc[iaxis] = 1;
inc[iaxis] = 1;
}
// User channel selection.
int startChan = cStartChan[iIF];
int endChan = cEndChan[iIF];
blc[cFreqAxis] = std::min(startChan, endChan);
trc[cFreqAxis] = std::max(startChan, endChan);
if (cALFA_CIMA > 1) {
// CIMAFITS 2.x has a legitimate STOKES axis...
blc[cStokesAxis] = iPol+1;
trc[cStokesAxis] = iPol+1;
} else {
// ...older ALFA data does not.
blc[cStokesAxis] = 1;
trc[cStokesAxis] = 1;
}
if (cTimeAxis >= 0) {
blc[cTimeAxis] = cTimeIdx;
trc[cTimeAxis] = cTimeIdx;
}
blc[cNAxes] = cRow;
trc[cNAxes] = cRow;
float spectrum[endChan];
int anynul;
if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
logMsg();
delete [] blc;
delete [] trc;
delete [] inc;
return 1;
}
// Factor to rescale according to the number of unblanked accumulations.
float factor = 1.0f;
if (cALFA_CIMA > 1) {
int colnum, naccum;
findCol("STAT", &colnum);
fits_read_col(cSDptr, TINT, colnum, cRow, 2, 1, 0, &naccum, &anynul,
&cStatus);
factor = cALFAacc / naccum;
}
// Average the spectrum.
float mean = 1e9f;
for (int k = 0; k < 2; k++) {
float discrim = 2.0f * mean;
int nChan = 0;
float sum = 0.0f;
float *chanN = spectrum + abs(endChan - startChan) + 1;
for (float *chan = spectrum; chan < chanN; chan++) {
// Simple discriminant that eliminates strong radar interference.
if (*chan < discrim) {
nChan++;
sum += *chan * factor;
}
}
mean = sum / nChan;
}
if (calOn) {
sALFAcalOn[iBeam][iPol] *= sALFAcalNon[iBeam][iPol];
sALFAcalOn[iBeam][iPol] += mean;
sALFAcalOn[iBeam][iPol] /= ++sALFAcalNon[iBeam][iPol];
} else {
sALFAcalOff[iBeam][iPol] *= sALFAcalNoff[iBeam][iPol];
sALFAcalOff[iBeam][iPol] += mean;
sALFAcalOff[iBeam][iPol] /= ++sALFAcalNoff[iBeam][iPol];
}
if (sALFAcalNon[iBeam][iPol] && sALFAcalNoff[iBeam][iPol]) {
// Tcal should come from the TCAL table, it varies weakly with beam,
// polarization, and frequency. However, TCAL is not written properly.
float Tcal = 12.0f;
sALFAcal[iBeam][iPol] = Tcal / (sALFAcalOn[iBeam][iPol] -
sALFAcalOff[iBeam][iPol]);
// Scale from K to Jy; the gain also varies weakly with beam,
// polarization, frequency, and zenith angle.
float fluxCal = 10.0f;
sALFAcal[iBeam][iPol] /= fluxCal;
}
return 0;
}
//----------------------------------------------------- SDFITSreader::alfaGain
// ALFA gain factor.
float SDFITSreader::alfaGain(
float zd)
{
// Gain vs zenith distance table from Robert Minchin, 2008/12/08.
const int nZD = 37;
const float zdLim[] = {1.5f, 19.5f};
const float zdInc = (nZD - 1) / (zdLim[1] - zdLim[0]);
float zdGain[] = { 1.00723708,
1.16644573, 1.15003645, 1.07117307, 1.02532673,
1.01788402, 1.01369524, 1.00000000, 0.989855111,
0.990888834, 0.993996620, 0.989964068, 0.982213855,
0.978662670, 0.979349494, 0.978478372, 0.974631131,
0.972126007, 0.972835243, 0.972742677, 0.968671739,
0.963891327, 0.963452935, 0.966831207, 0.969585896,
0.970700860, 0.972644389, 0.973754644, 0.967344403,
0.952168941, 0.937160134, 0.927843094, 0.914048433,
0.886700928, 0.864701211, 0.869126320, 0.854309499};
float gain;
// Do table lookup by linear interpolation.
float lambda = zdInc * (zd - zdLim[0]);
int j = int(lambda);
if (j < 0) {
gain = zdGain[0];
} else if (j >= nZD-1) {
gain = zdGain[nZD-1];
} else {
gain = zdGain[j] + (lambda - j) * (zdGain[j+1] - zdGain[j]);
}
return gain;
}