MINICERN/mathlib/gen/g/studin.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:42  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       FUNCTION STUDIN(Q,N)
  11 C    THIS ROUTINE COMPUTES THE INVERSE OF THE DISTRIBUTION
  12 C    FUNCTION OF THE STUDENT DISTRIBUTION WITH N DEGREES OF FREEDOM.
  13 C    Q=PROBABILITY AT WHICH THE FUNCTION IS INVERTED,
  14 C    N=INTEGER GREATER THAN ZERO.
  15 C     THE ROUTINE WAS WRITTEN BY G.W.HILL IN ALGOL
  16 C     C.A.C.M. ALGORITHM 396
  17 C    AT LEAST 6 SIGNIFICANT FIGURES ARE CORRECT.
  18       DATA HP/1.5707963268/
  19       IF(N. LT. 1) GO TO 10
  20       RL=1.
  21       IF(Q .GE. 0.5) GO TO 1
  22       RL=-1.
  23       P=2.*Q
  24       GO TO 2
  25     1 P=2.*(1.-Q)
  26     2 IF(P .LE. 0. .OR. P .GT. 1.) GO TO 20
  27       IF(N .GT. 1) GO TO 3
  28       PP=COS(HP*P)
  29       STUDIN=PP/SQRT(1.-PP*PP)*RL
  30       RETURN
  31     3 IF(N .GT. 2) GO TO 4
  32       STUDIN=SQRT(2./(P*(2.-P))-2.)*RL
  33       RETURN
  34     4 RN=N
  35       A=1./(RN-0.5)
  36       B=48./(A*A)
  37       C=((20700.*A/B-98.)*A-16.)*A+96.36
  38       D=((94.5/(B+C)-3.)/B+1.)*SQRT(A*HP)*RN
  39       X=D*P
  40       Y=X**(2./RN)
  41       IF(Y .LE. 0.05+A) GO TO 5
  42       PP=0.5*P
  43       X=GAUSIN(PP)
  44       Y=X*X
  45       IF(N .GE. 5) C=C+0.3*(RN-4.5)*(X+0.6)
  46       C=(((0.05*D*X-5.)*X-7.)*X-2.)*X+B+C
  47       Y=(((((0.4*Y+6.3)*Y+36.)*Y+94.5)/C-Y-3.)/B+1.)*X
  48       Y=A*Y*Y
  49       IF(Y .LE. 0.002) Y=0.5*Y*Y+Y
  50       IF(Y .GT. 0.002) Y=EXP(Y)-1.
  51       GO TO 6
  52     5 Y=((1./(((RN+6.)/(RN*Y)-0.089*D-0.822)*(RN+2.)*3.)+0.5/(RN+4.))*Y-
  53      11.)*(RN+1.)/(RN+2.)+1./Y
  54     6 STUDIN=SQRT(RN*Y)*RL
  55       RETURN
  56    10 WRITE(6,7) N
  57    20 WRITE(6,8) Q
  58       STOP
  59     7 FORMAT(/10X,'DEGREE OF FREEDOM N=',I5,' IN STUDIN ILLEGAL'/)
  60     8 FORMAT(/10X,'ARGUMENT  Q=',E15.5,' IN STUDIN ILLEGAL'/)
  61       END