/* 3 and 9 level van vleck correction based on Fred Schwab's excellent paper: DRAFT, Van Vleck Correction for the GBT correlator, Feb 11, 2002 Fred Schwab (NRAO) These routines were proven by Bill Sisk (NAIC) bsisk@naic.edu These routines were converted to C by Jeff Hagen (NAIC) jeffh@naic.edu to use these routines: vanvleck3lev( double *, npts ); vanvleck9lev( double *, npts ); vanvleck3levf( float *, npts ); vanvleck9levf( float *, npts ); double pow3lev(zho); - compute 3 level power based on zero lag (1 optimium) double pow9lev(zho); - compute 9 level power based on zero lag (1 optimium) double attndb3lev(zho); - compute db correction to optimium power double attndb9lev(zho); - compute db correction to optimium power a few other things left over: double cal2mjd(int iy, int im, int id) compute mjd */ #include #include #define NAT_SPLINE 1.0e30 #define FITPTS 65 double r3(double, double, double, double, double ); double r9(double, double, double, double, double ); double thresh(int, double); double inverse_cerf(double); double inverse_erf(double); int spline( int, double *, double *, double, double, double *); double splint( double *, double *, double *, int, double ); vanvleck3lev( acf, n ) double *acf; int n; { int i; double zerolag, r0, rtr0, alpha3, alpha1, mu, v, rt2; double *rho, *r, *spl; zerolag = *acf; r0 = 0.0; /* skip offset removal for now avglen = 1000; for(i = n - avglen; i< n; i++) r0 += data[i]; r0 /= avglen; if( r0< 0.0 ) rtr0 = 0.0; else rtr0 = sqrt(r0); */ rtr0 = 0.0; rt2 = sqrt(2.0); alpha3 = rt2 * inverse_cerf( zerolag + rtr0); alpha1 = rt2 * inverse_cerf( zerolag - rtr0); mu = 0.5 * (alpha1 - alpha3); v = 0.5 * (alpha1 + alpha3); rho = (double *)malloc( sizeof(double)*FITPTS); r = (double *)malloc( sizeof(double)*FITPTS); spl = (double *)malloc( sizeof(double)*FITPTS); for( i=0; i #define L(x,y,z) bvnd(x,y,z) #define MAX(x,y) ((x)>(y)?(x):(y)) #define MIN(x,y) ((x)<(y)?(x):(y)) double PHID(double); double sign(double,double); double bvnd(double, double, double ); #define STAR2(x) ((x)*(x)) double r3(mux,muy,v1x,v1y,rho) double mux,muy,v1x,v1y,rho; { double rt2, ret; rt2=sqrt(2.0); ret=0.5*(erf((-mux+v1x)/rt2)-erf((mux+v1x)/rt2)+ erf((-muy+v1y)/rt2)-erf((muy+v1y)/rt2))+ L(-mux-v1x,-muy-v1y,rho)+L(-mux-v1x,-muy+v1y,rho)+ L(-mux+v1x,-muy-v1y,rho)+L(-mux+v1x,-muy+v1y,rho)-1.0; return (ret); } double r9(mux,muy,v1x,v1y,rho) double mux,muy,v1x,v1y,rho; { double rt2, ret; rt2=sqrt(2.0); ret=-16.0+2.0*(-erf((mux-7*v1x)/rt2)-erf((mux-5*v1x)/rt2)- erf((mux-3*v1x)/rt2)+erf((-mux+v1x)/rt2)- erf((mux+v1x)/rt2)-erf((mux+3*v1x)/rt2)- erf((mux+5*v1x)/rt2)-erf((mux+7*v1x)/rt2)- erf((muy-7*v1y)/rt2)-erf((muy-5*v1y)/rt2)- erf((muy-3*v1y)/rt2)+erf((-muy+v1y)/rt2)- erf((muy+v1y)/rt2)-erf((muy+3*v1y)/rt2)- erf((muy+5*v1y)/rt2)-erf((muy+7*v1y)/rt2))+ L(-mux-7*v1x,-muy-7*v1y,rho)+L(-mux-7*v1x,-muy-5*v1y,rho)+ L(-mux-7*v1x,-muy-3*v1y,rho)+L(-mux-7*v1x,-muy-v1y,rho)+ L(-mux-7*v1x,-muy+v1y,rho)+L(-mux-7*v1x,-muy+3*v1y,rho)+ L(-mux-7*v1x,-muy+5*v1y,rho)+L(-mux-7*v1x,-muy+7*v1y,rho)+ L(-mux-5*v1x,-muy-7*v1y,rho)+L(-mux-5*v1x,-muy-5*v1y,rho)+ L(-mux-5*v1x,-muy-3*v1y,rho)+L(-mux-5*v1x,-muy-v1y,rho)+ L(-mux-5*v1x,-muy+v1y,rho)+L(-mux-5*v1x,-muy+3*v1y,rho)+ L(-mux-5*v1x,-muy+5*v1y,rho)+L(-mux-5*v1x,-muy+7*v1y,rho)+ L(-mux-3*v1x,-muy-7*v1y,rho)+L(-mux-3*v1x,-muy-5*v1y,rho)+ L(-mux-3*v1x,-muy-3*v1y,rho)+L(-mux-3*v1x,-muy-v1y,rho)+ L(-mux-3*v1x,-muy+v1y,rho)+L(-mux-3*v1x,-muy+3*v1y,rho)+ L(-mux-3*v1x,-muy+5*v1y,rho)+L(-mux-3*v1x,-muy+7*v1y,rho)+ L(-mux-v1x,-muy-7*v1y,rho)+L(-mux-v1x,-muy-5*v1y,rho)+ L(-mux-v1x,-muy-3*v1y,rho)+L(-mux-v1x,-muy-v1y,rho)+ L(-mux-v1x,-muy+v1y,rho)+L(-mux-v1x,-muy+3*v1y,rho)+ L(-mux-v1x,-muy+5*v1y,rho)+L(-mux-v1x,-muy+7*v1y,rho)+ L(-mux+v1x,-muy-7*v1y,rho)+L(-mux+v1x,-muy-5*v1y,rho)+ L(-mux+v1x,-muy-3*v1y,rho)+L(-mux+v1x,-muy-v1y,rho)+ L(-mux+v1x,-muy+v1y,rho)+L(-mux+v1x,-muy+3*v1y,rho)+ L(-mux+v1x,-muy+5*v1y,rho)+L(-mux+v1x,-muy+7*v1y,rho)+ L(-mux+3*v1x,-muy-7*v1y,rho)+L(-mux+3*v1x,-muy-5*v1y,rho)+ L(-mux+3*v1x,-muy-3*v1y,rho)+L(-mux+3*v1x,-muy-v1y,rho)+ L(-mux+3*v1x,-muy+v1y,rho)+L(-mux+3*v1x,-muy+3*v1y,rho)+ L(-mux+3*v1x,-muy+5*v1y,rho)+L(-mux+3*v1x,-muy+7*v1y,rho)+ L(-mux+5*v1x,-muy-7*v1y,rho)+L(-mux+5*v1x,-muy-5*v1y,rho)+ L(-mux+5*v1x,-muy-3*v1y,rho)+L(-mux+5*v1x,-muy-v1y,rho)+ L(-mux+5*v1x,-muy+v1y,rho)+L(-mux+5*v1x,-muy+3*v1y,rho)+ L(-mux+5*v1x,-muy+5*v1y,rho)+L(-mux+5*v1x,-muy+7*v1y,rho)+ L(-mux+7*v1x,-muy-7*v1y,rho)+L(-mux+7*v1x,-muy-5*v1y,rho)+ L(-mux+7*v1x,-muy-3*v1y,rho)+L(-mux+7*v1x,-muy-v1y,rho)+ L(-mux+7*v1x,-muy+v1y,rho)+L(-mux+7*v1x,-muy+3*v1y,rho)+ L(-mux+7*v1x,-muy+5*v1y,rho)+L(-mux+7*v1x,-muy+7*v1y,rho); return(ret); } /* * A function for computing bivariate normal probabilities. * * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenzwsu.edu * * This function is based on the method described by * Drezner, Z and G.O. Wesolowsky, (1989), * On the computation of the bivariate normal inegral, * Journal of Statist. Comput. Simul. 35, pp. 101-107, * with major modifications for double precision, and for |R| close to 1. * * BVND - calculate the probability that X is larger than DH and Y is * larger than DK. * * Parameters * * DH DOUBLE PRECISION, integration limit * DK DOUBLE PRECISION, integration limit * R DOUBLE PRECISION, correlation coefficient */ double bvnd( dh, dk, r ) double dh, dk, r; { static double ZERO, TWOPI; static int flag = 1; static double X[10][3], W[10][3]; int I, IS, LG, NG; double AS, A, B, C, D, RS, XS, BVN; double SN, ASR, H, K, BS, HS, HK; /* Gauss Legendre Points and Weights, N = 6 */ if( flag ) { ZERO = 0.0; TWOPI = 6.283185307179586e0; /* DATA ( W(I,1), X(I,1), I = 1,3) */ W[0][0] = 0.1713244923791705e+00; X[0][0] = -0.9324695142031522e+00; W[1][0] = 0.3607615730481384e+00; X[1][0] = -0.6612093864662647e+00; W[2][0] = 0.4679139345726904e+00; X[2][0] = -0.2386191860831970e+00; /* Gauss Legendre Points and Weights, N = 12 */ /* DATA ( W(I,2), X(I,2), I = 1,6) */ W[0][1] = 0.4717533638651177e-01; X[0][1] = -0.9815606342467191e+00; W[1][1] = 0.1069393259953183e+00; X[1][1] = -0.9041172563704750e+00; W[2][1] = 0.1600783285433464e+00; X[2][1] = -0.7699026741943050e+00; W[3][1] = 0.2031674267230659e+00; X[3][1] = -0.5873179542866171e+00; W[4][1] = 0.2334925365383547e+00; X[4][1] = -0.3678314989981802e+00; W[5][1] = 0.2491470458134029e+00; X[5][1] = -0.1252334085114692e+00; /* Gauss Legendre Points and Weights, N = 20 */ /* DATA ( W(I,3), X(I,3), I = 1, 10 ) */ W[0][2] = 0.1761400713915212e-01; X[0][2] = -0.9931285991850949e+00; W[1][2] = 0.4060142980038694e-01; X[1][2] = -0.9639719272779138e+00; W[2][2] = 0.6267204833410906e-01; X[2][2] = -0.9122344282513259e+00; W[3][2] = 0.8327674157670475e-01; X[3][2] = -0.8391169718222188e+00; W[4][2] = 0.1019301198172404e+00; X[4][2] = -0.7463319064601508e+00; W[5][2] = 0.1181945319615184e+00; X[5][2] = -0.6360536807265150e+00; W[6][2] = 0.1316886384491766e+00; X[6][2] = -0.5108670019508271e+00; W[7][2] = 0.1420961093183821e+00; X[7][2] = -0.3737060887154196e+00; W[8][2] = 0.1491729864726037e+00; X[8][2] = -0.2277858511416451e+00; W[9][2] = 0.1527533871307259e+00; X[9][2] = -0.7652652113349733e-01; flag = 0; } if( fabs(r) < 0.3 ) { NG = 0; LG = 3; } else if ( fabs(r) < 0.75 ) { NG = 1; LG = 6; } else { NG = 2; LG = 10; } H = dh; K = dk; HK = H*K; BVN = 0; if( fabs(r) < 0.925 ) { HS = ( H*H + K*K )/2; ASR = asin(r); for( I=0; I -100 ) { BVN = A*exp(ASR)*( 1.0 - C*( BS - AS )*( 1.0 - D*BS/5 )/3.0 + C*D*AS*AS/5 ); } if ( -HK < 100 ) { B = sqrt(BS); BVN -= exp( -HK/2 )*sqrt(TWOPI)*PHID(-B/A)*B *( 1 - C*BS*( 1 - D*BS/5.0 )/3.0 ); } A = A/2.0; for( I=0; I -100 ) { BVN = BVN + A*W[I][NG]*exp( ASR ) *( exp( -HK*( 1 - RS )/( 2*( 1 + RS ) ) )/RS - ( 1 + C*XS*( 1 + D*XS ) ) ); } } } BVN = -BVN/TWOPI; } if ( r > 0.0) BVN = BVN + PHID( -MAX( H, K ) ); if ( r < 0.0 ) BVN = -BVN + MAX( ZERO, PHID(-H) - PHID(-K) ); } return BVN; } /* * Normal distribution probabilities accurate to 1.e-15. * Z = no. of standard deviations from the mean. * * Based upon algorithm 5666 for the error function, from: * Hart, J.F. et al, 'Computer Approximations', Wiley 1968 * * Programmer: Alan Miller * * Latest revision - 30 March 1986 */ double PHID(Z) double Z; { static double P0, P1, P2, P3, P4, P5, P6; static double Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7; static double CUTOFF, ROOTPI; static flag = 1; double EXPNTL, ZABS, P; if( flag ) { P0 = 220.2068679123761; P1 = 221.2135961699311; P2 = 112.0792914978709; P3 = 33.91286607838300; P4 = 6.373962203531650; P5 = .7003830644436881; P6 = .03526249659989109; Q0 = 440.4137358247522; Q1 = 793.8265125199484; Q2 = 637.3336333788311; Q3 = 296.5642487796737; Q4 = 86.78073220294608; Q5 = 16.06417757920695; Q6 = 1.755667163182642; Q7 = .08838834764831844; ROOTPI = 2.506628274631001; CUTOFF = 7.071067811865475; flag = 0; } ZABS = fabs(Z); /* * |Z| > 37 */ if ( ZABS > 37 ) { P = 0; } else { /* * |Z| <= 37 */ EXPNTL = exp( -STAR2(ZABS)/2.0 ); /* * |Z| < CUTOFF = 10/sqrt(2) */ if ( ZABS < CUTOFF ) { P = EXPNTL*( ( ( ( ( ( P6*ZABS + P5 )*ZABS + P4 )*ZABS + P3 )*ZABS + P2 )*ZABS + P1 )*ZABS + P0 ) /( ( ( ( ( ( ( Q7*ZABS + Q6 )*ZABS + Q5 )*ZABS + Q4 )*ZABS + Q3 )*ZABS + Q2 )*ZABS + Q1 )*ZABS + Q0 ); /* * |Z| >= CUTOFF. */ } else { P = EXPNTL/( ZABS + 1/( ZABS + 2/( ZABS + 3/( ZABS + 4/( ZABS + 0.65 ) ) ) ) )/ROOTPI; } } if ( Z > 0 ) P = 1 - P; return(P); } double thresh(n,zerolag) int n; double zerolag; { double pi = M_PI; double x, tol, f, fp, deltax; int itmax, k, i; x=0.0; if( n%2 == 0) x=1.0; itmax=30; tol=1.0e-8; for( i=1; i<=itmax; i++ ) { f = zerolag; fp = 0.0; if( n%2 == 1) { for( k=1; k<=(n-1)/2; k++ ) { f = f-(2.0*k-1.0)*erfc((2.0*k-1.0)*x/sqrt(2.0)); fp = fp+sqrt(2.0/pi)*STAR2(2.0*k-1.0)*exp(-0.5*STAR2((2.0*k-1.0)*x)); } } else { f = f-1.0; for( k=1; k<=(n-1)/2; k++ ) { f=f-8*k*erfc(k*x/sqrt(2.0)); fp=fp+8*STAR2(k)*sqrt(2.0/pi)*exp(-0.5*STAR2(k*x)); } } deltax=-f/fp; deltax= sign(1.0,deltax)*MIN(0.5,fabs(deltax)); x=x+deltax; if (n%2 == 1) x = MAX(0.0,x); if( fabs(deltax/x) < tol ) break; } return(x); } /* This subroutine calculates the second derivatives of the interpolating cubic spline function. It was taken from "Numerical Recipes," by Press et al. 1986, p. 88. X, Y are the arrays of input interpolation values; N is the number of interpolation points; YP1, YPN are the first derivatives at the first and last points, resp. They can be set to 1.E30 for a natural spline; Y2 is the output array of second derivatives to be used by SPLINT. removed fortran stupidity - now 0 zero based arrays */ int spline(n, x, y, yp1, ypn, y2) int n; double x[], y[], yp1, ypn, y2[]; { double sig, p, qn, un; double *u; int i, k; u = (double *)malloc( sizeof(double) *n ); if(yp1 > 0.99e30) { y2[0] =0.0; u[0] =0.0; } else { y2[0] = -0.5; u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1); } for(i=1; i <= n-2; i++) { sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]); p = sig*y2[i-1]+2.0; y2[i] = (sig-1.0)/p; u[i] = (6.0*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1]) /(x[i]-x[i-1]))/(x[i+1]-x[i-1])-sig*u[i-1])/p; } if (ypn > 0.99e30) { qn = 0.0; un = 0.0; } else { qn = 0.5; un = (3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); } y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0); for(k = n-2; k >= 0; k--) y2[k] = y2[k]*y2[k+1] + u[k]; free(u); } /* This subroutine performs a cubic spline interpolation. It was taken from "Numerical Recipes," by Press et al. 1986, p. 88. XA, YA are the arrays of input interpolation data; Y2A is the array of second derivatives produced by SPLINE; N is the number of input interpolation points; X is the X value for which we wish an interpolation; Y is the output interpolation value. Modified 20 Oct 89 to change x**3 to x*x*x Mod. 25 Oct 89 to take out PAUSE and replace with call to WARNING remove fortran stupidity: return double as answer and make arrays zero based */ double splint(xa, ya, y2a, n, x) double xa[], ya[], y2a[], x; int n; { double h, a, b, y; int k, klo, khi, errlev; klo = 0; khi = n-1; while(khi - klo > 1) { k = (int)((khi+klo)/2); if(xa[k] > x) khi = k; else klo = k; } h = xa[khi] - xa[klo]; if(h == 0.0) { fprintf(stderr, "WARN Bad interpolation in routine splint"); return(x); } a = (xa[khi]-x)/h; b = (x-xa[klo])/h; y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; return(y); } double sign( double a, double b ) { if( b < 0.0 ) return( -fabs(a)); else return( fabs(a)); } /*---------------------------------------------------------------------* * return modified julian date given a gregorian calendar - based on slalib * routine sla_cldj which itself is based on Hatcher (1984) QJRAS 25, 53-55 * creation date 1999/07/10 [mjd=51639] (dunc@naic.edu) *---------------------------------------------------------------------*/ double cal2mjd(int iy, int im, int id) { int leap; /* month lengths in days for normal and leap years */ static int mtab[2][13] = { {0,31,28,31,30,31,30,31,31,30,31,30,31}, {0,31,29,31,30,31,30,31,31,30,31,30,31} }; /*validate year*/ if (iy<-4699) { fprintf(stderr,"Invalid year passed to cal2mjd\n"); return(0.0); } else { /* validate month */ if (im<1 || im>12) { fprintf(stderr,"Invalid month passed to cal2mjd\n"); return(0.0); } else { /* allow for leap year */ leap = iy%4 == 0 && iy%100 != 00 || iy%400 == 0; /* validate day */ if (id<1 || id>mtab[leap][im]) { fprintf(stderr,"Invalid day passed to cal2mjd\n"); return(0.0); } } } return ((1461*(iy-(12-im)/10+4712))/4+(5+306*((im+9)%12))/10 -(3*((iy-(12-im)/10+4900)/100))/4+id-2399904); } /* round an integer UP to the next factor of two this used to ensure an fft is applied on factors of 2, even if number of lags is an odd value Thus 1000 -> 1024 or 33 -> 64, etc */ int makefactor2(int in) { double tmp; int if2; tmp = ceil( log10((double)in)/log10(2.0)); if2 = rint(pow(2.0,tmp)); return(if2); } /* Hamming Window Function - Standard time domain version: + + + + + + + + ---------------------------------------- 0 1 2 3 4 5 6 7 npts = 8 Return the hanning window factor for index x, note index is "C" type i.e. index = 0,1,2, ... N-1 Reference "Signals and Systems" by Ziemer, et al */ double hamming(int index, int npts) { return( 0.54 - 0.46*cos(2*M_PI*((double)index + 0.5)/((double)npts))); }