C***********************************************************************
C   THIS SUBROUTINE COMPUTES THE PROBABILITY INTEGRAL TRANSFORMATION
C   FOR A PRODUCT MOMENT CORRELATION COEFFICIENT. IT IS ASSUMED THAT
C   THE CORRELATION WAS COMPUTED FROM A SAMPLE OF A BIVARIATE
C   NORMAL POPULATION WITH TRUE CORRELATION EQUAL ZERO.
C
C   CALLING ARGUMENTS:
C   N    THE SAMPLE SIZE FROM WHICH R IS COMPUTED.
C   R    THE CORRELATION COEFFICIENT.
C   P    THE PROBABILITY OF A CORRELATION VALUE LESS THAN OR EQUAL
C        TO THE COMPUTED ONE. IF THERE IS ANY PROBLEM WITH COMPUTING P
C        THEN P IS RETURNED EQUAL TO -1.
C   PSQ  THE PROBABILITY OF A CORRELATION COFFICIENT LESS THAN OR EQUAL
C        TO THE SQUARED R VALUE. ALSO RETURNED AS -1 IN CASE OF TROUBLE.
C
C   FOR SAMPLE SIZES GREATER THAN 100 R IS TRANSFORMED TO FISHERS Z:
C        Z   (102)*LN((1+R)/(1-R))
C   THE DISTRIBUTION OF Y Z*Z*(N-3) IS NEARLY CENTRAL CHI-SQUARE ON
C   ONE DEGREE OF FREEDOM.
C   RSQ HAS A BETA(1/2,(N-2)/2) DISTRIBUTION, THE DENSITY FUNCTION IS
C   B(N)*(1-RSQ)**((N-4)/2)/SQRT(RSQ).  THE ARRAY OF NORMALIZING
C   CONSTANTS B(N) HAS BEEN PRE SET FOR VALUES OF N 3,100.
C   THE DENSITY FORM HAS BEEN EXPANDED IN A SERIES THEN INTEGRATED
C   TERM WISE TO YEILD THE P. I. T. AS A SERIES. WHEN N IS EVEN THE
C   SERIES IS FINITE. FOR N ODD THE SERIES IS INFINITE AND ONLY FINITE
C   MANY TERMS ARE USED AS AN APPROXIMATION.
C   ALSO FOR RSQ VALUES VERY NEAR 1, SO THAT PSQ WOULD BE .999 OR MORE,
C   THE PROGRAM JUST SETS PSQ EQUAL TO 1.0  AND RETURNS.
C
C     SUBROUTINES CALLED: ONE
C**************************************
      SUBROUTINE COR (N,R,P,PSQ)
      DIMENSION B(100)
      DATA B /0.,0.,.318309886184,.5,0.63661977,0.75000000,0.84882636,0.
     1 93750000,1.01859164,1.09375000,1.16410473,1.23046875,1.29344970,1
     2 .35351562,1.41103603,1.46630859,1.51957727,1.57104492,1.62088242,
     3 1.66923523,1.71622844,1.76197052,1.80655625,1.85006905,1.89258274
     4 ,1.93416309,1.97486895,2.01475322,2.05386371,2.09224373,2.1299327
     5 3,2.16696672,2.20337869,2.23919895,2.27445542,2.30917391,2.343378
     6 31,2.37709079,2.41033198,2.44312109,2.47547609,2.50741375,2.53894
     7 983,2.57009910,2.60087544,2.63129193,2.66136091,2.69109402,2.7205
     8 0227,2.74959607,2.77838529,2.80687932,2.83508703,2.86301690,2.890
     9 67698,2.91807492,2.94521805,2.97211334,2.99876747,3.02518680,3.05
     $ 137742,3.07734519,3.10309569,3.12863428,3.15396611,3.17909612,3.2
     $ 0402906,3.22876950,3.25332182,3.27769025,3.30187886,3.32589157,3.
     $ 34973218,3.37340431,3.39691150,3.42025715,3.44344454,3.46647684,3
     $ .48935713,3.51208838,3.53467346,3.55711515,3.57941616,3.60157909,
     $ 3.62360648,3.64550079,3.66726439,3.68889961,3.71040868,3.73179379
     $ ,3.75305705,3.77420053,3.79522623,3.81613610,3.83693202,3.8576158
     $ 4,3.87818935,3.89865430,3.91901240,3.93926528/
C
C   CHECK INPUT
C
      P=-1.
      PSQ=-1.
      RSQ=R*R
      IF ((N.LT.3).OR.(RSQ.GT.1.)) RETURN
C
C   CHECK IF N IS LESS THAN 100
C
      IF (N.LE.100) GO TO 10
      SIGN=+1.
      IF (R.LT.0.) SIGN=-1.
      Z=.5*ALOG((1.+R)/(1.-R))
      Y=Z*Z*FLOAT(N-3)
      CALL ONE (Y,PSQ)
      P=(1.+SIGN*PSQ)*.5
      RETURN
C
C   N IS LESS THAN 100, CHECK THAT RSQ IS NOT CLOSE TO 1.
C
   10 DIV=N-3
      IF (N.EQ.3) DIV=1.
      TEST=EXP(SQRT(60./DIV))
      TEST=(TEST-1.)/(TEST+1.)
      PSQ=1.
      P=1.
      IF (R.LT.0.) P=0.
      IF (RSQ.GE.TEST) RETURN
C
C   COMPUTATION OF PSQ BEGINS HERE.
C
      A=2.*ABS(R)
      S=A
      IF (N.EQ.4) GO TO 50
      X=FLOAT(N)*.5
      M=X+.001
      X=X-2.
      TEST=(1.-RSQ)*.00000001
      IF ((2*M).NE.N) GO TO 30
C
C   N IS EVEN
C
      M=M-2
      DO 20 J=1,M
      XI=FLOAT(J-1)
      W=(RSQ*(X-XI)*(XI+.5))/((XI+1.5)*(XI+1))
      A=-A*W
      S=S+A
      IF (ABS(A).LT.TEST) GO TO 50
   20 CONTINUE
      GO TO 50
C
C   N IS ODD
C
   30 I=0
   40 XI=FLOAT(I)
      W=(RSQ*(X-XI)*(XI+.5))/((XI+1.5)*(XI+1))
      A=-A*W
      S=S+A
      IF (ABS(A).LT.TEST) GO TO 50
      I=I+1
      GO TO 40
C
C   FINISH COMPUTING PSQ AND P.
C
   50 PSQ=B(N)*S
      IF(PSQ.GT.1.)PSQ=1.
      SIGN=1.
      IF (R.LT.0.) SIGN=-1.
      P=(1.+SIGN*PSQ)*.5
      RETURN
      END