*-- Author : S.Burke / J.V. Morris SUBROUTINE FKCOVR(CPRO,H,GMES,CFIL,WT,IERR) *-----------------------------------------Updates 24/01/92------- **: FKCOVR 30205.SB. Overflow error trapped. *-----------------------------------------Updates---------------- ********************************************************************** * * * Calculate filtered covariance for a `radial' (i.e. two * * measurements) * * * * ERROR CONDITIONS; * * IERR = 0 ; normal termination * * -> IERR = 111 ; filtered covariance not positive definite * * -> IERR = 113 ; determinant .GT. 10**20 * * * * -> Fatal error * * * ********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IUTIL=1,IROUT=2) *KEEP,FKECODE. PARAMETER (IWARN=0,IFATAL=1,IFPRO=2,IFFLT=3,IFSMO=4,IFPRS=5, & IFPAS=6,IFPAF=7) PARAMETER (IINF1=1,IINF2=2,IINF3=3,IINV=4,IDONE=5) PARAMETER (IICV=6,IMCV=7,IOCV=11,IRCV=12,IOVCV=13, & ITHGP2=16,ITHG1=17) PARAMETER (IFREE=20,IFREE1=30,IFREE2=40,IFREE3=50) *KEND. PARAMETER (GCMAX=1.0D10) DIMENSION CPRO(5,5),H(2,2),GMES(2,2),CFIL(5,5),WT(2,5) ********************************************************************** IERR = 0 * First rotate the measurement covariance ... CR11 = H(1,1)*(GMES(1,1)*H(1,1) + 2.*GMES(2,1)*H(2,1)) + & H(2,1)*GMES(2,2)*H(2,1) HG21 = H(1,2)*GMES(1,1) + H(2,2)*GMES(2,1) HG22 = H(1,2)*GMES(2,1) + H(2,2)*GMES(2,2) CR21 = HG21*H(1,1) + HG22*H(2,1) CR22 = HG21*H(1,2) + HG22*H(2,2) * ... then calculate the top two rows of G.CPRO + 1 ... GC11 = CR11*CPRO(1,1) + CR21*CPRO(2,1) + 1.D0 GC12 = CR11*CPRO(2,1) + CR21*CPRO(2,2) GC13 = CR11*CPRO(3,1) + CR21*CPRO(3,2) GC14 = CR11*CPRO(4,1) + CR21*CPRO(4,2) GC15 = CR11*CPRO(5,1) + CR21*CPRO(5,2) GC21 = CR21*CPRO(1,1) + CR22*CPRO(2,1) GC22 = CR21*CPRO(2,1) + CR22*CPRO(2,2) + 1.D0 GC23 = CR21*CPRO(3,1) + CR22*CPRO(3,2) GC24 = CR21*CPRO(4,1) + CR22*CPRO(4,2) GC25 = CR21*CPRO(5,1) + CR22*CPRO(5,2) IF (GC11.GT.GCMAX .OR. GC22.GT.GCMAX .OR. & GC12.GT.GCMAX .OR. GC21.GT.GCMAX) THEN CALL FKERR(IUTIL,IROUT,IFATAL,IOVCV,IERR) RETURN ENDIF * ... and invert (only the top two rows are non-trivial) ... DET = GC11*GC22 - GC12*GC21 IF (DET.LE.0.) THEN CALL FKERR(IUTIL,IROUT,IFATAL,IOCV,IERR) RETURN ENDIF DET = 1.D0/DET GDET22 = GC22*DET GDET12 = GC12*DET WT(1,1) = GDET22 WT(1,2) = -GDET12 WT(1,3) = GDET12*GC23 - GC13*GDET22 WT(1,4) = GDET12*GC24 - GC14*GDET22 WT(1,5) = GDET12*GC25 - GC15*GDET22 GDET11 = GC11*DET GDET21 = GC21*DET WT(2,1) = -GDET21 WT(2,2) = GDET11 WT(2,3) = GC13*GDET21 - GDET11*GC23 WT(2,4) = GC14*GDET21 - GDET11*GC24 WT(2,5) = GC15*GDET21 - GDET11*GC25 * ... and finally, premultiply by CPRO CFIL(1,1) = CPRO(1,1)*WT(1,1) + CPRO(2,1)*WT(2,1) CFIL(2,1) = CPRO(2,1)*WT(1,1) + CPRO(2,2)*WT(2,1) CFIL(2,2) = CPRO(2,1)*WT(1,2) + CPRO(2,2)*WT(2,2) CFIL(3,1) = CPRO(3,1)*WT(1,1) + CPRO(3,2)*WT(2,1) CFIL(3,2) = CPRO(3,1)*WT(1,2) + CPRO(3,2)*WT(2,2) CFIL(3,3) = CPRO(3,1)*WT(1,3) + CPRO(3,2)*WT(2,3) + CPRO(3,3) CFIL(4,1) = CPRO(4,1)*WT(1,1) + CPRO(4,2)*WT(2,1) CFIL(4,2) = CPRO(4,1)*WT(1,2) + CPRO(4,2)*WT(2,2) CFIL(4,3) = CPRO(4,1)*WT(1,3) + CPRO(4,2)*WT(2,3) + CPRO(4,3) CFIL(4,4) = CPRO(4,1)*WT(1,4) + CPRO(4,2)*WT(2,4) + CPRO(4,4) CFIL(5,1) = CPRO(5,1)*WT(1,1) + CPRO(5,2)*WT(2,1) CFIL(5,2) = CPRO(5,1)*WT(1,2) + CPRO(5,2)*WT(2,2) CFIL(5,3) = CPRO(5,1)*WT(1,3) + CPRO(5,2)*WT(2,3) + CPRO(5,3) CFIL(5,4) = CPRO(5,1)*WT(1,4) + CPRO(5,2)*WT(2,4) + CPRO(5,4) CFIL(5,5) = CPRO(5,1)*WT(1,5) + CPRO(5,2)*WT(2,5) + CPRO(5,5) RETURN END *