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
39namespace 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
95{
96public:
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
167private:
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
219
220 // inverse complementary err fn - used by threshN3
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