C @(#)pre.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:51 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 *LSL* : CALCULATES THE (REWEIGHTED) LS WHEN 'NVAR' EQUALS 1 CC AND 'JCST'=0. CC SUBROUTINE LSL(NCAS,NDXY,X,Y,RESDU,A,FCKW,H,NMXV,NDXX) INCLUDE 'MID_REL_INCL:implicit.inc' DOUBLE PRECISION H(NMXV,NDXX) DIMENSION A(NMXV),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY) AL=NCAS SXY=0.0 SX2=0.0 DO 10 JNC=1,NCAS SXY=SXY+X(1,JNC)*Y(JNC)*RESDU(JNC) SX2=SX2+(X(1,JNC)**2)*RESDU(JNC) 10 CONTINUE A(1)=SXY/SX2 FCKW=0.0 DO 20 JNC=1,NCAS 20 FCKW=FCKW+( (Y(JNC)-X(1,JNC)*A(1))**2 )*RESDU(JNC) H(1,1)=DBLE((FCKW/(AL-1.0))/SX2) RETURN END CC CC *FCN* : PUTS A ROW OF THE MATRIX X IN A VECTOR. CC SUBROUTINE FCN(K,M,F,X,JFLAG) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION F(K),X(M) DO 10 J=1,K F(J)=X(J) 10 CONTINUE RETURN END CC CC *LSREG* : CALCULATES THE LEAST SQUARES REGRESSION ESTIMATES. CC SUBROUTINE LSREG(NDXX,NDXY,NMXV,K,N,M,F,X,Y,W,DA, 1 H,FCKW,HVEC,JDMB,JNDEX) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),F(K),Y(N),W(N),H(NMXV,NDXX),DA(K) DIMENSION HVEC(JDMB),JNDEX(NDXX) DOUBLE PRECISION H,HVEC,DTEM,DFCKW,DFACT DOUBLE PRECISION DWJNC,DYJ,DFKA KPLUS=K+1 DO 10 JNC=1,K DO 20 J=1,KPLUS H(JNC,J) = 0.D0 20 CONTINUE 10 CONTINUE JFLAG=1 ANUL=0.0 DO 30 JNC=1,N CALL FCN(K,M,F,X(1,JNC),JFLAG) JFLAG=4 DWJNC = DBLE(W(JNC)) ANUL=ANUL+W(JNC) DYJ = DBLE(Y(JNC)) DO 40 KA=1,K DFKA = DBLE(F(KA)) H(KA,K+1) = H(KA,K+1) + DWJNC*DFKA*DYJ DO 50 L=1,KA H(KA,L) = H(KA,L) + DWJNC*DFKA*DBLE(F(L)) 50 CONTINUE 40 CONTINUE 30 CONTINUE DO 60 J=1,K DO 70 JNC=1,J H(JNC,J)=H(J,JNC) 70 CONTINUE 60 CONTINUE CALL MATNV(H,NMXV,NDXX,HVEC,JDMB,K,1,NERR,DTEM,JNDEX) MM=K+1 FCKW = QLSRG(K,N,M,NDXX,NDXY,NMXV,F,X,Y,W,DA,H,MM) DO 80 JNC=1,K F(JNC)=H(JNC,K+1) 80 CONTINUE DFCKW=DBLE(FCKW) ANK=ANUL-K DFACT=DBLE(ANK) DFACT=DFCKW/DFACT DO 90 JNC=1,K DO 100 J=1,K H(JNC,J)=H(JNC,J)*DFACT 100 CONTINUE 90 CONTINUE DO 110 JNC=1,K HDA=H(JNC,JNC) DA(JNC)=SQRT(HDA) 110 CONTINUE RETURN END CC CC *QLSRG* : EVALUATES THE OBJECTIVE FUNCTION FOR LS REGRESSION. CC FUNCTION QLSRG(K,N,M,NDXX,NDXY,NMXV,F,X,Y,W,DA,H,MM) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION F(K),X(NDXX,NDXY),DA(K),H(NMXV,NDXX),Y(N),W(N) DOUBLE PRECISION H,Q,HSUM JFLAG = 4 Q=0. DO 30 JNC=1,N CALL FCN(K,M,F,X(1,JNC),JFLAG) HSUM=0. DO 20 JNCB=1,K 20 HSUM=H(JNCB,MM)*F(JNCB)+HSUM 30 Q=(HSUM-Y(JNC))*(HSUM-Y(JNC))*W(JNC)+Q QLSRG = Q RETURN END CC CC *MATNV* : PERFORMS A MATRIX INVERSION. CC SUBROUTINE MATNV(AM,NMXV,NDXX,HVEC,JDMB,NA,NB,NERR,DTEM,JNDEX) INCLUDE 'MID_REL_INCL:implicit.inc' DOUBLE PRECISION AM,HVEC,DTEM,DETER,TURN,SWAP DIMENSION HVEC(JDMB),JNDEX(NDXX),AM(NMXV,NDXX) DETER=1.0D0 N=NA NPNB=N+NB JNK=0 DO 10 J=1,NPNB JNK=(J-1)*NMXV DO 10 NC=1,NMXV JNK=JNK+1 HVEC(JNK)=AM(NC,J) 10 CONTINUE JDM=NMXV NMA=N-1 JDELC=1-JDM DO 130 JHFD=1,N TURN=0.0D0 JDELC=JDELC+JDM JDLA=JDELC+JHFD-1 JDLB=JDELC +NMA DO 40 JNCB=JDLA,JDLB IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30 30 TURN=HVEC(JNCB) LDEL=JNCB 40 CONTINUE IF(TURN) 50,170,50 50 JPAAL=LDEL-JDELC+1 JNDEX(JHFD)=JPAAL IF(JPAAL-JHFD) 80,80,60 60 DETER=-DETER JPAAL=JPAAL-JDM JNCD=JHFD-JDM DO 70 JNC=1,NPNB JPAAL=JPAAL+JDM JNCD=JNCD+JDM SWAP=HVEC(JNCD) HVEC(JNCD)=HVEC(JPAAL) 70 HVEC(JPAAL)=SWAP 80 DETER=DETER*TURN TURN=1./TURN JNCD=JDELC+NMA DO 90 JNC=JDELC,JNCD 90 HVEC(JNC)=-HVEC(JNC)*TURN HVEC(JDLA)=TURN JNCB=JHFD-JDM JPAAL=1-JDM DO 120 JNC=1,NPNB JPAAL=JPAAL+JDM JNCB=JNCB+JDM IF(JNC-JHFD) 100,120,100 100 JCL=JPAAL+NMA SWAP=HVEC(JNCB) JNCD=JDELC-1 DO 110 JNCC=JPAAL,JCL JNCD=JNCD+1 110 HVEC(JNCC)=HVEC(JNCC)+SWAP*HVEC(JNCD) HVEC(JNCB)=SWAP*TURN 120 CONTINUE 130 CONTINUE DO 160 JNCB=1,N JHFD=N+1-JNCB LDEL=JNDEX(JHFD) IF(LDEL-JHFD) 140,160,140 140 JPAAL=(LDEL-1)*JDM+1 JCL=JPAAL+NMA JDELC=(JHFD-1)*JDM+1-JPAAL DO 150 JNCC=JPAAL,JCL JNCD=JNCC+JDELC SWAP=HVEC(JNCC) HVEC(JNCC)=HVEC(JNCD) 150 HVEC(JNCD)=SWAP 160 CONTINUE DTEM=DETER NERR=0 GOTO 180 170 NERR=JHFD DTEM=DETER 180 JNK=0 DO 190 J=1,NPNB DO 190 NC=1,NMXV JNK=JNK+1 AM(NC,J)=HVEC(JNK) 190 CONTINUE RETURN END CC CC *LCAT* : CALCULATES LOCATION ESTIMATES. CC SUBROUTINE LCAT(NCAS,NVAR,JCST,JPRT,NVAD,X,Y,RESDU,WEIGHTS,PREC, 1 XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,NDXX,NDXY,MVAL,LUA,LUB,LUC, 1 JREG,JHEAD,FNAMEA,FNAMEB,FNAMEC,YNSAVE,LAB,JFMT,JVARS,YN,JPLACE) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION A(11),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY) DIMENSION XMED(NDXX),XMAD(NDXX),AW(NDXY),NMVAL(NDXY),JPLACE(NDXX) CHARACTER YN,YNSAVE CHARACTER*10 LAB(NDXX) CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*60 JFMT,JHEAD CHARACTER*80 CBUF INTEGER ISTAT LOGICAL NULL INTEGER STATUS 10 WRITE(LUB,8000) NCAS WRITE(LUB,8010) JHEAD IF (JPRT.EQ.0) THEN WRITE(LUB,8020) ELSEIF (JPRT.EQ.1) THEN WRITE(LUB,8030) ELSEIF (JPRT.EQ.2) THEN WRITE(LUB,8040) ENDIF IF (JPLT.NE.0) WRITE(LUB,8050) DO 5 J=1,11 5 A(J)=0.0 IF (NCAS.LE.2) THEN WRITE(CBUF,8070) NCAS CALL STTPUT(CBUF,ISTAT) STOP ENDIF IF (FNAMEA.EQ.'CON') THEN WRITE(CBUF,8060) CALL STTPUT(CBUF,ISTAT) ENDIF DO 40 JNC=1,NCAS NMVAL(JNC)=JNC IF (FNAMEA.EQ.'CON'.AND.JNC.EQ.1) THEN WRITE(CBUF,8080) JNC CALL STTPUT(CBUF,ISTAT) ENDIF IF (FNAMEA.EQ.'CON'.AND.JNC.NE.1) THEN WRITE(CBUF,8090) JNC CALL STTPUT(CBUF,ISTAT) ENDIF IF (YN.NE.'Y') GOTO 30 C READ(LUA,*) (AW(J),J=1,JVARS) JH=JPLACE(NVAD) CALL TBERDR(LUA,JNC,JH,AW(JH),NULL,STATUS) CALL STTPUT('Module PRE.F/LCAT: Beware the selection flag', 1 ISTAT) Y(JNC)=AW(JH) IF (YNSAVE.EQ.'Y') WRITE(LUC,*) Y(JNC) GOTO 40 30 READ(LUA,JFMT) (AW(J),J=1,JVARS) JH=JPLACE(NVAD) Y(JNC)=AW(JH) IF (YNSAVE.EQ.'Y') WRITE(LUC,JFMT) Y(JNC) 40 X(1,JNC)=Y(JNC) IF (YNSAVE.EQ.'Y') CLOSE(LUC,STATUS='KEEP') IF (MVAL.NE.0) CALL SMISLOC(NCAS,NDXX,NDXY,X,Y,NMVAL,NSTOP) IF (NCAS.LE.2) THEN WRITE(CBUF,8070) NCAS CALL STTPUT(CBUF,ISTAT) STOP ENDIF IF (NSTOP.EQ.1) STOP AL=NCAS ANVAR=NVAR JQU=INT(AL/2.0)+INT((ANVAR+1.0)/2.0) IF (JPRT.EQ.0) GOTO 60 WRITE(LUB,8100) DO 50 JNC=1,NCAS 50 WRITE(LUB,8110) JNC,Y(JNC) 60 XMED(1)=AMDAN(AW,NDXY,Y,NCAS) DO 70 JNC=1,NCAS RESDU(JNC)=ABS(Y(JNC)-XMED(1)) 70 A(1)=A(1)+Y(JNC) XMAD(1)=AMDAN(AW,NDXY,RESDU,NCAS) XMAD(2)=XMAD(1)*1.4826 WRITE(LUB,8120) XMED(1),XMAD(2) IF (ABS(XMAD(1)).LE.1.0E-12) THEN WRITE(LUB,8125) STOP ENDIF WRITE(LUB,8130) A(1)=A(1)/AL DO 80 JNC=1,NCAS 80 A(2)=A(2)+(Y(JNC)-A(1))**2 A(2)=SQRT(A(2)/(AL-1.0)) WRITE(LUB,8140) A(1),A(2) JREG=1 CALL RDUAL(A,0,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL, 1 NDXX,NDXY,JHEAD,LAB) 90 WRITE(LUB,8130) CALL RANGS(Y,NCAS) CALL SHHLF(Y,NCAS,JQU,SLUTN,BSTD,PREC) BSTD=BSTD*1.4826*(5.0/(AL-ANVAR)+1.0) WRITE(LUB,8150) SLUTN,BSTD A(1)=SLUTN A(2)=BSTD CC-----DETERMINATION OF THE WEIGHTS NZWE=0 DO 100 JNC=1,NCAS Y(JNC)=X(1,JNC) RESDU(JNC)=(Y(JNC)-A(1))/A(2) IF (ABS(RESDU(JNC)).LT.2.5) THEN WEIGHTS(JNC)=1.0 NZWE=NZWE+1 ELSE WEIGHTS(JNC)=0.0 ENDIF write(CBUF,*) jnc,y(jnc),resdu(jnc),nzwe,weights(jnc) CALL STTPUT(CBUF,ISTAT) 100 CONTINUE AVW=(NZWE*1.0)/AL AL=NZWE JREG=2 CALL RDUAL(A,1,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL, 1 NDXX,NDXY,JHEAD,LAB) WRITE(LUB,8130) WRITE(LUB,8160) IF (NZWE.LE.2) THEN WRITE(CBUF,8170) NZWE CALL STTPUT(CBUF,ISTAT) STOP ENDIF A(1)=0.0 A(2)=0.0 DO 110 JNC=1,NCAS 110 A(1)=A(1)+Y(JNC)*WEIGHTS(JNC) A(1)=A(1)/AL DO 120 JNC=1,NCAS 120 A(2)=A(2)+((Y(JNC)-A(1))**2)*WEIGHTS(JNC) A(2)=SQRT(A(2)/(AL-1.0)) WRITE(LUB,8180) A(1),A(2) WRITE(LUB,8190) NZWE WRITE(LUB,8200) AVW JREG=3 CALL RDUAL(A,2,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL, 1 NDXX,NDXY,JHEAD,LAB) 130 WRITE (LUB,8130) WRITE(CBUF,8210) CALL STTPUT(CBUF,ISTAT) IF (YNSAVE.EQ.'Y') THEN WRITE(CBUF,8220) FNAMEC CALL STTPUT(CBUF,ISTAT) ENDIF IF (FNAMEA.NE.'CON') THEN WRITE(CBUF,8230) FNAMEA CALL STTPUT(CBUF,ISTAT) ENDIF IF (FNAMEB.NE.'CON'.AND.FNAMEB.NE.'PRN') THEN WRITE(CBUF,8240) FNAMEB CALL STTPUT(CBUF,ISTAT) ENDIF IF (JPRT.NE.10) STOP 8000 FORMAT(1X,34('*')/' * ROBUST ESTIMATION OF LOCATION. *'/ 1 1X,34('*')///' NUMBER OF CASES = ',I5//) 8010 FORMAT(' DATA SET = ',A60/) 8020 FORMAT(/' THE DESIRED OUTPUT IS SMALL.'/) 8030 FORMAT(/' THE DESIRED OUTPUT IS MEDIUM-SIZED.'/) 8040 FORMAT(/' THE DESIRED OUTPUT IS LARGE.'/) 8050 FORMAT(/' AN INDEX PLOT WILL BE DRAWN.'/) 8070 FORMAT(/' THERE ARE ONLY ',I4,' CASES. THE ANALYSIS MUST', 1 ' BE STOPPED.'/) 8060 FORMAT(//' PLEASE ENTER YOUR DATA.'/) 8080 FORMAT(1X,' THE DATA FOR CASE NUMBER ',I4,' : ',$) 8090 FORMAT(' THE DATA FOR CASE NUMBER ',I4,' : ',$) 8100 FORMAT(' THE OBSERVATIONS ARE:'/) 8110 FORMAT(I5,2X,F10.4) 8120 FORMAT(//' THE MEDIAN',12X,'= ',F11.4,2X, 1 ' THE MAD (X 1.4826)',4X,'= ',F11.4/) 8125 FORMAT(/' MORE THAN HALF OF THE DATA ARE EQUAL.'//) 8130 FORMAT(/1X,78('*')/) 8140 FORMAT(//' THE MEAN',14X,'= ',F11.4,2X, 1 ' STANDARD DEVIATION',4X,'= ',F11.4//) 8150 FORMAT(//' LMS',19X,'= ',F11.4,3X, 1 'CORRESPONDING SCALE',3X,'= ',F11.4//) 8160 FORMAT(/' REWEIGHTING USING THE WEIGHTS BASED ON THE LMS.'/ 1 1X,47('*')/) 8170 FORMAT(/' THERE ARE ONLY ',I4,' CASES WITH NON-ZERO WEIGHT.'/ 1 ' THE ANALYSIS MUST BE STOPPED.'/) 8180 FORMAT(//' WEIGHTED MEAN',9X,'= ',F11.4,2X, 1 ' WEIGHTED STAND. DEV. = ',F11.4//) 8190 FORMAT(/' THERE ARE',2X,I4,' POINTS WITH NON-ZERO WEIGHT. ') 8200 FORMAT(/' AVERAGE WEIGHT',13X,'= ',F20.6//) 8210 FORMAT(////' THE RUN HAS SUCCESSFULLY BEEN EXECUTED . '//) 8220 FORMAT(' THE DATA IS SAVED IN FILE : ',A30) 8230 FORMAT(' THE DATA HAS BEEN READ FROM FILE : ',A30) 8240 FORMAT(' THE OUTPUT HAS BEEN WRITTEN IN FILE : ',A30) RETURN END CC CC *SMISLOC* : HANDLING OF MISSING VALUES IN THE CASE OF LOCATION CC SUBROUTINE SMISLOC(NCAS,NDXX,NDXY,X,Y,NMVAL,NSTOP) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),Y(NDXY),NMVAL(NDXY) CHARACTER*80 CBUF INTEGER ISTAT MAXM=(NCAS*4)/5 JND=0 JCAS=0 CC-----GIVE THE MISSING VALUE CODE WRITE(CBUF,8000) CALL STTPUT(CBUF,ISTAT) 5 READ(*,*,ERR=5)CODE DO 10 JNC=1,NCAS IF (Y(JNC).EQ.CODE) THEN JND=JND+1 ELSE JCAS=JCAS+1 X(1,JCAS)=Y(JNC) NMVAL(JCAS)=JNC ENDIF 10 CONTINUE NCAS=JCAS WRITE(CBUF,8010) NCAS CALL STTPUT(CBUF,ISTAT) IF (JND.GT.MAXM) THEN CALL STTPUT(' MORE THAN 80 PERCENT OF THE CASES HAD TO 1 BE DELETED BECAUSE OF THE MISSING VALUES.',ISTAT) CALL STTPUT('THE ANALYSIS WILL BE STOPPED.',ISTAT) NSTOP=1 RETURN ENDIF DO 20 JNC=1,NCAS Y(JNC)=X(1,JNC) 20 CONTINUE 8000 FORMAT(/' PLEASE ENTER THE VALUE OF THIS VARIABLE WHICH'/ 1 ' HAS TO BE INTERPRETED AS THE MISSING VALUE CODE : ',$) 8010 FORMAT(/' THERE ARE ',I4,' CASES STAYING IN THE ANALYSIS.'/) RETURN END