LCOV - code coverage report
Current view: top level - components/ComponentModels - FluxCalcVQS.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 104 120 86.7 %
Date: 2024-11-06 17:42:47 Functions: 11 13 84.6 %

          Line data    Source code
       1             : //# FluxCalcVQS.cc: Implementation of FluxCalcVQS.h
       2             : //# Copyright (C) 2013
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : #include <components/ComponentModels/FluxCalcVQS.h>
      27             : #include <casacore/casa/Arrays/Vector.h>
      28             : #include <casacore/casa/BasicSL/String.h>
      29             : #include <casacore/casa/Logging/LogIO.h>
      30             : #include <casacore/casa/Quanta/MVTime.h>
      31             : #include <casacore/measures/Measures/MDirection.h>
      32             : #include <casacore/measures/Measures/MFrequency.h>
      33             : #include <casacore/measures/Measures/MEpoch.h>
      34             : #include <casacore/tables/Tables/ScalarColumn.h>
      35             : #include <casacore/tables/Tables/ArrayColumn.h>
      36             : #include <casacore/tables/Tables/TableRecord.h>
      37             : #include <casacore/tables/Tables/TableRecord.h>
      38             : #include <casacore/tables/Tables/TableProxy.h>
      39             : 
      40             : // Handy for passing anonymous arrays to functions.
      41             : #include <casacore/scimath/Mathematics/RigidVector.h>
      42             : #include <casacore/scimath/Functionals/ScalarSampledFunctional.h>
      43             : 
      44             : #include <components/ComponentModels/FluxCalcLogFreqPolynomial.h>
      45             : 
      46             : #include <map>
      47             : 
      48             : 
      49             : using namespace casacore;
      50             : namespace casa { //# NAMESPACE CASA - BEGIN
      51         103 : FluxCalcVQS::FluxCalcVQS() :
      52         103 :   srcEnum_p(FluxStdSrcs::UNKNOWN_SOURCE),
      53         103 :   validfreqrange_p(2),
      54         206 :   istimevar_p(false)
      55         103 : { }
      56             : 
      57             : // Defined even though it's pure virtual; see http://www.gotw.ca/gotw/031.htm
      58         103 : FluxCalcVQS::~FluxCalcVQS() { }
      59          10 : Bool FluxCalcVQS::operator()(Vector<Flux<Double> >& values,
      60             :                             Vector<Flux<Double> >& errors,
      61             :                             const Vector<MFrequency>& mfreqs)
      62             : {
      63          10 :   uInt nfreqs = mfreqs.nelements();
      64             :   
      65          10 :   values.resize(nfreqs);
      66          10 :   errors.resize(nfreqs);
      67             :   
      68          10 :   if (coeffsmat_p.nelements()) {
      69             :       // coeffs have been read from a table.
      70           0 :       uInt nep=0; // should be a single epoch (single row)
      71           0 :       setSourceCoeffsfromVec(nep);    
      72             :   }
      73             : 
      74          10 :   Bool success = true;
      75          20 :   for(uInt f = 0; f < nfreqs; ++f)
      76          10 :     success &= (*this)(values[f], errors[f], mfreqs[f], false);
      77             :   
      78          10 :   return success;
      79             : }
      80             : 
      81             : //with interpolation over epochs
      82          93 : Bool FluxCalcVQS::operator()(Vector<Flux<Double> >& values,
      83             :                             Vector<Flux<Double> >& errors,
      84             :                             const Vector<MFrequency>& mfreqs,
      85             :                             const MEpoch& mtime,
      86             :                             const String& interpmethod)
      87             : 
      88             : {
      89          93 :   uInt nfreqs = mfreqs.nelements();
      90             : 
      91          93 :   values.resize(nfreqs);
      92          93 :   errors.resize(nfreqs);
      93             : 
      94             :   // for those considered to be variable, for each set of coefficients
      95             :   // at each epoch, call the follwoing to get values (fluxes)
      96             :   // and accumate in vector of vector to do interpolation.
      97             :   // istimevar_p to determine if the source is variable. If not
      98             :   // no interpolation is done, use first row of the coeffs table. 
      99             :   
     100          93 :   Bool success = true;
     101             : 
     102          93 :   Int nep = 1;
     103          93 :   if (istimevar_p) nep=epochvec_p.nelements();
     104             : 
     105          93 :   Flux<Double> tmpfluxes;
     106          93 :   Flux<Double> tmperrors;
     107          93 :   Vector<Double> tmpfluxvec;
     108          93 :   Vector<Double> tmperrvec;
     109             :   
     110          93 :   fluxes_p.resize(nep);
     111        2734 :   for(uInt f = 0; f < nfreqs; ++f){
     112       14735 :     for(uInt iep=0; iep < (uInt)nep; ++iep) {
     113       12094 :       setSourceCoeffsfromVec(iep);    
     114       12094 :       success &= (*this)(tmpfluxes,tmperrors,mfreqs[f],true);
     115       12094 :       tmpfluxes.value(tmpfluxvec);
     116       12094 :       tmperrors.value(tmperrvec);
     117             :       // currently only I flux is returned...
     118             :       //cerr<<"tmpfluxvec[0]="<<tmpfluxvec[0]<<endl;
     119             :       //cerr<<"epochvec_p[iep]="<<epochvec_p[iep]<<endl;
     120       12094 :       fluxes_p[iep]=tmpfluxvec[0];
     121             :     }   
     122        2641 :     if (istimevar_p) {
     123             :       //setup interpolator 
     124         571 :       ScalarSampledFunctional<Double> dts(epochvec_p);
     125         571 :       ScalarSampledFunctional<Float> flxs(fluxes_p);
     126         571 :       Interpolate1D<Double, Float> interpolateFlux(dts,flxs);
     127         571 :       interpolateFlux.setMethod(getInterpMethod_p(interpmethod));
     128             : 
     129         571 :       values[f].setValue(interpolateFlux(mtime.get("d").getValue()));
     130         571 :     }
     131             :     else { // no interpolation for non-var source, use first row data
     132        2070 :       values[f].setValue(fluxes_p[0]);  
     133             :     }
     134             :   }
     135             : 
     136             : //      success &= (*this)(values[f], errors[f], mfreqs[f]);
     137             :     
     138          93 :   return success;
     139          93 : }
     140             : 
     141         103 : Bool FluxCalcVQS::setSource(const String& sourceName, const MDirection& sourceDir)
     142             : {
     143         103 :   srcEnum_p = srcNameToEnum(sourceName,sourceDir);
     144             :   //return srcEnum_p != FCQS::UNKNOWN_SOURCE;
     145         103 :   return srcEnum_p != FSS::UNKNOWN_SOURCE;
     146             : }
     147             : 
     148         103 : FluxCalcVQS::Source FluxCalcVQS::getSrcEnum()
     149             : {
     150         103 :   return srcEnum_p;
     151             : }
     152             : 
     153          93 : void FluxCalcVQS::readQSCoeffsTable(const Path& fileName)
     154             : {
     155             :   //table containing the coefficents has
     156             :   //epoch, vector of coefficients
     157          93 :   const String& fullName = fileName.absoluteName();
     158         186 :   LogIO os(LogOrigin("FluxCalcVQS", "readQSCoeffsTable", WHERE));
     159             :    os << LogIO::NORMAL1
     160             :       << "Reading the coefficient data from a table, "<<  fullName
     161          93 :       << LogIO::POST;
     162             :   
     163          93 :   AlwaysAssert(Table::isReadable(fullName), AipsError);
     164          93 :   Table_p = Table(fullName, Table::Old);
     165             :   ///TableProxy tab2(Table_p);
     166             :   //String srcName(names_p[srcEnum_p](0));
     167          93 :   String srcName(EnumToSrcName(srcEnum_p));
     168          93 :   String srcCoeffColName=srcName+"_coeffs";
     169          93 :   String srcCoeffErrorColName=srcName+"_coefferrs";
     170             : 
     171             :   //check the source data exist in the table
     172          93 :   const ColumnDescSet& cds=Table_p.tableDesc().columnDescSet();
     173          93 :   if (!cds.isDefined(srcCoeffColName)) 
     174           0 :     throw(AipsError(srcName+" does not appears to be in "+fullName));
     175          93 :   const ScalarColumn<Double> epochCol(Table_p, "Epoch");
     176          93 :   const ArrayColumn<Float> CoeffCol(Table_p, srcCoeffColName);
     177          93 :   const ArrayColumn<Float> CoeffErrorCol(Table_p, srcCoeffErrorColName);
     178             :   // check if col contains a valid freq range info in a keyword
     179          93 :   if (CoeffCol.keywordSet().isDefined("ValidFreqRange")) {
     180          87 :     Vector<Double> validfreqRange;
     181          87 :     String freqRangeUnit;
     182          87 :     CoeffCol.keywordSet().asRecord("ValidFreqRange").get("freqRange",validfreqRange);
     183          87 :     CoeffCol.keywordSet().asRecord("ValidFreqRange").get("freqRangeUnit",freqRangeUnit);
     184          87 :     Vector<MFrequency> mvalidfreqs;
     185          87 :     validfreqrange_p(0) = MFrequency(Quantity(validfreqRange(0),freqRangeUnit), MFrequency::TOPO);
     186          87 :     validfreqrange_p(1) = MFrequency(Quantity(validfreqRange(1),freqRangeUnit), MFrequency::TOPO);
     187          87 :   } 
     188             :   else {
     189           6 :     validfreqrange_p(0) = MFrequency(Quantity(0.0,"GHz"), MFrequency::TOPO);
     190           6 :     validfreqrange_p(1) = MFrequency(Quantity(0.0,"GHz"), MFrequency::TOPO);
     191             :   }
     192          93 :   Vector<Double> tempEpochs;
     193          93 :   epochCol.getColumn(tempEpochs,true);
     194          93 :   CoeffCol.getColumn(coeffsmat_p,true);
     195          93 :   CoeffErrorCol.getColumn(coefferrsmat_p,true);
     196             :   //convert the epoch (year + fraction) to mjd
     197          93 :   convertYearFracToMjd(tempEpochs,epochvec_p);
     198             :   os << LogIO::DEBUG1
     199             :      << "nepoch="<<epochvec_p.nelements()
     200         186 :      << "coeff_0 for the first epoch (coeffsmat_p(0,0))="<<coeffsmat_p(0,0) 
     201          93 :      << LogIO::POST;
     202             : 
     203          93 : }
     204             : 
     205           0 : void FluxCalcVQS::interpolate(const String& interpmethod)
     206             : {
     207           0 :   ScalarSampledFunctional<Double> dts(epochvec_p);
     208           0 :   ScalarSampledFunctional<Float> flxs(fluxes_p);
     209           0 :   Interpolate1D<Double, Float> interpolateFlux(dts,flxs);
     210           0 :   interpolateFlux.setMethod(getInterpMethod_p(interpmethod));
     211           0 : }
     212             : 
     213         571 : Interpolate1D<Double,Float>::Method FluxCalcVQS::getInterpMethod_p(const String& interpmethod)
     214             : {
     215         571 :   if (interpmethod.contains("nearest"))
     216         571 :     return Interpolate1D<Double,Float>::nearestNeighbour;
     217           0 :   if (interpmethod.contains("linear"))
     218           0 :     return Interpolate1D<Double,Float>::linear;
     219           0 :   if (interpmethod.contains("cubic")) 
     220           0 :     return Interpolate1D<Double,Float>::cubic;
     221           0 :   if (interpmethod.contains("spline"))
     222           0 :     return Interpolate1D<Double,Float>::spline;
     223           0 :   throw(AipsError("Unknown interpolation method: "+interpmethod));
     224             : }
     225             : 
     226          93 : void FluxCalcVQS::convertYearFracToMjd(const Vector<Double>& yearfrac, Vector<Double>& mjds)
     227             : {
     228          93 :   uInt daysintheyr=365;
     229          93 :   uInt n = yearfrac.nelements();
     230          93 :   mjds.resize(n);
     231        1761 :   for (uInt i=0; i < n; i++) {
     232        1668 :     if (Time::isLeapYear(uInt(yearfrac(i)))) daysintheyr=366;
     233             :     Double fintyr, fyrfrac;
     234        1668 :     fyrfrac = modf((Double)yearfrac(i),&fintyr);
     235        1668 :     Float days=fyrfrac*daysintheyr;
     236        1668 :     Time time0(Int(yearfrac(i)),1,1);
     237        1668 :     Double mjdtime0 = time0.modifiedJulianDay();
     238             : 
     239        1668 :     mjds[i] = mjdtime0 + days;
     240             :   }
     241          93 : }
     242             :     
     243       12094 : void FluxCalcVQS::setSourceCoeffsfromVec(uInt& i) 
     244             : {
     245             :   //Vector<Float> err(4,0.0);
     246       12094 :   tvcoeffs_p(0)=coeffsmat_p.column(i); 
     247             :   //cerr<<"coeffsmat_p.column="<<coeffsmat_p.column(i)(0)<<endl;
     248       12094 :   tvcoeffs_p(1)=coefferrsmat_p.column(i); 
     249       12094 : }
     250             : 
     251          93 : void FluxCalcVQS::isTimeVar(Bool istimevar) 
     252             : {
     253          93 :   istimevar_p=istimevar;
     254          93 : }
     255             : 
     256             : } //# NAMESPACE CASA - END

Generated by: LCOV version 1.16