C @(#)prf.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:52 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