C @(#)prd.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 *RANGS* : SORTS A VECTOR OF LENGTH 'NCAS'. CC SUBROUTINE RANGS(AM,LAME) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION JLV(30),JRV(30),AM(LAME) INTEGER J,JNC,JNDL,JR,JSS,JTWE,LAME REAL AMM,XX INTEGER JLV(30),JRV(30) REAL AM(LAME) JSS=1 JLV(1)=1 JRV(1)=LAME 10 JNDL=JLV(JSS) JR=JRV(JSS) JSS=JSS-1 20 JNC=JNDL J=JR JTWE=(JNDL+JR)/2 XX=AM(JTWE) 30 IF(AM(JNC).GE.XX)GOTO 40 JNC=JNC+1 GOTO 30 40 IF(XX.GE.AM(J))GOTO 50 J=J-1 GOTO 40 50 IF(JNC.GT.J)GOTO 60 AMM=AM(JNC) AM(JNC)=AM(J) AM(J)=AMM JNC=JNC+1 J=J-1 60 IF (JNC.LE.J) GOTO 30 IF ((J-JNDL).LT.(JR-JNC)) GOTO 80 IF (JNDL.GE.J) GOTO 70 JSS=JSS+1 JLV(JSS)=JNDL JRV(JSS)=J 70 JNDL=JNC GOTO 100 80 IF (JNC.GE.JR) GOTO 90 JSS=JSS+1 JLV(JSS)=JNC JRV(JSS)=JR 90 JR=J 100 IF (JNDL.LT.JR) GOTO 20 IF (JSS.NE.0) GOTO 10 RETURN END CC CC *RDUAL* : CALCULATES THE RESIDUALS OF ALL CASES CC SUBROUTINE RDUAL(AA,JKEUS,JDA,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC, 1 JREG,X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL, 1 NDXX,NDXY,JHEAD,LAB) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY) C DIMENSION XMED(NDXX),XMAD(NDXX),AW(NDXY),NMVAL(NDXY) INTEGER J,JCST,JDA,JKEUS,JNC,JPL2,JPLT,JPRT,JREG INTEGER NCAS,NDXX,NDXY,NVAD,NVAR,NZWE,LUB REAL AL,AVW,BBB,PREC,YJ REAL AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY) REAL XMED(NDXX),XMAD(NDXX),AW(NDXY) INTEGER NMVAL(NDXY) CHARACTER*10 LAB(NDXX) CHARACTER*60 JHEAD JPL2=1 AL=NCAS IF(JPRT.EQ.0) GOTO 10 IF (NVAR.EQ.1.AND.JCST.EQ.1) THEN IF (JKEUS.NE.2) WRITE(LUB,8000) LAB(NVAD) IF (JKEUS.EQ.2) WRITE(LUB,8010) LAB(NVAD) ELSE IF (JKEUS.NE.2) WRITE(LUB,8020) LAB(NVAD),LAB(NVAD) IF (JKEUS.EQ.2) WRITE(LUB,8030) LAB(NVAD),LAB(NVAD) ENDIF 10 DO 20 JNC=1,NCAS IF (.NOT.(NVAR.EQ.1.AND.JCST.EQ.1)) GOTO 40 BBB=Y(JNC)-AA(1) GOTO 60 40 AW(JNC)=0.0 DO 30 J=1,NVAR 30 AW(JNC)=AW(JNC)+AA(J)*(X(J,JNC)*XMAD(J)+XMED(J)) YJ=Y(JNC)*XMAD(NVAD)+XMED(NVAD) BBB=YJ-AW(JNC) 60 IF (AA(NVAD).GT.PREC) GOTO 70 IF (ABS(BBB).LT.PREC) BBB=0.0 RESDU(JNC)=BBB GOTO 50 70 RESDU(JNC)=BBB/AA(NVAD) 50 IF (JPRT.EQ.0) GOTO 20 IF (NVAR.EQ.1.AND.JCST.EQ.1) GOTO 80 IF (AA(NVAD).GT.PREC.AND.JKEUS.NE.2) 1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC) IF (AA(NVAD).LE.PREC) 1 WRITE(LUB,8050) YJ,AW(JNC),BBB,NMVAL(JNC) IF (AA(NVAD).GT.PREC.AND.JKEUS.EQ.2) 1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC), 1 WEIGHTS(JNC) GOTO 20 80 IF (JKEUS.NE.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC) IF (JKEUS.EQ.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC), 1 WEIGHTS(JNC) 20 CONTINUE IF(JPLT.LT.1) GOTO 90 IF (AA(NVAD).LE.PREC) JPL2=0 CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG, 1 JHEAD,NDXX,NVAD,LAB) IF (JPLT.LT.3) GOTO 90 JPLT=2 CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG, 1 JHEAD,NDXX,NVAD,LAB) JPLT=3 90 CONTINUE 8000 FORMAT(//11X,'OBSERVED',12X,'RESIDUAL', 1 ' NO',' RES/SC'/9X,A10/) 8010 FORMAT(//11X,'OBSERVED',12X,'RESIDUAL', 1 ' NO',' RES/SC',' WEIGHT'/9X,A10/) 8020 FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL', 1 3X,'NO',' RES/SC'/9X,A10,10X,A10/) 8030 FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL', 1 3X,'NO',' RES/SC',' WEIGHT'/9X,A10,10X,A10/) 8040 FORMAT(1X,F18.5,2(1X,F19.5),I5,1X,F7.2,F6.1) 8050 FORMAT(1X,F18.5,2(1X,F19.5),I5,' *****') 8060 FORMAT(1X,F18.5,1X,F19.5,I5,1X,F7.2,F6.1) RETURN END CC CC *SHHLF* : CALCULATES THE LMS IN THE ONE-DIMENSIONAL CASE. CC SUBROUTINE SHHLF(W,NCAS,JQU,SLUTN,BSTD,PREC) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION W(NCAS) INTEGER JNC,JPK,JQU,JTEL,K,LNK,NCAS REAL ABSW,BSTD,PREC,SLUTN,TEL,VERSC REAL W(NCAS) K=JQU-1 VERSC=W(K+1)-W(1) SLUTN=(W(1)+W(K+1))/2.0 BSTD=((W(K+1)-W(1))/2.0) IF (NCAS.EQ.2) RETURN JTEL=1 LNK=NCAS-K DO 10 JNC=2,LNK JPK=JNC+K ABSW=W(JPK)-W(JNC) IF ( ABS(ABSW-VERSC).GT.PREC) GOTO 20 JTEL=JTEL+1 SLUTN=SLUTN+(W(JNC)+W(JPK))/2.0 GOTO 10 20 IF (ABSW.LT.VERSC) THEN VERSC=ABSW JTEL=1 BSTD=((W(JPK)-W(JNC))/2.0) SLUTN=(W(JNC)+W(JPK))/2.0 ENDIF 10 CONTINUE TEL=JTEL SLUTN=SLUTN/TEL RETURN END CC CC *TRC* : RETRANSFORMATION OF THE VARIANCE-COVARIANCE MATRIX CC SUBROUTINE TRC(H,DA,NMXV,NDXX,NDXY,NVAR,JCST,NVSB,NVAD, 1 XMED,XMAD) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION H(NMXV,NDXX),DA(NMXV) C DOUBLE PRECISION H,XMP2,HNN C DIMENSION XMED(NDXX),XMAD(NDXX) INTEGER J,JCST,JU,K,K2,NDXX,NDXY,NMXV,NVAD,NVAR,NVSB REAL DA(NMXV) DOUBLE PRECISION H(NMXV,NDXX),XMP2,HNN REAL XMED(NDXX),XMAD(NDXX) XMP2=DBLE(XMAD(NVAD))*DBLE(XMAD(NVAD)) IF (JCST.EQ.0) THEN DO 10 J=1,NVAR DO 20 K=1,J H(J,K)=H(J,K)*(XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K)))) 20 CONTINUE DA(J)=DSQRT(H(J,J)) 10 CONTINUE ELSE DO 25 J=1,NVAR H(J,NVAD)=H(J,J) 25 CONTINUE DO 30 J=1,NVAR DO 40 K=1,J H(J,K)=H(J,K)*XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K))) 40 CONTINUE DA(J)=DSQRT(H(J,J)) 30 CONTINUE DO 50 K=1,NVSB H(NVAR,K)=H(K,NVAR)*XMP2/DBLE(XMAD(K)) DO 60 K2=1,NVAR IF (K.EQ.K2) THEN H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K))*XMP2/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,NVAD) GOTO 60 ENDIF IF (K.LT.K2) THEN H(NVAR,K)=H(NVAR,K)-(DBLE(XMED(K2))*XMP2)/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K,K2) ELSE H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K2))*XMP2/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,K) ENDIF 60 CONTINUE 50 CONTINUE H(NVAR,NVAR)=H(NVAR,NVAD)*XMP2 DO 70 K=1,NVAR H(NVAR,NVAR)=H(NVAR,NVAR)+ 1 (DBLE(XMED(K))*DBLE(XMED(K)))*XMP2/ 1 (DBLE(XMAD(K))*DBLE(XMAD(K)))*H(K,NVAD) 70 CONTINUE DO 80 K=1,NVAR IF (K.NE.NVAR) THEN H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/ 1 (DBLE(XMAD(K)))*H(K,NVAR) ELSE H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/ 1 (DBLE(XMAD(K)))*H(NVAR,NVAD) ENDIF 80 CONTINUE DO 90 J=1,NVSB JU=J+1 DO 90 K=JU,NVAR HNN=2.0D0*DBLE(XMED(J))*DBLE(XMED(K))*XMP2 H(NVAR,NVAR)=H(NVAR,NVAR)+HNN/ 1 (DBLE(XMAD(J))*DBLE(XMAD(K)))*H(J,K) 90 CONTINUE DA(NVAR)=DSQRT(H(NVAR,NVAR)) ENDIF RETURN END CC CC *SCHCV* : PRINTS THE VARIANCE-COVARIANCE MATRIX OF THE CC REGRESSION COEFFICIENTS. CC SUBROUTINE SCHCV(NVAR,NDXX,NMXV,H,LUB) C INCLUDE 'MID_REL_INCL:implicit.inc' INTEGER J,K,NDXX,NMXV,NVAR,LUB DOUBLE PRECISION H(NMXV,NDXX) WRITE(LUB,9000) DO 10 J=1,NVAR WRITE(LUB,9010) (H(J,K),K=1,J) 10 CONTINUE 9000 FORMAT(//' VARIANCE - COVARIANCE MATRIX = '/) 9010 FORMAT(5(5X,D10.4)) RETURN END CC CC *RTRAN* : RETRANSFORMATION OF THE REGRESSION COEFFICIENTS. CC SUBROUTINE RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,AA,JAL,FCKW) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION AA(JAL),XMED(NDXX),XMAD(NDXX) INTEGER J,JAL,JCST,NDXX,NVAD,NVAR,NVSB REAL FCKW REAL AA(JAL),XMED(NDXX),XMAD(NDXX) IF (NVAR.LE.1) THEN AA(1)=AA(1)*XMAD(NVAD)/XMAD(1) GOTO 30 ENDIF DO 10 J=1,NVSB 10 AA(J)=AA(J)*XMAD(NVAD)/XMAD(J) IF (JCST.EQ.0) THEN AA(NVAR)=AA(NVAR)*XMAD(NVAD)/XMAD(NVAR) ELSE AA(NVAR)=AA(NVAR)*XMAD(NVAD) DO 20 J=1,NVSB 20 AA(NVAR)=AA(NVAR)-AA(J)*XMED(J) AA(NVAR)=AA(NVAR)+XMED(NVAD) ENDIF 30 FCKW=FCKW*(XMAD(NVAD)*XMAD(NVAD)) RETURN END CC CC *RANPN* : RANDOM NUMBER GENERATOR . CC SUBROUTINE RANPN(NCAS,NVAR,JNDVC,NDXX,JEY,JLP,MAXRP) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION JNDVC(NDXX) INTEGER J,JEY,JEYQT,JLP,JNC,JNCM,JRAN,MAXRP,NCAS,NDXX,NVAR REAL AL,RNULE,RY INTEGER JNDVC(NDXX) CC WE PROGRAMMED THIS RANDOM NUMBER GENERATOR OURSELVES BECAUSE CC WE WANTED IT TO BE MACHINE INDEPENDENT. IT SHOULD RUN ON MOST CC COMPUTERS BECAUSE THE LARGEST INTEGER USED IS LESS THAN 2**30. CC THE PERIOD IS 2**16=65536, WHICH IS SUFFICIENT FOR THIS CC APPLICATION. AL=NCAS JLP=JLP+1 IF (JLP.LE.MAXRP) GOTO 10 RETURN 10 DO 20 JNC=1,NVAR 30 JEY=JEY*5761+999 JEYQT=JEY/65536 JEY=JEY-JEYQT*65536 RY=JEY RNULE=RY/65536.0 JRAN=INT(RNULE*AL)+1 IF (JNC.EQ.1) GOTO 50 JNCM=JNC-1 DO 40 J=1,JNCM IF (JNDVC(J).EQ.JRAN) GOTO 30 40 CONTINUE 50 JNDVC(JNC)=JRAN 20 CONTINUE RETURN END CC CC *GENPN* : GENERATES ALL COMBINATIONS OF 'NVAR' INTEGERS CC BETWEEN 1 AND 'NCAS', WHERE NVAR IS AT MOST FIVE. CC SUBROUTINE GENPN(NCAS,NVAR,JNDVC,NDXX,JJJ) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION JNDVC(NDXX) INTEGER J,JJJ,NDXX,NCAS,NVAR INTEGER JNDVC(NDXX) IF (JJJ.LE.1) THEN DO 5 J=1,NVAR 5 JNDVC(J)=J RETURN ENDIF IF (JNDVC(NVAR).EQ.NCAS) GOTO 10 JNDVC(NVAR)=JNDVC(NVAR)+1 RETURN 10 IF (JNDVC(NVAR-1).EQ.NCAS-1) GOTO 20 JNDVC(NVAR-1)=JNDVC(NVAR-1)+1 JNDVC(NVAR)=JNDVC(NVAR-1)+1 RETURN 20 IF(JNDVC(NVAR-2).EQ.NCAS-2) GOTO 30 JNDVC(NVAR-2)=JNDVC(NVAR-2)+1 JNDVC(NVAR-1)=JNDVC(NVAR-2)+1 JNDVC(NVAR)=JNDVC(NVAR-1)+1 RETURN 30 IF(JNDVC(NVAR-3).EQ.NCAS-3) GOTO 40 JNDVC(NVAR-3)=JNDVC(NVAR-3)+1 JNDVC(NVAR-2)=JNDVC(NVAR-3)+1 JNDVC(NVAR-1)=JNDVC(NVAR-2)+1 JNDVC(NVAR)=JNDVC(NVAR-1)+1 RETURN 40 IF(JNDVC(NVAR-4).EQ.NCAS-4) GOTO 50 JNDVC(NVAR-4)=JNDVC(NVAR-4)+1 JNDVC(NVAR-3)=JNDVC(NVAR-4)+1 JNDVC(NVAR-2)=JNDVC(NVAR-3)+1 JNDVC(NVAR-1)=JNDVC(NVAR-2)+1 JNDVC(NVAR)=JNDVC(NVAR-1)+1 RETURN 50 IF(JNDVC(NVAR-5).EQ.NCAS-5) RETURN JNDVC(NVAR-5)=JNDVC(NVAR-5)+1 JNDVC(NVAR-4)=JNDVC(NVAR-5)+1 JNDVC(NVAR-3)=JNDVC(NVAR-4)+1 JNDVC(NVAR-2)=JNDVC(NVAR-3)+1 JNDVC(NVAR-1)=JNDVC(NVAR-2)+1 JNDVC(NVAR)=JNDVC(NVAR-1)+1 RETURN END CC CC *EQUAT* : SOLVES A SYSTEM OF LINEAR EQUATIONS. CC SUBROUTINE EQUAT(AM,NMXV,NDXX,HVEC,JDMB,NA,NB,NERR,DTEM) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION HVEC(JDMB),AM(NMXV,NDXX) C DOUBLE PRECISION HVEC,TURN,SWAP,DTEM,DETER INTEGER J,JBEGC,JBEGX,JDEL,JENDC,JENDX,JDM,JDMB,JHFD INTEGER JMAT,JNC,JNCB,JNCC,JNCD,JNCE,JNCF,JNK,JROW INTEGER LCLPL,LDEL,N,NA,NB,NC,NDXX,NEQA,NERR,NMXV,NZNDE REAL AM(NMXV,NDXX) DOUBLE PRECISION HVEC(JDMB),TURN,SWAP,DTEM,DETER JDM=NMXV DETER=1.0 N=NA JMAT=N+NB JNK=0 DO 10 J=1,JMAT JNK=(J-1)*NMXV DO 10 NC=1,NMXV JNK=JNK+1 HVEC(JNK)=AM(NC,J) 10 CONTINUE NZNDE=N-1 LCLPL=-JDM DO 120 JHFD=1,N TURN=0. LCLPL=LCLPL+JDM+1 JDEL=LCLPL+N-JHFD DO 40 JNCB=LCLPL,JDEL IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30 30 TURN=HVEC(JNCB) LDEL=JNCB 40 CONTINUE IF(TURN) 50,170,50 50 IF(LDEL-LCLPL) 60,80,60 60 DETER=-DETER LDEL=LDEL-JDM JNCB=LCLPL-JDM DO 70 JNCC=JHFD,JMAT LDEL=LDEL+JDM JNCB=JNCB+JDM SWAP=HVEC(JNCB) HVEC(JNCB)=HVEC(LDEL) 70 HVEC(LDEL)=SWAP 80 DETER=DETER*TURN IF (JHFD.EQ.N) GOTO 120 TURN=1./TURN JNCB=LCLPL+1 DO 90 JNCC=JNCB,JDEL 90 HVEC(JNCC)=HVEC(JNCC)*TURN JNCD=LCLPL JROW=JHFD+1 DO 110 JNCB=JROW,N JNCD=JNCD+1 JNCE=LCLPL JNCF=JNCD DO 100 JNCC=JROW,JMAT JNCE=JNCE+JDM JNCF=JNCF+JDM 100 HVEC(JNCF)=HVEC(JNCF)-HVEC(JNCE)*HVEC(JNCD) 110 CONTINUE 120 CONTINUE DTEM=DETER NERR=0 NEQA=N+1 JBEGX=NZNDE*JDM+1 DO 150 JNC=NEQA,JMAT JBEGX=JBEGX+JDM JENDX=JBEGX+N JBEGC=N*JDM+1 JENDC=JBEGC+NZNDE DO 140 JNCB=1,NZNDE JENDX=JENDX-1 JBEGC=JBEGC-JDM JENDC=JENDC-JDM-1 HVEC(JENDX)=HVEC(JENDX)/HVEC(JENDC+1) SWAP=HVEC(JENDX) JNCD=JBEGX-1 DO 130 JNCC=JBEGC,JENDC JNCD=JNCD+1 130 HVEC(JNCD)=HVEC(JNCD)-HVEC(JNCC)*SWAP 140 CONTINUE HVEC(JBEGX)=HVEC(JBEGX)/HVEC(1) 150 CONTINUE JNC=-JDM JBEGX=NZNDE*JDM+1 JENDX=JBEGX+NZNDE DO 160 JNCB=NEQA,JMAT JBEGX=JBEGX+JDM JENDX=JENDX+JDM JNC=JNC+JDM JNCD=JNC DO 160 JNCC=JBEGX,JENDX JNCD=JNCD+1 160 HVEC(JNCD)=HVEC(JNCC) GOTO 180 170 NERR=-1 DTEM=DETER 180 JNK=0 DO 190 J=1,JMAT DO 190 NC=1,NMXV JNK=JNK+1 AM(NC,J)=HVEC(JNK) 190 CONTINUE RETURN END