C @(#)prf.for	19.1 (ES0-DMD) 02/25/03 14:29:02
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
CC
CC   *MOVE* : MOVES A CHARACTER FROM LEFT TO RIGHT
CC              (USED FOR ADJUSTING VARIABLE LABELS TO THE RIGHT)
CC
       SUBROUTINE MOVE(KARAK)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       CHARACTER KARAK(10)
       NB=0
       DO 10 J=1,10
       JPL=11-J
       IF (KARAK(JPL).EQ.' ') NB=NB+1
C       write(6,*) 'JPL, KARAK,NB:',JPL,':',KARAK(JPL),':',NB
       IF((JPL.NE.1.AND.KARAK(JPL-1).NE.' ').OR.(NB.EQ.10))GOTO 20
 10    CONTINUE
 20    CONTINUE
C       write(6,*) 'MOVE/NB = ',NB
       IF(NB.EQ.10.OR.NB.EQ.0) RETURN
       NBM=10-NB
       DO 30 J=1,NBM
       JPB=NBM+1-J
       JPL=11-J
       KARAK(JPL)=KARAK(JPB)
       write(6,*) 'JPL, JPB, KARAK(JPL), KARAK(JPB), J',
     +             JPL, JPB, KARAK(JPL), KARAK(JPB), J
       KARAK(JPB)=' '
 30    CONTINUE
       RETURN
       END
CC
CC    *PVAL* : THIS IS A FORTRAN VERSION OF THE LACKRITZ ALGORITHM
CC             (AMERICAN STATISTICIAN 1984, VOL. 38, P.312-314)
CC             FOR CALCULATING THE ONE-SIDED TAIL AREA OF AN F-DISTRIBUTION
CC             WITH N1 AND N2 DEGREES OF FREEDOM. IMPLEMENTED BY
CC             J. VAN SOEST AND B. VAN ZOMEREN, JUNE 1986.
CC
       FUNCTION PVAL(F,N1,N2)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       EXTERNAL SUM
       DOUBLE PRECISION F,PI
       DOUBLE PRECISION PVAL,SUM,T,TINV,T1,TINV1,DSUM,T5,A,B
       LOGICAL ODD
       DATA PI/3.14159265/
       IF(F.GT.(0.000001))GOTO 5
       PVAL=1.0
       RETURN
  5    T=N1*F/N2
       TINV=1/T
       T1=T+1
       TINV1=TINV+1
       IF (.NOT.ODD(N1) .AND. .NOT. ODD(N2)) THEN
          M1=N1/2
          M2=N2/2
          PVAL=SUM(M1-1,DBLE(M2-1),T/T1)/DEXP(DBLE(M2)*DLOG(T1))
       ENDIF
       IF (.NOT.ODD(N1).AND.ODD(N2)) THEN
           M1=N1/2
           M2=(N2-1)/2
           A=SUM(M1-1,DBLE(M2-0.5),T/T1)
           PVAL=A/DEXP(DBLE(M2+0.5)*DLOG(T1))
       ENDIF
       IF (ODD(N1).AND..NOT.ODD(N2)) THEN
           M1=(N1-1)/2
           M2=N2/2
           A=SUM(M2-1,DBLE(M1-0.5),TINV/TINV1)
           PVAL=1-A/DEXP(DBLE(M1+0.5)*DLOG(TINV1))
       ENDIF
       IF (ODD(N1).AND.ODD(N2)) THEN
           M1=(N1-1)/2
           M2=(N2-1)/2
           A=0.0
           IF (M1.GT.0) THEN
              DSUM=2*(DSQRT(T))/(PI*T1)
              DO 10 I=1,M2
                 DSUM=DSUM*DBLE(I)/(T1*DBLE(I-0.5))
 10                 CONTINUE
              A=DSUM
              DO 20 I=2,M1
                 DSUM=DSUM*T*2*DBLE(M2+I-1)/(T1*DBLE(2*I-1))
                 A=A+DSUM
 20                 CONTINUE
           ENDIF
           B=0.0
           IF (M2.GT.0) THEN
              DSUM=DSQRT(TINV)*2/(PI*TINV1)
              B=DSUM
              T5=2*TINV/TINV1
              DO 30 I=2,M2
              DSUM=DSUM*T5*DBLE(I-1)/DBLE(2*I-1)
              B=B+DSUM
 30              CONTINUE
           ENDIF
           PVAL=A-B+(2/PI)*DATAN(DSQRT(TINV))
       ENDIF
       RETURN
       END
CC
CC   *PTVAL* : FUNCTION (BASED ON THE FUNCTION PVAL) FOR CALCULATING
CC             THE TWO-SIDED TAIL AREA OF A STUDENT DISTRIBUTION
CC             WITH N DEGREES OF FREEDOM
CC
       FUNCTION PTVAL(TVAL,N)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       DOUBLE PRECISION PTVAL
       DOUBLE PRECISION PVAL
       DOUBLE PRECISION TVAL
       PTVAL=PVAL(TVAL*TVAL,1,N)
       RETURN
       END
CC
CC   *SUM* : AUXILIARY FUNCTION (FOR FUNCTION PVAL)
CC
       FUNCTION SUM(N,F1,F2)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       DOUBLE PRECISION SUM
       DOUBLE PRECISION F1,F2,DSUM,TSUM
       DSUM=1
       TSUM=1
       IF (N.NE.0) THEN
       DO 10 I=1,N
          DSUM=DSUM*(F1+I)*F2/I
          TSUM=TSUM+DSUM
 10    CONTINUE
       ENDIF
       SUM=TSUM
       RETURN
       END
CC
CC   *ODD* : AUXILIARY FUNCTION (FOR FUNCTION PVAL)
CC
       FUNCTION ODD(N)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       LOGICAL ODD
       ODD=.TRUE.
       IF(2*(N/2).EQ.N) ODD=.FALSE.
       RETURN
       END
CC
CC  *GRAF* : PERFORMS RESIDUAL PLOTS AND/OR A PLOT OF THE
CC           OBSERVATIONS IN THE TWO-DIMENSIONAL CASE.
CC
       SUBROUTINE GRAF(X,Y,NCAS,JPLXY,JPL2,JPLT,NMVAL,NDXY,LUB,JREG,
     1 JHEAD,NDXX,NVAD,LAB)

       INCLUDE 'MID_REL_INCL:implicit.inc'

       DIMENSION X(NCAS),Y(NCAS),MGRAF(41,51),NMVAL(NDXY)
       CHARACTER JDRAW(51),NUM(14)
       CHARACTER*10 LAB(NDXX)
       CHARACTER*60 JHEAD
       LOGICAL ODD
       NUM(1)=' '
       NUM(2)='1'
       NUM(3)='2'
       NUM(4)='3'
       NUM(5)='4'
       NUM(6)='5'
       NUM(7)='6'
       NUM(8)='7'
       NUM(9)='8'
       NUM(10)='9'
       NUM(11)='*'
       NUM(12)='-'
       NUM(13)='+'
       NUM(14)=' '
       JXSCA=51
       JYSCA=41
       AFR=10.0E-02
       DO 10 J=1,JXSCA
       DO 10 K=1,JYSCA
 10    MGRAF(K,J)=0
       IF (JPLT.EQ.2.AND.JPLXY.NE.1) X(1)=NMVAL(1)
       XMIN=X(1)
       XMAX=X(1)
       YMAX2=2.5
       YKLE2=-2.5
       YKLE=Y(1)
       YMAX=Y(1)
       DO 20 JNC=1,NCAS
       IF (JPLT.EQ.2.AND.JPLXY.NE.1) X(JNC)=NMVAL(JNC)
       IF (X(JNC).GT.XMAX) XMAX=X(JNC)
       IF (X(JNC).LT.XMIN) XMIN=X(JNC)
       IF (Y(JNC).GT.YMAX) YMAX=Y(JNC)
       IF (Y(JNC).LT.YKLE) YKLE=Y(JNC)
 20    CONTINUE
       IF((YMAX-YKLE).LE.AFR.AND.JPLXY.EQ.1)THEN
         YKLE=YKLE-1.0
         YMAX=YMAX+1.0
       ENDIF
       IF (NVAD.EQ.2.AND.JPLXY.EQ.1) THEN
       IF(0.0.GT.XMAX)XMAX=0.0
       IF(0.0.LT.XMIN)XMIN=0.0
       IF(0.0.GT.YMAX)YMAX=0.0
       IF(0.0.LT.YKLE)YKLE=0.0
       ENDIF
 30    IF (YMAX.GT.2.5.OR.JPLXY.EQ.1.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX))
     1 YMAX2=YMAX
       IF (YKLE.LT.(-2.5).OR.JPLXY.EQ.1.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX))
     1 YKLE2=YKLE
       XSCA=JXSCA
       YSCA=JYSCA
       IF((XMAX-XMIN).LE.AFR)THEN
         XMIN=XMIN-1.0
         XMAX=XMAX+1.0
       ENDIF
 40    A=(XSCA-1.0)/(XMAX-XMIN)
       B=1.0-A*XMIN
       C=(YSCA-1.0)/(YMAX2-YKLE2)
       D=1.0-C*YKLE2
       IF (JPLXY.EQ.1) GOTO 50
       JNUL=INT(YSCA-D)
       IF ((JNUL.LE.1.OR.ABS(YMAX2).LT.AFR).AND.JPL2.NE.0) THEN
          J2P=1
          JNUL=2
          J2M=3
          GOTO 50
       ENDIF
       IF ((JNUL.GE.JYSCA.OR.ABS(YKLE2).LT.AFR).AND.JPL2.NE.0) THEN
          J2P=JYSCA-2
          JNUL=JYSCA-1
          J2M=JYSCA
          GOTO 50
       ENDIF
       IF (JPL2.EQ.0) THEN
          IF (JNUL.LE.1.OR.ABS(YMAX2).LE.AFR) JNUL=1
          IF (JNUL.GT.JYSCA.OR.ABS(YKLE2).LE.AFR) JNUL=JYSCA
       ENDIF
       IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 50
       J2P=INT(YSCA-(2.5*C+D))
       IF (JNUL.EQ.J2P) J2P=JNUL-1
       IF (ABS(YMAX2-2.5).LT.AFR) J2P=1
       IF (J2P.LE.0) THEN
          IF (ABS(YMAX2-2.5).LT.AFR) THEN
             J2P=1
                                     ELSE
             J2P=2
          ENDIF
       ENDIF
       IF (J2P.GT.JYSCA) J2P=JYSCA-2
       IF (ABS(YKLE2-2.5).LT.AFR) J2P=JYSCA-2
       J2M=2*JNUL-J2P
       IF (J2M.LE.0.OR.ABS(YMAX2+2.5).LT.AFR) THEN
          J2P=1
          JNUL=2
          J2M=3
          GOTO 50
       ENDIF
       IF (J2M.GT.JYSCA.OR.ABS(YKLE2+2.5).LT.AFR) THEN
          J2M=JYSCA
          JS=J2M+J2P
          IF (ODD(JS)) THEN
             J2P=J2P+1
          ENDIF
          JNUL=(J2M+J2P)/2
       ENDIF
 50    DO 110 JNC=1,NCAS
       XJNC=X(JNC)
       YJNC=Y(JNC)
 60    JX=INT(XJNC*A+B)
       IF (JX.LE.0) JX=1
       IF (JX.GT.JXSCA) JX=JXSCA
       IF(.NOT.(ABS(YJNC-YMAX2).LT.AFR.OR.ABS(YJNC-YKLE2).LT.AFR))
     1 GOTO 70
       IF (ABS(YJNC-YMAX2).LT.AFR) JEY=1
       IF (ABS(YJNC-YKLE2).LT.AFR) JEY=JYSCA
       IF (ABS(YJNC).LT.AFR.AND.JPLXY.EQ.0) JEY=JNUL
       GOTO 90
 70    JEY=INT(YSCA-(YJNC*C+D))
       IF (JPLXY.NE.1) THEN
         IF (ABS(YJNC).LT.AFR) JEY=JNUL
         IF (JEY.EQ.J2P.AND.YJNC.LT.(2.5-AFR)) JEY=J2P+1
         IF (JEY.EQ.J2P.AND.YJNC.GT.(2.5+AFR)) JEY=J2P-1
         IF (JEY.EQ.J2M.AND.YJNC.LT.(-2.5-AFR)) JEY=J2M+1
         IF (JEY.EQ.J2M.AND.YJNC.GT.(-2.5+AFR)) JEY=J2M-1
         IF (ABS(YJNC-2.5).LT.AFR.OR.(JEY.GT.J2P.AND.YJNC.GT.2.5)
     1            .OR.(JEY.LT.J2P.AND.YJNC.LT.2.5)) JEY=J2P
         IF (ABS(YJNC+2.5).LT.AFR.OR.(JEY.GT.J2M.AND.YJNC.GT.-2.5)
     1         .OR.(JEY.LT.J2M.AND.YJNC.LT.-2.5)) JEY=J2M
       ENDIF
 80    IF (JEY.LE.0) JEY=1
       IF (JEY.GT.JYSCA) JEY=JYSCA
 90    MGRAF(JEY,JX)=MGRAF(JEY,JX)+1
       IF(NUM(14).NE.'0') GOTO 100
       MGRAF(JEY,JX)=-20
       GOTO 120
 100   IF (.NOT.(JNC.EQ.NCAS.AND.JPLXY.EQ.1.AND.NVAD.EQ.2)) GOTO 110
       XJNC=0.0
       YJNC=0.0
       NUM(14)='0'
       GOTO 60
 110   CONTINUE
 120   IF (JREG.EQ.0) WRITE(LUB,8000) JHEAD
       IF (JREG.EQ.1) WRITE(LUB,8010) JHEAD
       IF (JREG.EQ.2) WRITE(LUB,8020) JHEAD
       IF (JREG.EQ.3) WRITE(LUB,8030) JHEAD
       IF (JPLXY.EQ.0.AND.JPL2.NE.0) WRITE(LUB,8050)
       IF (JPLXY.EQ.0.AND.JPL2.EQ.0) WRITE(LUB,8060)
       IF (JPLXY.EQ.1) WRITE(LUB,8070) LAB(NVAD)
       WRITE(LUB,8075)
       DO 240 J=1,JYSCA
       IF (JPLXY.EQ.1) GOTO 160
       IF(J.NE.JNUL) GOTO 140
       DO 130 K=1,JXSCA
       MJK=MGRAF(J,K)+1
       IF (MJK.EQ.1) JDRAW(K)=NUM(12)
       IF (MJK.GT.10) JDRAW(K)=NUM(11)
       IF (MJK.GE.2.AND.MJK.LE.10) JDRAW(K)=NUM(MJK)
 130   CONTINUE
       GOTO 180
 140   IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 160
       IF (J.NE.J2P.AND.J.NE.J2M) GOTO 160
       DO 150 K=1,JXSCA
       MJK=MGRAF(J,K)+1
       IF (MJK.EQ.1) JDRAW(K)=NUM(13)
       IF (MJK.GT.10) JDRAW(K)=NUM(11)
       IF (MJK.GE.2.AND.MJK.LE.10) JDRAW(K)=NUM(MJK)
 150   CONTINUE
       GOTO 180
 160   DO 170 K=1,JXSCA
       MJK=MGRAF(J,K)+1
       IF(MJK.GT.10) JDRAW(K)=NUM(11)
       IF(MJK.LE.10) JDRAW(K)=NUM(MJK)
       IF(MJK.LT.0) JDRAW(K)=NUM(14)
 170   CONTINUE
 180   JREST=MOD(J,5)
       IF (JPLXY.EQ.1) GOTO 190
       IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 190
       IF (J.NE.J2P) GOTO 190
       IF (JREST.EQ.1) WRITE(LUB,8080)(JDRAW(K),K=1,JXSCA)
       IF (JREST.NE.1) WRITE(LUB,8085)(JDRAW(K),K=1,JXSCA)
       GOTO 240
 190   IF (J.NE.1) GOTO 200
       WRITE(LUB,8130) YMAX2,(JDRAW(K),K=1,JXSCA)
       GOTO 240
 200   IF (JPLXY.EQ.1) GOTO 220
       IF (J.NE.JNUL) GOTO 210
       IF (JREST.EQ.1) WRITE(LUB,8090)(JDRAW(K),K=1,JXSCA)
       IF (JREST.NE.1) WRITE(LUB,8100)(JDRAW(K),K=1,JXSCA)
       GOTO 240
 210   IF (J.NE.J2M.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 220
       IF (JREST.EQ.1) WRITE(LUB,8110) (JDRAW(K),K=1,JXSCA)
       IF (JREST.NE.1) WRITE(LUB,8120) (JDRAW(K),K=1,JXSCA)
       GOTO 240
 220   IF(J.NE.JYSCA) GOTO 230
       WRITE(LUB,8130) YKLE2,(JDRAW(K),K=1,JXSCA)
       GOTO 240
 230   IF (JREST.EQ.1) WRITE(LUB,8140)(JDRAW(K),K=1,JXSCA)
       IF (JREST.NE.1) WRITE(LUB,8150)(JDRAW(K),K=1,JXSCA)
 240   CONTINUE
       WRITE(LUB,8075)
       IF (JPLT.NE.2.OR.JPLXY.EQ.1) WRITE(LUB,8160) XMIN,XMAX
       IF ((JPLT.EQ.1.OR.JPLT.EQ.3).AND.JPLXY.EQ.0)
     1 WRITE(LUB,8170) LAB(NVAD)
       IF (JPLXY.EQ.1) WRITE(LUB,8180) LAB(1)
       IF (JPLT.NE.2.OR.JPLXY.EQ.1) GOTO 250
       J1=XMIN
       JN=XMAX
       WRITE(LUB,8190) J1,JN
       WRITE(LUB,8200)
 8000  FORMAT(//16X,A60)
 8010  FORMAT(//16X,A60//29X,'--- L E A S T  S Q U A R E S ---')
 8020  FORMAT(//16X,A60//20X,'--- L E A S T  M E D I A N  O F ',
     1 ' S Q U A R E S ---')
 8030  FORMAT(//16X,A60//13X,'--- R E W E I G H T E D',3X,'L S',6X,
     1 '( B A S E D',3X,'O N',3X,'L M S ) ---')
 8050  FORMAT(/'  STAND. RESIDUAL',1X,'I-+',10('----+'),'-I')
 8060  FORMAT(/7X,'RESIDUAL',3X,'I-+',10('----+'),'-I')
 8070  FORMAT(/6X,'OBSERVED'/6X,A10,2X,'I-+',10('----+'),'-I')
 8075  FORMAT(18X,'I',53X,'I')
 8080  FORMAT(14X,'2.5',1X,'++',51A1,'++')
 8085  FORMAT(14X,'2.5',1X,'I+',51A1,'+I')
 8090  FORMAT(14X,'0.0',1X,'+-',51A1,'-+')
 8100  FORMAT(14X,'0.0',1X,'I-',51A1,'-I')
 8110  FORMAT(13X,'-2.5',1X,'++',51A1,'++')
 8120  FORMAT(13X,'-2.5',1X,'I+',51A1,'+I')
 8130  FORMAT(6X,E11.4,1X,'+',1X,51A1,1X,'+')
 8140  FORMAT(18X,'+',1X,51A1,1X,'+')
 8150  FORMAT(18X,'I',1X,51A1,1X,'I')
 8160  FORMAT(18X,'I-+',10('----+'),'-I'/10X,E11.4,46X,E11.4)
 8170  FORMAT(/49X,'ESTIMATED ',A10//)
 8180  FORMAT(/54X,'OBSERVED ',A10//)
 8190  FORMAT(18X,'I-+',10('----+'),'-I'/17X,I4,46X,I4)
 8200  FORMAT(/48X,' INDEX OF THE OBSERVATION'//)
 250   RETURN
       END
