FUNCTION PCORRE, P1, P2, CORRELATION=CORRELATION, ERRMSG=ERRMSG
;+
; Project     :	SOHO - CDS
;
; Name        :	PCORRE()
;
; Purpose     :	Calculate significance of correlation coefficient.
;
; Category    :	Error-analysis, Class3
;
; Explanation :	Calculates the probability that a given cross-correlation
;		coefficient could be generated by chance.
;
; Syntax      :	Result = PCORRE(R, NPTS)
;		Result = PCORRE(P1, P2, CORRELATION=CORR)
;
; Examples    :	IDL> PRINT, PCORRE(0.1, 30)
;		      0.59904802
;		IDL> PRINT, PCORRE(0.5, 30)
;		    0.0048999337
;
;		This means that, for two random populations of 30 points each,
;		there's 59.9% probability that the cross-correlation
;		coefficient between them will be 0.1 or better.  However,
;		there's only a 0.49% probability of generating a coefficient of
;		0.5 or higher.
;
; Inputs      :	One of two sets of input parameters must be passed, either
;
;		  R	= The cross correlation coefficient to test.  This must
;			  be between -1 and 1.
;
;		  NPTS	= The number of data points in each of the two
;			  populations.  This must be at least 3.
;
;		or
;
;		  P1,P2	= Two arrays to take the cross-correlation of.  These
;			  must both have the same number of elements.
;
; Opt. Inputs :	None.
;
; Outputs     :	The result of the function is the probability of generating a
;		coefficient greater than or equal to R strictly by chance with
;		two arrays of NPTS values each.  This will be between 0 and 1.
;		If an error occurs, then -1 will be returned.
;
; Opt. Outputs:	None.
;
; Keywords    :	CORRELATION = If P1 and P2 are passed, then the correlation
;			      between then will be returned in this parameter.
;
;		ERRMSG = If defined and passed, then any error messages will be
;			 returned to the user in this parameter rather than
;			 depending on the MESSAGE routine in IDL.  If no errors
;			 are encountered, then a null string is returned.  In
;			 order to use this feature, ERRMSG must be defined
;			 first, e.g.
;
;				ERRMSG = ''
;				RESULT = PCORRE(ERRMSG=ERRMSG, ...)
;				IF ERRMSG NE '' THEN ...
;
; Calls       :	None.
;
; Common      :	None.
;
; Restrictions:	The input parameters must be scalar.
;
; Side effects:	Different algorithms are used for even and odd values of NPTS.
;		Above NPTS=70, only the odd algorithm is used, due to problems
;		with the even algorithm for large numbers.  When NPTS is even
;		and >70, then NPTS-1 is substituted.  This will result in a
;		slight overestimate of the probability.
;
; Prev. Hist. :	Adapted from P. R. Bevington, 1969, "Data Reduction and Error
;		Analysis for the Physical Sciences".
;
; History     :	Version 1, 10-Sep-1998, W. Thompson, GSFC
;		Version 2, 11-Apr-2000, W. Thompson, GSFC
;			Added P1,P2 option, and keyword CORRELATION
;		Version 3, 17-Nov-2000, W. Thompson, GSFC
;			Corrected error with large numbers of points.
; Contact     :	WTHOMPSON
;-
;
;
	ON_ERROR, 2
;
;  Check the input parameters.
;
	IF N_PARAMS() NE 2 THEN BEGIN
	    MESSAGE = 'Syntax:  Result=PCORRE(RR,NPTS) or PCORRE(P1,P2)'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  Check to see if the alternate P1,P2 calling sequence was used.
;
	NPTS = N_ELEMENTS(P2)
	IF NPTS GT 1 THEN BEGIN
	    IF N_ELEMENTS(P1) NE NPTS THEN BEGIN
		MESSAGE = 'The two arrays must be the same size'
		GOTO, HANDLE_ERROR
	    ENDIF
	    RR = CORRELATE(P1(*),P2(*))
	END ELSE BEGIN
	    RR   = P1
	    NPTS = P2
	ENDELSE
	CORRELATION = RR
;
	IF N_ELEMENTS(RR) NE 1 THEN BEGIN
	    MESSAGE = 'The correlation coefficient R must be scalar'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  Evaluate the number of degrees of freedom.  If NFREE is even and greater
;  than 70, then subtract 1.  The algorithm for even numbers cannot calculate
;  beyond 70 points, but the algorithm for odd numbers works for all values.
;  Subtracting 1 will overestimate the probability slightly.  If there aren't
;  enough points, then signal an error.
;
	NFREE = LONG(NPTS) - 2
	IF ((NFREE/2.) EQ (NFREE/2)) AND (NFREE GT 70) THEN NFREE = NFREE-1
	IF NFREE LE 0 THEN BEGIN
	    MESSAGE = 'NPTS must be at least 3.'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  Convert RR to double precision.  All calculations will be done in double
;  precision.  If R is outside the range -1 to 1, then signal an error.
;
	R = DOUBLE(RR)
	R2 = R^2
	IF R2 EQ 1 THEN RETURN, 1.D0 ELSE IF R2 GT 1 THEN BEGIN
	    MESSAGE = 'R must be between -1 and 1.'
	    GOTO, HANDLE_ERROR
	ENDIF
;
;  The number of degrees of freedom is even.
;
	IF (NFREE/2.) EQ (NFREE/2) THEN BEGIN
	    IMAX = (NFREE-2)/2
	    TERM = ABS(R)
	    SUM = TERM
	    IF IMAX LT 0 THEN BEGIN
		PCORRE = 1.D0
	    END ELSE IF IMAX EQ 0 THEN BEGIN
		PCORRE = 1.D0 - TERM
	    END ELSE BEGIN
		FOR I = 1L, IMAX DO BEGIN
		    TERM = -TERM * R2 * (IMAX - I + 1.D0)/I
		    SUM = SUM + TERM / (2.D0*I + 1.D0)
		ENDFOR
		PCORRE = (2.D0/SQRT(!DPI)) * (GAMMA((NFREE+1.D0)/2.D0) / $
			GAMMA(NFREE/2.D0))
		PCORRE = 1.D0 - PCORRE * SUM
	    ENDELSE
;
;  Number of degrees of freedom is odd.
;
	END ELSE BEGIN
	    IMAX = (NFREE-3)/2
	    TERM = ABS(R) * SQRT(1.D0-R2)
	    SUM = ATAN(R2/TERM)
	    IF IMAX EQ 0 THEN BEGIN
		SUM = SUM + TERM
	    END ELSE IF IMAX GT 0 THEN BEGIN
		SUM = SUM + TERM
		FOR I=1L,IMAX DO BEGIN
		    TERM = TERM * (1.D0 - R2) * (2.D0*I) / (2.D0*I + 1.D0)
		    SUM = SUM + TERM
		ENDFOR
	    ENDIF
	    PCORRE = 1.D0 - (2.D0/!DPI) * SUM
	ENDELSE
;
;  Return the result.
;
	IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = ''
	RETURN, PCORRE
;
;  Error handling point.
;
HANDLE_ERROR:
	IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = 'PCORRE: ' + MESSAGE $
		ELSE MESSAGE, MESSAGE, /CONTINUE
	RETURN, -1.D0
	END