C===========================================================================
C Copyright (C) 1995-2005 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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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=50)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
       RETURN
 50    CALL STETER('Error reading input!')
       RETURN

 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.'/)
       END
