casacore
VanVleck.h
Go to the documentation of this file.
1 //# VanVleck.h: Class of static functions to aid with vanVleck corrections.
2 //# Copyright (C) 2002
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: aips2-request@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 //# $Id$
27 
28 #ifndef SCIMATH_VANVLECK_H
29 #define SCIMATH_VANVLECK_H
30 
31 //#! Includes go here
32 #include <casacore/casa/aips.h>
33 #include <casacore/casa/Arrays/Matrix.h>
34 #include <casacore/scimath/Functionals/Interpolate1D.h>
35 #include <casacore/casa/BasicSL/Constants.h>
36 
37 #include <mutex>
38 
39 namespace casacore { //# NAMESPACE CASACORE - BEGIN
40 
41 //# Forward Declarations
42 
43 // <summary>
44 // A class of static functions to aid with vanVleck corrections of lag data.
45 // </summary>
46 
47 // <use visibility=export>
48 
49 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
50 // </reviewed>
51 
52 // <prerequisite>
53 // <li> Familiarity with the issues involved in turning digitally
54 // sampled lag data from a correlator into spectral data.
55 // </prerequisite>
56 //
57 // <etymology>
58 // This provides the functions necessary to determine the van Vleck correction
59 // for a general n-level by m-level correlator.
60 // </etymology>
61 //
62 // <synopsis>
63 // This provides the functions necessary to determine the van Vleck correction
64 // for a general n-level by m-level correlator.
65 // </synopsis>
66 //
67 // <example>
68 // </example>
69 //
70 // <motivation>
71 // The GBT spectrometer provides the measured auto-correlation and
72 // cross-correlation lags. The gbt MeasurementSet filler (gbtmsfiller)
73 // needs to convert those lags to the spectral domain. These functions
74 // allow the filler to calculate the van Vleck correction appropriate
75 // for each measured set of lags. They are of general and hence are
76 // not specific to the GBT spectrometer.
77 //
78 // The functions here are static because of the nature of the underlying
79 // numerical quadrature fortran code used to integrate the
80 // drbyrho function.
81 // </motivation>
82 //
83 // <thrown>
84 // <li>
85 // <li>
86 // </thrown>
87 //
88 // <todo asof="2002/07/19">
89 // <li> The inverse error functions may be more generally useful.
90 // It exists here only as a private member function to be
91 // used internally.
92 // </todo>
93 
94 class VanVleck
95 {
96 public:
97  // Set the interpolation table size.
98  // Should be an odd number. The default size is 65.
99  static void size(uInt npts);
100 
101  // get the current size.
102  static uInt getsize();
103 
104  // Set the x and y quantization functions.
105  // Each matrix should have dimensions (2,n)
106  // where n is the number of levels. The first
107  // row (0,...) is the (n-1) threshold levels and
108  // the second row is the n quantizations based
109  // on those thresholds. The thresholds may
110  // include a DC offset. The (0,(n-1)) element is
111  // never used and need not be set.
112  static void setQuantization(const Matrix<Double> &qx,
113  const Matrix<Double> &qy);
114 
115  // Set the x and y quantization levels for the case
116  // of equi-spaced levels with a possible non-zero
117  // offset. The total number of levels is given by n,
118  // which must be 3 or 9. If n is not 3 or 9, False
119  // will be returned and no quantization will have been
120  // set. For the 3- and 9- level cases a bivarate normal
121  // integral calculation will be used. That is much faster
122  // than the more general numerical integration used
123  // by setQuantization.
124  static Bool setEquiSpaced(Double xlev, Double ylev,
125  Double xmean, Double ymean,
126  Int n);
127 
128  // Get the data used in setting up the interpolation
129  static void getTable(Vector<Double> &rs, Vector<Double> &rhos);
130 
131  // Given a rho return the corresponding corrected r
132  // Returns 0.0 if no quantization has been set yet.
133  static Double r(const Double rho);
134 
135  // Given a measured zero-lag autocorrelation and number of
136  // levels (n>=3) return the first positive quantizer input
137  // threshold level. This can be used to set the up the
138  // matrix arguments used in setQuantization.
139  static Double thresh(Int n, Double zerolag)
140  { return ( (n>3) ? threshNgt3(n,zerolag) : threshN3(zerolag) ); }
141 
142  // Predict a given zero-lag given n and a threshold. This
143  // is included here to be used as a check against the output
144  // of thresh.
145  static Double predict(Int n, Double threshhold)
146  { return ( (n>3) ? predictNgt3(n,threshhold) : predictN3(threshhold));}
147 
148  // Compute an approximation to the mean signal level (DC offset)
149  // and quantizer threshold setting (both in terms of the r.m.s.
150  // signal input level) given the observed positive bias (the
151  // asymptotic limit of the measured autocorrelation at large
152  // lags) and the zero-lag autocorrelation.
153  // dcoffset is the mean signal level, threshold is the quantizer
154  // setting, n is the number of levels, zerolag is the zero-lag
155  // value and bias is the asymptotic bias.
156  // Currently, this is only available for the n==3 level case,
157  // all other cases set the returned dcoffset to 0 and use thresh()
158  // to set the returned value of threshold. A return value of F
159  // indicates that the zerolag and bias values are inconsistent
160  // and the dcoffset can not be determined. In that case,
161  // the returned dcoffset is 0 and thresh() is used to set
162  // the threshold level.
163  static Bool dcoff(Double &dcoffset, Double &threshold,
164  Int n, Double zerolag, Double bias);
165 
166 
167 private:
168  // the number of points to use in setting up the interpolator
170 
172 
174 
175  // The interpolator
177 
178  // the quantization functions
180 
181  // Useful combinations of the above - to speed up drbydrho
182  // these are -1/2*(Qx0*Qx0) and -1/2*(Qy0*Qy0)
183  // These are only used for i up to (itsQx0.nelements() and
184  // for j up to (itsQy0.nelements()).
186  // This is Qx0[i]*Qy0[j]
188  // This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j])
190  // The mutex to make the functions thread-safe.
191  static std::mutex theirMutex;
192 
193 
194  // The fortran numerical integration function will call this.
195  // For a given rho and quantization functions, this computes,
196  // via Price's theorem, the value dr/drho of the derivative,
197  // with respect to rho, of the expected value of the correlator
198  // output.
199  static Double drbydrho(Double *rho);
200 
201  // For a given rhoi, rhof, this produces a high-accuracy numerical
202  // approximation to the integral of drbydrho over the range
203  // rhoi to rhof. It calls the standard QUADPACK adaptive Gaussian quadrature
204  // procedure, dqags, to do the numerical integration.
205  static Double rinc(Double &rhoi, Double &rhof);
206 
207  // initialize the interpolator
208  static void initInterpolator();
209 
210  // compute first threshhold for a given zerolag for n>3
211  static Double threshNgt3(Int n, Double zerolag);
212 
213  // compute first threshhold for a given zerolag for n==3
214  static Double threshN3(Double zerolag)
215  { return sqrt(2.0)*invErfc(zerolag);}
216 
217  // inverse err fn - used by invErfc
218  static Double invErf(Double x);
219 
220  // inverse complementary err fn - used by threshN3
221  static Double invErfc(Double x);
222 
223  // Predict a zero-lag value given the indicated first threshold level
224  // for n>3.
225  static Double predictNgt3(Int n, Double threshhold);
226 
227  // Predict a zero-lag value given the indicated first threshold level
228  // for n=3.
229  static Double predictN3(Double threshhold)
230  { return ::erfc(threshhold/sqrt(2.0));}
231 
232  // implementation of dcoff for the 3-level case
233  static Bool dcoff3(Double &dcoffset, Double &threshold,
234  Double zerolag, Double bias);
235 };
236 
237 
238 } //# NAMESPACE CASACORE - END
239 
240 #endif
static Double predictNgt3(Int n, Double threshhold)
Predict a zero-lag value given the indicated first threshold level for n>3.
static Double threshNgt3(Int n, Double zerolag)
compute first threshhold for a given zerolag for n>3
static Double itsXmean
Definition: VanVleck.h:173
static Vector< Double > itsQy1
Definition: VanVleck.h:179
static Matrix< Double > itsQx1Qy1diffs
This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j])
Definition: VanVleck.h:189
static Vector< Double > itsQx1
Definition: VanVleck.h:179
static Bool dcoff3(Double &dcoffset, Double &threshold, Double zerolag, Double bias)
implementation of dcoff for the 3-level case
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:185
static void getTable(Vector< Double > &rs, Vector< Double > &rhos)
Get the data used in setting up the interpolation.
static Double threshN3(Double zerolag)
compute first threshhold for a given zerolag for n==3
Definition: VanVleck.h:214
static Double predict(Int n, Double threshhold)
Predict a given zero-lag given n and a threshold.
Definition: VanVleck.h:145
static Vector< Double > itsQx0
the quantization functions
Definition: VanVleck.h:179
static void size(uInt npts)
Set the interpolation table size.
static Double itsYlev
Definition: VanVleck.h:173
static uInt itsNx
Definition: VanVleck.h:169
static Interpolate1D< Double, Double > * itsInterp
The interpolator.
Definition: VanVleck.h:176
static Matrix< Double > itsQx0Qy0
This is Qx0[i]*Qy0[j].
Definition: VanVleck.h:187
static Double itsXlev
Definition: VanVleck.h:173
static Double rinc(Double &rhoi, Double &rhof)
For a given rhoi, rhof, this produces a high-accuracy numerical approximation to the integral of drby...
static Bool itsEquiSpaced
Definition: VanVleck.h:171
static void setQuantization(const Matrix< Double > &qx, const Matrix< Double > &qy)
Set the x and y quantization functions.
static uInt itsNy
Definition: VanVleck.h:169
static uInt getsize()
get the current size.
static uInt itsSize
the number of points to use in setting up the interpolator
Definition: VanVleck.h:169
static Double invErf(Double x)
inverse err fn - used by invErfc
static Double invErfc(Double x)
inverse complementary err fn - used by threshN3
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:139
static Double drbydrho(Double *rho)
The fortran numerical integration function will call this.
static void initInterpolator()
initialize the interpolator
static Double r(const Double rho)
Given a rho return the corresponding corrected r Returns 0.0 if no quantization has been set yet.
static Double predictN3(Double threshhold)
Predict a zero-lag value given the indicated first threshold level for n=3.
Definition: VanVleck.h:229
static std::mutex theirMutex
The mutex to make the functions thread-safe.
Definition: VanVleck.h:191
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...
static Double itsYmean
Definition: VanVleck.h:173
static Vector< Double > itsQy0Qy0
Definition: VanVleck.h:185
static Vector< Double > itsQy0
Definition: VanVleck.h:179
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...
this file contains all the compiler specific defines
Definition: mainpage.dox:28
unsigned int uInt
Definition: aipstype.h:51
LatticeExprNode sqrt(const LatticeExprNode &expr)
int Int
Definition: aipstype.h:50
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
double Double
Definition: aipstype.h:55