#include #include #include #include #include #include #include #include #include #include #include #include #include "ASDMFiller.h" using namespace std ; using namespace casa ; using namespace asap ; ASDMFiller::ASDMFiller( CountedPtr stable ) : FillerBase( stable ), antennaId_( -1 ), antennaName_( "" ), className_("ASDMFiller") { reader_ = new ASDMReader() ; } ASDMFiller::~ASDMFiller() { // nothing to do? logsink_ = 0 ; } void ASDMFiller::setLogger( CountedPtr &logsink ) { logsink_ = logsink ; if ( !(reader_.null()) ) { reader_->setLogger( logsink ) ; } } bool ASDMFiller::open( const string &filename, const Record &rec ) { String funcName = "open" ; bool status = reader_->open( filename, rec ) ; antennaId_ = reader_->getAntennaId() ; antennaName_ = reader_->getAntennaName() ; //logsink->postLocally( LogMessage("antennaId_ = "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ; //logsink->postLocally( LogMessage("antennaName_ = "+antennaName_,LogOrigin(className_,funcName,WHERE)) ) ; return status ; } void ASDMFiller::fill() { String funcName = "fill" ; // header fillHeader() ; Vector antpos = table_->getHeader().antennaposition ; //STHeader hdr = table_->getHeader() ; // data selection reader_->select() ; // pick up valid configDescriptionId Vector configDescIdList = reader_->getConfigDescriptionIdList() ; uInt numConfigDescId = configDescIdList.size() ; //logsink->postLocally( LogMessage("configDescIdList = "+String::toString(configDescIdList),LogOrigin(className_,funcName,WHERE)) ) ; // get field list Vector fieldIdList = reader_->getFieldIdList() ; uInt numFieldId = fieldIdList.size() ; //logsink->postLocally( LogMessage("fieldIdList = "+String::toString(fieldIdList),LogOrigin(className_,funcName,WHERE)) ) ; // BEAMNO is always 0 since ALMA antenna is single beam uInt beamno = 0 ; // REFBEAMNO is -1 setReferenceBeam() ; // fill FOCUS_ID and add FOCUS row if necessary setFocus() ; for ( uInt icon = 0 ; icon < numConfigDescId ; icon++ ) { for ( unsigned int ifield = 0 ; ifield < numFieldId ; ifield++ ) { //logsink_->postLocally( LogMessage("start configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ; //Bool status = reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] ) ; if ( !(reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] )) ) { //logsink_->postLocally( LogMessage("skip configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ; continue ; } // number of rows uInt nrow = reader_->getNumMainRow() ; //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows in Main table corresponding to configDescId "+String::toString(configDescIdList[icon]+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ; // CYCLENO unsigned int cycleno = 0 ; for ( uInt irow = 0 ; irow < nrow ; irow++ ) { // set main row if ( !(reader_->setMainRow( irow )) ) { // skip row since the row doesn't have valid configDescId //logsink_->postLocally( LogMessage("skip "+String::toString(irow),LogOrigin(className_,funcName,WHERE)) ) ; continue ; } // scan and subscan unsigned int scanno = reader_->getScanNo() ; //uInt subscanno = reader_->getSubscanNo() ; // set data if ( !(reader_->setData()) ) { // skip row since reader failed to retrieve data //logsink_->postLocally( LogMessage("skip "+String::toString(irow),LogOrigin(className_,funcName,WHRER)) ) ; continue ; } unsigned int numData = reader_->getNumData() ; double refpix = 0.0 ; double refval = 0.0 ; double incr = 0.0 ; for ( unsigned int idata = 0 ; idata < numData ; idata++ ) { // subscan number unsigned int subscanno = reader_->getSubscanNo( idata ) ; // IFNO uInt ifno = reader_->getIFNo( idata ) ; // REFPIX, REFVAL, INCREMENT String ifkey = getIFKey( ifno ) ; if ( ifrec_.isDefined( ifkey ) ) { getFrequencyRec( ifkey, refpix, refval, incr ) ; } else { reader_->getFrequency( idata, refpix, refval, incr ) ; setFrequencyRec( ifkey, refpix, refval, incr ) ; } // fill FREQ_ID and add FREQUENCIES row if necessary setFrequency( (casa::Double)refpix, (casa::Double)refval, (casa::Double)incr ) ; // rest frequency vector rf = reader_->getRestFrequency( idata ) ; // fill MOLECULE_ID and add MOLECULES row if necessary Vector restFreqs( rf.size() ) ; for ( uInt i = 0 ; i < rf.size() ; i++ ) restFreqs[i] = (casa::Double)(rf[i]) ; setMolecule( restFreqs ) ; // time and interval casa::Double mjd = (casa::Double)(reader_->getTime( idata )) ; casa::Double interval = (casa::Double)(reader_->getInterval( idata )) ; // fill TIME and INTERVAL setTime( mjd, interval ) ; // source spec string srcname = reader_->getSourceName( idata ) ; string fieldname = reader_->getFieldName( idata ) ; int srctype = reader_->getSrcType( scanno, subscanno ) ; vector srcDirection = reader_->getSourceDirection( idata ) ; vector srcProperMotion = reader_->getSourceProperMotion( idata ) ; double sysVel = reader_->getSysVel( idata ) ; // fill SRCNAME, SRCTYPE, FIELDNAME, SRCDIRECTION, SRCPROPERMOTION, and SRCVELOCITY Vector srcDir( 2 ) ; srcDir[0] = (casa::Double)(srcDirection[0]) ; srcDir[1] = (casa::Double)(srcDirection[1]) ; Vector srcPM( 2 ) ; srcPM[0] = (casa::Double)(srcProperMotion[0]) ; srcPM[1] = (casa::Double)(srcProperMotion[1]) ; setSource( srcname, srctype, fieldname, srcDir, srcPM, (casa::Double)sysVel ) ; // fill FLAGROW unsigned int flagrow = reader_->getFlagRow( idata ) ; setFlagrow( (uInt)flagrow ) ; // fill WEATHER_ID and add WEATHER row if necessary float temperature ; float pressure ; float humidity ; float windspeed ; float windaz ; reader_->getWeatherInfo( idata, temperature, pressure, humidity, windspeed, windaz ) ; setWeather2( (casa::Float)temperature, (casa::Float)pressure, (casa::Float)humidity, (casa::Float)windspeed, (casa::Float)windaz ) ; // fill AZIMUTH, ELEVATION, DIRECTION and SCANRATE vector dir ; double az ; double el ; vector srate ; reader_->getPointingInfo( idata, dir, az, el, srate ) ; Vector scanRate( 2, 0.0 ) ; Vector direction( 2, 0.0 ) ; if ( srate.size() > 0 ) { scanRate[0] = (casa::Double)(srate[0]) ; scanRate[1] = (casa::Double)(srate[1]) ; } setScanRate( scanRate ) ; if ( dir.size() > 0 ) { direction[0] = (casa::Double)(dir[0]) ; direction[1] = (casa::Double)(dir[1]) ; } else { toJ2000( direction, az, el, mjd, antpos ) ; } //logsink_->postLocally( LogMessage("direction = "+String::toString(direction),LogOrigin(className_,funcName,WHERE)) ) ; setDirection( direction, (casa::Float)az, (casa::Float)el ) ; // loop on polarization vector dataShape = reader_->getDataShape( idata ) ; // ostringstream oss ; // for ( unsigned int i = 0 ; i < dataShape.size() ; i++ ) { // if ( i == 0 ) // oss << "dataShape=[" << dataShape[i] << ", " ; // else if ( i == dataShape.size()-1 ) // oss << dataShape[i] << "]" ; // else // oss << dataShape[i] << ", " ; // } // logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ; //int numPol = reader_->getNumPol( idata ) ; unsigned int numPol = dataShape[0] ; unsigned int numChan = dataShape[1] ; //logsink_->postLocally( LogMessage("numPol = "+String::toString(numPol),LogOrigin(className_,funcName,WHERE)) ) ; // OPACITY vector tau = reader_->getOpacity( idata ) ; Vector opacity = toVector( tau, numPol ) ; // SPECTRA, FLAGTRA, TSYS, TCAL float *sp = reader_->getSpectrum( idata ) ; vector< vector > ts = reader_->getTsys( idata ) ; vector< vector > tc = reader_->getTcal( idata ) ; Matrix spectra = toMatrix( sp, numPol, numChan ) ; Vector flagtra( numChan, 0 ) ; Matrix tsys = toMatrix( ts, numPol, numChan ) ; Matrix tcal = toMatrix( tc, numPol, numChan ) ; // String caltime = "" ; // if ( anyNE( tcal, (casa::Float)1.0 ) ) // caltime = toTcalTime( mjd ) ; String caltime = toTcalTime( mjd ) ; for ( unsigned int ipol = 0 ; ipol < numPol ; ipol++ ) { // fill SCANNO, CYCLENO, IFNO, POLNO, and BEAMNO setIndex( (uInt)scanno-1, (uInt)cycleno, ifno, ipol, beamno ) ; // fill SPECTRA, FLAGTRA, TSYS setSpectrum( spectra.row(ipol), flagtra, tsys.row(ipol) ) ; // fill TCAL_ID and add TCAL row if necessary setTcal2( caltime, tcal.row(ipol) ) ; // fill OPACITY setOpacity( opacity[ipol] ) ; // commit row commitRow() ; } // increment CYCLENO cycleno++ ; } } } } return ; } void ASDMFiller::close() { reader_->close() ; reader_ = 0 ; return ; } void ASDMFiller::fillHeader() { STHeader hdr ; reader_->fillHeader( hdr.nchan, hdr.npol, hdr.nif, hdr.nbeam, hdr.observer, hdr.project, hdr.obstype, hdr.antennaname, hdr.antennaposition, hdr.equinox, hdr.freqref, hdr.reffreq, hdr.bandwidth, hdr.utc, hdr.fluxunit, hdr.epoch, hdr.poltype ) ; setHeader( hdr ) ; } String ASDMFiller::getIFKey( uInt ifno ) { return "IFNO"+String::toString( ifno ) ; } void ASDMFiller::getFrequencyRec( String key, double &refpix, double &refval, double &incr ) { Record frec = ifrec_.asRecord( key ) ; refpix = frec.asdouble( "REFPIX" ) ; refval = frec.asdouble( "REFVAL" ) ; incr = frec.asdouble( "INCREMENT" ) ; } void ASDMFiller::setFrequencyRec( String key, double refpix, double refval, double incr ) { Record frec ; frec.define( "REFPIX", refpix ) ; frec.define( "REFVAL", refval ) ; frec.define( "INCREMENT", incr ) ; ifrec_.defineRecord( key, frec ) ; } Matrix ASDMFiller::toMatrix( float *sp, unsigned int npol, unsigned int nchan ) { Matrix mSp( npol, nchan ) ; if ( npol <= 2 ) { // 1 or 2 polarization case for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) { for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) { mSp(ipol,ich) = (casa::Float)(sp[npol*ich+ipol]) ; } } } else { // 4 polarization case for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) { mSp(0,ich) = (casa::Float)(sp[4*ich]) ; // Re(XX) mSp(1,ich) = (casa::Float)(sp[4*ich+4]) ; // Re(YY) mSp(2,ich) = (casa::Float)(sp[4*ich+2]) ; // Re(XY) mSp(3,ich) = (casa::Float)(sp[4*ich+3]) ; // Im(XY) } } return mSp ; } Matrix ASDMFiller::toMatrix( vector< vector > &tsys, unsigned int npol, unsigned int nchan ) { unsigned int numRec = tsys.size() ; unsigned int numChan = tsys[0].size() ; Matrix ret ; if ( npol == numRec && nchan == numChan ) { ret.resize( npol, nchan ) ; for ( unsigned int ip = 0 ; ip < npol ; ip++ ) for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) ret( ip, ic ) = (casa::Float)(tsys[ip][ic]) ; } else if ( npol == numRec && numChan == 1 ) { ret.resize( npol, 1 ) ; for ( unsigned int ip = 0 ; ip < npol ; ip++ ) ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ; } else if ( numRec == 1 && nchan == numChan ) { ret.resize( npol, nchan ) ; for ( unsigned int ip = 0 ; ip < npol ; ip++ ) for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ; } else if ( numRec == 1 && numChan == 1 ) { ret.resize( npol, 1 ) ; for ( unsigned int ip = 0 ; ip < npol ; ip++ ) ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ; } else if ( numRec == 2 && npol == 4 && numChan == nchan ) { // TODO: How to determine Tsys for XY? // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y]) ret.resize( npol, nchan ) ; for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) { casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][ic]+tsys[1][ic])) ; ret( 0, ic ) = (casa::Float)(tsys[0][ic]) ; ret( 1, ic ) = (casa::Float)(tsys[1][ic]) ; ret( 2, ic ) = tsysxy ; ret( 3, ic ) = tsysxy ; } } else if ( numRec == 2 && npol == 4 && numChan == 1 ) { // TODO: How to determine Tsys for XY? // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y]) ret.resize( npol, 1 ) ; casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][0]+tsys[1][0])) ; ret( 0, 0 ) = (casa::Float)(tsys[0][0]) ; ret( 1, 0 ) = (casa::Float)(tsys[1][0]) ; ret( 2, 0 ) = tsysxy ; ret( 3, 0 ) = tsysxy ; } else { // I don't know how to handle ... for ( unsigned int ip = 0 ; ip < npol ; ip++ ) for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ; } return ret ; } Vector ASDMFiller::toVector( vector &tau, unsigned int npol ) { Vector ret( npol ) ; if ( npol == 4 ) { ret[0] = (casa::Float)tau[0] ; ret[1] = (casa::Float)tau[1] ; ret[2] = 0.5 * ( ret[0] + ret[1] ) ; ret[3] = ret[2] ; } else if ( npol == tau.size() ) { for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) ret[ipol] = (casa::Float)tau[ipol] ; } else { // I don't know how to handle... for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) ret[ipol] = (casa::Float)tau[0] ; } return ret ; } String ASDMFiller::toTcalTime( casa::Double mjd ) { return MVTime( mjd ).string( MVTime::YMD ) ; } void ASDMFiller::toJ2000( Vector &dir, double az, double el, casa::Double mjd, Vector antpos ) { String funcName = "toJ2000" ; Vector azel( 2 ) ; azel[0] = az ; azel[1] = el ; MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ; Vector qantpos( 3 ) ; qantpos[0] = Quantity( antpos[0], "m" ) ; qantpos[1] = Quantity( antpos[1], "m" ) ; qantpos[2] = Quantity( antpos[2], "m" ) ; MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ; // mp.print( os_.output() ) ; MeasFrame mf( me, mp ) ; MDirection::Convert toj2000( MDirection::AZELGEO, MDirection::Ref( MDirection::J2000, mf ) ) ; dir = toj2000( azel ).getAngle( "rad" ).getValue() ; //logsink->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ; }