pro fisher_r,rho,n,r,prob_r,plot=plot,oplot=oplot,check=check,$
             gaussian=gaussian,dr0=dr0,mark=mark,$
             student0=student0,dxc0=dxc0

if(n_params() le 0) then begin
    print,'-----------------------------------------'
    print,'fisher_r, rho,n,  r,prob_r'
    print,'-----------------------------------------'
    print,'Calculate the Fisher pdf (WJ4.8) of Pearson r'
    print,'for a Bivariate Gaussian distribution of
    print,'N pairs with a given rho'
    print,'  input: rho,n
    print,'  output: r,prob_r'
    print,'NOTE: works for N <100'
    print,''
    print,'/plot  ->  plot distribution'
    print,'/oplot -> oplot distribution'
    print,' check=ntrial -> check the distribution by creating'
    print,'               ntrial bivariate gaussians with given rho'
    print,'               and sampling the distribution of r'
    print,''
    print,' gaussian -> overplot Gaussian approximation for large N'
    print,' mark= value of cumulative pdf'
    print,' student=value --> in checking use t-distribution with df=value
    print,'                   SLOW!
    print,' dxc = tabulation interval in check (def=0.02)'
    print,'-----------------------------------------'
    return
endif


student=0
if(keyword_set(student0)) then student=student0


r=dindgen(1001)/1000.*2.-1.
if(keyword_set(dr0)) then begin
    nr=2./dr0
    r=dindgen(nr)/nr*2.-1.
endif
apu=1.d0-rho^2

prob_r=(1.d0-rho^2)^(N/2.d0-0.5d0)*(1-r^2)^(N/2.d0-2.d0)/$
  (1.d0-rho*r)^(N-1.5d0)*(1.d0+1.d0/(N-0.5d0)*(1.d0+r*rho)/8.d0)

;normalize
prob_r=prob_r/mean(prob_r)*.5



;---------------------------------------------
;plot

if(keyword_set(plot)) then begin
    nwin
    plot,r,prob_r,$
      xtitle='r',ytitle='Fisher Prob(r)',$
      title='rho='+string(rho,'(f8.4)')+'  N='+string(n,'(i6)')
endif

if(keyword_set(oplot)) then begin
    oplot,r,prob_r     
endif

;check the distribution
if(keyword_set(check)) then begin
    ntrial=check

    rtab=fltarr(ntrial)
    for i=0l,ntrial-1 do begin
        bivariate_gaussian_f,x,y,n,rho,sx=sx,sy=sy,student=student
        rtab(i)=correlate(x,y)
    endfor
    
    dxc=0.01
    if(keyword_set(dxc0)) then dxc=dxc0
    histo_f,rtab,-1.,1.,dxc,xt,yt,/noscale
    prob_r_num=1.*yt/ntrial/dxc    
    oplot,xt,prob_r_num,psym=10,col=2
endif


if(keyword_set(gaussian)) then begin
    sigma=(1.-rho^2)/sqrt(N-1)
    oplot,r,1./sqrt(2.*!pi)/sigma*exp(-0.5*(r-rho)^2/sigma^2),col=2
endif


;-------------------------------------------
;mark regions corresponding to tails?
;denote F=cumulative probability function
; mark > 0.5 -> mark region where F>mark
; mark < 0.5 -> mark region where F<mark

if(keyword_set(mark)) then begin
    f=prob_r*0.
    for i=1l,n_elements(r)-1 do begin
        f(i)=f(i-1)+prob_r(i)*(r(i)-r(i-1))
    endfor

    if(mark gt 0.5) then begin
        ind=where(f ge mark)
        print,'probability of r <'+string(r(ind(0)))+' = '+string(mark)
        xapu=[r(ind),reverse(r(ind))]
        yapu=[prob_r(ind)*0,reverse(prob_r(ind))]
        polyfill,xapu,yapu
    endif

    if(mark lt 0.5) then begin
        ind=reverse(where(f le mark))
        print,'probability of r >'+string(r(ind(0)))+' = '+string(mark)
        xapu=[r(ind),reverse(r(ind))]
        yapu=[prob_r(ind)*0,reverse(prob_r(ind))]
        polyfill,xapu,yapu

    endif
endif


end