Changeset 1452 for trunk/external/atnf/PKSIO/PKSMS2reader.cc
- Timestamp:
- 11/19/08 20:41:16 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/external/atnf/PKSIO/PKSMS2reader.cc
r1427 r1452 26 26 //# Charlottesville, VA 22903-2475 USA 27 27 //# 28 //# $Id: PKSMS2reader.cc,v 19. 13 2008-06-26 02:08:02cal103 Exp $28 //# $Id: PKSMS2reader.cc,v 19.21 2008-11-17 06:55:18 cal103 Exp $ 29 29 //#--------------------------------------------------------------------------- 30 30 //# Original: 2000/08/03, Mark Calabretta, ATNF 31 31 //#--------------------------------------------------------------------------- 32 32 33 34 // AIPS++ includes. 33 #include <atnf/pks/pks_maths.h> 34 #include <atnf/PKSIO/PKSmsg.h> 35 #include <atnf/PKSIO/PKSrecord.h> 36 #include <atnf/PKSIO/PKSMS2reader.h> 37 35 38 #include <casa/stdio.h> 36 39 #include <casa/Arrays/ArrayMath.h> … … 39 42 #include <tables/Tables.h> 40 43 41 // Parkes includes.42 #include <atnf/pks/pks_maths.h>43 #include <atnf/PKSIO/PKSMS2reader.h>44 45 46 44 //------------------------------------------------- PKSMS2reader::PKSMS2reader 47 45 … … 51 49 { 52 50 cMSopen = False; 51 52 // By default, messages are written to stderr. 53 initMsg(); 53 54 } 54 55 … … 353 354 const Bool getSpectra, 354 355 const Bool getXPol, 355 const Bool getFeedPos)356 const Int coordSys) 356 357 { 357 358 if (!cMSopen) { … … 432 433 cGetXPol = cGetXPol && getXPol; 433 434 434 // Get feed positions? (Notavailable.)435 c GetFeedPos = False;435 // Coordinate system? (Only equatorial available.) 436 cCoordSys = 0; 436 437 437 438 return maxNChan; … … 486 487 // Read the next data record. 487 488 488 Int PKSMS2reader::read( MBrecord &MBrec)489 Int PKSMS2reader::read(PKSrecord &pksrec) 489 490 { 490 491 if (!cMSopen) { … … 514 515 515 516 // Renumerate scan no. Here still is 1-based 516 MBrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;517 518 if ( MBrec.scanNo != cScanNo) {517 pksrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1; 518 519 if (pksrec.scanNo != cScanNo) { 519 520 // Start of new scan. 520 cScanNo = MBrec.scanNo;521 cScanNo = pksrec.scanNo; 521 522 cCycleNo = 1; 522 523 cTime = cTimeCol(cIdx); … … 524 525 525 526 Double time = cTimeCol(cIdx); 526 MBrec.mjd = time/86400.0;527 MBrec.interval = cIntervalCol(cIdx);527 pksrec.mjd = time/86400.0; 528 pksrec.interval = cIntervalCol(cIdx); 528 529 529 530 // Reconstruct the integration cycle number; due to small latencies the 530 531 // integration time is usually slightly less than the time between cycles, 531 532 // resetting cTime will prevent the difference from accumulating. 532 cCycleNo += nint((time - cTime)/ MBrec.interval);533 MBrec.cycleNo = cCycleNo;534 cTime 533 cCycleNo += nint((time - cTime)/pksrec.interval); 534 pksrec.cycleNo = cCycleNo; 535 cTime = time; 535 536 536 537 Int fieldId = cFieldIdCol(cIdx); 537 MBrec.fieldName = cFieldNameCol(fieldId);538 pksrec.fieldName = cFieldNameCol(fieldId); 538 539 539 540 Int srcId = cSrcIdCol(fieldId); 540 MBrec.srcName = cSrcNameCol(srcId);541 MBrec.srcDir = cSrcDirCol(srcId);542 MBrec.srcPM = cSrcPMCol(srcId);541 pksrec.srcName = cSrcNameCol(srcId); 542 pksrec.srcDir = cSrcDirCol(srcId); 543 pksrec.srcPM = cSrcPMCol(srcId); 543 544 544 545 // Systemic velocity. 545 546 if (!cHaveSrcVel) { 546 MBrec.srcVel = 0.0f;547 pksrec.srcVel = 0.0f; 547 548 } else { 548 MBrec.srcVel= cSrcVelCol(srcId)(IPosition(1,0));549 pksrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0)); 549 550 } 550 551 551 552 // Observation type. 552 553 Int stateId = cStateIdCol(cIdx); 553 MBrec.obsType = cObsModeCol(stateId);554 555 MBrec.IFno = iIF + 1;554 pksrec.obsType = cObsModeCol(stateId); 555 556 pksrec.IFno = iIF + 1; 556 557 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1; 557 558 … … 560 561 if (nChan == 1) { 561 562 cout << "The input is continuum data. "<< endl; 562 MBrec.freqInc = chanFreq(0);563 MBrec.refFreq = chanFreq(0);564 MBrec.restFreq = 0.0f;563 pksrec.freqInc = chanFreq(0); 564 pksrec.refFreq = chanFreq(0); 565 pksrec.restFreq = 0.0f; 565 566 } else { 566 567 if (cStartChan(iIF) <= cEndChan(iIF)) { 567 MBrec.freqInc = chanFreq(1) - chanFreq(0);568 pksrec.freqInc = chanFreq(1) - chanFreq(0); 568 569 } else { 569 MBrec.freqInc = chanFreq(0) - chanFreq(1); 570 } 571 572 MBrec.refFreq = chanFreq(cRefChan(iIF)-1); 573 MBrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0)); 574 } 575 MBrec.bandwidth = abs(MBrec.freqInc * nChan); 576 577 MBrec.tcal.resize(cNPol(iIF)); 578 MBrec.tcal = 0.0f; 579 MBrec.tcalTime = ""; 580 MBrec.azimuth = 0.0f; 581 MBrec.elevation = 0.0f; 582 MBrec.parAngle = 0.0f; 583 MBrec.focusAxi = 0.0f; 584 MBrec.focusTan = 0.0f; 585 MBrec.focusRot = 0.0f; 570 pksrec.freqInc = chanFreq(0) - chanFreq(1); 571 } 572 573 pksrec.refFreq = chanFreq(cRefChan(iIF)-1); 574 pksrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0)); 575 } 576 pksrec.bandwidth = abs(pksrec.freqInc * nChan); 577 578 pksrec.tcal.resize(cNPol(iIF)); 579 pksrec.tcal = 0.0f; 580 pksrec.tcalTime = ""; 581 pksrec.azimuth = 0.0f; 582 pksrec.elevation = 0.0f; 583 pksrec.parAngle = 0.0f; 584 585 pksrec.focusAxi = 0.0f; 586 pksrec.focusTan = 0.0f; 587 pksrec.focusRot = 0.0f; 586 588 587 589 // Find the appropriate entry in the WEATHER subtable. … … 596 598 if (weatherIdx < 0) { 597 599 // No appropriate WEATHER entry. 598 MBrec.temperature = 0.0f;599 MBrec.pressure = 0.0f;600 MBrec.humidity = 0.0f;600 pksrec.temperature = 0.0f; 601 pksrec.pressure = 0.0f; 602 pksrec.humidity = 0.0f; 601 603 } else { 602 MBrec.temperature = cTemperatureCol(weatherIdx);603 MBrec.pressure = cPressureCol(weatherIdx);604 MBrec.humidity = cHumidityCol(weatherIdx);605 } 606 607 MBrec.windSpeed = 0.0f;608 MBrec.windAz = 0.0f;609 610 MBrec.refBeam = 0;611 MBrec.beamNo = ibeam + 1;604 pksrec.temperature = cTemperatureCol(weatherIdx); 605 pksrec.pressure = cPressureCol(weatherIdx); 606 pksrec.humidity = cHumidityCol(weatherIdx); 607 } 608 609 pksrec.windSpeed = 0.0f; 610 pksrec.windAz = 0.0f; 611 612 pksrec.refBeam = 0; 613 pksrec.beamNo = ibeam + 1; 612 614 613 615 Matrix<Double> pointingDir = cPointingCol(fieldId); 614 MBrec.direction = pointingDir.column(0); 616 pksrec.direction = pointingDir.column(0); 617 pksrec.pCode = 0; 618 pksrec.rateAge = 0.0f; 615 619 uInt ncols = pointingDir.ncolumn(); 616 620 if (ncols == 1) { 617 MBrec.scanRate = 0.0f;621 pksrec.scanRate = 0.0f; 618 622 } else { 619 MBrec.scanRate = pointingDir.column(1);620 }621 MBrec.rateAge = 0;622 MBrec.rateson = 0;623 pksrec.scanRate(0) = pointingDir.column(1)(0); 624 pksrec.scanRate(1) = pointingDir.column(1)(1); 625 } 626 pksrec.paRate = 0.0f; 623 627 624 628 // Get Tsys assuming that entries in the SYSCAL table match the main table. … … 630 634 } 631 635 if (cHaveTsys) { 632 cTsysCol.get(cIdx, MBrec.tsys, True);636 cTsysCol.get(cIdx, pksrec.tsys, True); 633 637 } else { 634 638 Int numReceptor; 635 639 cNumReceptorCol.get(0, numReceptor); 636 MBrec.tsys.resize(numReceptor);637 MBrec.tsys = 1.0f;638 } 639 cSigmaCol.get(cIdx, MBrec.sigma, True);640 pksrec.tsys.resize(numReceptor); 641 pksrec.tsys = 1.0f; 642 } 643 cSigmaCol.get(cIdx, pksrec.sigma, True); 640 644 641 645 // Calibration factors (if available). 642 MBrec.calFctr.resize(cNPol(iIF));646 pksrec.calFctr.resize(cNPol(iIF)); 643 647 if (cHaveCalFctr) { 644 cCalFctrCol.get(cIdx, MBrec.calFctr);648 cCalFctrCol.get(cIdx, pksrec.calFctr); 645 649 } else { 646 MBrec.calFctr = 0.0f;650 pksrec.calFctr = 0.0f; 647 651 } 648 652 649 653 // Baseline parameters (if available). 650 654 if (cHaveBaseLin) { 651 MBrec.baseLin.resize(2,cNPol(iIF));652 cBaseLinCol.get(cIdx, MBrec.baseLin);653 654 MBrec.baseSub.resize(9,cNPol(iIF));655 cBaseSubCol.get(cIdx, MBrec.baseSub);655 pksrec.baseLin.resize(2,cNPol(iIF)); 656 cBaseLinCol.get(cIdx, pksrec.baseLin); 657 658 pksrec.baseSub.resize(9,cNPol(iIF)); 659 cBaseSubCol.get(cIdx, pksrec.baseSub); 656 660 657 661 } else { 658 MBrec.baseLin.resize(0,0);659 MBrec.baseSub.resize(0,0);662 pksrec.baseLin.resize(0,0); 663 pksrec.baseSub.resize(0,0); 660 664 } 661 665 … … 670 674 // Transpose spectra. 671 675 Int nPol = tmpData.nrow(); 672 MBrec.spectra.resize(nChan, nPol);673 MBrec.flagged.resize(nChan, nPol);676 pksrec.spectra.resize(nChan, nPol); 677 pksrec.flagged.resize(nChan, nPol); 674 678 if (cEndChan(iIF) >= cStartChan(iIF)) { 675 679 // Simple transposition. 676 680 for (Int ipol = 0; ipol < nPol; ipol++) { 677 681 for (Int ichan = 0; ichan < nChan; ichan++) { 678 MBrec.spectra(ichan,ipol) = tmpData(ipol,ichan);679 MBrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);682 pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan); 683 pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan); 680 684 } 681 685 } … … 686 690 for (Int ipol = 0; ipol < nPol; ipol++) { 687 691 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) { 688 MBrec.spectra(ichan,ipol) = tmpData(ipol,jchan);689 MBrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);692 pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan); 693 pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan); 690 694 } 691 695 } … … 696 700 if (cGetXPol) { 697 701 if (cHaveXCalFctr) { 698 cXCalFctrCol.get(cIdx, MBrec.xCalFctr);702 cXCalFctrCol.get(cIdx, pksrec.xCalFctr); 699 703 } else { 700 MBrec.xCalFctr = Complex(0.0f, 0.0f);701 } 702 703 cDataCol.get(cIdx, MBrec.xPol, True);704 pksrec.xCalFctr = Complex(0.0f, 0.0f); 705 } 706 707 cDataCol.get(cIdx, pksrec.xPol, True); 704 708 705 709 if (cEndChan(iIF) < cStartChan(iIF)) { … … 707 711 Int jchan = nChan - 1; 708 712 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) { 709 ctmp = MBrec.xPol(ichan); 710 MBrec.xPol(ichan) = MBrec.xPol(jchan); 711 MBrec.xPol(jchan) = ctmp; 712 } 713 } 714 } 715 716 cIdx++; 717 718 return 0; 719 } 720 721 //--------------------------------------------------------- PKSMS2reader::read 722 723 // Read the next data record, just the basics. 724 725 Int PKSMS2reader::read( 726 Int &IFno, 727 Vector<Float> &tsys, 728 Vector<Float> &calFctr, 729 Matrix<Float> &baseLin, 730 Matrix<Float> &baseSub, 731 Matrix<Float> &spectra, 732 Matrix<uChar> &flagged) 733 { 734 if (!cMSopen) { 735 return 1; 736 } 737 738 // Check for EOF. 739 if (cIdx >= cNRow) { 740 return -1; 741 } 742 743 // Find the next selected beam and IF. 744 Int ibeam; 745 Int iIF; 746 while (True) { 747 ibeam = cBeamNoCol(cIdx); 748 iIF = cDataDescIdCol(cIdx); 749 if (cBeams(ibeam) && cIFs(iIF)) { 750 break; 751 } 752 753 // Check for EOF. 754 if (++cIdx >= cNRow) { 755 return -1; 756 } 757 } 758 759 IFno = iIF + 1; 760 761 // Get Tsys assuming that entries in the SYSCAL table match the main table. 762 cTsysCol.get(cIdx, tsys, True); 763 764 // Calibration factors (if available). 765 if (cHaveCalFctr) { 766 cCalFctrCol.get(cIdx, calFctr, True); 767 } else { 768 calFctr.resize(cNPol(iIF)); 769 calFctr = 0.0f; 770 } 771 772 // Baseline parameters (if available). 773 if (cHaveBaseLin) { 774 baseLin.resize(2,cNPol(iIF)); 775 cBaseLinCol.get(cIdx, baseLin); 776 777 baseSub.resize(9,cNPol(iIF)); 778 cBaseSubCol.get(cIdx, baseSub); 779 780 } else { 781 baseLin.resize(0,0); 782 baseSub.resize(0,0); 783 } 784 785 if (cGetSpectra) { 786 // Get spectral data. 787 Matrix<Float> tmpData; 788 Matrix<Bool> tmpFlag; 789 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True); 790 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True); 791 792 // Transpose spectra. 793 Int nChan = tmpData.ncolumn(); 794 Int nPol = tmpData.nrow(); 795 spectra.resize(nChan, nPol); 796 flagged.resize(nChan, nPol); 797 if (cEndChan(iIF) >= cStartChan(iIF)) { 798 // Simple transposition. 799 for (Int ipol = 0; ipol < nPol; ipol++) { 800 for (Int ichan = 0; ichan < nChan; ichan++) { 801 spectra(ichan,ipol) = tmpData(ipol,ichan); 802 flagged(ichan,ipol) = tmpFlag(ipol,ichan); 803 } 804 } 805 806 } else { 807 // Transpose with inversion. 808 Int jchan = nChan - 1; 809 for (Int ipol = 0; ipol < nPol; ipol++) { 810 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) { 811 spectra(ichan,ipol) = tmpData(ipol,jchan); 812 flagged(ichan,ipol) = tmpFlag(ipol,jchan); 813 } 713 ctmp = pksrec.xPol(ichan); 714 pksrec.xPol(ichan) = pksrec.xPol(jchan); 715 pksrec.xPol(jchan) = ctmp; 814 716 } 815 717 }
Note: See TracChangeset
for help on using the changeset viewer.