Changeset 1635 for trunk/external/atnf/PKSIO
- Timestamp:
- 09/25/09 11:47:11 (15 years ago)
- Location:
- trunk/external/atnf/PKSIO
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/external/atnf/PKSIO/FITSreader.h
r1452 r1635 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: FITSreader.h,v 19. 9 2008-11-17 06:28:04 cal103 Exp $29 //# $Id: FITSreader.h,v 19.10 2008-11-27 04:28:24 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The FITSreader class is an abstract base class for the Parkes Multibeam … … 96 96 // Coordinate systems are 97 97 // 0: equatorial (RA,Dec), 98 // 1: vertical (Az,El), 99 // 2: feed-plane. 98 // 1: horizontal (Az,El), 99 // 2: feed-plane, 100 // 3: zenithal position angle of feed and elevation, (ZPA,El). 100 101 int select( 101 102 const int startChan[], -
trunk/external/atnf/PKSIO/MBFITSreader.cc
r1452 r1635 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: MBFITSreader.cc,v 19.5 4 2008-11-17 06:51:55cal103 Exp $29 //# $Id: MBFITSreader.cc,v 19.55 2009-01-20 06:45:33 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes … … 141 141 int jstat = -3; 142 142 if (rpfitsin(jstat)) { 143 sprintf(cMsg, "ERROR: failed to open MBFITS file\n %s", rpname);143 sprintf(cMsg, "ERROR: Failed to open MBFITS file\n %s", rpname); 144 144 logMsg(cMsg); 145 145 return 1; … … 159 159 jstat = -1; 160 160 if (rpfitsin(jstat)) { 161 sprintf(cMsg, "ERROR: failed to read MBFITS header in file\n"161 sprintf(cMsg, "ERROR: Failed to read MBFITS header in file\n" 162 162 " %s", rpname); 163 163 logMsg(cMsg); … … 205 205 206 206 if (cNBeam <= 0) { 207 logMsg("ERROR , couldn't determine number of beams.");207 logMsg("ERROR: Couldn't determine number of beams."); 208 208 close(); 209 209 return 1; … … 337 337 // Read the first syscal record. 338 338 if (rpget(1, cEOS)) { 339 logMsg("ERROR , failed to read first syscal record.");339 logMsg("ERROR: Failed to read first syscal record."); 340 340 close(); 341 341 return 1; … … 373 373 { 374 374 if (!cMBopen) { 375 logMsg("ERROR , an MBFITS file has not been opened.");375 logMsg("ERROR: An MBFITS file has not been opened."); 376 376 return 1; 377 377 } … … 515 515 516 516 if (!cMBopen) { 517 logMsg("ERROR , an MBFITS file has not been opened.");517 logMsg("ERROR: An MBFITS file has not been opened."); 518 518 return 1; 519 519 } … … 588 588 cXpolOff = new int[cNIF]; 589 589 590 int simulIF = 0;591 590 int maxChan = 0; 592 591 int maxXpol = 0; 593 592 593 cSimulIF = 0; 594 594 for (int iIF = 0; iIF < cNIF; iIF++) { 595 595 if (cIFs[iIF]) { … … 617 617 618 618 // Maximum number of selected IFs in any simultaneous set. 619 simulIF = max(simulIF, cIFSel[iIF]+1);619 cSimulIF = max(cSimulIF, cIFSel[iIF]+1); 620 620 621 621 // Maximum memory required for any simultaneous set. … … 644 644 645 645 if (cNBin > 1 && cNBeamSel > 1) { 646 logMsg("ERROR, cannot handle binning mode for multiple beams."); 646 logMsg("ERROR: Cannot handle binning mode for multiple beams.\n" 647 " Select a single beam for input."); 647 648 close(); 648 649 return 1; … … 657 658 // Allocate memory for spectral arrays. 658 659 for (int ibuff = 0; ibuff < nBuff; ibuff++) { 659 cBuffer[ibuff].setNIFs( simulIF);660 cBuffer[ibuff].setNIFs(cSimulIF); 660 661 cBuffer[ibuff].allocate(0, maxChan, maxXpol); 661 662 662 663 // Signal that this IF in this buffer has been flushed. 663 for (int iIF = 0; iIF < simulIF; iIF++) {664 for (int iIF = 0; iIF < cSimulIF; iIF++) { 664 665 cBuffer[ibuff].IFno[iIF] = 0; 665 666 } … … 679 680 cCycleNo = 0; 680 681 cPrevUTC = -1.0; 682 } 683 684 // Apply beam and IF selection before the change-of-day test to allow 685 // a single selected beam and IF to be handled in binning-mode. 686 beamNo = int(cBaseline / 256.0); 687 if (beamNo == 1) { 688 // Store the position of beam 1 for grid convergence corrections. 689 cRA0 = cU; 690 cDec0 = cV; 691 } 692 iBeamSel = cBeamSel[beamNo-1]; 693 if (iBeamSel < 0) continue; 694 695 // Sanity check (mainly for MOPS). 696 if (cIFno > cNIF) continue; 697 698 // Apply IF selection. 699 iIFSel = cIFSel[cIFno - 1]; 700 if (iIFSel < 0) continue; 701 702 703 if (cNBin > 1) { 704 // Binning mode: correct the time. 705 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0); 681 706 } 682 707 … … 689 714 // the start of the next. 690 715 #ifdef PKSIO_DEBUG 691 fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f ",716 fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f\n", 692 717 cPrevUTC, cUTC); 693 718 #endif … … 707 732 } 708 733 709 if (cNBin > 1) {710 // Binning mode: correct the time.711 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);712 }713 714 734 // New integration cycle? 715 735 if ((cUTC+cod) > cPrevUTC) { … … 717 737 cPrevUTC = cUTC + 0.0001; 718 738 } 719 720 // Apply beam selection.721 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 }727 iBeamSel = cBeamSel[beamNo-1];728 if (iBeamSel < 0) continue;729 730 // Sanity check (mainly for MOPS).731 if (cIFno > cNIF) continue;732 733 // Apply IF selection.734 iIFSel = cIFSel[cIFno - 1];735 if (iIFSel < 0) continue;736 739 737 740 sprintf(cDateObs, "%-10.10s", names_.datobs); … … 784 787 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin; 785 788 786 // iMBuff->nIF is set to zero (below) to signal that all IFs in787 // an integration have been flushed.789 // iMBuff->nIF is decremented (below) and if zero signals that all 790 // IFs in an integration have been flushed. 788 791 if (iMBuff->nIF) { 789 792 if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) { … … 805 808 806 809 // Find the IF to flush. 807 for (; cFlushIF < iMBuff->nIF; cFlushIF++) {810 for (; cFlushIF < cSimulIF; cFlushIF++) { 808 811 if (iMBuff->IFno[cFlushIF]) break; 809 812 } … … 827 830 828 831 829 if (c Flushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {832 if (cInterp && cFlushing == 1) { 830 833 // Start of flush cycle, interpolate the beam position. 831 834 // … … 1098 1101 #endif 1099 1102 } 1103 1104 cFlushing = 2; 1100 1105 } 1101 1106 … … 1114 1119 iMBuff->IFno[cFlushIF] = 0; 1115 1120 1116 if (cFlushIF == iMBuff->nIF - 1) { 1117 // Signal that all IFs in this buffer location have been flushed. 1118 iMBuff->nIF = 0; 1119 1120 // Stop cEOS being set when the next integration is read. 1121 iMBuff->nIF--; 1122 if (iMBuff->nIF == 0) { 1123 // All IFs in this buffer location have been flushed. Stop cEOS 1124 // being set when the next integration is read. 1121 1125 iMBuff->cycleNo = 0; 1122 1126 -
trunk/external/atnf/PKSIO/MBFITSreader.h
r1452 r1635 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: MBFITSreader.h,v 19.2 1 2008-11-17 06:33:10cal103 Exp $29 //# $Id: MBFITSreader.h,v 19.22 2009-01-20 06:36:03 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes … … 116 116 char cDateObs[12]; 117 117 int *cBeamSel, *cChanOff, cFirst, *cIFSel, cInterp, cIntTime, cMBopen, 118 cMopra, cNBeamSel, cNBin, cRetry, cS Upos, *cXpolOff;118 cMopra, cNBeamSel, cNBin, cRetry, cSimulIF, cSUpos, *cXpolOff; 119 119 120 120 // The data has to be bufferred to allow positions to be interpolated. -
trunk/external/atnf/PKSIO/MBrecord.cc
r1452 r1635 2 2 //# MBrecord.cc: Class to store an MBFITS single-dish data record. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Mark Calabretta, ATNF 6 6 //# … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: MBrecord.cc,v 19.1 2 2008-11-17 06:52:35cal103 Exp $29 //# $Id: MBrecord.cc,v 19.13 2009-03-24 06:15:33 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBrecord class stores an MBFITS single-dish data record. … … 93 93 xcalfctr = new float[nif][2]; 94 94 baseLin = new float[nif][2][2]; 95 baseSub = new float[nif][2][ 9];95 baseSub = new float[nif][2][24]; 96 96 spectra = new float*[nif]; 97 97 flagged = new unsigned char*[nif]; … … 260 260 baseLin[iIF][ipol][1] = other.baseLin[iIF][ipol][1]; 261 261 262 for (int j = 0; j < 9; j++) {262 for (int j = 0; j < 24; j++) { 263 263 baseSub[iIF][ipol][j] = other.baseSub[iIF][ipol][j]; 264 264 } … … 380 380 baseLin[0][ipol][1] = other.baseLin[iIF][ipol][1]; 381 381 382 for (int j = 0; j < 9; j++) {382 for (int j = 0; j < 24; j++) { 383 383 baseSub[0][ipol][j] = other.baseSub[iIF][ipol][j]; 384 384 } -
trunk/external/atnf/PKSIO/MBrecord.h
r1452 r1635 2 2 //# MBrecord.h: Class to store an MBFITS single-dish data record. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Mark Calabretta, ATNF 6 6 //# … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: MBrecord.h,v 19.1 4 2008-11-17 06:36:12cal103 Exp $29 //# $Id: MBrecord.h,v 19.15 2009-03-24 06:15:33 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBrecord class stores an MBFITS single-dish data record. … … 134 134 int haveBase; // Are baseline parameters present? 135 135 float (*baseLin)[2][2]; // Linear baseline fit for each polarization. 136 float (*baseSub)[2][ 9]; // Polynomial baseline subtracted.136 float (*baseSub)[2][24]; // Polynomial baseline subtracted. 137 137 int haveSpectra; // Is spectral data present? 138 138 float* *spectra; // Spectra for each polarization, Jy. -
trunk/external/atnf/PKSIO/PKSFITSreader.cc
r1452 r1635 1 1 //# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file. 2 2 //#--------------------------------------------------------------------------- 3 //# Copyright (C) 2000-200 83 //# Copyright (C) 2000-2009 4 4 //# Associated Universities, Inc. Washington DC, USA. 5 5 //# … … 25 25 //# Charlottesville, VA 22903-2475 USA 26 26 //# 27 //# $Id: PKSFITSreader.cc,v 19.2 1 2008-11-17 06:54:25cal103 Exp $27 //# $Id: PKSFITSreader.cc,v 19.22 2009-03-24 06:15:33 cal103 Exp $ 28 28 //#--------------------------------------------------------------------------- 29 29 //# Original: 2000/08/02, Mark Calabretta, ATNF … … 446 446 if (cMBrec.haveBase) { 447 447 pksrec.baseLin.resize(2,nPol); 448 pksrec.baseSub.resize( 9,nPol);448 pksrec.baseSub.resize(24,nPol); 449 449 450 450 for (uInt ipol = 0; ipol < nPol; ipol++) { … … 452 452 pksrec.baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1]; 453 453 454 for (uInt j = 0; j < 9; j++) {454 for (uInt j = 0; j < 24; j++) { 455 455 pksrec.baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j]; 456 456 } -
trunk/external/atnf/PKSIO/PKSMS2reader.cc
r1452 r1635 2 2 //# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Associated Universities, Inc. Washington DC, USA. 6 6 //# … … 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: PKSMS2reader.cc,v 19.2 1 2008-11-17 06:55:18cal103 Exp $28 //# $Id: PKSMS2reader.cc,v 19.22 2009-03-24 06:15:33 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 //# Original: 2000/08/03, Mark Calabretta, ATNF … … 656 656 cBaseLinCol.get(cIdx, pksrec.baseLin); 657 657 658 pksrec.baseSub.resize( 9,cNPol(iIF));658 pksrec.baseSub.resize(24,cNPol(iIF)); 659 659 cBaseSubCol.get(cIdx, pksrec.baseSub); 660 660 -
trunk/external/atnf/PKSIO/PKSMS2writer.cc
r1452 r1635 2 2 //# PKSMS2writer.cc: Class to write Parkes multibeam data to a measurementset. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Associated Universities, Inc. Washington DC, USA. 6 6 //# … … 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: PKSMS2writer.cc,v 19.1 4 2008-11-17 06:56:13 cal103 Exp $28 //# $Id: PKSMS2writer.cc,v 19.15 2009-03-24 06:15:33 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 … … 123 123 IPosition(2,2,maxNPol), ColumnDesc::Direct)); 124 124 pksDesc.addColumn(ArrayColumnDesc<Float>("BASESUB", "Baseline subtracted", 125 IPosition(2, 9,maxNPol), ColumnDesc::Direct));125 IPosition(2,24,maxNPol), ColumnDesc::Direct)); 126 126 } 127 127 … … 469 469 cBaseLinCol->put(irow, pksrec.baseLin); 470 470 471 if (pksrec.baseSub.nrow() == 9) {471 if (pksrec.baseSub.nrow() == 24) { 472 472 cBaseSubCol->put(irow, pksrec.baseSub); 473 473 474 474 } else { 475 Matrix<Float> tmp( 9, 2, 0.0f);475 Matrix<Float> tmp(24, 2, 0.0f); 476 476 for (Int ipol = 0; ipol < nPol; ipol++) { 477 477 for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) { -
trunk/external/atnf/PKSIO/PKSSDwriter.cc
r1452 r1635 2 2 //# PKSSDwriter.cc: Class to write Parkes multibeam data to an SDFITS file. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Associated Universities, Inc. Washington DC, USA. 6 6 //# … … 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: PKSSDwriter.cc,v 19.1 5 2008-11-17 06:56:50cal103 Exp $28 //# $Id: PKSSDwriter.cc,v 19.16 2009-03-24 06:15:33 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 … … 239 239 mbrec.baseSub[0][ipol][j] = pksrec.baseSub(j,ipol); 240 240 } 241 for (uInt j = pksrec.baseSub.nrow(); j < 9; j++) {241 for (uInt j = pksrec.baseSub.nrow(); j < 24; j++) { 242 242 mbrec.baseSub[0][ipol][j] = 0.0f; 243 243 } -
trunk/external/atnf/PKSIO/PKSmsg.cc
r1466 r1635 36 36 37 37 #include <string> 38 #include <cstring>39 38 40 39 //------------------------------------------------------------- PKSmsg::PKSmsg -
trunk/external/atnf/PKSIO/PKSreader.h
r1452 r1635 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: PKSreader.h,v 19.2 2 2008-11-17 06:44:34 cal103 Exp $28 //# $Id: PKSreader.h,v 19.23 2008-11-27 04:28:24 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 //# Original: 2000/08/02, Mark Calabretta, ATNF … … 134 134 // Coordinate system selection (only supported for SDFITS input): 135 135 // 0: equatorial (RA,Dec), 136 // 1: vertical (Az,El), 137 // 2: feed-plane. 136 // 1: horizontal (Az,El), 137 // 2: feed-plane, 138 // 3: zenithal position angle of feed and elevation, (ZPA,El). 138 139 virtual uInt select( 139 140 const Vector<Bool> beamSel, -
trunk/external/atnf/PKSIO/SDFITSreader.cc
r1509 r1635 1 1 //#--------------------------------------------------------------------------- 2 //# SDFITSreader.cc: ATNF CFITSIO interface class for SDFITS input.2 //# SDFITSreader.cc: ATNF interface class for SDFITS input using CFITSIO. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Associated Universities, Inc. Washington DC, USA. 6 6 //# … … 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: SDFITSreader.cc,v 19. 33 2008-11-17 06:58:34cal103 Exp $28 //# $Id: SDFITSreader.cc,v 19.43 2009-05-06 03:28:17 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 //# The SDFITSreader class reads single dish FITS files such as those written … … 41 41 #include <casa/math.h> 42 42 #include <casa/stdio.h> 43 #include <cstring>44 43 45 44 #include <algorithm> … … 64 63 const double D2R = PI / 180.0; 65 64 65 //---------------------------------------------------- SDFITSreader::(statics) 66 67 int SDFITSreader::sInit = 1; 68 int SDFITSreader::sReset = 0; 69 int (*SDFITSreader::sALFAcalNon)[2] = (int (*)[2])(new float[16]); 70 int (*SDFITSreader::sALFAcalNoff)[2] = (int (*)[2])(new float[16]); 71 float (*SDFITSreader::sALFAcalOn)[2] = (float (*)[2])(new float[16]); 72 float (*SDFITSreader::sALFAcalOff)[2] = (float (*)[2])(new float[16]); 73 float (*SDFITSreader::sALFAcal)[2] = (float (*)[2])(new float[16]); 74 66 75 //------------------------------------------------- SDFITSreader::SDFITSreader 67 76 … … 69 78 { 70 79 // Default constructor. 71 cSDptr = 0 ;80 cSDptr = 0x0; 72 81 73 82 // Allocate space for data descriptors. … … 157 166 // Arecibo ALFA data of some kind. 158 167 cALFA = 1; 159 for (int iBeam = 0; iBeam < 8; iBeam++) { 160 for (int iPol = 0; iPol < 2; iPol++) { 161 cALFAcalOn[iBeam][iPol] = 0.0f; 162 cALFAcalOff[iBeam][iPol] = 0.0f; 163 164 // Nominal factor to calibrate spectra in Jy. 165 cALFAcal[iBeam][iPol] = 3.0f; 166 } 168 if (sInit) { 169 for (int iBeam = 0; iBeam < 8; iBeam++) { 170 for (int iPol = 0; iPol < 2; iPol++) { 171 sALFAcalOn[iBeam][iPol] = 0.0f; 172 sALFAcalOff[iBeam][iPol] = 0.0f; 173 174 // Nominal factor to calibrate spectra in Jy. 175 sALFAcal[iBeam][iPol] = 3.0f; 176 } 177 } 178 179 sInit = 0; 167 180 } 168 181 } … … 174 187 strncmp(telescope, "NRAO_GBT", 8) == 0; 175 188 176 cRow = 0;177 178 189 179 190 // Check that the DATA array column is present. … … 181 192 haveSpectra = cHaveSpectra = cData[DATA].colnum > 0; 182 193 194 cNAxisTime = 0; 183 195 if (cHaveSpectra) { 184 196 // Find the number of data axes (must be the same for each IF). 185 cNAx is = 5;186 if (readDim(DATA, 1, &cNAx is, cNAxes)) {197 cNAxes = 5; 198 if (readDim(DATA, 1, &cNAxes, cNAxis)) { 187 199 logMsg(); 188 200 close(); … … 193 205 // ALFA BDFITS: variable length arrays don't actually vary and there is 194 206 // no TDIM (or MAXISn) card; use the LAGS_IN value. 195 cNAx is = 5;196 readParm("LAGS_IN", TLONG, cNAx es);197 cNAx es[1] = 1;198 cNAx es[2] = 1;199 cNAx es[3] = 1;200 cNAx es[4] = 1;201 cData[DATA].nelem = cNAx es[0];202 } 203 204 if (cNAx is < 4) {207 cNAxes = 5; 208 readParm("LAGS_IN", TLONG, cNAxis); 209 cNAxis[1] = 1; 210 cNAxis[2] = 1; 211 cNAxis[3] = 1; 212 cNAxis[4] = 1; 213 cData[DATA].nelem = cNAxis[0]; 214 } 215 216 if (cNAxes < 4) { 205 217 // Need at least four axes (for now). 206 218 logMsg("ERROR: DATA array contains fewer than four axes."); 207 219 close(); 208 220 return 1; 209 } else if (cNAx is > 5) {221 } else if (cNAxes > 5) { 210 222 // We support up to five axes. 211 223 logMsg("ERROR: DATA array contains more than five axes."); … … 229 241 readParm("DATAXED", TSTRING, dataxed); 230 242 231 for (int iaxis = 0; iaxis < 5; iaxis++) cNAx es[iaxis] = 0;232 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAx es, cNAxes+1, cNAxes+2,233 cNAx es+3, cNAxes+4);243 for (int iaxis = 0; iaxis < 5; iaxis++) cNAxis[iaxis] = 0; 244 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1, cNAxis+2, 245 cNAxis+3, cNAxis+4); 234 246 for (int iaxis = 4; iaxis > -1; iaxis--) { 235 if (cNAx es[iaxis] == 0) cNAxis = iaxis;247 if (cNAxis[iaxis] == 0) cNAxes = iaxis; 236 248 } 237 249 } … … 244 256 // Find required DATA array axes. 245 257 char ctype[5][72]; 246 for (int iaxis = 0; iaxis < cNAx is; iaxis++) {258 for (int iaxis = 0; iaxis < cNAxes; iaxis++) { 247 259 strcpy(ctype[iaxis], ""); 248 260 readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]); // Core. … … 255 267 } 256 268 257 char *fqCRPIX = 0;258 269 char *fqCRVAL = 0; 259 270 char *fqCDELT = 0; 271 char *fqCRPIX = 0; 260 272 char *raCRVAL = 0; 261 273 char *decCRVAL = 0; 262 274 char *timeCRVAL = 0; 275 char *timeCDELT = 0; 276 char *timeCRPIX = 0; 263 277 char *beamCRVAL = 0; 264 278 265 for (int iaxis = 0; iaxis < cNAxis; iaxis++) { 279 cFreqAxis = -1; 280 cStokesAxis = -1; 281 cRaAxis = -1; 282 cDecAxis = -1; 283 cTimeAxis = -1; 284 cBeamAxis = -1; 285 286 for (int iaxis = 0; iaxis < cNAxes; iaxis++) { 266 287 if (strncmp(ctype[iaxis], "FREQ", 4) == 0) { 267 c Reqax[0]= iaxis;268 fqCR PIX = CRPIX[iaxis];269 fqC RVAL = CRVAL[iaxis];270 fqC DELT = CDELT[iaxis];288 cFreqAxis = iaxis; 289 fqCRVAL = CRVAL[iaxis]; 290 fqCDELT = CDELT[iaxis]; 291 fqCRPIX = CRPIX[iaxis]; 271 292 272 293 } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) { 273 c Reqax[1]= iaxis;294 cStokesAxis = iaxis; 274 295 275 296 } else if (strncmp(ctype[iaxis], "RA", 2) == 0) { 276 cR eqax[2]= iaxis;277 raCRVAL = CRVAL[iaxis];297 cRaAxis = iaxis; 298 raCRVAL = CRVAL[iaxis]; 278 299 279 300 } else if (strncmp(ctype[iaxis], "DEC", 3) == 0) { 280 c Reqax[3]= iaxis;281 decCRVAL = CRVAL[iaxis];301 cDecAxis = iaxis; 302 decCRVAL = CRVAL[iaxis]; 282 303 283 304 } else if (strcmp(ctype[iaxis], "TIME") == 0) { 284 // TIME (UTC seconds since midnight) can be a keyword or axis type. 305 // TIME (UTC seconds since midnight); axis type, if present, takes 306 // precedence over keyword. 307 cTimeAxis = iaxis; 285 308 timeCRVAL = CRVAL[iaxis]; 309 310 // Check for non-degeneracy. 311 if ((cNAxisTime = cNAxis[iaxis]) > 1) { 312 timeCDELT = CDELT[iaxis]; 313 timeCRPIX = CRPIX[iaxis]; 314 sprintf(cMsg, "DATA array contains a TIME axis of length %ld.", 315 cNAxisTime); 316 logMsg(cMsg); 317 } 286 318 287 319 } else if (strcmp(ctype[iaxis], "BEAM") == 0) { 288 320 // BEAM can be a keyword or axis type. 321 cBeamAxis = iaxis; 289 322 beamCRVAL = CRVAL[iaxis]; 290 323 } … … 293 326 if (cALFA_BD) { 294 327 // Fixed in ALFA CIMAFITS. 295 cR eqax[2]= 2;328 cRaAxis = 2; 296 329 raCRVAL = "CRVAL2A"; 297 330 298 c Reqax[3]= 3;331 cDecAxis = 3; 299 332 decCRVAL = "CRVAL3A"; 300 333 } 301 334 302 // Check that all are present. 303 for (int iaxis = 0; iaxis < 4; iaxis++) { 304 if (cReqax[iaxis] < 0) { 305 logMsg("ERROR: Could not find required DATA array axes."); 306 close(); 307 return 1; 308 } 335 // Check that required axes are present. 336 if (cFreqAxis < 0 || cStokesAxis < 0 || cRaAxis < 0 || cDecAxis < 0) { 337 logMsg("ERROR: Could not find required DATA array axes."); 338 close(); 339 return 1; 309 340 } 310 341 … … 313 344 findData(CYCLE, "CYCLE", TINT); // Additional. 314 345 findData(DATE_OBS, "DATE-OBS", TSTRING); // Core. 315 findData(TIME, "TIME", TDOUBLE); // Core. 346 347 if (cTimeAxis >= 0) { 348 // The DATA array has a TIME axis. 349 if (cNAxisTime > 1) { 350 // Non-degenerate. 351 findData(TimeRefVal, timeCRVAL, TDOUBLE); // Time reference value. 352 findData(TimeDelt, timeCDELT, TDOUBLE); // Time increment. 353 findData(TimeRefPix, timeCRPIX, TFLOAT); // Time reference pixel. 354 } else { 355 // Degenerate, treat its like a simple TIME keyword. 356 findData(TIME, timeCRVAL, TDOUBLE); 357 } 358 359 } else { 360 findData(TIME, "TIME", TDOUBLE); // Core. 361 } 362 316 363 findData(EXPOSURE, "EXPOSURE", TFLOAT); // Core. 317 364 findData(OBJECT, "OBJECT", TSTRING); // Core. … … 323 370 findData(BEAM, "BEAM", TSHORT); // Additional. 324 371 findData(IF, "IF", TSHORT); // Additional. 325 findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel.326 372 findData(FqRefVal, fqCRVAL, TDOUBLE); // Frequency reference value. 327 373 findData(FqDelt, fqCDELT, TDOUBLE); // Frequency increment. 374 findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel. 328 375 findData(RA, raCRVAL, TDOUBLE); // Right ascension. 329 376 findData(DEC, decCRVAL, TDOUBLE); // Declination. … … 367 414 findData(SCAN, "SCAN_ID", TINT); 368 415 if (cALFA_CIMA > 1) { 416 // Note that RECNUM increases by cNAxisTime per row. 369 417 findData(CYCLE, "RECNUM", TINT); 370 418 } else { … … 410 458 } 411 459 cIF_1rel = 0; 412 }413 414 if (cData[TIME].colnum < 0) {415 if (timeCRVAL) {416 // There is a TIME axis.417 findData(TIME, timeCRVAL, TDOUBLE);418 }419 460 } 420 461 … … 607 648 if (cData[DATA].nelem < 0) { 608 649 // Variable dimension array. 609 if (readDim(DATA, irow+1, &cNAx is, cNAxes)) {650 if (readDim(DATA, irow+1, &cNAxes, cNAxis)) { 610 651 logMsg(); 611 652 close(); … … 619 660 readParm("DATAXED", TSTRING, dataxed); 620 661 621 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAx es, cNAxes+1,622 cNAx es+2, cNAxes+3, cNAxes+4);662 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1, 663 cNAxis+2, cNAxis+3, cNAxis+4); 623 664 } 624 665 } 625 666 626 667 // Number of channels and polarizations. 627 cNChan[iIF] = cNAx es[cReqax[0]];628 cNPol[iIF] = cNAx es[cReqax[1]];668 cNChan[iIF] = cNAxis[cFreqAxis]; 669 cNPol[iIF] = cNAxis[cStokesAxis]; 629 670 cHaveXPol[iIF] = 0; 630 671 … … 665 706 666 707 // Number of channels and polarizations. 667 cNChan[0] = cNAx es[cReqax[0]];668 cNPol[0] = cNAx es[cReqax[1]];708 cNChan[0] = cNAxis[cFreqAxis]; 709 cNPol[0] = cNAxis[cStokesAxis]; 669 710 cHaveXPol[0] = 0; 670 711 } … … 726 767 727 768 extraSysCal = cExtraSysCal; 769 770 771 // Extras for ALFA data. 772 cALFAacc = 0.0f; 773 if (cALFA_CIMA > 1) { 774 // FFTs per second when the Mock correlator operates in RFI blanking mode. 775 readData("PHFFTACC", TFLOAT, 0, &cALFAacc); 776 } 777 778 779 cRow = 0; 780 cTimeIdx = cNAxisTime; 728 781 729 782 return 0; … … 825 878 // No, try digging it out of the CTYPE card (AIPS convention). 826 879 char keyw[9], ctype[9]; 827 sprintf(keyw, "CTYPE%ld", c Reqax[0]+1);880 sprintf(keyw, "CTYPE%ld", cFreqAxis+1); 828 881 readParm(keyw, TSTRING, ctype); 829 882 … … 860 913 861 914 // Get parameters from first row of table. 862 readData(DATE_OBS, 1, datobs); 863 readData(TIME, 1, &utc); 915 readTime(1, 1, datobs, utc); 864 916 readData(FqRefVal, 1, &refFreq); 865 917 readParm("BANDWID", TDOUBLE, &bandwidth); // Core. 866 867 if (cALFA_BD) utc *= 3600.0;868 918 869 919 if (cStatus) { 870 920 logMsg(); 871 921 return 1; 872 }873 874 // Check DATE-OBS format.875 if (datobs[2] == '/') {876 // Translate an old-format DATE-OBS.877 datobs[9] = datobs[1];878 datobs[8] = datobs[0];879 datobs[2] = datobs[6];880 datobs[5] = datobs[3];881 datobs[3] = datobs[7];882 datobs[6] = datobs[4];883 datobs[7] = '-';884 datobs[4] = '-';885 datobs[1] = '9';886 datobs[0] = '1';887 datobs[10] = '\0';888 889 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {890 // Dig UTC out of a new-format DATE-OBS.891 int hh, mm;892 float ss;893 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);894 utc = (hh*60 + mm)*60 + ss;895 datobs[10] = '\0';896 922 } 897 923 … … 993 1019 994 1020 // Find the number of rows selected. 995 short *sel = new short[ nRow];996 for (int irow = 0; irow < nRow; irow++) {1021 short *sel = new short[cNRow]; 1022 for (int irow = 0; irow < cNRow; irow++) { 997 1023 sel[irow] = 1; 998 1024 } … … 1010 1036 } 1011 1037 1012 for (int irow = 0; irow < nRow; irow++) {1038 for (int irow = 0; irow < cNRow; irow++) { 1013 1039 if (!cBeams[beamCol[irow]-cBeam_1rel]) { 1014 1040 sel[irow] = 0; … … 1030 1056 } 1031 1057 1032 for (int irow = 0; irow < nRow; irow++) {1058 for (int irow = 0; irow < cNRow; irow++) { 1033 1059 if (!cIFs[IFCol[irow]-cIF_1rel]) { 1034 1060 sel[irow] = 0; … … 1040 1066 1041 1067 nSel = 0; 1042 for (int irow = 0; irow < nRow; irow++) {1068 for (int irow = 0; irow < cNRow; irow++) { 1043 1069 nSel += sel[irow]; 1044 1070 } … … 1046 1072 1047 1073 // Find the time range assuming the data is in chronological order. 1048 readData(DATE_OBS, 1, dateSpan[0]); 1049 readData(DATE_OBS, nRow, dateSpan[1]); 1050 readData(TIME, 1, utcSpan); 1051 readData(TIME, nRow, utcSpan+1); 1052 1053 if (cALFA_BD) { 1054 utcSpan[0] *= 3600.0; 1055 utcSpan[1] *= 3600.0; 1056 } 1057 1058 // Check DATE-OBS format. 1059 for (int i = 0; i < 2; i++) { 1060 if (dateSpan[0][2] == '/') { 1061 // Translate an old-format DATE-OBS. 1062 dateSpan[i][9] = dateSpan[i][1]; 1063 dateSpan[i][8] = dateSpan[i][0]; 1064 dateSpan[i][2] = dateSpan[i][6]; 1065 dateSpan[i][5] = dateSpan[i][3]; 1066 dateSpan[i][3] = dateSpan[i][7]; 1067 dateSpan[i][6] = dateSpan[i][4]; 1068 dateSpan[i][7] = '-'; 1069 dateSpan[i][4] = '-'; 1070 dateSpan[i][1] = '9'; 1071 dateSpan[i][0] = '1'; 1072 dateSpan[i][10] = '\0'; 1073 } 1074 1075 if (dateSpan[i][10] == 'T' && cData[TIME].colnum < 0) { 1076 // Dig UTC out of a new-format DATE-OBS. 1077 int hh, mm; 1078 float ss; 1079 sscanf(dateSpan[i]+11, "%d:%d:%f", &hh, &mm, &ss); 1080 utcSpan[i] = (hh*60 + mm)*60 + ss; 1081 } 1082 } 1074 readTime(1, 1, dateSpan[0], utcSpan[0]); 1075 readTime(cNRow, cNAxisTime, dateSpan[1], utcSpan[1]); 1083 1076 1084 1077 … … 1088 1081 1089 1082 if (cCoordSys == 1) { 1090 // Vertical (Az,El).1091 if (cData[AZIMUTH].colnum < 1||1092 cData[ELEVATIO].colnum < 1) {1083 // Horizontal (Az,El). 1084 if (cData[AZIMUTH].colnum < 0 || 1085 cData[ELEVATIO].colnum < 0) { 1093 1086 logMsg("WARNING: Azimuth/elevation information absent."); 1094 1087 cStatus = -1; … … 1097 1090 float *az = new float[cNRow]; 1098 1091 float *el = new float[cNRow]; 1099 fits_read_col(cSDptr, TFLOAT, cData[AZIMUTH].colnum, 1, 1, nRow, 0, az, 1100 &anynul, &cStatus); 1101 fits_read_col(cSDptr, TFLOAT, cData[ELEVATIO].colnum, 1, 1, nRow, 0, el, 1102 &anynul, &cStatus); 1092 readCol(AZIMUTH, az); 1093 readCol(ELEVATIO, el); 1103 1094 1104 1095 if (!cStatus) { 1105 for (int irow = 0; irow < nRow; irow++) {1096 for (int irow = 0; irow < cNRow; irow++) { 1106 1097 if (sel[irow]) { 1107 1098 positions[isel++] = az[irow] * D2R; … … 1115 1106 } 1116 1107 1108 } else if (cCoordSys == 3) { 1109 // ZPA-EL. 1110 if (cData[BEAM].colnum < 0 || 1111 cData[FOCUSROT].colnum < 0 || 1112 cData[ELEVATIO].colnum < 0) { 1113 logMsg("WARNING: ZPA/elevation information absent."); 1114 cStatus = -1; 1115 1116 } else { 1117 short *beam = new short[cNRow]; 1118 float *rot = new float[cNRow]; 1119 float *el = new float[cNRow]; 1120 readCol(BEAM, beam); 1121 readCol(FOCUSROT, rot); 1122 readCol(ELEVATIO, el); 1123 1124 if (!cStatus) { 1125 for (int irow = 0; irow < cNRow; irow++) { 1126 if (sel[irow]) { 1127 Int beamNo = beam[irow]; 1128 Double zpa = rot[irow]; 1129 if (beamNo > 1) { 1130 // Beam geometry for the Parkes multibeam. 1131 if (beamNo < 8) { 1132 zpa += -60.0 + 60.0*(beamNo-2); 1133 } else { 1134 zpa += -90.0 + 60.0*(beamNo-8); 1135 } 1136 1137 if (zpa < -180.0) { 1138 zpa += 360.0; 1139 } else if (zpa > 180.0) { 1140 zpa -= 360.0; 1141 } 1142 } 1143 1144 positions[isel++] = zpa * D2R; 1145 positions[isel++] = el[irow] * D2R; 1146 } 1147 } 1148 } 1149 1150 delete [] beam; 1151 delete [] rot; 1152 delete [] el; 1153 } 1154 1117 1155 } else { 1118 1156 double *ra = new double[cNRow]; 1119 1157 double *dec = new double[cNRow]; 1120 fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum, 1, 1, nRow, 0, ra, 1121 &anynul, &cStatus); 1122 fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec, 1123 &anynul, &cStatus); 1158 readCol(RA, ra); 1159 readCol(DEC, dec); 1160 1124 1161 if (cStatus) { 1125 1162 delete [] ra; … … 1129 1166 1130 1167 if (cALFA_BD) { 1131 for (int irow = 0; irow < nRow; irow++) {1168 for (int irow = 0; irow < cNRow; irow++) { 1132 1169 // Convert hours to degrees. 1133 1170 ra[irow] *= 15.0; … … 1137 1174 if (cCoordSys == 0) { 1138 1175 // Equatorial (RA,Dec). 1139 for (int irow = 0; irow < nRow; irow++) {1176 for (int irow = 0; irow < cNRow; irow++) { 1140 1177 if (sel[irow]) { 1141 1178 positions[isel++] = ra[irow] * D2R; … … 1148 1185 if (cData[OBJ_RA].colnum < 0 || 1149 1186 cData[OBJ_DEC].colnum < 0 || 1150 cData[PARANGLE].colnum < 1||1151 cData[FOCUSROT].colnum < 1) {1187 cData[PARANGLE].colnum < 0 || 1188 cData[FOCUSROT].colnum < 0) { 1152 1189 logMsg("WARNING: Insufficient information to compute feed-plane\n" 1153 1190 " coordinates."); … … 1160 1197 float *rot = new float[cNRow]; 1161 1198 1162 if (cData[OBJ_RA].colnum == 0) { 1163 // Header keyword. 1164 readData(OBJ_RA, 0, srcRA); 1165 for (int irow = 1; irow < nRow; irow++) { 1166 srcRA[irow] = *srcRA; 1167 } 1168 } else { 1169 // Table column. 1170 fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum, 1, 1, nRow, 1171 0, srcRA, &anynul, &cStatus); 1172 } 1173 1174 if (cData[OBJ_DEC].colnum == 0) { 1175 // Header keyword. 1176 readData(OBJ_DEC, 0, srcDec); 1177 for (int irow = 1; irow < nRow; irow++) { 1178 srcDec[irow] = *srcDec; 1179 } 1180 } else { 1181 // Table column. 1182 fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum, 1, 1, nRow, 1183 0, srcDec, &anynul, &cStatus); 1184 } 1185 1186 fits_read_col(cSDptr, TFLOAT, cData[PARANGLE].colnum, 1, 1, nRow, 0, 1187 par, &anynul, &cStatus); 1188 fits_read_col(cSDptr, TFLOAT, cData[FOCUSROT].colnum, 1, 1, nRow, 0, 1189 rot, &anynul, &cStatus); 1199 readCol(OBJ_RA, srcRA); 1200 readCol(OBJ_DEC, srcDec); 1201 readCol(PARANGLE, par); 1202 readCol(FOCUSROT, rot); 1190 1203 1191 1204 if (!cStatus) { 1192 for (int irow = 0; irow < nRow; irow++) {1205 for (int irow = 0; irow < cNRow; irow++) { 1193 1206 if (sel[irow]) { 1194 1207 // Convert to feed-plane coordinates. … … 1247 1260 // Find the next selected beam and IF. 1248 1261 short iBeam = 0, iIF = 0; 1249 while (++cRow <= cNRow) { 1262 while (1) { 1263 if (++cTimeIdx > cNAxisTime) { 1264 if (++cRow > cNRow) break; 1265 cTimeIdx = 1; 1266 } 1267 1250 1268 if (cData[BEAM].colnum > 0) { 1251 1269 readData(BEAM, cRow, &iBeam); … … 1269 1287 char chars[32]; 1270 1288 readData(OBSMODE, cRow, chars); 1271 if (strcmp(chars, "CAL") == 0) { 1289 if (strcmp(chars, "DROP") == 0) { 1290 // Completely flagged integration. 1291 continue; 1292 1293 } else if (strcmp(chars, "CAL") == 0) { 1294 sReset = 1; 1272 1295 if (cALFA_CIMA > 1) { 1273 1296 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) { … … 1280 1303 continue; 1281 1304 } 1305 1306 } else { 1307 // Reset for the next CAL record. 1308 if (sReset) { 1309 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) { 1310 sALFAcalNon[iBeam][iPol] = 0; 1311 sALFAcalNoff[iBeam][iPol] = 0; 1312 sALFAcalOn[iBeam][iPol] = 0.0f; 1313 sALFAcalOff[iBeam][iPol] = 0.0f; 1314 } 1315 sReset = 0; 1316 1317 sprintf(cMsg, "ALFA cal factors for beam %d: %.3e, %.3e", 1318 iBeam+1, sALFAcal[iBeam][0], sALFAcal[iBeam][1]); 1319 logMsg(cMsg); 1320 } 1282 1321 } 1283 1322 } … … 1312 1351 // Times. 1313 1352 char datobs[32]; 1314 readData(DATE_OBS, cRow, datobs); 1315 readData(TIME, cRow, &mbrec.utc); 1316 if (cALFA_BD) mbrec.utc *= 3600.0; 1317 1318 if (datobs[2] == '/') { 1319 // Translate an old-format DATE-OBS. 1320 datobs[9] = datobs[1]; 1321 datobs[8] = datobs[0]; 1322 datobs[2] = datobs[6]; 1323 datobs[5] = datobs[3]; 1324 datobs[3] = datobs[7]; 1325 datobs[6] = datobs[4]; 1326 datobs[7] = '-'; 1327 datobs[4] = '-'; 1328 datobs[1] = '9'; 1329 datobs[0] = '1'; 1330 1331 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) { 1332 // Dig UTC out of a new-format DATE-OBS. 1333 int hh, mm; 1334 float ss; 1335 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss); 1336 mbrec.utc = (hh*60 + mm)*60 + ss; 1337 } 1338 1339 datobs[10] = '\0'; 1353 readTime(cRow, cTimeIdx, datobs, mbrec.utc); 1340 1354 strcpy(mbrec.datobs, datobs); 1341 1355 1342 1356 if (cData[CYCLE].colnum > 0) { 1343 1357 readData(CYCLE, cRow, &mbrec.cycleNo); 1358 mbrec.cycleNo += cTimeIdx - 1; 1344 1359 if (cALFA_BD) mbrec.cycleNo++; 1345 1360 } else { … … 1471 1486 1472 1487 // Read data, sectioning and transposing it in the process. 1473 long *blc = new long[cNAx is+1];1474 long *trc = new long[cNAx is+1];1475 long *inc = new long[cNAx is+1];1476 for (int iaxis = 0; iaxis <= cNAx is; iaxis++) {1488 long *blc = new long[cNAxes+1]; 1489 long *trc = new long[cNAxes+1]; 1490 long *inc = new long[cNAxes+1]; 1491 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) { 1477 1492 blc[iaxis] = 1; 1478 1493 trc[iaxis] = 1; … … 1480 1495 } 1481 1496 1482 blc[cReqax[0]] = std::min(startChan, endChan); 1483 trc[cReqax[0]] = std::max(startChan, endChan); 1484 blc[cNAxis] = cRow; 1485 trc[cNAxis] = cRow; 1497 blc[cFreqAxis] = std::min(startChan, endChan); 1498 trc[cFreqAxis] = std::max(startChan, endChan); 1499 if (cTimeAxis >= 0) { 1500 blc[cTimeAxis] = cTimeIdx; 1501 trc[cTimeAxis] = cTimeIdx; 1502 } 1503 blc[cNAxes] = cRow; 1504 trc[cNAxes] = cRow; 1486 1505 1487 1506 mbrec.haveSpectra = cGetSpectra; … … 1489 1508 int anynul; 1490 1509 1491 for (int i pol = 0; ipol < nPol; ipol++) {1492 blc[c Reqax[1]] = ipol+1;1493 trc[c Reqax[1]] = ipol+1;1510 for (int iPol = 0; iPol < nPol; iPol++) { 1511 blc[cStokesAxis] = iPol+1; 1512 trc[cStokesAxis] = iPol+1; 1494 1513 1495 1514 if (cALFA && cALFA_CIMA < 2) { 1496 1515 // ALFA data: polarizations are stored in successive rows. 1497 blc[c Reqax[1]] = 1;1498 trc[c Reqax[1]] = 1;1499 1500 if (i pol) {1516 blc[cStokesAxis] = 1; 1517 trc[cStokesAxis] = 1; 1518 1519 if (iPol) { 1501 1520 if (++cRow > cNRow) { 1502 1521 return -1; 1503 1522 } 1504 1523 1505 blc[cNAx is] = cRow;1506 trc[cNAx is] = cRow;1524 blc[cNAxes] = cRow; 1525 trc[cNAxes] = cRow; 1507 1526 } 1508 1527 1509 1528 } else if (cData[DATA].nelem < 0) { 1510 1529 // Variable dimension array; get axis lengths. 1511 int naxis = 5, status;1512 1513 if ((status = readDim(DATA, cRow, &nax is, cNAxes))) {1530 int naxes = 5, status; 1531 1532 if ((status = readDim(DATA, cRow, &naxes, cNAxis))) { 1514 1533 logMsg(); 1515 1534 1516 } else if ((status = (nax is != cNAxis))) {1535 } else if ((status = (naxes != cNAxes))) { 1517 1536 logMsg("ERROR: DATA array dimensions changed."); 1518 1537 } … … 1526 1545 } 1527 1546 1528 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAx is, cNAxes,1529 blc, trc, inc, 0, mbrec.spectra[0] + i pol*nChan, &anynul,1547 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis, 1548 blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul, 1530 1549 &cStatus)) { 1531 1550 logMsg(); … … 1538 1557 if (endChan < startChan) { 1539 1558 // Reverse the spectrum. 1540 float *iptr = mbrec.spectra[0] + i pol*nChan;1559 float *iptr = mbrec.spectra[0] + iPol*nChan; 1541 1560 float *jptr = iptr + nChan - 1; 1542 1561 float *mid = iptr + nChan/2; … … 1550 1569 if (cALFA) { 1551 1570 // ALFA data, rescale the spectrum. 1552 float *chan = mbrec.spectra[0] + ipol*nChan; 1571 float el, zd; 1572 readData(ELEVATIO, cRow, &el); 1573 zd = 90.0f - el; 1574 1575 float factor = sALFAcal[iBeam][iPol] / alfaGain(zd); 1576 1577 if (cALFA_CIMA > 1) { 1578 // Rescale according to the number of unblanked accumulations. 1579 int colnum, naccum; 1580 findCol("STAT", &colnum); 1581 fits_read_col(cSDptr, TINT, colnum, cRow, 10*(cTimeIdx-1)+2, 1, 0, 1582 &naccum, &anynul, &cStatus); 1583 factor *= cALFAacc / naccum; 1584 } 1585 1586 float *chan = mbrec.spectra[0] + iPol*nChan; 1553 1587 float *chanN = chan + nChan; 1554 1588 while (chan < chanN) { 1555 1589 // Approximate conversion to Jy. 1556 *(chan++) *= cALFAcal[iBeam][iIF];1590 *(chan++) *= factor; 1557 1591 } 1558 1592 } 1559 1593 1560 if (mbrec.tsys[0][i pol] == 0.0) {1594 if (mbrec.tsys[0][iPol] == 0.0) { 1561 1595 // Compute Tsys as the average across the spectrum. 1562 float *chan = mbrec.spectra[0] + i pol*nChan;1596 float *chan = mbrec.spectra[0] + iPol*nChan; 1563 1597 float *chanN = chan + nChan; 1564 float *tsys = mbrec.tsys[0] + i pol;1598 float *tsys = mbrec.tsys[0] + iPol; 1565 1599 while (chan < chanN) { 1566 1600 *tsys += *(chan++); … … 1572 1606 // Read data flags. 1573 1607 if (cData[FLAGGED].colnum > 0) { 1574 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAx is,1575 cNAx es, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,1608 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxes, 1609 cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul, 1576 1610 &cStatus)) { 1577 1611 logMsg(); … … 1584 1618 if (endChan < startChan) { 1585 1619 // Reverse the flag vector. 1586 unsigned char *iptr = mbrec.flagged[0] + i pol*nChan;1620 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan; 1587 1621 unsigned char *jptr = iptr + nChan - 1; 1588 1622 for (int ichan = 0; ichan < nChan/2; ichan++) { … … 1595 1629 } else { 1596 1630 // All channels are unflagged by default. 1597 unsigned char *iptr = mbrec.flagged[0] + i pol*nChan;1631 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan; 1598 1632 for (int ichan = 0; ichan < nChan; ichan++) { 1599 1633 *(iptr++) = 0; … … 1708 1742 int status = 0; 1709 1743 fits_close_file(cSDptr, &status); 1710 cSDptr = 0 ;1744 cSDptr = 0x0; 1711 1745 1712 1746 if (cBeams) delete [] cBeams; … … 1791 1825 } 1792 1826 1827 //------------------------------------------------------ SDFITSreader::findCol 1828 1829 // Locate a parameter in the SDFITS file. 1830 1831 void SDFITSreader::findCol( 1832 char *name, 1833 int *colnum) 1834 { 1835 *colnum = 0; 1836 int status = 0; 1837 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status); 1838 1839 if (status) { 1840 // Not a real column - maybe it's virtual. 1841 char card[81]; 1842 1843 status = 0; 1844 fits_read_card(cSDptr, name, card, &status); 1845 if (status) { 1846 // Not virtual either. 1847 *colnum = -1; 1848 } 1849 1850 // Clear error messages. 1851 fits_clear_errmsg(); 1852 } 1853 } 1854 1793 1855 //------------------------------------------------------ SDFITSreader::readDim 1794 1856 … … 1798 1860 int iData, 1799 1861 long iRow, 1800 int *nax is,1801 long nax es[])1862 int *naxes, 1863 long naxis[]) 1802 1864 { 1803 1865 int colnum = cData[iData].colnum; … … 1806 1868 } 1807 1869 1808 int maxdim = *nax is;1870 int maxdim = *naxes; 1809 1871 if (cData[iData].tdimcol < 0) { 1810 1872 // No TDIMnnn column for this array. 1811 1873 if (cData[iData].nelem < 0) { 1812 1874 // Variable length array; read the array descriptor. 1813 *nax is = 1;1875 *naxes = 1; 1814 1876 long dummy; 1815 if (fits_read_descript(cSDptr, colnum, iRow, nax es, &dummy, &cStatus)) {1877 if (fits_read_descript(cSDptr, colnum, iRow, naxis, &dummy, &cStatus)) { 1816 1878 return 1; 1817 1879 } … … 1819 1881 } else { 1820 1882 // Read the repeat count from TFORMnnn. 1821 if (fits_read_tdim(cSDptr, colnum, maxdim, nax is, naxes, &cStatus)) {1883 if (fits_read_tdim(cSDptr, colnum, maxdim, naxes, naxis, &cStatus)) { 1822 1884 return 1; 1823 1885 } … … 1837 1899 1838 1900 tp++; 1839 *nax is = 0;1901 *naxes = 0; 1840 1902 for (size_t j = 1; j < strlen(tdimval); j++) { 1841 1903 if (tdimval[j] == ',' || tdimval[j] == ')') { 1842 sscanf(tp, "%ld", nax es + (*naxis)++);1904 sscanf(tp, "%ld", naxis + (*naxes)++); 1843 1905 if (tdimval[j] == ')') break; 1844 1906 tp = tdimval + j + 1; … … 1875 1937 findCol(name, &colnum); 1876 1938 1877 if (colnum > 0 ) {1939 if (colnum > 0 && iRow > 0) { 1878 1940 // Read the first value from the specified row of the table. 1879 1941 int coltype; … … 1938 2000 void *value) 1939 2001 { 1940 char *name = cData[iData].name;1941 2002 int type = cData[iData].type; 1942 2003 int colnum = cData[iData].colnum; 1943 long nelem = cData[iData].nelem; 1944 1945 if (colnum > 0) { 2004 2005 if (colnum > 0 && iRow > 0) { 1946 2006 // Read the required number of values from the specified row of the table. 2007 long nelem = cData[iData].nelem; 1947 2008 int anynul; 1948 2009 if (type == TSTRING) { … … 1973 2034 } else if (colnum == 0) { 1974 2035 // Read keyword value. 2036 char *name = cData[iData].name; 1975 2037 fits_read_key(cSDptr, type, name, value, 0, &cStatus); 1976 2038 … … 1993 2055 } 1994 2056 1995 //------------------------------------------------------ SDFITSreader:: findCol1996 1997 // Locate a parameter inthe SDFITS file.1998 1999 void SDFITSreader::findCol(2000 char *name,2001 int *colnum)2057 //------------------------------------------------------ SDFITSreader::readCol 2058 2059 // Read a scalar column from the SDFITS file. 2060 2061 int SDFITSreader::readCol( 2062 int iData, 2063 void *value) 2002 2064 { 2003 *colnum = 0; 2004 int status = 0; 2005 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status); 2006 2007 if (status) { 2008 // Not a real column - maybe it's virtual. 2009 char card[81]; 2010 2011 status = 0; 2012 fits_read_card(cSDptr, name, card, &status); 2013 if (status) { 2014 // Not virtual either. 2015 *colnum = -1; 2016 } 2017 2018 // Clear error messages. 2019 fits_clear_errmsg(); 2020 } 2065 int type = cData[iData].type; 2066 2067 if (cData[iData].colnum > 0) { 2068 // Table column. 2069 int anynul; 2070 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRow, 0, 2071 value, &anynul, &cStatus); 2072 2073 } else { 2074 // Header keyword. 2075 readData(iData, 0, value); 2076 for (int irow = 1; irow < cNRow; irow++) { 2077 if (type == TSHORT) { 2078 ((short *)value)[irow] = *((short *)value); 2079 } else if (type == TINT) { 2080 ((int *)value)[irow] = *((int *)value); 2081 } else if (type == TFLOAT) { 2082 ((float *)value)[irow] = *((float *)value); 2083 } else if (type == TDOUBLE) { 2084 ((double *)value)[irow] = *((double *)value); 2085 } 2086 } 2087 } 2088 2089 return cData[iData].colnum < 0; 2090 } 2091 2092 //----------------------------------------------------- SDFITSreader::readTime 2093 2094 // Read the time from the SDFITS file. 2095 2096 int SDFITSreader::readTime( 2097 long iRow, 2098 int iPix, 2099 char *datobs, 2100 double &utc) 2101 { 2102 readData(DATE_OBS, iRow, datobs); 2103 if (cData[TIME].colnum >= 0) { 2104 readData(TIME, iRow, &utc); 2105 } else if (cNAxisTime > 1) { 2106 double timeDelt, timeRefPix, timeRefVal; 2107 readData(TimeRefVal, iRow, &timeRefVal); 2108 readData(TimeDelt, iRow, &timeDelt); 2109 readData(TimeRefPix, iRow, &timeRefPix); 2110 utc = timeRefVal + (iPix - timeRefPix) * timeDelt; 2111 } 2112 2113 if (cALFA_BD) utc *= 3600.0; 2114 2115 // Check DATE-OBS format. 2116 if (datobs[2] == '/') { 2117 // Translate an old-format DATE-OBS. 2118 datobs[9] = datobs[1]; 2119 datobs[8] = datobs[0]; 2120 datobs[2] = datobs[6]; 2121 datobs[5] = datobs[3]; 2122 datobs[3] = datobs[7]; 2123 datobs[6] = datobs[4]; 2124 datobs[7] = '-'; 2125 datobs[4] = '-'; 2126 datobs[1] = '9'; 2127 datobs[0] = '1'; 2128 2129 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) { 2130 // Dig UTC out of a new-format DATE-OBS. 2131 int hh, mm; 2132 float ss; 2133 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss); 2134 utc = (hh*60 + mm)*60 + ss; 2135 } 2136 2137 datobs[10] = '\0'; 2138 2139 return 0; 2021 2140 } 2022 2141 … … 2047 2166 2048 2167 // Read cal data. 2049 long *blc = new long[cNAx is+1];2050 long *trc = new long[cNAx is+1];2051 long *inc = new long[cNAx is+1];2052 for (int iaxis = 0; iaxis <= cNAx is; iaxis++) {2168 long *blc = new long[cNAxes+1]; 2169 long *trc = new long[cNAxes+1]; 2170 long *inc = new long[cNAxes+1]; 2171 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) { 2053 2172 blc[iaxis] = 1; 2054 2173 trc[iaxis] = 1; … … 2060 2179 int endChan = cEndChan[iIF]; 2061 2180 2062 blc[cNAxis] = cRow; 2063 trc[cNAxis] = cRow; 2064 blc[cReqax[0]] = std::min(startChan, endChan); 2065 trc[cReqax[0]] = std::max(startChan, endChan); 2181 blc[cFreqAxis] = std::min(startChan, endChan); 2182 trc[cFreqAxis] = std::max(startChan, endChan); 2066 2183 if (cALFA_CIMA > 1) { 2067 2184 // CIMAFITS 2.x has a legitimate STOKES axis... 2068 blc[c Reqax[1]] = iPol+1;2069 trc[c Reqax[1]] = iPol+1;2185 blc[cStokesAxis] = iPol+1; 2186 trc[cStokesAxis] = iPol+1; 2070 2187 } else { 2071 2188 // ...older ALFA data does not. 2072 blc[cReqax[1]] = 1; 2073 trc[cReqax[1]] = 1; 2074 } 2189 blc[cStokesAxis] = 1; 2190 trc[cStokesAxis] = 1; 2191 } 2192 if (cTimeAxis >= 0) { 2193 blc[cTimeAxis] = cTimeIdx; 2194 trc[cTimeAxis] = cTimeIdx; 2195 } 2196 blc[cNAxes] = cRow; 2197 trc[cNAxes] = cRow; 2075 2198 2076 2199 float spectrum[endChan]; 2077 2200 int anynul; 2078 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAx is, cNAxes,2201 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis, 2079 2202 blc, trc, inc, 0, spectrum, &anynul, &cStatus)) { 2080 2203 logMsg(); … … 2085 2208 } 2086 2209 2210 // Factor to rescale according to the number of unblanked accumulations. 2211 float factor = 1.0f; 2212 if (cALFA_CIMA > 1) { 2213 int colnum, naccum; 2214 findCol("STAT", &colnum); 2215 fits_read_col(cSDptr, TINT, colnum, cRow, 2, 1, 0, &naccum, &anynul, 2216 &cStatus); 2217 factor = cALFAacc / naccum; 2218 } 2219 2087 2220 // Average the spectrum. 2088 2221 float mean = 1e9f; … … 2098 2231 if (*chan < discrim) { 2099 2232 nChan++; 2100 sum += *chan ;2233 sum += *chan * factor; 2101 2234 } 2102 2235 } … … 2106 2239 2107 2240 if (calOn) { 2108 cALFAcalOn[iBeam][iPol] += mean; 2241 sALFAcalOn[iBeam][iPol] *= sALFAcalNon[iBeam][iPol]; 2242 sALFAcalOn[iBeam][iPol] += mean; 2243 sALFAcalOn[iBeam][iPol] /= ++sALFAcalNon[iBeam][iPol]; 2109 2244 } else { 2110 cALFAcalOff[iBeam][iPol] += mean; 2111 } 2112 2113 if (cALFAcalOn[iBeam][iPol] != 0.0f && 2114 cALFAcalOff[iBeam][iPol] != 0.0f) { 2245 sALFAcalOff[iBeam][iPol] *= sALFAcalNoff[iBeam][iPol]; 2246 sALFAcalOff[iBeam][iPol] += mean; 2247 sALFAcalOff[iBeam][iPol] /= ++sALFAcalNoff[iBeam][iPol]; 2248 } 2249 2250 if (sALFAcalNon[iBeam][iPol] && sALFAcalNoff[iBeam][iPol]) { 2115 2251 // Tcal should come from the TCAL table, it varies weakly with beam, 2116 2252 // polarization, and frequency. However, TCAL is not written properly. 2117 2253 float Tcal = 12.0f; 2118 cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -2119 cALFAcalOff[iBeam][iPol]);2254 sALFAcal[iBeam][iPol] = Tcal / (sALFAcalOn[iBeam][iPol] - 2255 sALFAcalOff[iBeam][iPol]); 2120 2256 2121 2257 // Scale from K to Jy; the gain also varies weakly with beam, 2122 2258 // polarization, frequency, and zenith angle. 2123 2259 float fluxCal = 10.0f; 2124 cALFAcal[iBeam][iPol] /= fluxCal; 2125 2126 cALFAcalOn[iBeam][iPol] = 0.0f; 2127 cALFAcalOff[iBeam][iPol] = 0.0f; 2260 sALFAcal[iBeam][iPol] /= fluxCal; 2128 2261 } 2129 2262 2130 2263 return 0; 2131 2264 } 2265 2266 //----------------------------------------------------- SDFITSreader::alfaGain 2267 2268 // ALFA gain factor. 2269 2270 float SDFITSreader::alfaGain( 2271 float zd) 2272 { 2273 // Gain vs zenith distance table from Robert Minchin, 2008/12/08. 2274 const int nZD = 37; 2275 const float zdLim[] = {1.5f, 19.5f}; 2276 const float zdInc = (nZD - 1) / (zdLim[1] - zdLim[0]); 2277 float zdGain[] = { 1.00723708, 2278 1.16644573, 1.15003645, 1.07117307, 1.02532673, 2279 1.01788402, 1.01369524, 1.00000000, 0.989855111, 2280 0.990888834, 0.993996620, 0.989964068, 0.982213855, 2281 0.978662670, 0.979349494, 0.978478372, 0.974631131, 2282 0.972126007, 0.972835243, 0.972742677, 0.968671739, 2283 0.963891327, 0.963452935, 0.966831207, 0.969585896, 2284 0.970700860, 0.972644389, 0.973754644, 0.967344403, 2285 0.952168941, 0.937160134, 0.927843094, 0.914048433, 2286 0.886700928, 0.864701211, 0.869126320, 0.854309499}; 2287 2288 float gain; 2289 // Do table lookup by linear interpolation. 2290 float lambda = zdInc * (zd - zdLim[0]); 2291 int j = int(lambda); 2292 if (j < 0) { 2293 gain = zdGain[0]; 2294 } else if (j >= nZD-1) { 2295 gain = zdGain[nZD-1]; 2296 } else { 2297 gain = zdGain[j] + (lambda - j) * (zdGain[j+1] - zdGain[j]); 2298 } 2299 2300 return gain; 2301 } 2302 -
trunk/external/atnf/PKSIO/SDFITSreader.h
r1452 r1635 2 2 //# SDFITSreader.h: ATNF CFITSIO interface class for SDFITS input. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Associated Universities, Inc. Washington DC, USA. 6 6 //# … … 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: SDFITSreader.h,v 19. 16 2008-11-17 06:47:05cal103 Exp $28 //# $Id: SDFITSreader.h,v 19.21 2009-03-18 07:11:51 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 //# The SDFITSreader class reads single dish FITS files such as those written … … 108 108 109 109 private: 110 int cCycleNo, cExtraSysCal, cNAxis, cStatus; 111 long cNAxes[5], cNRow, cReqax[4], cRow; 110 int cCycleNo, cExtraSysCal, cNAxes, cStatus; 111 long cBeamAxis, cDecAxis, cFreqAxis, cNAxis[5], cNAxisTime, cNRow, 112 cRaAxis, cRow, cStokesAxis, cTimeAxis, cTimeIdx; 112 113 double cLastUTC; 113 114 fitsfile *cSDptr; … … 118 119 119 120 enum {SCAN, CYCLE, DATE_OBS, TIME, EXPOSURE, OBJECT, OBJ_RA, OBJ_DEC, 120 RESTFRQ, OBSMODE, BEAM, IF, FqRef Pix, FqRefVal, FqDelt, RA, DEC,121 SCANRATE, TSYS, CALFCTR, XCALFCTR, BASELIN, BASESUB, DATA, FLAGGED,122 DATAXED, XPOLDATA, REFBEAM, TCAL, TCALTIME, AZIMUTH, ELEVATIO,123 PARANGLE, FOCUSAXI, FOCUSTAN, FOCUSROT, TAMBIENT, PRESSURE,124 HUMIDITY, WINDSPEE, WINDDIRE, NDATA};121 RESTFRQ, OBSMODE, BEAM, IF, FqRefVal, FqDelt, FqRefPix, RA, DEC, 122 TimeRefVal, TimeDelt, TimeRefPix, SCANRATE, TSYS, CALFCTR, XCALFCTR, 123 BASELIN, BASESUB, DATA, FLAGGED, DATAXED, XPOLDATA, REFBEAM, TCAL, 124 TCALTIME, AZIMUTH, ELEVATIO, PARANGLE, FOCUSAXI, FOCUSTAN, FOCUSROT, 125 TAMBIENT, PRESSURE, HUMIDITY, WINDSPEE, WINDDIRE, NDATA}; 125 126 126 127 // Message handling. … … 128 129 129 130 void findData(int iData, char *name, int type); 131 void findCol(char *name, int *colnum); 130 132 int readDim(int iData, long iRow, int *naxis, long naxes[]); 131 133 int readParm(char *name, int type, void *value); 132 134 int readData(char *name, int type, long iRow, void *value); 133 135 int readData(int iData, long iRow, void *value); 134 void findCol(char *name, int *colnum); 136 int readCol(int iData, void *value); 137 int readTime(long iRow, int iPix, char *datobs, double &utc); 135 138 136 // These are for ALFA data: "BDFITS" or "CIMAFITS". 139 // These are for ALFA data: "BDFITS" or "CIMAFITS". Statics are required 140 // for CIMAFITS v2.0 because CAL ON/OFF data is split into separate files. 141 static int sInit, sReset; 142 static int (*sALFAcalNon)[2], (*sALFAcalNoff)[2]; 143 static float (*sALFAcal)[2], (*sALFAcalOn)[2], (*sALFAcalOff)[2]; 144 137 145 int cALFA, cALFA_BD, cALFA_CIMA, cALFAscan, cScanNo; 138 float cALFA cal[8][2], cALFAcalOn[8][2], cALFAcalOff[8][2];146 float cALFAacc; 139 147 int alfaCal(short iBeam, short iIF, short iPol); 148 float alfaGain(float zd); 140 149 141 150 // These are for GBT data. -
trunk/external/atnf/PKSIO/SDFITSwriter.cc
r1466 r1635 2 2 //# SDFITSwriter.cc: ATNF CFITSIO interface class for SDFITS output. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 84 //# Copyright (C) 2000-2009 5 5 //# Mark Calabretta, ATNF 6 6 //# … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id: SDFITSwriter.cc,v 19.1 5 2008-11-17 06:58:58cal103 Exp $29 //# $Id: SDFITSwriter.cc,v 19.17 2009-09-18 07:29:26 cal103 Exp $ 30 30 //#--------------------------------------------------------------------------- 31 31 //# Original: 2000/07/24, Mark Calabretta, ATNF … … 37 37 38 38 #include <casa/iostream.h> 39 #include <cstring>40 39 41 40 #include <algorithm> … … 174 173 char version[7]; 175 174 char date[11]; 176 sscanf("$Revision: 19.1 5$", "%*s%s", version);177 sscanf("$Date: 200 8-11-17 06:58:58$", "%*s%s", date);175 sscanf("$Revision: 19.17 $", "%*s%s", version); 176 sscanf("$Date: 2009-09-18 07:29:26 $", "%*s%s", date); 178 177 sprintf(text, "SDFITSwriter (v%s, %s)", version, date); 179 178 fits_write_key_str(cSDptr, "ORIGIN", text, "output class", &cStatus); … … 231 230 232 231 // CYCLE (additional, real). 233 fits_insert_col(cSDptr, ++ncol, "CYCLE", "1 I", &cStatus);232 fits_insert_col(cSDptr, ++ncol, "CYCLE", "1J", &cStatus); 234 233 235 234 // DATE-OBS (core, real). … … 391 390 392 391 // BASESUB (additional, real). 393 sprintf(tform, "%dE", 9*maxNPol);392 sprintf(tform, "%dE", 24*maxNPol); 394 393 fits_insert_col(cSDptr, ++ncol, "BASESUB", tform, &cStatus); 395 tdim[0] = 9;394 tdim[0] = 24; 396 395 fits_write_tdim(cSDptr, ncol, 2, tdim, &cStatus); 397 396 } … … 681 680 682 681 // BASESUB. 683 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 9*nPol, mbrec.baseSub[0][0],682 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 24*nPol, mbrec.baseSub[0][0], 684 683 &cStatus); 685 684 }
Note:
See TracChangeset
for help on using the changeset viewer.