Loading...
Searching...
No Matches
VanVleck.h
Go to the documentation of this file.
//# VanVleck.h: Class of static functions to aid with vanVleck corrections.
//# Copyright (C) 2002
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//# Internet email: casa-feedback@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
#ifndef SCIMATH_VANVLECK_H
#define SCIMATH_VANVLECK_H
//#! Includes go here
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/scimath/Functionals/Interpolate1D.h>
#include <casacore/casa/BasicSL/Constants.h>
#include <mutex>
namespace casacore { //# NAMESPACE CASACORE - BEGIN
//# Forward Declarations
// <summary>
// A class of static functions to aid with vanVleck corrections of lag data.
// </summary>
// <use visibility=export>
// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>
// <prerequisite>
// <li> Familiarity with the issues involved in turning digitally
// sampled lag data from a correlator into spectral data.
// </prerequisite>
//
// <etymology>
// This provides the functions necessary to determine the van Vleck correction
// for a general n-level by m-level correlator.
// </etymology>
//
// <synopsis>
// This provides the functions necessary to determine the van Vleck correction
// for a general n-level by m-level correlator.
// </synopsis>
//
// <example>
// </example>
//
// <motivation>
// The GBT spectrometer provides the measured auto-correlation and
// cross-correlation lags. The gbt MeasurementSet filler (gbtmsfiller)
// needs to convert those lags to the spectral domain. These functions
// allow the filler to calculate the van Vleck correction appropriate
// for each measured set of lags. They are of general and hence are
// not specific to the GBT spectrometer.
//
// The functions here are static because of the nature of the underlying
// numerical quadrature fortran code used to integrate the
// drbyrho function.
// </motivation>
//
// <thrown>
// <li>
// <li>
// </thrown>
//
// <todo asof="2002/07/19">
// <li> The inverse error functions may be more generally useful.
// It exists here only as a private member function to be
// used internally.
// </todo>
class VanVleck
{
public:
// Set the interpolation table size.
// Should be an odd number. The default size is 65.
static void size(uInt npts);
// get the current size.
static uInt getsize();
// Set the x and y quantization functions.
// Each matrix should have dimensions (2,n)
// where n is the number of levels. The first
// row (0,...) is the (n-1) threshold levels and
// the second row is the n quantizations based
// on those thresholds. The thresholds may
// include a DC offset. The (0,(n-1)) element is
// never used and need not be set.
static void setQuantization(const Matrix<Double> &qx,
const Matrix<Double> &qy);
// Set the x and y quantization levels for the case
// of equi-spaced levels with a possible non-zero
// offset. The total number of levels is given by n,
// which must be 3 or 9. If n is not 3 or 9, False
// will be returned and no quantization will have been
// set. For the 3- and 9- level cases a bivarate normal
// integral calculation will be used. That is much faster
// than the more general numerical integration used
// by setQuantization.
static Bool setEquiSpaced(Double xlev, Double ylev,
Double xmean, Double ymean,
Int n);
// Get the data used in setting up the interpolation
static void getTable(Vector<Double> &rs, Vector<Double> &rhos);
// Given a rho return the corresponding corrected r
// Returns 0.0 if no quantization has been set yet.
static Double r(const Double rho);
// Given a measured zero-lag autocorrelation and number of
// levels (n>=3) return the first positive quantizer input
// threshold level. This can be used to set the up the
// matrix arguments used in setQuantization.
static Double thresh(Int n, Double zerolag)
{ return ( (n>3) ? threshNgt3(n,zerolag) : threshN3(zerolag) ); }
// Predict a given zero-lag given n and a threshold. This
// is included here to be used as a check against the output
// of thresh.
static Double predict(Int n, Double threshhold)
{ return ( (n>3) ? predictNgt3(n,threshhold) : predictN3(threshhold));}
// Compute an approximation to the mean signal level (DC offset)
// and quantizer threshold setting (both in terms of the r.m.s.
// signal input level) given the observed positive bias (the
// asymptotic limit of the measured autocorrelation at large
// lags) and the zero-lag autocorrelation.
// dcoffset is the mean signal level, threshold is the quantizer
// setting, n is the number of levels, zerolag is the zero-lag
// value and bias is the asymptotic bias.
// Currently, this is only available for the n==3 level case,
// all other cases set the returned dcoffset to 0 and use thresh()
// to set the returned value of threshold. A return value of F
// indicates that the zerolag and bias values are inconsistent
// and the dcoffset can not be determined. In that case,
// the returned dcoffset is 0 and thresh() is used to set
// the threshold level.
static Bool dcoff(Double &dcoffset, Double &threshold,
Int n, Double zerolag, Double bias);
private:
// the number of points to use in setting up the interpolator
static uInt itsSize, itsNx, itsNy;
static Bool itsEquiSpaced;
static Double itsXlev, itsYlev, itsXmean, itsYmean;
// The interpolator
static Interpolate1D<Double, Double> *itsInterp;
// the quantization functions
static Vector<Double> itsQx0, itsQx1, itsQy0, itsQy1;
// Useful combinations of the above - to speed up drbydrho
// these are -1/2*(Qx0*Qx0) and -1/2*(Qy0*Qy0)
// These are only used for i up to (itsQx0.nelements() and
// for j up to (itsQy0.nelements()).
static Vector<Double> itsQx0Qx0, itsQy0Qy0;
// This is Qx0[i]*Qy0[j]
static Matrix<Double> itsQx0Qy0;
// This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j])
static Matrix<Double> itsQx1Qy1diffs;
// The mutex to make the functions thread-safe.
static std::mutex theirMutex;
// The fortran numerical integration function will call this.
// For a given rho and quantization functions, this computes,
// via Price's theorem, the value dr/drho of the derivative,
// with respect to rho, of the expected value of the correlator
// output.
static Double drbydrho(Double *rho);
// For a given rhoi, rhof, this produces a high-accuracy numerical
// approximation to the integral of drbydrho over the range
// rhoi to rhof. It calls the standard QUADPACK adaptive Gaussian quadrature
// procedure, dqags, to do the numerical integration.
static Double rinc(Double &rhoi, Double &rhof);
// initialize the interpolator
static void initInterpolator();
// compute first threshhold for a given zerolag for n>3
static Double threshNgt3(Int n, Double zerolag);
// compute first threshhold for a given zerolag for n==3
static Double threshN3(Double zerolag)
{ return sqrt(2.0)*invErfc(zerolag);}
// inverse err fn - used by invErfc
static Double invErf(Double x);
// inverse complementary err fn - used by threshN3
static Double invErfc(Double x);
// Predict a zero-lag value given the indicated first threshold level
// for n>3.
static Double predictNgt3(Int n, Double threshhold);
// Predict a zero-lag value given the indicated first threshold level
// for n=3.
static Double predictN3(Double threshhold)
{ return ::erfc(threshhold/sqrt(2.0));}
// implementation of dcoff for the 3-level case
static Bool dcoff3(Double &dcoffset, Double &threshold,
Double zerolag, Double bias);
};
} //# NAMESPACE CASACORE - END
#endif
casacore::Interpolate1D
Definition Interpolate1D.h:140
casacore::Matrix
Definition Matrix.h:87
casacore::VanVleck
Definition VanVleck.h:93
casacore::VanVleck::predictNgt3
static Double predictNgt3(Int n, Double threshhold)
Predict a zero-lag value given the indicated first threshold level for n>3.
casacore::VanVleck::threshNgt3
static Double threshNgt3(Int n, Double zerolag)
compute first threshhold for a given zerolag for n>3
casacore::VanVleck::itsXmean
static Double itsXmean
Definition VanVleck.h:171
casacore::VanVleck::itsQy1
static Vector< Double > itsQy1
Definition VanVleck.h:177
casacore::VanVleck::itsQx1Qy1diffs
static Matrix< Double > itsQx1Qy1diffs
This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j]).
Definition VanVleck.h:187
casacore::VanVleck::itsQx1
static Vector< Double > itsQx1
Definition VanVleck.h:177
casacore::VanVleck::dcoff3
static Bool dcoff3(Double &dcoffset, Double &threshold, Double zerolag, Double bias)
implementation of dcoff for the 3-level case
casacore::VanVleck::itsQx0Qx0
static Vector< Double > itsQx0Qx0
Useful combinations of the above - to speed up drbydrho these are -1/2*(Qx0*Qx0) and -1/2*(Qy0*Qy0) T...
Definition VanVleck.h:183
casacore::VanVleck::getTable
static void getTable(Vector< Double > &rs, Vector< Double > &rhos)
Get the data used in setting up the interpolation.
casacore::VanVleck::threshN3
static Double threshN3(Double zerolag)
compute first threshhold for a given zerolag for n==3
Definition VanVleck.h:212
casacore::VanVleck::predict
static Double predict(Int n, Double threshhold)
Predict a given zero-lag given n and a threshold.
Definition VanVleck.h:143
casacore::VanVleck::itsQx0
static Vector< Double > itsQx0
the quantization functions
Definition VanVleck.h:177
casacore::VanVleck::size
static void size(uInt npts)
Set the interpolation table size.
casacore::VanVleck::itsYlev
static Double itsYlev
Definition VanVleck.h:171
casacore::VanVleck::itsNx
static uInt itsNx
Definition VanVleck.h:167
casacore::VanVleck::itsInterp
static Interpolate1D< Double, Double > * itsInterp
The interpolator.
Definition VanVleck.h:174
casacore::VanVleck::itsQx0Qy0
static Matrix< Double > itsQx0Qy0
This is Qx0[i]*Qy0[j].
Definition VanVleck.h:185
casacore::VanVleck::itsXlev
static Double itsXlev
Definition VanVleck.h:171
casacore::VanVleck::rinc
static Double rinc(Double &rhoi, Double &rhof)
For a given rhoi, rhof, this produces a high-accuracy numerical approximation to the integral of drby...
casacore::VanVleck::itsEquiSpaced
static Bool itsEquiSpaced
Definition VanVleck.h:169
casacore::VanVleck::setQuantization
static void setQuantization(const Matrix< Double > &qx, const Matrix< Double > &qy)
Set the x and y quantization functions.
casacore::VanVleck::itsNy
static uInt itsNy
Definition VanVleck.h:167
casacore::VanVleck::getsize
static uInt getsize()
get the current size.
casacore::VanVleck::itsSize
static uInt itsSize
the number of points to use in setting up the interpolator
Definition VanVleck.h:167
casacore::VanVleck::invErf
static Double invErf(Double x)
inverse err fn - used by invErfc
casacore::VanVleck::invErfc
static Double invErfc(Double x)
inverse complementary err fn - used by threshN3
casacore::VanVleck::thresh
static Double thresh(Int n, Double zerolag)
Given a measured zero-lag autocorrelation and number of levels (n>=3) return the first positive quant...
Definition VanVleck.h:137
casacore::VanVleck::drbydrho
static Double drbydrho(Double *rho)
The fortran numerical integration function will call this.
casacore::VanVleck::initInterpolator
static void initInterpolator()
initialize the interpolator
casacore::VanVleck::r
static Double r(const Double rho)
Given a rho return the corresponding corrected r Returns 0.0 if no quantization has been set yet.
casacore::VanVleck::predictN3
static Double predictN3(Double threshhold)
Predict a zero-lag value given the indicated first threshold level for n=3.
Definition VanVleck.h:227
casacore::VanVleck::theirMutex
static std::mutex theirMutex
The mutex to make the functions thread-safe.
Definition VanVleck.h:189
casacore::VanVleck::dcoff
static Bool dcoff(Double &dcoffset, Double &threshold, Int n, Double zerolag, Double bias)
Compute an approximation to the mean signal level (DC offset) and quantizer threshold setting (both i...
casacore::VanVleck::itsYmean
static Double itsYmean
Definition VanVleck.h:171
casacore::VanVleck::itsQy0Qy0
static Vector< Double > itsQy0Qy0
Definition VanVleck.h:183
casacore::VanVleck::itsQy0
static Vector< Double > itsQy0
Definition VanVleck.h:177
casacore::VanVleck::setEquiSpaced
static Bool setEquiSpaced(Double xlev, Double ylev, Double xmean, Double ymean, Int n)
Set the x and y quantization levels for the case of equi-spaced levels with a possible non-zero offse...
casacore::Vector
Definition Vector.h:82
casacore
this file contains all the compiler specific defines
Definition mainpage.dox:28
casacore::uInt
unsigned int uInt
Definition aipstype.h:49
casacore::sqrt
LatticeExprNode sqrt(const LatticeExprNode &expr)
casacore::Int
int Int
Definition aipstype.h:48
casacore::Bool
bool Bool
Define the standard types used by Casacore.
Definition aipstype.h:40
casacore::Double
double Double
Definition aipstype.h:53