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 37 : FluxCalcVQS::FluxCalcVQS() :
52 37 : srcEnum_p(FluxStdSrcs::UNKNOWN_SOURCE),
53 37 : validfreqrange_p(2),
54 74 : istimevar_p(false)
55 37 : { }
56 :
57 : // Defined even though it's pure virtual; see http://www.gotw.ca/gotw/031.htm
58 37 : FluxCalcVQS::~FluxCalcVQS() { }
59 5 : Bool FluxCalcVQS::operator()(Vector<Flux<Double> >& values,
60 : Vector<Flux<Double> >& errors,
61 : const Vector<MFrequency>& mfreqs)
62 : {
63 5 : uInt nfreqs = mfreqs.nelements();
64 :
65 5 : values.resize(nfreqs);
66 5 : errors.resize(nfreqs);
67 :
68 5 : 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 5 : Bool success = true;
75 10 : for(uInt f = 0; f < nfreqs; ++f)
76 5 : success &= (*this)(values[f], errors[f], mfreqs[f], false);
77 :
78 5 : return success;
79 : }
80 :
81 : //with interpolation over epochs
82 32 : 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 32 : uInt nfreqs = mfreqs.nelements();
90 :
91 32 : values.resize(nfreqs);
92 32 : 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 32 : Bool success = true;
101 :
102 32 : Int nep = 1;
103 32 : if (istimevar_p) nep=epochvec_p.nelements();
104 :
105 32 : Flux<Double> tmpfluxes;
106 32 : Flux<Double> tmperrors;
107 32 : Vector<Double> tmpfluxvec;
108 32 : Vector<Double> tmperrvec;
109 :
110 32 : fluxes_p.resize(nep);
111 2080 : for(uInt f = 0; f < nfreqs; ++f){
112 4096 : for(uInt iep=0; iep < (uInt)nep; ++iep) {
113 2048 : setSourceCoeffsfromVec(iep);
114 2048 : success &= (*this)(tmpfluxes,tmperrors,mfreqs[f],true);
115 2048 : tmpfluxes.value(tmpfluxvec);
116 2048 : 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 2048 : fluxes_p[iep]=tmpfluxvec[0];
121 : }
122 2048 : if (istimevar_p) {
123 : //setup interpolator
124 0 : ScalarSampledFunctional<Double> dts(epochvec_p);
125 0 : ScalarSampledFunctional<Float> flxs(fluxes_p);
126 0 : Interpolate1D<Double, Float> interpolateFlux(dts,flxs);
127 0 : interpolateFlux.setMethod(getInterpMethod_p(interpmethod));
128 :
129 0 : values[f].setValue(interpolateFlux(mtime.get("d").getValue()));
130 0 : }
131 : else { // no interpolation for non-var source, use first row data
132 2048 : values[f].setValue(fluxes_p[0]);
133 : }
134 : }
135 :
136 : // success &= (*this)(values[f], errors[f], mfreqs[f]);
137 :
138 32 : return success;
139 32 : }
140 :
141 37 : Bool FluxCalcVQS::setSource(const String& sourceName, const MDirection& sourceDir)
142 : {
143 37 : srcEnum_p = srcNameToEnum(sourceName,sourceDir);
144 : //return srcEnum_p != FCQS::UNKNOWN_SOURCE;
145 37 : return srcEnum_p != FSS::UNKNOWN_SOURCE;
146 : }
147 :
148 37 : FluxCalcVQS::Source FluxCalcVQS::getSrcEnum()
149 : {
150 37 : return srcEnum_p;
151 : }
152 :
153 32 : void FluxCalcVQS::readQSCoeffsTable(const Path& fileName)
154 : {
155 : //table containing the coefficents has
156 : //epoch, vector of coefficients
157 32 : const String& fullName = fileName.absoluteName();
158 64 : LogIO os(LogOrigin("FluxCalcVQS", "readQSCoeffsTable", WHERE));
159 : os << LogIO::NORMAL1
160 : << "Reading the coefficient data from a table, "<< fullName
161 32 : << LogIO::POST;
162 :
163 32 : AlwaysAssert(Table::isReadable(fullName), AipsError);
164 32 : Table_p = Table(fullName, Table::Old);
165 : ///TableProxy tab2(Table_p);
166 : //String srcName(names_p[srcEnum_p](0));
167 32 : String srcName(EnumToSrcName(srcEnum_p));
168 32 : String srcCoeffColName=srcName+"_coeffs";
169 32 : String srcCoeffErrorColName=srcName+"_coefferrs";
170 :
171 : //check the source data exist in the table
172 32 : const ColumnDescSet& cds=Table_p.tableDesc().columnDescSet();
173 32 : if (!cds.isDefined(srcCoeffColName))
174 0 : throw(AipsError(srcName+" does not appears to be in "+fullName));
175 32 : const ScalarColumn<Double> epochCol(Table_p, "Epoch");
176 32 : const ArrayColumn<Float> CoeffCol(Table_p, srcCoeffColName);
177 32 : const ArrayColumn<Float> CoeffErrorCol(Table_p, srcCoeffErrorColName);
178 : // check if col contains a valid freq range info in a keyword
179 32 : if (CoeffCol.keywordSet().isDefined("ValidFreqRange")) {
180 32 : Vector<Double> validfreqRange;
181 32 : String freqRangeUnit;
182 32 : CoeffCol.keywordSet().asRecord("ValidFreqRange").get("freqRange",validfreqRange);
183 32 : CoeffCol.keywordSet().asRecord("ValidFreqRange").get("freqRangeUnit",freqRangeUnit);
184 32 : Vector<MFrequency> mvalidfreqs;
185 32 : validfreqrange_p(0) = MFrequency(Quantity(validfreqRange(0),freqRangeUnit), MFrequency::TOPO);
186 32 : validfreqrange_p(1) = MFrequency(Quantity(validfreqRange(1),freqRangeUnit), MFrequency::TOPO);
187 32 : }
188 : else {
189 0 : validfreqrange_p(0) = MFrequency(Quantity(0.0,"GHz"), MFrequency::TOPO);
190 0 : validfreqrange_p(1) = MFrequency(Quantity(0.0,"GHz"), MFrequency::TOPO);
191 : }
192 32 : Vector<Double> tempEpochs;
193 32 : epochCol.getColumn(tempEpochs,true);
194 32 : CoeffCol.getColumn(coeffsmat_p,true);
195 32 : CoeffErrorCol.getColumn(coefferrsmat_p,true);
196 : //convert the epoch (year + fraction) to mjd
197 32 : convertYearFracToMjd(tempEpochs,epochvec_p);
198 : os << LogIO::DEBUG1
199 : << "nepoch="<<epochvec_p.nelements()
200 64 : << "coeff_0 for the first epoch (coeffsmat_p(0,0))="<<coeffsmat_p(0,0)
201 32 : << LogIO::POST;
202 :
203 32 : }
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 0 : Interpolate1D<Double,Float>::Method FluxCalcVQS::getInterpMethod_p(const String& interpmethod)
214 : {
215 0 : if (interpmethod.contains("nearest"))
216 0 : 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 32 : void FluxCalcVQS::convertYearFracToMjd(const Vector<Double>& yearfrac, Vector<Double>& mjds)
227 : {
228 32 : uInt daysintheyr=365;
229 32 : uInt n = yearfrac.nelements();
230 32 : mjds.resize(n);
231 608 : for (uInt i=0; i < n; i++) {
232 576 : if (Time::isLeapYear(uInt(yearfrac(i)))) daysintheyr=366;
233 : Double fintyr, fyrfrac;
234 576 : fyrfrac = modf((Double)yearfrac(i),&fintyr);
235 576 : Float days=fyrfrac*daysintheyr;
236 576 : Time time0(Int(yearfrac(i)),1,1);
237 576 : Double mjdtime0 = time0.modifiedJulianDay();
238 :
239 576 : mjds[i] = mjdtime0 + days;
240 : }
241 32 : }
242 :
243 2048 : void FluxCalcVQS::setSourceCoeffsfromVec(uInt& i)
244 : {
245 : //Vector<Float> err(4,0.0);
246 2048 : tvcoeffs_p(0)=coeffsmat_p.column(i);
247 : //cerr<<"coeffsmat_p.column="<<coeffsmat_p.column(i)(0)<<endl;
248 2048 : tvcoeffs_p(1)=coefferrsmat_p.column(i);
249 2048 : }
250 :
251 32 : void FluxCalcVQS::isTimeVar(Bool istimevar)
252 : {
253 32 : istimevar_p=istimevar;
254 32 : }
255 :
256 : } //# NAMESPACE CASA - END
|
