Changeset 1453 for branches/alma/external
- Timestamp:
- 11/27/08 12:47:15 (16 years ago)
- Location:
- branches/alma/external/atnf/PKSIO
- Files:
-
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/alma/external/atnf/PKSIO/FITSreader.cc
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : FITSreader.cc,v 19.2 2006/05/19 02:14:50 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The FITSreader class is an abstract base class for the Parkes Multibeam -
branches/alma/external/atnf/PKSIO/FITSreader.h
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : FITSreader.h,v 19.5 2006/05/19 02:18:02 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The FITSreader class is an abstract base class for the Parkes Multibeam -
branches/alma/external/atnf/PKSIO/MBFITSreader.cc
r1372 r1453 2 2 //# MBFITSreader.cc: ATNF single-dish RPFITS reader. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 74 //# Copyright (C) 2000-2006 5 5 //# Mark Calabretta, ATNF 6 6 //# … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : MBFITSreader.cc,v 19.34 2007/07/02 06:12:18 cal103 Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes … … 81 81 cRefChan = 0x0; 82 82 83 cVis = 0x0; 84 cWgt = 0x0; 83 cVis = new float[2*4*8163]; 85 84 86 85 cBeamSel = 0x0; … … 92 91 93 92 cMBopen = 0; 93 jstat = -3; 94 94 } 95 95 … … 127 127 128 128 // Open the RPFITS file. 129 int jstat = -3; 130 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, 131 &cBin, &cIFno, &cSrcNo); 129 rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin, 130 &if_no, &sourceno); 132 131 133 132 if (jstat) { 134 fprintf(stderr, " ERROR, failed to open MBFITS file: %s\n", rpname);133 fprintf(stderr, "Failed to open MBFITS file: %s\n", rpname); 135 134 return 1; 136 135 } … … 148 147 // Read the first header. 149 148 jstat = -1; 150 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,151 & cBin, &cIFno, &cSrcNo);149 rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin, 150 &if_no, &sourceno); 152 151 153 152 if (jstat) { 154 fprintf(stderr, " ERROR, failed to read MBFITS header: %s\n", rpname);153 fprintf(stderr, "Failed to read MBFITS header: %s\n", rpname); 155 154 close(); 156 155 return 1; … … 160 159 cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0; 161 160 162 // Non-ATNF data may not store the position in (u,v,w). 163 if (strncmp(names_.sta, "tid", 3) == 0) { 164 fprintf(stderr, "WARNING, found Tidbinbilla data"); 165 cSUpos = 1; 166 } else if (strncmp(names_.sta, "HOB", 3) == 0) { 167 fprintf(stderr, "WARNING, found Hobart data"); 168 cSUpos = 1; 169 } else if (strncmp(names_.sta, "CED", 3) == 0) { 170 fprintf(stderr, "WARNING, found Ceduna data"); 171 cSUpos = 1; 172 } else { 173 cSUpos = 0; 174 } 175 176 if (cSUpos) { 177 fprintf(stderr, ", using telescope position from SU table.\n"); 161 // Tidbinbilla data has some more. 162 cTid = strncmp(names_.sta, "tid", 3) == 0; 163 if (cTid) { 164 // Telescope position is stored in the source table. 178 165 cInterp = 0; 179 166 } … … 189 176 190 177 if (cNBeam <= 0) { 191 fprintf(stderr, " ERROR, couldn't determine number of beams.\n");178 fprintf(stderr, "Couldn't determine number of beams.\n"); 192 179 close(); 193 180 return 1; … … 250 237 } 251 238 252 // Allocate memory for RPFITSIN subroutine arguments.253 if ( cVis) delete [] cVis;254 if (cWgt) delete [] cWgt;255 cVis = new float[2*maxProd];256 cWgt = new float[maxProd];239 // Is the vis array declared by RPFITS.h large enough? 240 if (8*8193 < maxProd) { 241 // Need to allocate more memory for RPFITSIN. 242 cVis = new float[2*maxProd]; 243 } 257 244 258 245 nChan = cNChan; … … 291 278 // Read the first syscal record. 292 279 if (rpget(1, cEOS)) { 293 fprintf(stderr, "E RROR, failed to read first syscal record.\n");280 fprintf(stderr, "Error: Failed to read first syscal record.\n"); 294 281 close(); 295 282 return 1; … … 326 313 { 327 314 if (!cMBopen) { 328 fprintf(stderr, " ERROR, an MBFITS file has not been opened.\n");315 fprintf(stderr, "An MBFITS file has not been opened.\n"); 329 316 return 1; 330 317 } … … 399 386 // Time at start of observation. 400 387 sprintf(datobs, "%-10.10s", names_.datobs); 401 utc = cUTC;388 utc = ut; 402 389 403 390 // Spectral parameters. … … 448 435 449 436 if (!cMBopen) { 450 fprintf(stderr, " ERROR, an MBFITS file has not been opened.\n");437 fprintf(stderr, "An MBFITS file has not been opened.\n"); 451 438 return 1; 452 439 } … … 576 563 577 564 if (cNBin > 1 && cNBeamSel > 1) { 578 fprintf(stderr, " ERROR, cannot handle binning mode for multiple "565 fprintf(stderr, "Cannot handle binning mode for multiple " 579 566 "beams.\n"); 580 567 close(); … … 582 569 } 583 570 584 // Allocate buffer data storage; the PKSMBrecord constructor zeroes 585 // class members such as cycleNo that are tested in the first pass 586 // below. 571 // Allocate buffer data storage. 587 572 int nBuff = cNBeamSel * cNBin; 588 573 cBuffer = new PKSMBrecord[nBuff]; … … 599 584 cScanNo = 1; 600 585 cCycleNo = 0; 601 c PrevUTC = 0.0;586 cUTC = 0.0; 602 587 cStaleness = new int[cNBeamSel]; 603 588 for (int iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) { … … 610 595 cScanNo++; 611 596 cCycleNo = 0; 612 c PrevUTC = 0.0;597 cUTC = 0.0; 613 598 } 614 599 615 600 // Apply beam selection. 616 beamNo = int( cBaseline / 256.0);601 beamNo = int(baseline / 256.0); 617 602 iBeamSel = cBeamSel[beamNo-1]; 618 603 if (iBeamSel < 0) continue; 619 604 620 605 // Sanity check (mainly for MOPS). 621 if ( cIFno > cNIF) continue;606 if (if_no > cNIF) continue; 622 607 623 608 // Apply IF selection. 624 iIFSel = cIFSel[ cIFno - 1];609 iIFSel = cIFSel[if_no - 1]; 625 610 if (iIFSel < 0) continue; 626 611 627 612 sprintf(cDateObs, "%-10.10s", names_.datobs); 628 613 629 // Check for change-of-day. 630 if (cUTC < cPrevUTC - 85800.0) { 631 cUTC += 86400.0; 614 // Change-of-day; note that the ut variable from RPFITS.h is global 615 // and will be preserved between calls to this function. 616 if (ut < cUTC - 85800.0) { 617 ut += 86400.0; 618 } 619 620 // New integration cycle? 621 if (ut > cUTC) { 622 cCycleNo++; 623 cUTC = ut + 0.0001; 632 624 } 633 625 634 626 if (cNBin > 1) { 635 627 // Binning mode: correct the time. 636 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0); 637 } 638 639 // New integration cycle? 640 if (cUTC > cPrevUTC) { 641 cCycleNo++; 642 cPrevUTC = cUTC + 0.0001; 628 ut += param_.intbase * (bin - (cNBin + 1)/2.0); 643 629 } 644 630 645 631 // Compute buffer number. 646 632 iMBuff = cBuffer + iBeamSel; 647 if (cNBin > 1) iMBuff += cNBeamSel*( cBin-1);633 if (cNBin > 1) iMBuff += cNBeamSel*(bin-1); 648 634 649 635 if (cCycleNo < iMBuff->cycleNo) { … … 654 640 655 641 // Begin flush cycle? 656 if (cEOS || (iMBuff->nIF && cUTC> iMBuff->utc + 0.0001)) {642 if (cEOS || (iMBuff->nIF && ut > iMBuff->utc + 0.0001)) { 657 643 cFlushing = 1; 658 644 cFlushBin = 0; … … 661 647 662 648 #ifdef PKSIO_DEBUG 663 printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);649 printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, if_no); 664 650 if (cEOS) printf("Start of new scan, flushing previous scan.\n"); 665 651 #endif … … 710 696 711 697 // The last record read must have been the first of a new cycle. 712 beamNo = int( cBaseline / 256.0);698 beamNo = int(baseline / 256.0); 713 699 iBeamSel = cBeamSel[beamNo-1]; 714 700 715 701 // Compute buffer number. 716 702 iMBuff = cBuffer + iBeamSel; 717 if (cNBin > 1) iMBuff += cNBeamSel*( cBin-1);703 if (cNBin > 1) iMBuff += cNBeamSel*(bin-1); 718 704 } 719 705 } … … 735 721 // returned as the 'u' and 'v' visibility coordinates, must be 736 722 // interpolated to the integration time which RPFITSIN returns as 737 // ' cUTC', this usually being a second or two later.723 // 'ut', this usually being a second or two later. 738 724 // 739 725 // Note that the time recorded as the 'w' visibility coordinate 740 // cycles through 86400 back to 0 at midnight, whereas that in ' cUTC'726 // cycles through 86400 back to 0 at midnight, whereas that in 'ut' 741 727 // continues to increase past 86400. 742 728 743 double thisRA = cU;744 double thisDec = cV;745 double thisUTC = cW;729 double thisRA = u; 730 double thisDec = v; 731 double thisUTC = w; 746 732 747 733 if (thisUTC < prevUTC) { … … 865 851 // Signal that all IFs in this buffer location have been flushed. 866 852 iMBuff->nIF = 0; 867 868 // Stop cEOS being set when the next integration is read.869 iMBuff->cycleNo = 0;870 871 853 } else { 872 854 // Carry on flushing the other IFs. … … 877 859 if (cFlushBin == cNBin - 1) { 878 860 if (cEOS || cEOF) { 861 // Stop cEOS being set when the next integration is read. 862 iMBuff->cycleNo = 0; 863 879 864 // Carry on flushing other buffers. 880 865 cFlushIF = 0; … … 884 869 cFlushing = 0; 885 870 886 beamNo = int( cBaseline / 256.0);871 beamNo = int(baseline / 256.0); 887 872 iBeamSel = cBeamSel[beamNo-1]; 888 873 889 874 // Compute buffer number. 890 875 iMBuff = cBuffer + iBeamSel; 891 if (cNBin > 1) iMBuff += cNBeamSel*( cBin-1);876 if (cNBin > 1) iMBuff += cNBeamSel*(bin-1); 892 877 } 893 878 } … … 903 888 } 904 889 905 // Sanity check on incomplete integrations within a scan.906 if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {907 // Force the incomplete integration to be flushed before proceeding.908 cFlushing = 1;909 continue;910 }911 912 890 iMBuff->scanNo = cScanNo; 913 891 iMBuff->cycleNo = cCycleNo; … … 915 893 // Times. 916 894 strncpy(iMBuff->datobs, cDateObs, 10); 917 iMBuff->utc = cUTC;895 iMBuff->utc = ut; 918 896 iMBuff->exposure = param_.intbase; 919 897 920 898 // Source identification. 921 899 sprintf(iMBuff->srcName, "%-16.16s", 922 names_.su_name + ( cSrcNo-1)*16);923 iMBuff->srcRA = doubles_.su_ra[ cSrcNo-1];924 iMBuff->srcDec = doubles_.su_dec[ cSrcNo-1];900 names_.su_name + (sourceno-1)*16); 901 iMBuff->srcRA = doubles_.su_ra[sourceno-1]; 902 iMBuff->srcDec = doubles_.su_dec[sourceno-1]; 925 903 926 904 // Rest frequency of the line of interest. … … 948 926 949 927 // Beam position at the specified time. 950 if (c SUpos) {951 // Non-ATNF data that does not store the position in (u,v,w).952 iMBuff->ra = doubles_.su_ra[ cSrcNo-1];953 iMBuff->dec = doubles_.su_dec[ cSrcNo-1];928 if (cTid) { 929 // Tidbinbilla data. 930 iMBuff->ra = doubles_.su_ra[sourceno-1]; 931 iMBuff->dec = doubles_.su_dec[sourceno-1]; 954 932 } else { 955 iMBuff->ra = cU;956 iMBuff->dec = cV;957 } 958 cPosUTC[iBeamSel] = cW;933 iMBuff->ra = u; 934 iMBuff->dec = v; 935 } 936 cPosUTC[iBeamSel] = w; 959 937 960 938 // IF-dependent parameters. 961 int iIF = cIFno - 1;939 int iIF = if_no - 1; 962 940 int startChan = cStartChan[iIF]; 963 941 int endChan = cEndChan[iIF]; … … 968 946 iIFSel = cIFSel[iIF]; 969 947 iMBuff->nIF++; 970 iMBuff->IFno[iIFSel] = cIFno;948 iMBuff->IFno[iIFSel] = if_no; 971 949 iMBuff->nChan[iIFSel] = nChan; 972 950 iMBuff->nPol[iIFSel] = cNPol[iIF]; … … 1027 1005 1028 1006 // The baseline flag may be set independently. 1029 if (rpflag == 0) rpflag = cFlag;1007 if (rpflag == 0) rpflag = flag; 1030 1008 1031 1009 // Copy and scale data. … … 1157 1135 int numErr = 0; 1158 1136 1159 intjstat = 0;1137 jstat = 0; 1160 1138 while (numErr < 10) { 1161 1139 int lastjstat = jstat; 1162 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,1163 & cBin, &cIFno, &cSrcNo);1140 rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin, 1141 &if_no, &sourceno); 1164 1142 1165 1143 switch(jstat) { … … 1174 1152 // Successful read. 1175 1153 if (lastjstat == 0) { 1176 if ( cBaseline == -1) {1154 if (baseline == -1) { 1177 1155 // Syscal data. 1178 1156 if (syscalonly) { … … 1230 1208 } 1231 1209 1232 fprintf(stderr, " ERROR,RPFITS read failed too many times.\n");1210 fprintf(stderr, "RPFITS read failed too many times.\n"); 1233 1211 return 2; 1234 1212 } … … 1241 1219 { 1242 1220 if (cMBopen) { 1243 intjstat = 1;1244 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,1245 & cBin, &cIFno, &cSrcNo);1221 jstat = 1; 1222 rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin, 1223 &if_no, &sourceno); 1246 1224 1247 1225 if (cBeams) delete [] cBeams; … … 1254 1232 if (cRefChan) delete [] cRefChan; 1255 1233 1256 if (cVis) delete [] cVis; 1257 if (cWgt) delete [] cWgt; 1234 if (cVis) delete [] cVis; 1258 1235 1259 1236 if (cBeamSel) delete [] cBeamSel; -
branches/alma/external/atnf/PKSIO/MBFITSreader.h
r1372 r1453 2 2 //# MBFITSreader.h: ATNF single-dish RPFITS reader. 3 3 //#--------------------------------------------------------------------------- 4 //# Copyright (C) 2000-200 74 //# Copyright (C) 2000-2006 5 5 //# Mark Calabretta, ATNF 6 6 //# … … 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : MBFITSreader.h,v 19.13 2007/07/02 06:11:22 cal103 Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes … … 107 107 108 108 private: 109 // RPFITSIN subroutine arguments.110 int cBaseline, cFlag, cBin, cIFno, cSrcNo;111 float cUTC, cU, cV, cW, *cVis, *cWgt;112 113 109 char cDateObs[10]; 114 110 int *cBeamSel, *cChanOff, cFirst, *cIFSel, cInterp, cIntTime, cMBopen, 115 cMopra, cNBeamSel, cNBin, cRetry, *cStaleness, cSUpos, *cXpolOff; 111 cMopra, cNBeamSel, cNBin, cRetry, *cStaleness, cTid, *cXpolOff; 112 float *cVis; 116 113 117 114 // The data has to be bufferred to allow positions to be interpolated. … … 121 118 122 119 int cCycleNo, cScanNo; 123 double c PrevUTC;120 double cUTC; 124 121 125 122 // Read the next data record from the RPFITS file. 126 123 int rpget(int syscalonly, int &EOS); 124 // These are no longer define in the RPFITS.h file so moved here 125 int jstat; 126 float *weight; 127 int baseline; 128 int flag; 129 int bin; 130 int if_no; 131 int sourceno; 132 float *vis; 133 float u; 134 float v; 135 float w; 136 float ut; 137 127 138 }; 128 139 -
branches/alma/external/atnf/PKSIO/PKSFITSreader.cc
r1393 r1453 337 337 Double &bandwidth, 338 338 Double &freqInc, 339 Double&restFreq,339 Vector<Double> &restFreq, 340 340 Vector<Float> &tcal, 341 341 String &tcalTime, … … 403 403 bandwidth = chanWidth * nChan; 404 404 freqInc = cMBrec.fqDelt[0]; 405 restFreq = cMBrec.restFreq; 405 //restFreq = cMBrec.restFreq; 406 restFreq.resize(1); 407 restFreq(0) = cMBrec.restFreq; 406 408 407 409 tcal.resize(nPol); -
branches/alma/external/atnf/PKSIO/PKSFITSreader.h
r1393 r1453 129 129 Double &bandwidth, 130 130 Double &freqInc, 131 Double&restFreq,131 Vector<Double> &restFreq, 132 132 Vector<Float> &tcal, 133 133 String &tcalTime, -
branches/alma/external/atnf/PKSIO/PKSMBrecord.cc
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : PKSMBrecord.cc,v 19.6 2006/07/05 04:52:07 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The PKSMBrecord class stores an MBFITS single-dish data record. -
branches/alma/external/atnf/PKSIO/PKSMBrecord.h
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : PKSMBrecord.h,v 19.8 2006/07/05 04:52:07 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# The PKSMBrecord class stores an MBFITS single-dish data record. -
branches/alma/external/atnf/PKSIO/PKSMS2reader.cc
r1393 r1453 92 92 cPKSMS = MeasurementSet(msName); 93 93 // taql access to the syscal table 94 String tmpcalName = msName; 95 tmpcalName.append("/SYSCAL"); 96 cSysCalTab = Table(tmpcalName); 94 cHaveSysCal = False; 95 if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) { 96 cSysCalTab = Table(cPKSMS.sysCalTableName()); 97 } 97 98 98 99 cIdx = 0; … … 116 117 ROMSSysCalColumns sysCalCols(cPKSMS.sysCal()); 117 118 ROMSWeatherColumns weatherCols(cPKSMS.weather()); 119 ROMSAntennaColumns antennaCols(cPKSMS.antenna()); 118 120 119 121 // Column accessors for required columns. … … 153 155 154 156 // Optional columns. 157 cHaveTsys = False; 158 cHaveTcal = False; 155 159 if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) { 156 160 cSrcVelCol.attach(cPKSMS.source(), "SYSVEL"); 157 161 } 158 162 159 if ( (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {163 if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) { 160 164 cTsysCol.attach(cPKSMS.sysCal(), "TSYS"); 161 165 } 162 166 163 if ( (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {167 if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) { 164 168 cTcalCol.attach(cPKSMS.sysCal(), "TCAL"); 165 169 } … … 176 180 // Spectral data should always be present. 177 181 haveSpectra = True; 178 cFloatDataCol.reference(msCols.floatData()); 182 cHaveDataCol = False; 183 cHaveCorrectedDataCol = False; 184 //String telName = antennaCols.name()(0); 185 ROMSObservationColumns observationCols(cPKSMS.observation()); 186 String telName = observationCols.telescopeName()(0); 187 //cATF = (telName.contains("DA41") || telName.contains("DV01")); 188 cATF = telName.contains("ATF"); 189 190 if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) { 191 if (cATF) { 192 //try to read a single baseline interferometeric data 193 //and treat it as single dish data 194 //maybe extended for ALMA commissioning later 195 cDataCol.reference(msCols.data()); 196 if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) { 197 //cerr<<"Do have CORRECTED_DATA column"<<endl; 198 cCorrectedDataCol.reference(msCols.correctedData()); 199 } 200 } 201 } 202 else { 203 cFloatDataCol.reference(msCols.floatData()); 204 } 179 205 cFlagCol.reference(msCols.flag()); 180 206 181 207 182 if ( (cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {208 if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cATF))) { 183 209 if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) { 184 210 cXCalFctrCol.attach(cPKSMS, "XCALFCTR"); … … 223 249 haveXPol = False; 224 250 225 if (cGetXPol ) {251 if (cGetXPol && !(cATF)) { 226 252 for (Int irow = 0; irow < cNRow; irow++) { 227 253 if (cDataCol.isDefined(irow)) { … … 303 329 304 330 fluxunit = ""; 305 const TableRecord& keywordSet 306 = cFloatDataCol.columnDesc().keywordSet(); 331 if (cHaveDataCol) { 332 const TableRecord& keywordSet2 333 = cDataCol.columnDesc().keywordSet(); 334 if(keywordSet2.isDefined("UNIT")) { 335 fluxunit = keywordSet2.asString("UNIT"); 336 } 337 } else { 338 const TableRecord& keywordSet 339 = cFloatDataCol.columnDesc().keywordSet(); 340 if(keywordSet.isDefined("UNIT")) { 341 fluxunit = keywordSet.asString("UNIT"); 342 } 343 } 344 345 /*** 346 const TableRecord& keywordSet 347 = cFloatDataCol.columnDesc().keywordSet(); 307 348 if(keywordSet.isDefined("UNIT")) { 308 349 fluxunit = keywordSet.asString("UNIT"); 309 350 } 310 351 ***/ 311 352 // Coordinate equinox. 312 353 ROMSPointingColumns pointingCols(cPKSMS.pointing()); … … 529 570 Double &bandwidth, 530 571 Double &freqInc, 531 Double&restFreq,572 Vector<Double> &restFreq, 532 573 Vector<Float> &tcal, 533 574 String &tcalTime, … … 584 625 } 585 626 } 586 587 627 // Renumerate scan no. Here still is 1-based 588 628 //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1; … … 624 664 625 665 // Systemic velocity. 626 if (!cHaveSrcVel ) {666 if (!cHaveSrcVel || cATF) { 627 667 srcVel = 0.0f; 628 668 } else { … … 634 674 Bool cGBT = telescope.contains("GBT"); 635 675 // Observation type. 636 Int stateId = cStateIdCol(cIdx); 637 if (stateId == -1) { 638 //obsMode = " "; 639 } else { 640 obsMode = cObsModeCol(stateId); 641 Bool sigState =cSigStateCol(stateId); 642 Bool refState =cRefStateCol(stateId); 643 //DEBUG 644 //cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl; 645 if (cGBT) { 646 // split the obsMode string and append a proper label 647 // (these are GBT specific) 648 int epos = obsMode.find_first_of(':'); 649 int nextpos = obsMode.find_first_of(':',epos+1); 650 string obsMode1 = obsMode.substr(0,epos); 651 string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1); 676 // check if State Table exist 677 //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName()); 678 Int stateId = 0; 679 Int StateNRow = 0; 680 StateNRow=cObsModeCol.nrow(); 681 if (Table::isReadable(cPKSMS.stateTableName())) { 682 obsMode = " "; 683 if (StateNRow > 0) { 684 stateId = cStateIdCol(cIdx); 685 if (stateId == -1) { 686 //obsMode = " "; 687 } else { 688 obsMode = cObsModeCol(stateId); 689 Bool sigState =cSigStateCol(stateId); 690 Bool refState =cRefStateCol(stateId); 691 //DEBUG 692 //cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl; 693 if (cGBT) { 694 // split the obsMode string and append a proper label 695 // (these are GBT specific) 696 int epos = obsMode.find_first_of(':'); 697 int nextpos = obsMode.find_first_of(':',epos+1); 698 string obsMode1 = obsMode.substr(0,epos); 699 string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1); 652 700 653 //cerr <<"obsMode2= "<<obsMode2<<endl; 654 if (!srcName.contains("_ps") 655 &&!srcName.contains("_psr") 656 &&!srcName.contains("_nod") 657 &&!srcName.contains("_fs") 658 &&!srcName.contains("_fsr")) { 659 // if Nod mode observation , append '_nod' 660 if (obsMode1 == "Nod") { 661 srcName.append("_nod"); 662 } else if (obsMode1 == "OffOn") { 663 // for GBT position switch observations (OffOn or OnOff) 664 if (obsMode2 == "PSWITCHON") srcName.append("_ps"); 665 if (obsMode2 == "PSWITCHOFF") srcName.append("_psr"); 666 } else { 667 if (obsMode2 == "FSWITCH") { 668 // for GBT frequency switch mode 669 if (sigState) srcName.append("_fs"); 670 if (refState) srcName.append("_fsr"); 701 //cerr <<"obsMode2= "<<obsMode2<<endl; 702 if (!srcName.contains("_ps") 703 &&!srcName.contains("_psr") 704 &&!srcName.contains("_nod") 705 &&!srcName.contains("_fs") 706 &&!srcName.contains("_fsr")) { 707 // if Nod mode observation , append '_nod' 708 if (obsMode1 == "Nod") { 709 srcName.append("_nod"); 710 } else if (obsMode1 == "OffOn") { 711 // for GBT position switch observations (OffOn or OnOff) 712 if (obsMode2 == "PSWITCHON") srcName.append("_ps"); 713 if (obsMode2 == "PSWITCHOFF") srcName.append("_psr"); 714 } else { 715 if (obsMode2 == "FSWITCH") { 716 // for GBT frequency switch mode 717 if (sigState) srcName.append("_fs"); 718 if (refState) srcName.append("_fsr"); 719 } 720 } 671 721 } 672 } 673 } 674 } 675 } 722 } 723 } 724 } 725 } 676 726 // CAL state 677 727 // this should be apply just for GBT data? 678 728 Double Cal; 679 if (stateId==-1 ) {729 if (stateId==-1 || StateNRow==0) { 680 730 Cal = 0; 681 731 } else { … … 694 744 Vector<Double> chanFreq = cChanFreqCol(iIF); 695 745 if (nChan == 1) { 696 cout << "The input is continuum data. "<< endl;697 746 freqInc = chanFreq(0); 698 747 refFreq = chanFreq(0); 699 748 restFreq = 0.0f; 700 749 } else { 750 701 751 if (cStartChan(iIF) <= cEndChan(iIF)) { 702 752 freqInc = chanFreq(1) - chanFreq(0); … … 705 755 } 706 756 refFreq = chanFreq(cRefChan(iIF)-1); 707 restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0)); 757 Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId); 758 if (HaveSrcRestFreq) { 759 //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0)); 760 restFreq = cSrcRestFrqCol(srcId); 761 } else { 762 restFreq = 0.0f; 763 } 708 764 } 709 765 bandwidth = abs(freqInc * nChan); … … 720 776 721 777 // Find the appropriate entry in the WEATHER subtable. 722 Vector<Double> wTimes = cWeatherTimeCol.getColumn(); 723 Int weatherIdx; 724 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) { 725 if (cWeatherTimeCol(weatherIdx) <= time) { 726 break; 727 } 728 } 729 730 if (weatherIdx < 0) { 778 //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName()); 779 Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName()); 780 Int weatherIdx=-1; 781 if (cHaveWeatherTab) { 782 Vector<Double> wTimes = cWeatherTimeCol.getColumn(); 783 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) { 784 if (cWeatherTimeCol(weatherIdx) <= time) { 785 break; 786 } 787 } 788 } 789 if (weatherIdx < 0 || !cHaveWeatherTab) { 731 790 // No appropriate WEATHER entry. 732 791 pressure = 0.0f; … … 874 933 baseSub.resize(0,0); 875 934 } 876 877 935 // Get spectral data. 878 936 if (cGetSpectra) { 879 937 Matrix<Float> tmpData; 880 938 Matrix<Bool> tmpFlag; 881 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True); 939 if (cHaveDataCol) { 940 Matrix<Complex> tmpCmplxData; 941 Matrix<Float> tmpReData; 942 Matrix<Float> tmpImData; 943 //cerr<<"reading spectra..."<<endl; 944 //# TODO - should have a flag to user to select DATA or CORRECTED_DATA 945 //# currently just automatically determined, --- read CORRECTED one 946 //# if the column exist. 947 if (cHaveCorrectedDataCol) { 948 cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True); 949 } else { 950 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True); 951 } 952 tmpReData = real(tmpCmplxData); 953 tmpImData = imag(tmpCmplxData); 954 tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData); 955 } else { 956 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True); 957 } 882 958 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True); 883 959 … … 909 985 // Get cross-polarization data. 910 986 if (cGetXPol) { 987 //cerr<<"cGetXPol="<<cGetXPol<<endl; 988 //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl; 989 911 990 if (cHaveXCalFctr) { 912 991 cXCalFctrCol.get(cIdx, xCalFctr); … … 915 994 } 916 995 917 cDataCol.get(cIdx, xPol, True); 918 919 if (cEndChan(iIF) < cStartChan(iIF)) { 920 Complex ctmp; 921 Int jchan = nChan - 1; 922 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) { 923 ctmp = xPol(ichan); 924 xPol(ichan) = xPol(jchan); 925 xPol(jchan) = ctmp; 926 } 927 } 928 } 929 996 if(!cATF) { 997 cDataCol.get(cIdx, xPol, True); 998 999 if (cEndChan(iIF) < cStartChan(iIF)) { 1000 Complex ctmp; 1001 Int jchan = nChan - 1; 1002 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) { 1003 ctmp = xPol(ichan); 1004 xPol(ichan) = xPol(jchan); 1005 xPol(jchan) = ctmp; 1006 } 1007 } 1008 } 1009 } 1010 /** 1011 cerr<<"scanNo="<<scanNo<<endl; 1012 cerr<<"cycleNo="<<cycleNo<<endl; 1013 cerr<<"mjd="<<mjd<<endl; 1014 cerr<<"interval="<<interval<<endl; 1015 cerr<<"fieldName="<<fieldName<<endl; 1016 cerr<<"srcNmae="<<srcName<<endl; 1017 cerr<<"srcDir="<<srcDir<<endl; 1018 cerr<<"srcPM="<<srcPM<<endl; 1019 cerr<<"srcVel="<<srcVel<<endl; 1020 cerr<<"obsMode="<<obsMode<<endl; 1021 cerr<<"IFno="<<IFno<<endl; 1022 cerr<<"refFreq="<<refFreq<<endl; 1023 cerr<<"tcal="<<tcal<<endl; 1024 cerr<<"direction="<<direction<<endl; 1025 cerr<<"scanRate="<<scanRate<<endl; 1026 cerr<<"tsys="<<tsys<<endl; 1027 cerr<<"sigma="<<sigma<<endl; 1028 cerr<<"calFctr="<<calFctr<<endl; 1029 cerr<<"baseLin="<<baseLin<<endl; 1030 cerr<<"baseSub="<<baseSub<<endl; 1031 cerr<<"spectra="<<spectra.shape()<<endl; 1032 cerr<<"flagged="<<flagged.shape()<<endl; 1033 cerr<<"xCalFctr="<<xCalFctr<<endl; 1034 cerr<<"xPol="<<xPol<<endl; 1035 **/ 930 1036 cIdx++; 931 1037 … … 1003 1109 Matrix<Float> tmpData; 1004 1110 Matrix<Bool> tmpFlag; 1005 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True); 1111 if (cHaveDataCol) { 1112 Matrix<Complex> tmpCmplxData; 1113 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True); 1114 tmpData = real(tmpCmplxData); 1115 } else { 1116 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True); 1117 } 1006 1118 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True); 1007 1119 -
branches/alma/external/atnf/PKSIO/PKSMS2reader.h
r1393 r1453 126 126 Double &bandwidth, 127 127 Double &freqInc, 128 Double&restFreq,128 Vector<Double> &restFreq, 129 129 Vector<Float> &tcal, 130 130 String &tcalTime, … … 169 169 private: 170 170 Bool cHaveBaseLin, cHaveCalFctr, cHaveSrcVel, cHaveTsys, cHaveXCalFctr, 171 cMSopen, cHaveTcal ;171 cMSopen, cHaveTcal, cHaveDataCol, cATF, cHaveSysCal, cHaveCorrectedDataCol; 172 172 Int cCycleNo, cIdx, cNRow, cScanNo; 173 173 Double cTime, lastmjd; … … 213 213 ROScalarColumn<Complex> cXCalFctrCol; 214 214 ROArrayColumn<Complex> cDataCol; 215 ROArrayColumn<Complex> cCorrectedDataCol; 215 216 ROScalarColumn<Int> cNumReceptorCol; 216 217 ROScalarColumn<Bool> cSigStateCol; -
branches/alma/external/atnf/PKSIO/PKSMS2writer.cc
r1393 r1453 93 93 94 94 Int maxNPol = max(cNPol); 95 cGBT = cAPEX = cSMT = cALMA = False; 95 96 96 97 // check if it is GBT data 97 cGBT = antName.contains("GBT"); 98 if (antName.contains("GBT")) { 99 cGBT = True; 100 } 101 else if (antName.contains("APEX")) { 102 cAPEX = True; 103 } 104 else if (antName.contains("HHT") || antName.contains("SMT")) { 105 cSMT = True; 106 } 107 else if (antName.contains("ALMA")) { 108 cALMA = True; 109 } 110 111 112 113 //cGBT = antName.contains("GBT"); 114 //cAPEX = antName.contains("APEX"); 115 //cSMT = antName.contains("HHT"); 116 //cALMA = antName.contains("ALMA"); 98 117 99 118 // Add the non-standard CALFCTR column. … … 125 144 126 145 MS::addColumnToDesc(pksDesc, MS::DATA, 2); 146 //pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet(). 147 // define("UNIT", "Jy"); 127 148 pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet(). 128 define("UNIT", "Jy");149 define("UNIT", fluxUnit); 129 150 pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet(). 130 151 define("MEASURE_TYPE", ""); … … 391 412 const Double bandwidth, 392 413 const Double freqInc, 393 const Double restFreq, 414 //const Double restFreq, 415 const Vector<Double> restFreq, 394 416 const Vector<Float> tcal, 395 417 const String tcalTime, … … 533 555 Vector<Float> weight(nPol, 1.0f); 534 556 cMSCols->weight().put(irow, weight); 557 //imaging weight 558 //Vector<Float> imagingWeight(nChan); 559 //cMSCols->imagingWeight().put(irow, imagingWeight); 535 560 536 561 // Flag information. … … 612 637 // do specific things for GBT 613 638 // Data. 639 // plus some more telescopes 614 640 cAntennaCols->name().put(n, antName); 615 641 //cAntennaCols->station().put(n, "ATNF_PARKES"); 616 642 if (cGBT) { 617 643 cAntennaCols->station().put(n, "GREENBANK"); 644 cAntennaCols->dishDiameter().put(n, 110.0); 645 } 646 else if (cAPEX) { 647 cAntennaCols->station().put(n, "CHAJNANTOR"); 648 cAntennaCols->dishDiameter().put(n, 12.0); 649 } 650 else if (cALMA) { 651 cAntennaCols->station().put(n, "CHAJNANTOR"); 652 cAntennaCols->dishDiameter().put(n, 12.0); 653 } 654 else if (cSMT) { 655 cAntennaCols->station().put(n, "MT_GRAHAM"); 656 cAntennaCols->dishDiameter().put(n, 10.0); 618 657 } 619 658 else { 620 659 cAntennaCols->station().put(n, "ATNF_PARKES"); 660 cAntennaCols->dishDiameter().put(n, 64.0); 621 661 } 622 662 cAntennaCols->type().put(n, "GROUND-BASED"); … … 626 666 cAntennaCols->offset().put(n, antOffset); 627 667 //cAntennaCols->dishDiameter().put(n, 64.0); 628 if (cGBT) {629 cAntennaCols->dishDiameter().put(n, 110.0);630 }631 else {632 cAntennaCols->dishDiameter().put(n, 64.0);633 }668 //if (cGBT) { 669 // cAntennaCols->dishDiameter().put(n, 110.0); 670 //} 671 //else { 672 // cAntennaCols->dishDiameter().put(n, 64.0); 673 //} 634 674 // Flags. 635 675 cAntennaCols->flagRow().put(n, False); … … 950 990 const Vector<Double> direction, 951 991 const Vector<Double> properMotion, 952 const Double restFreq, 992 //const Double restFreq, 993 const Vector<Double> restFreq, 953 994 const Double radialVelocity) 954 995 { … … 984 1025 // cSourceCols->position().put(n, position); 985 1026 cSourceCols->properMotion().put(n, properMotion); 986 Vector<Double> restFrequency(1, restFreq); 987 cSourceCols->restFrequency().put(n, restFrequency); 1027 // Vector<Double> restFrequency(1, restFreq); 1028 // cSourceCols->restFrequency().put(n, restFrequency); 1029 cSourceCols->restFrequency().put(n, restFreq); 988 1030 Vector<Double> sysvel(1, radialVelocity); 989 1031 cSourceCols->sysvel().put(n, sysvel); -
branches/alma/external/atnf/PKSIO/PKSMS2writer.h
r1393 r1453 88 88 const Double bandwidth, 89 89 const Double freqInc, 90 const Double restFreq, 90 //const Double restFreq, 91 const Vector<Double> restFreq, 91 92 const Vector<Float> tcal, 92 93 const String tcalTime, … … 166 167 ScalarColumn<Complex> *cXCalFctrCol; 167 168 168 // for GBT specific data handling169 // for handling parameters specific to GBT and other telescopes 169 170 Bool cGBT; 171 Bool cSMT; 172 Bool cAPEX; 173 Bool cALMA; 170 174 171 175 // Add an entry to the ANTENNA subtable. … … 225 229 const Vector<Double> direction, 226 230 const Vector<Double> properMotion, 227 const Double restFreq, 231 //const Double restFreq, 232 const Vector<Double> restFreq, 228 233 const Double radialVelocity); 229 234 -
branches/alma/external/atnf/PKSIO/PKSSDwriter.cc
r1393 r1453 135 135 const Double bandwidth, 136 136 const Double freqInc, 137 const Double restFreq, 137 //const Double restFreq, 138 const Vector<Double> restFreq, 138 139 const Vector<Float> tcal, 139 140 const String tcalTime, … … 211 212 mbrec.srcDec = srcDir(1); 212 213 213 mbrec.restFreq = restFreq; 214 //mbrec.restFreq = restFreq; 215 mbrec.restFreq = restFreq(0); 214 216 215 217 strncpy(mbrec.obsType, (char *)obsMode.chars(), 16); -
branches/alma/external/atnf/PKSIO/PKSSDwriter.h
r1393 r1453 88 88 const Double bandwidth, 89 89 const Double freqInc, 90 const Double restFreq, 90 //const Double restFreq, 91 const Vector<Double> restFreq, 91 92 const Vector<Float> tcal, 92 93 const String tcalTime, -
branches/alma/external/atnf/PKSIO/PKSreader.cc
r1325 r1453 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id : PKSreader.cc,v 19.4 2006/05/19 02:14:50 mcalabre Exp$28 //# $Id$ 29 29 //#--------------------------------------------------------------------------- 30 30 //# Original: 2000/08/23, Mark Calabretta, ATNF -
branches/alma/external/atnf/PKSIO/PKSreader.h
r1393 r1453 148 148 Double &bandwidth, 149 149 Double &freqInc, 150 Double&restFreq,150 Vector<Double> &restFreq, 151 151 Vector<Float> &tcal, 152 152 String &tcalTime, -
branches/alma/external/atnf/PKSIO/PKSwriter.h
r1393 r1453 81 81 const Double bandwidth, 82 82 const Double freqInc, 83 const Double restFreq, 83 //const Double restFreq, 84 const Vector<Double> restFreq, 84 85 const Vector<Float> tcal, 85 86 const String tcalTime, -
branches/alma/external/atnf/PKSIO/SDFITSreader.cc
r1325 r1453 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id : SDFITSreader.cc,v 19.22 2006/07/12 00:14:26 mcalabre Exp$28 //# $Id$ 29 29 //#--------------------------------------------------------------------------- 30 30 //# The SDFITSreader class reads single dish FITS files such as those written -
branches/alma/external/atnf/PKSIO/SDFITSreader.h
r1325 r1453 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id : SDFITSreader.h,v 19.12 2006/07/05 05:24:35 mcalabre Exp$28 //# $Id$ 29 29 //#--------------------------------------------------------------------------- 30 30 //# The SDFITSreader class reads single dish FITS files such as those written -
branches/alma/external/atnf/PKSIO/SDFITSwriter.cc
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : SDFITSwriter.cc,v 19.10 2006/07/05 05:44:52 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# Original: 2000/07/24, Mark Calabretta, ATNF -
branches/alma/external/atnf/PKSIO/SDFITSwriter.h
r1325 r1453 27 27 //# AUSTRALIA 28 28 //# 29 //# $Id : SDFITSwriter.h,v 19.5 2006/05/19 02:19:58 mcalabre Exp$29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 31 //# Original: 2000/07/24, Mark Calabretta, ATNF
Note:
See TracChangeset
for help on using the changeset viewer.