- Timestamp:
- 12/07/05 14:14:45 (19 years ago)
- Location:
- branches/Release12/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/Release12/src/SDMath.cc
r716 r779 101 101 } 102 102 103 SDMath& SDMath::operator=(const SDMath& other) 103 SDMath& SDMath::operator=(const SDMath& other) 104 104 { 105 105 if (this != &other) { … … 113 113 114 114 115 SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in, 116 const String& refTime, 117 const String& method, 118 115 SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in, 116 const String& refTime, 117 const String& method, 118 Bool perFreqID) 119 119 { 120 120 // Get frame info from Table … … 136 136 137 137 138 CountedPtr<SDMemTable> 138 CountedPtr<SDMemTable> 139 139 SDMath::average(const std::vector<CountedPtr<SDMemTable> >& in, 140 141 140 const Vector<Bool>& mask, Bool scanAv, 141 const String& weightStr, Bool alignFreq) 142 142 // Weighted averaging of spectra from one or more Tables. 143 143 { 144 // Convert weight type 144 // Convert weight type 145 145 WeightType wtType = NONE; 146 146 convertWeightString(wtType, weightStr, True); … … 242 242 // Should check that the frequency tables don't change if doing 243 243 // FreqAlignment 244 244 245 245 // Attach columns to Table 246 246 const Table& tabIn = in[iTab]->table(); … … 258 258 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape(); 259 259 if (!shp.isEqual(shp2)) { 260 260 delete pTabOut; 261 261 throw (AipsError("Shapes for all rows must be the same")); 262 262 } 263 263 264 265 264 // If we are not doing scan averages, make checks for source 265 // and frequency setup and warn if averaging across them 266 266 scanIDCol.getScalar(iRow, scanID); 267 267 268 268 // Get quantities from columns 269 269 srcNameCol.getScalar(iRow, sourceName); 270 270 mjdCol.get(iRow, time); … … 273 273 fqIDCol.get(iRow, freqID); 274 274 275 275 // Initialize first source and freqID 276 276 if (iRow==0 && iTab==0) { 277 277 sourceNameStart = sourceName; … … 279 279 } 280 280 281 282 283 284 285 281 // If we are doing scan averages, see if we are at the end of 282 // an accumulation period (scan). We must check soutce names 283 // too, since we might have two tables with one scan each but 284 // different source names; we shouldn't average different 285 // sources together 286 286 if (scanAv && ( (scanID != oldScanID) || 287 287 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) { 288 288 289 290 291 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 292 293 294 289 // Normalize data in 'sum' accumulation array according to 290 // weighting scheme 291 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 292 asap::ChanAxis, nAxesSub); 293 294 // Get ScanContainer for the first row of this averaged Scan 295 295 SDContainer scOut = in[iTab]->getSDContainer(rowStart); 296 296 297 298 299 297 // Fill scan container. The source and freqID come from the 298 // first row of the first table that went into this average 299 // ( should be the same for all rows in the scan average) 300 300 301 301 Float nR(nAccum); 302 302 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID, 303 304 305 303 timeSum/nR, intSum, sourceNameStart, freqIDStart); 304 305 // Write container out to Table 306 306 pTabOut->putSDContainer(scOut); 307 308 // Reset accumulators 307 308 // Reset accumulators 309 309 sum = 0.0; 310 310 sumSq = 0.0; … … 315 315 timeSum = 0.0; 316 316 intSum = 0.0; 317 318 319 317 nPts = 0.0; 318 319 // Increment 320 320 rowStart = iRow; // First row for next accumulation 321 321 tableStart = iTab; // First table for next accumulation 322 322 sourceNameStart = sourceName; // First source name for next 323 323 // accumulation 324 324 freqIDStart = freqID; // First FreqID for next accumulation 325 325 326 326 oldScanID = scanID; 327 327 outScanID += 1; // Scan ID for next 328 328 // accumulation period 329 329 } 330 330 331 332 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, 333 tSysSum, tSysSqSum, tSys, 334 nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis, 331 // Accumulate 332 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, 333 tSysSum, tSysSqSum, tSys, 334 nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis, 335 335 nAxesSub, useMask, wtType); 336 337 336 oldSourceName = sourceName; 337 oldFreqID = freqID; 338 338 } 339 339 } … … 347 347 // scheme 348 348 349 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 350 349 normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, 350 asap::ChanAxis, nAxesSub); 351 351 352 352 // Create and fill container. The container we clone will be from … … 356 356 SDContainer scOut = in[tableStart]->getSDContainer(rowStart); 357 357 fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID, 358 358 timeSum/nR, intSum, sourceNameStart, freqIDStart); 359 359 pTabOut->putSDContainer(scOut); 360 360 pTabOut->resetCursor(); … … 365 365 366 366 367 CountedPtr<SDMemTable> 367 CountedPtr<SDMemTable> 368 368 SDMath::binaryOperate(const CountedPtr<SDMemTable>& left, 369 370 369 const CountedPtr<SDMemTable>& right, 370 const String& op, Bool preserve, Bool doTSys) 371 371 { 372 372 … … 399 399 } 400 400 401 // Input Tables 401 // Input Tables 402 402 const Table& tLeft = left->table(); 403 403 const Table& tRight = right->table(); … … 413 413 Array<Float> tSysLeftArr, tSysRightArr; 414 414 if (doTSys) tSysRightCol.get(0, tSysRightArr); 415 MaskedArray<Float>* pMRight = 415 MaskedArray<Float>* pMRight = 416 416 new MaskedArray<Float>(right->rowAsMaskedArray(0)); 417 417 … … 424 424 // Loop over rows 425 425 for (uInt i=0; i<nRowLeft; i++) { 426 427 // Get data 426 427 // Get data 428 428 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i)); 429 429 IPosition shpLeft = mLeft.shape(); 430 430 if (doTSys) tSysLeftCol.get(i, tSysLeftArr); 431 431 432 432 if (nRowRight>1) { 433 433 delete pMRight; … … 444 444 if (doTSys) { 445 445 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) { 446 447 448 446 delete pTabOut; 447 delete pMRight; 448 throw(AipsError("left and right Tsys data are not conformant")); 449 449 } 450 450 if (!shpRight.isEqual(tSysRightArr.shape())) { 451 452 453 451 delete pTabOut; 452 delete pMRight; 453 throw(AipsError("left and right scan tables are not conformant")); 454 454 } 455 455 } … … 459 459 460 460 // Operate on data and TSys 461 if (what==0) { 461 if (what==0) { 462 462 MaskedArray<Float> tmp = mLeft + *pMRight; 463 463 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); … … 476 476 if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr); 477 477 } else if (what==4) { 478 if (preserve) { 479 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 480 481 478 if (preserve) { 479 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 480 tSysRightArr; 481 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 482 482 } else { 483 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 484 485 483 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - 484 tSysLeftArr; 485 putDataInSDC(sc, tmp.getArray(), tmp.getMask()); 486 486 } 487 487 sc.putTsys(tSysRightArr); … … 493 493 if (pMRight) delete pMRight; 494 494 pTabOut->resetCursor(); 495 495 496 496 return CountedPtr<SDMemTable>(pTabOut); 497 497 } … … 499 499 500 500 std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in, 501 502 501 const Vector<Bool>& mask, 502 const String& which, Int row) const 503 503 // 504 504 // Perhaps iteration over pol/beam/if should be in here … … 519 519 uInt iStart = 0; 520 520 uInt iEnd = in->nRow()-1; 521 // 521 // 522 522 if (row>=0) { 523 523 iStart = row; … … 569 569 CoordinateSystem cSys; 570 570 cSys.addCoordinate(in.getSpectralCoordinate(j)); 571 CoordinateSystem cSysBin = 571 CoordinateSystem cSysBin = 572 572 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False); 573 573 // … … 584 584 585 585 // Loop over rows and bin along channel axis 586 586 587 587 for (uInt i=0; i < in.nRow(); ++i) { 588 588 SDContainer sc = in.getSDContainer(i); … … 603 603 putDataInSDC(sc, marrout.getArray(), marrout.getMask()); 604 604 605 // Bin up Tsys. 605 // Bin up Tsys. 606 606 607 607 Array<Bool> allGood(tSys.shape(),True); 608 608 MaskedArray<Float> tSysIn(tSys, allGood, True); 609 609 // 610 MaskedArray<Float> tSysOut; 610 MaskedArray<Float> tSysOut; 611 611 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width); 612 612 sc.putTsys(tSysOut.getArray()); … … 618 618 619 619 SDMemTable* SDMath::resample(const SDMemTable& in, const String& methodStr, 620 620 Float width) 621 621 // 622 622 // Should add the possibility of width being specified in km/s. This means 623 // that for each freqID (SpectralCoordinate) we will need to convert to an 624 // average channel width (say at the reference pixel). Then we would need 625 // to be careful to make sure each spectrum (of different freqID) 623 // that for each freqID (SpectralCoordinate) we will need to convert to an 624 // average channel width (say at the reference pixel). Then we would need 625 // to be careful to make sure each spectrum (of different freqID) 626 626 // is the same length. 627 627 // … … 632 632 SpectralCoordinate sC = in.getSpectralCoordinate(j); 633 633 } 634 } 634 } 635 635 636 636 // Interpolation method … … 691 691 692 692 Array<Float> valuesOut; 693 Array<Bool> maskOut; 693 Array<Bool> maskOut; 694 694 Array<Float> tSysOut; 695 695 Array<Bool> tSysMaskIn(shapeIn,True); … … 702 702 703 703 // Get data and Tsys 704 704 705 705 const Array<Float>& tSysIn = sc.getTsys(); 706 706 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i)); … … 710 710 // Interpolate data 711 711 712 InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut, 712 InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut, 713 713 xIn, valuesIn, maskIn, 714 714 interpMethod, True, True); … … 718 718 // Interpolate TSys 719 719 720 InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut, 720 InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut, 721 721 xIn, tSysIn, tSysMaskIn, 722 722 interpMethod, True, True); … … 739 739 SDMemTable* pOut = new SDMemTable(in,False); 740 740 const Table& tOut = pOut->table(); 741 ArrayColumn<Float> specCol(tOut,"SPECTRA"); 742 ArrayColumn<Float> tSysCol(tOut,"TSYS"); 741 ArrayColumn<Float> specCol(tOut,"SPECTRA"); 742 ArrayColumn<Float> tSysCol(tOut,"TSYS"); 743 743 Array<Float> tSysArr; 744 744 … … 797 797 const Table& tabIn = in.table(); 798 798 799 // Shape of input and output data 799 // Shape of input and output data 800 800 801 801 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape(); … … 826 826 tSysCol.attach(tabIn,"TSYS"); 827 827 } 828 // 828 // 829 829 const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis)); 830 830 … … 838 838 Array<Float>& arr = marr.getRWArray(); 839 839 const Array<Bool>& barr = marr.getMask(); 840 840 841 841 // Get Tsys 842 842 … … 852 852 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes); 853 853 ReadOnlyArrayIterator<Float>* pItTsysPlane = 0; 854 if (wtType==TSYS) 855 854 if (wtType==TSYS) 855 pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes); 856 856 857 857 // Accumulations … … 868 868 Vector<Float> t1(nChan); t1 = 0.0; 869 869 Vector<Bool> t2(nChan); t2 = True; 870 870 Float tSys = 0.0; 871 871 MaskedArray<Float> vecSum(t1,t2); 872 872 Float norm = 0.0; … … 876 876 // 877 877 ReadOnlyVectorIterator<Float>* pItTsysVec = 0; 878 879 pItTsysVec = 880 881 } 882 // 883 while (!itDataVec.pastEnd()) { 878 if (wtType==TSYS) { 879 pItTsysVec = 880 new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1); 881 } 882 // 883 while (!itDataVec.pastEnd()) { 884 884 885 885 // Create MA of data & mask (optionally including OTF mask) and get variance for this spectrum … … 887 887 if (useMask) { 888 888 const MaskedArray<Float> spec(itDataVec.vector(), 889 889 mask&&itMaskVec.vector()); 890 890 if (wtType==VAR) { 891 891 fac = 1.0 / variance(spec); 892 892 } else if (wtType==TSYS) { 893 893 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency 894 894 fac = 1.0 / tSys / tSys; 895 } 895 } 896 896 } else { 897 897 const MaskedArray<Float> spec(itDataVec.vector(), 898 898 itMaskVec.vector()); 899 899 if (wtType==VAR) { 900 900 fac = 1.0 / variance(spec); 901 901 } else if (wtType==TSYS) { 902 902 tSys = pItTsysVec->vector()[0]; // Drop pseudo channel dependency 903 903 fac = 1.0 / tSys / tSys; 904 904 } … … 908 908 909 909 const MaskedArray<Float> spec(fac*itDataVec.vector(), 910 910 itMaskVec.vector()); 911 911 vecSum += spec; 912 912 norm += fac; … … 916 916 itDataVec.next(); 917 917 itMaskVec.next(); 918 918 if (wtType==TSYS) pItTsysVec->next(); 919 919 } 920 920 921 921 // Clean up 922 922 923 924 delete pItTsysVec; 923 if (pItTsysVec) { 924 delete pItTsysVec; 925 925 pItTsysVec = 0; 926 } 926 } 927 927 } 928 928 … … 936 936 IPosition pos = itDataPlane.pos(); 937 937 938 // Write out data. This is a bit messy. We have to reform the Vector 938 // Write out data. This is a bit messy. We have to reform the Vector 939 939 // accumulator into an Array of shape (1,1,1,nChan) 940 940 941 941 start = pos; 942 end = pos; 942 end = pos; 943 943 end(asap::ChanAxis) = nChan-1; 944 944 outData(start,end) = vecSum.getArray().reform(vecShapeOut); … … 949 949 itDataPlane.next(); 950 950 itMaskPlane.next(); 951 951 if (wtType==TSYS) pItTsysPlane->next(); 952 952 } 953 953 … … 974 974 975 975 976 SDMemTable* SDMath::smooth(const SDMemTable& in, 977 978 976 SDMemTable* SDMath::smooth(const SDMemTable& in, 977 const casa::String& kernelType, 978 casa::Float width, Bool doAll) 979 979 // 980 980 // Should smooth TSys as well … … 1014 1014 Array<Bool> maskIn = dataIn.getMask(); 1015 1015 1016 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn 1016 Array<Float> valuesIn2 = valuesIn(start,end); // ref to valuesIn 1017 1017 Array<Bool> maskIn2 = maskIn(start,end); 1018 1018 … … 1021 1021 VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis); 1022 1022 while (!itValues.pastEnd()) { 1023 1024 1025 1026 mathutil::hanning(valuesOut, maskOut, itValues.vector(), 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1023 1024 // Smooth 1025 if (kernelType==VectorKernel::HANNING) { 1026 mathutil::hanning(valuesOut, maskOut, itValues.vector(), 1027 itMask.vector()); 1028 itMask.vector() = maskOut; 1029 } else { 1030 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector()); 1031 conv.linearConv(valuesOut, itValues.vector()); 1032 } 1033 1034 itValues.vector() = valuesOut; 1035 itValues.next(); 1036 itMask.next(); 1037 1037 } 1038 1038 1039 1039 // Create and put back 1040 1040 SDContainer sc = in.getSDContainer(ri); … … 1049 1049 1050 1050 1051 SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp, 1052 1053 // 1051 SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp, 1052 Float JyPerK, Bool doAll) 1053 // 1054 1054 // etaAp = aperture efficiency (-1 means find) 1055 1055 // D = geometric diameter (m) (-1 means find) 1056 // JyPerK 1056 // JyPerK 1057 1057 // 1058 1058 { … … 1063 1063 // mJy or mK) Also automatically find out what we are converting to 1064 1064 // according to the flux unit 1065 Unit fluxUnit(sh.fluxunit); 1065 Unit fluxUnit(sh.fluxunit); 1066 1066 Unit K(String("K")); 1067 1067 Unit JY(String("Jy")); 1068 1068 1069 1069 Bool toKelvin = True; 1070 Double cFac = 1.0; 1070 Double cFac = 1.0; 1071 1071 1072 1072 if (fluxUnit==JY) { … … 1075 1075 Quantum<Double> t2 = t.get(JY); 1076 1076 cFac = (t2 / t).getValue(); // value to Jy 1077 1077 1078 1078 toKelvin = True; 1079 1079 sh.fluxunit = "K"; … … 1083 1083 Quantum<Double> t2 = t.get(K); 1084 1084 cFac = (t2 / t).getValue(); // value to K 1085 1085 1086 1086 toKelvin = False; 1087 1087 sh.fluxunit = "Jy"; … … 1091 1091 1092 1092 pTabOut->putSDHeader(sh); 1093 1093 1094 1094 // Make sure input values are converted to either Jy or K first... 1095 1095 Float factor = cFac; … … 1121 1121 scaleByVector(pTabOut, in, doAll, factors, False); 1122 1122 } else { 1123 1123 1124 1124 // OK now we must deal with automatic look up of values. 1125 1125 // We must also deal with the fact that the factors need 1126 1126 // to be computed per IF and may be different and may 1127 1127 // change per integration. 1128 1128 1129 1129 pushLog("Looking up conversion factors"); 1130 1130 convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll); … … 1134 1134 1135 1135 1136 SDMemTable* SDMath::gainElevation(const SDMemTable& in, 1137 1138 1139 1136 SDMemTable* SDMath::gainElevation(const SDMemTable& in, 1137 const Vector<Float>& coeffs, 1138 const String& fileName, 1139 const String& methodStr, Bool doAll) 1140 1140 { 1141 1141 … … 1149 1149 Vector<Float> x = elev.getColumn(); 1150 1150 x *= Float(180 / C::pi); // Degrees 1151 1151 1152 1152 const uInt nC = coeffs.nelements(); 1153 1153 if (fileName.length()>0 && nC>0) { 1154 1154 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both")); 1155 1155 } 1156 1156 1157 1157 // Correct 1158 1158 if (nC>0 || fileName.length()==0) { … … 1160 1160 Bool throwIt = True; 1161 1161 Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt); 1162 1162 1163 1163 // Set polynomial 1164 1164 Polynomial<Float>* pPoly = 0; … … 1175 1175 msg = String("built in"); 1176 1176 } 1177 1177 1178 1178 if (coeff.nelements()>0) { 1179 1179 pPoly->setCoefficients(coeff); … … 1182 1182 throw(AipsError("There is no known gain-elevation polynomial known for this instrument")); 1183 1183 } 1184 pushLog("Making polynomial correction with "+msg+" coefficients"); 1184 ostringstream oss; 1185 oss << "Making polynomial correction with " << msg << " coefficients:" << endl; 1186 oss << " " << coeff; 1187 pushLog(String(oss)); 1185 1188 const uInt nRow = in.nRow(); 1186 1189 Vector<Float> factor(nRow); … … 1192 1195 1193 1196 } else { 1194 1197 1195 1198 // Indicate which columns to read from ascii file 1196 1199 String col0("ELEVATION"); 1197 1200 String col1("FACTOR"); 1198 1201 1199 1202 // Read and correct 1200 1203 1201 1204 pushLog("Making correction from ascii Table"); 1202 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1, 1203 1205 scaleFromAsciiTable (pTabOut, in, fileName, col0, col1, 1206 methodStr, doAll, x, True); 1204 1207 } 1205 1208 1206 1209 return pTabOut; 1207 1210 } 1208 1211 1209 1212 1210 1213 SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll) … … 1256 1259 throw(AipsError("Only linear polarizations are supported")); 1257 1260 } 1258 // 1261 // 1259 1262 const Table& tabIn = in.table(); 1260 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1263 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1261 1264 IPosition start(asap::nAxes,0); 1262 1265 IPosition end(asap::nAxes); 1263 1266 1264 1267 // Set cursor slice. Assumes shape the same for all rows 1265 1268 1266 1269 setCursorSlice (start, end, doAll, in); 1267 1270 IPosition start3(start); 1268 1271 start3(asap::PolAxis) = 2; // Real(XY) 1269 1272 IPosition end3(end); 1270 end3(asap::PolAxis) = 2; 1273 end3(asap::PolAxis) = 2; 1271 1274 // 1272 1275 IPosition start4(start); … … 1274 1277 IPosition end4(end); 1275 1278 end4(asap::PolAxis) = 3; 1276 // 1279 // 1277 1280 uInt nRow = in.nRow(); 1278 1281 Array<Float> data; … … 1280 1283 specCol.get(i,data); 1281 1284 IPosition shape = data.shape(); 1282 1285 1283 1286 // Get polarization slice references 1284 1287 Array<Float> C3 = data(start3,end3); 1285 1288 Array<Float> C4 = data(start4,end4); 1286 1289 1287 1290 // Rotate 1288 1291 SDPolUtil::rotatePhase(C3, C4, value); 1289 1292 1290 1293 // Put 1291 1294 specCol.put(i,data); 1292 1295 } 1293 } 1296 } 1294 1297 1295 1298 … … 1310 1313 throw(AipsError("Only linear polarizations are supported")); 1311 1314 } 1312 // 1315 // 1313 1316 const Table& tabIn = in.table(); 1314 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1315 ROArrayColumn<Float> stokesCol(tabIn,"STOKES"); 1317 ArrayColumn<Float> specCol(tabIn,"SPECTRA"); 1318 ROArrayColumn<Float> stokesCol(tabIn,"STOKES"); 1316 1319 IPosition start(asap::nAxes,0); 1317 1320 IPosition end(asap::nAxes); 1318 1321 1319 1322 // Set cursor slice. Assumes shape the same for all rows 1320 1323 1321 1324 setCursorSlice (start, end, doAll, in); 1322 1325 // … … 1324 1327 start1(asap::PolAxis) = 0; // C1 (XX) 1325 1328 IPosition end1(end); 1326 end1(asap::PolAxis) = 0; 1329 end1(asap::PolAxis) = 0; 1327 1330 // 1328 1331 IPosition start2(start); 1329 1332 start2(asap::PolAxis) = 1; // C2 (YY) 1330 1333 IPosition end2(end); 1331 end2(asap::PolAxis) = 1; 1334 end2(asap::PolAxis) = 1; 1332 1335 // 1333 1336 IPosition start3(start); 1334 1337 start3(asap::PolAxis) = 2; // C3 ( Real(XY) ) 1335 1338 IPosition end3(end); 1336 end3(asap::PolAxis) = 2; 1339 end3(asap::PolAxis) = 2; 1337 1340 // 1338 1341 IPosition startI(start); 1339 1342 startI(asap::PolAxis) = 0; // I 1340 1343 IPosition endI(end); 1341 endI(asap::PolAxis) = 0; 1344 endI(asap::PolAxis) = 0; 1342 1345 // 1343 1346 IPosition startQ(start); 1344 1347 startQ(asap::PolAxis) = 1; // Q 1345 1348 IPosition endQ(end); 1346 endQ(asap::PolAxis) = 1; 1349 endQ(asap::PolAxis) = 1; 1347 1350 // 1348 1351 IPosition startU(start); 1349 1352 startU(asap::PolAxis) = 2; // U 1350 1353 IPosition endU(end); 1351 endU(asap::PolAxis) = 2; 1354 endU(asap::PolAxis) = 2; 1352 1355 1353 1356 // … … 1358 1361 stokesCol.get(i,stokes); 1359 1362 IPosition shape = data.shape(); 1360 1363 1361 1364 // Get linear polarization slice references 1362 1365 1363 1366 Array<Float> C1 = data(start1,end1); 1364 1367 Array<Float> C2 = data(start2,end2); … … 1370 1373 Array<Float> Q = stokes(startQ,endQ); 1371 1374 Array<Float> U = stokes(startU,endU); 1372 1375 1373 1376 // Rotate 1374 1377 1375 1378 SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value); 1376 1379 1377 1380 // Put 1378 1381 1379 1382 specCol.put(i,data); 1380 1383 } 1381 } 1384 } 1382 1385 1383 1386 // 'private' functions 1384 1387 1385 void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in, 1388 void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in, 1386 1389 Bool toKelvin, Float cFac, Bool doAll) 1387 1390 { … … 1438 1441 // Loop over rows and apply correction factor 1439 1442 1440 Float factor = 1.0; 1443 Float factor = 1.0; 1441 1444 const uInt axis = asap::ChanAxis; 1442 1445 for (uInt i=0; i < in.nRow(); ++i) { … … 1457 1460 1458 1461 // Now the conversion factor depends only upon frequency 1459 // So we need to iterate through by IF only giving 1462 // So we need to iterate through by IF only giving 1460 1463 // us BEAM/POL/CHAN cubes 1461 1464 … … 1479 1482 SDMemTable* SDMath::frequencyAlign(const SDMemTable& in, 1480 1483 MFrequency::Types freqSystem, 1481 const String& refTime, 1484 const String& refTime, 1482 1485 const String& methodStr, 1483 1486 Bool perFreqID) … … 1501 1504 Vector<Double> times = mjdCol.getColumn(); 1502 1505 1503 // Generate DataDesc table 1504 1506 // Generate DataDesc table 1507 1505 1508 Matrix<uInt> ddIdx; 1506 1509 SDDataDesc dDesc; 1507 generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol, 1508 1510 generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol, 1511 fqIDCol, perFreqID); 1509 1512 1510 1513 // Get reference Epoch to time of first row or given String 1511 1514 1512 1515 Unit DAY(String("d")); 1513 1516 MEpoch::Ref epochRef(in.getTimeReference()); … … 1518 1521 refEpoch = in.getEpoch(0); 1519 1522 } 1520 pushLog("Aligning at reference Epoch "+formatEpoch(refEpoch) 1521 +" in frame "+MFrequency::showType(freqSystem)); 1522 1523 ostringstream oss; 1524 oss << "Aligned at reference Epoch " << formatEpoch(refEpoch) << " in frame " << MFrequency::showType(freqSystem); 1525 pushLog(String(oss)); 1526 1523 1527 // Get Reference Position 1524 1528 1525 1529 MPosition refPos = in.getAntennaPosition(); 1526 1530 1527 1531 // Create FrequencyAligner Block. One FA for each possible 1528 1532 // source/freqID (perFreqID=True) or source/IF (perFreqID=False) 1529 1533 // combination 1530 1534 1531 1535 PtrBlock<FrequencyAligner<Float>* > a(dDesc.length()); 1532 generateFrequencyAligners (a, dDesc, in, nChan, freqSystem, refPos,1536 generateFrequencyAligners(a, dDesc, in, nChan, freqSystem, refPos, 1533 1537 refEpoch, perFreqID); 1534 1538 1535 1539 // Generate and fill output Frequency Table. WHen perFreqID=True, 1536 1540 // there is one output FreqID for each entry in the SDDataDesc 1537 1541 // table. However, in perFreqID=False mode, there may be some 1538 1542 // degeneracy, so we need a little translation map 1539 1543 1540 1544 SDFrequencyTable freqTabOut = in.getSDFreqTable(); 1541 1545 freqTabOut.setLength(0); … … 1548 1552 1549 1553 // Get Aligned SC in Hz 1550 1554 1551 1555 SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear); 1552 1556 sC.setWorldAxisUnits(units); 1553 1557 1554 1558 // Add FreqID 1555 1559 1556 1560 uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0], 1557 1558 1561 sC.referenceValue()[0], 1562 sC.increment()[0]); 1559 1563 // output FreqID = ddFQTrans(ddIdx) 1560 1564 ddFQTrans(i) = idx; 1561 1565 } 1562 1566 1563 1567 // Interpolation method 1564 1568 1565 1569 InterpolateArray1D<Double,Float>::InterpolationMethod interp; 1566 1570 convertInterpString(interp, methodStr); 1567 1571 1568 1572 // New output Table 1569 1570 pushLog("Create output table");1573 1574 //pushLog("Create output table"); 1571 1575 SDMemTable* pTabOut = new SDMemTable(in,True); 1572 1576 pTabOut->putSDFreqTable(freqTabOut); 1573 1577 1574 1578 // Loop over rows in Table 1575 1579 1576 1580 Bool extrapolate=False; 1577 1581 const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis); … … 1586 1590 // 1587 1591 for (uInt iRow=0; iRow<nRows; ++iRow) { 1588 if (iRow%10==0) { 1589 pushLog("Processing row "+iRow); 1590 } 1591 1592 // if (iRow%10==0) { 1593 // ostringstream oss; 1594 // oss << "Processing row " << iRow; 1595 // pushLog(String(oss)); 1596 // } 1597 1592 1598 // Get EPoch 1593 1599 1594 1600 Quantum<Double> tQ2(times[iRow],DAY); 1595 1601 MVEpoch mv2(tQ2); 1596 1602 MEpoch epoch(mv2, epochRef); 1597 1603 1598 1604 // Get copy of data 1599 1605 1600 1606 const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow)); 1601 1607 Array<Float> values = mArrIn.getArray(); 1602 Array<Bool> mask = mArrIn.getMask(); 1603 1608 Array<Bool> mask = mArrIn.getMask(); 1609 1604 1610 // For each row, the Frequency abcissa will be the same 1605 1611 // regardless of polarization. For all other axes (IF and BEAM) … … 1607 1613 // pol-chan planes to mimimize the work. Probably won't work for 1608 1614 // multiple beams at this point. 1609 1615 1610 1616 ArrayIterator<Float> itValuesPlane(values, polChanAxes); 1611 1617 ArrayIterator<Bool> itMaskPlane(mask, polChanAxes); … … 1613 1619 1614 1620 // Find the IF index and then the FA PtrBlock index 1615 1621 1616 1622 const IPosition& pos = itValuesPlane.pos(); 1617 1623 ifIdx = pos(asap::IFAxis); 1618 1624 faIdx = ddIdx(iRow,ifIdx); 1619 1625 1620 1626 // Generate abcissa for perIF. Could cache this in a Matrix on 1621 1627 // a per scan basis. Pretty expensive doing it for every row. 1622 1628 1623 1629 if (!perFreqID) { 1624 1625 1626 1630 xIn.resize(nChan); 1631 uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx)); 1632 SpectralCoordinate sC = in.getSpectralCoordinate(fqID); 1627 1633 Double w; 1628 1634 for (uInt i=0; i<nChan; i++) { 1629 1635 sC.toWorld(w,Double(i)); 1630 1636 xIn[i] = w; 1631 1637 } 1632 1638 } 1633 1634 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1); 1639 1640 VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1); 1635 1641 VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1); 1636 1642 1637 1643 // Iterate through the plane by vector and align 1638 1644 1639 1645 first = True; 1640 1646 useCachedAbcissa=False; 1641 while (!itValuesVec.pastEnd()) { 1642 1643 1644 1645 interp, extrapolate); 1646 1647 1648 1649 interp, extrapolate); 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1647 while (!itValuesVec.pastEnd()) { 1648 if (perFreqID) { 1649 ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(), 1650 itMaskVec.vector(), epoch, useCachedAbcissa, 1651 interp, extrapolate); 1652 } else { 1653 ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(), 1654 itMaskVec.vector(), epoch, useCachedAbcissa, 1655 interp, extrapolate); 1656 } 1657 // 1658 itValuesVec.vector() = yOut; 1659 itMaskVec.vector() = maskOut; 1660 // 1661 itValuesVec.next(); 1662 itMaskVec.next(); 1663 // 1664 if (first) { 1665 useCachedAbcissa = True; 1666 first = False; 1667 } 1662 1668 } 1663 1664 1665 1669 // 1670 itValuesPlane.next(); 1671 itMaskPlane.next(); 1666 1672 } 1667 1673 1668 1674 // Create SDContainer and put back 1669 1675 1670 1676 SDContainer sc = in.getSDContainer(iRow); 1671 1677 putDataInSDC(sc, values, mask); 1672 1678 1673 1679 // Set output FreqIDs 1674 1680 1675 1681 for (uInt i=0; i<nIF; i++) { 1676 1682 uInt idx = ddIdx(iRow,i); // Index into SDDataDesc table … … 1681 1687 pTabOut->putSDContainer(sc); 1682 1688 } 1683 1689 1684 1690 // Now we must set the base and extra frames to the input frame 1685 1691 std::vector<string> info = pTabOut->getCoordInfo(); … … 1688 1694 pTabOut->setCoordInfo(info); 1689 1695 1690 // Clean up PointerBlock 1696 // Clean up PointerBlock 1691 1697 for (uInt i=0; i<a.nelements(); i++) delete a[i]; 1692 1698 … … 1704 1710 const Table& tabIn = in.table(); 1705 1711 1706 // Shape of input and output data 1712 // Shape of input and output data 1707 1713 const IPosition& shapeIn = in.rowAsMaskedArray(0).shape(); 1708 1714 … … 1750 1756 1751 1757 void SDMath::fillSDC(SDContainer& sc, 1752 1753 1754 1755 1756 1757 1758 const Array<Bool>& mask, 1759 const Array<Float>& data, 1760 const Array<Float>& tSys, 1761 Int scanID, Double timeStamp, 1762 Double interval, const String& sourceName, 1763 const Vector<uInt>& freqID) 1758 1764 { 1759 1765 // Data and mask … … 1776 1782 1777 1783 void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum, 1778 MaskedArray<Float>& sum, Array<Float>& sumSq, 1779 Array<Float>& nPts, Array<Float>& tSysSum, 1784 MaskedArray<Float>& sum, Array<Float>& sumSq, 1785 Array<Float>& nPts, Array<Float>& tSysSum, 1780 1786 Array<Float>& tSysSqSum, 1781 const Array<Float>& tSys, const Array<Float>& nInc, 1782 1783 1784 uInt iTab, uInt iRow, uInt axis, 1785 1786 1787 const Array<Float>& tSys, const Array<Float>& nInc, 1788 const Vector<Bool>& mask, Double time, Double interval, 1789 const std::vector<CountedPtr<SDMemTable> >& in, 1790 uInt iTab, uInt iRow, uInt axis, 1791 uInt nAxesSub, Bool useMask, 1792 WeightType wtType) 1787 1793 { 1788 1794 … … 1810 1816 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis); 1811 1817 Float fac = 1.0; 1812 IPosition pos(nAxesSub,0); 1818 IPosition pos(nAxesSub,0); 1813 1819 // 1814 1820 while (!itData.pastEnd()) { … … 1846 1852 // We will take a short cut and just use the value from the first channel 1847 1853 // for now. 1848 // 1854 // 1849 1855 VectorIterator<Float> itData(valuesIn, axis); 1850 1856 ReadOnlyVectorIterator<Float> itTSys(tSys, axis); 1851 IPosition pos(nAxesSub,0); 1857 IPosition pos(nAxesSub,0); 1852 1858 // 1853 1859 Float fac = 1.0; … … 1923 1929 } 1924 1930 } else if (wtType==TSYS || wtType==TINTSYS) { 1925 1931 1926 1932 // Normalize each spectrum by sum(1/Tsys**2) (TSYS) or 1927 1933 // sum(Tint/Tsys**2) (TINTSYS) where the pseudo … … 2000 2006 2001 2007 2002 void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type, 2008 void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type, 2003 2009 const casa::String& interp) 2004 2010 { … … 2019 2025 2020 2026 void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data, 2021 2027 const Array<Bool>& mask) 2022 2028 { 2023 2029 sc.putSpectrum(data); … … 2051 2057 } 2052 2058 2053 void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in, 2054 const Table& tTable, const String& col0, 2059 void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in, 2060 const Table& tTable, const String& col0, 2055 2061 const String& col1, 2056 2062 const String& methodStr, Bool doAll, … … 2074 2080 Vector<Float> yOut; 2075 2081 Vector<Bool> maskOut; 2076 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut, 2082 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut, 2077 2083 xIn, yIn, maskIn, intMethod, 2078 2084 True, True); 2079 // Apply 2085 // Apply 2080 2086 2081 2087 scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys); … … 2083 2089 2084 2090 2085 void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in, 2091 void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in, 2086 2092 Bool doAll, const Vector<Float>& factor, 2087 2093 Bool doTSys) … … 2100 2106 2101 2107 // Loop over rows and apply correction factor 2102 2108 2103 2109 const uInt axis = asap::ChanAxis; 2104 2110 for (uInt i=0; i < in.nRow(); ++i) { … … 2193 2199 2194 2200 2195 void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a, 2196 2197 2198 2199 2200 2201 2201 void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a, 2202 const SDDataDesc& dDesc, 2203 const SDMemTable& in, uInt nChan, 2204 MFrequency::Types system, 2205 const MPosition& refPos, 2206 const MEpoch& refEpoch, 2207 Bool perFreqID) 2202 2208 { 2203 2209 for (uInt i=0; i<dDesc.length(); i++) { … … 2208 2214 if (perFreqID) { 2209 2215 2210 // One aligner per source/FreqID pair. 2216 // One aligner per source/FreqID pair. 2211 2217 2212 2218 SpectralCoordinate sC = in.getSpectralCoordinate(ID); … … 2230 2236 return range; 2231 2237 } 2232 2238 2233 2239 2234 2240 Bool SDMath::rowInRange(uInt i, const Vector<uInt>& range) const -
branches/Release12/src/SDMemTable.cc
r773 r779 1703 1703 } 1704 1704 ostringstream oss; 1705 oss << "Replac ingrest frequencies, one per IF, with given list : " << restFreqs;1705 oss << "Replaced rest frequencies, one per IF, with given list : " << restFreqs; 1706 1706 sdft.deleteRestFrequencies(); 1707 1707 for (uInt i=0; i<nRestFreqs; i++) { … … 1715 1715 idx = sdft.addRestFrequency(rf.getValue("Hz")); 1716 1716 if (whichIF>=0) { 1717 oss << "Select inggiven rest frequency (" << restFreqs[0] << ") for IF " << whichIF << endl;1717 oss << "Selected given rest frequency (" << restFreqs[0] << ") for IF " << whichIF << endl; 1718 1718 } else { 1719 oss << "Select inggiven rest frequency (" << restFreqs[0] << ") for all IFs" << endl;1719 oss << "Selected given rest frequency (" << restFreqs[0] << ") for all IFs" << endl; 1720 1720 } 1721 1721 } … … 1839 1839 { 1840 1840 MPosition mp = getAntennaPosition(); 1841 1841 ostringstream oss; 1842 oss << "Computed azimuth/elevation using " << endl 1843 << mp << endl; 1842 1844 for (uInt i=0; i<nRow();++i) { 1843 1845 MEpoch me = getEpoch(i); 1844 1846 MDirection md = getDirection(i,False); 1847 oss << " Time: " << getTime(i,False) << " Direction: " << formatDirection(md) 1848 << endl << " => "; 1845 1849 MeasFrame frame(mp, me); 1846 1850 Vector<Double> azel = … … 1850 1854 azCol_.put(i,azel[0]); 1851 1855 elCol_.put(i,azel[1]); 1852 } 1853 } 1856 oss << "azel: " << azel[0]/C::pi*180.0 << " " 1857 << azel[1]/C::pi*180.0 << " (deg)" << endl; 1858 } 1859 pushLog(String(oss)); 1860 }
Note:
See TracChangeset
for help on using the changeset viewer.