;+
; Contains the correl procedure
;
; :Author:
; Laurent Lamy or Philippe Zarka
;
; :History:
; 2009/11/24: Created
;
; 2009/11/24: Last Edit
;-
;+
; Calcule le coefficient de correlation C et son degre
; de confiance P de 2 tableaux X(N) et Y(N).
; $$C = \frac{\sigma_{XY}}{\sigma_{X} * \sigma_{Y}}$$
; avec $$\sigma_X^{2} = \frac{1}{N} \sum_{i=1}^{N} (X_{i}-X_{m})^{2}$$
; $$\sigma_{Y}^{2} = \frac{1}{N} \sum_{i=1}^{N} (Y_{i}-Y_{m})^{2}$$
; $$\sigma_{XY} = \frac{1}{N} \sum_{i=1}^{N} (X_{i}-X_{m})*(Y_{i}-Y_{m})$$
; et $$P = 100 * \left|\frac{2}{\sqrt{\pi}} \frac{\Gamma(\frac{N-1}{2})}{\Gamma(\frac{N-2}{2})} \int_0^C (1-u^{2})^{\frac{N-4}{2}} \mathrm{d}u \right|$$
;
; Pour le calcul de P, on opere de la maniere suivante:
;
; 1) $$\frac{\Gamma(\frac{N-1}{2})}{\Gamma(\frac{N-2}{2})}$$ est calculé avec la fonction
; LNGAMMA (erreur < 1%)
;
; 2) $$\int_0^C (1-u^{2})^{\frac{N-4}{2}} \mathrm{d}u $$ est calculée par recurrence:
; si $$I_{k} = \int_0^C (1-u^{2})^{k} \mathrm{d}u $$
; $$ I_{k} = \frac{1}{2k+1} C(1-C^{2})^{k} + \frac{2k}{2k+1} I_{k-1} $$
; $$I_{0} = C$$
; $$I_{1/2} = \frac{1}{2} \arcsin C + \frac{C}{2} \sqrt{1-C^{2}}$$
;
; Le calcul de P n'est pas effectue pour N > 23000.
;
; :Params:
; X: in, required, type=fltarr
; Tableau
; Y: in, required, type=fltarr
; Tableau
; C: out, required, type=float
; coefficient de correlation (-1.<=C<=1)
; P: out, required, type=float
; degre de confiance ( 0.<=P<=100.)
;-
pro CORREL, X, Y, C, P
on_error,1
C=0. & P=999.
N=n_elements(X)
Xm=float(total(X)/N) & Ym=float(total(Y)/N)
sigma_X=sqrt(total((X-Xm)^2)/N) & sigma_Y=sqrt(total((Y-Ym)^2)/N)
sigma_XY=total((X-Xm)*(Y-Ym))/N
if sigma_X eq 0. or sigma_Y eq 0. then return
C=float(sigma_XY/sigma_X/sigma_Y)
; print,n,xm,ym,sigma_x,sigma_y,sigma_xy,c
if N gt 100000 then return
R=2/sqrt(!Pi)*exp(lngamma(double(N-1)/2.)-lngamma(double(N-2)/2.))
if (N mod 2) eq 0 then begin
ii=double(C)
k=1.
endif
if (N mod 2) eq 1 then begin
ii=.5*asin(double(C))+C/2.d0*sqrt(1.d0-C*C)
k=1.5
endif
i0=1000.
while (k le (N-4.)/2.) and (abs(ii-i0)/i0 ge 1.d-5) do begin
i0=ii
ii=(C*(1-C*C)^k+2.*k*i0)/(2.*k+1.)
k=k+1.
endwhile
P=100.*abs(R*ii)
; print,'c,p,k,ii=',c,p,k,ii
return
end