*-- Author : S Burke / J.V. Morris
SUBROUTINE FKLFTR(IERR)
**********************************************************************
* *
* Kalman fit with removal of the initial state vector *
* *
* Calling sequence is as for FKLFIT *
* *
* ERROR CONDITIONS; *
* *
* IERR = 0 ; normal termination *
* -> IERR = 101 ; no starting point was provided *
* -> IERR = 102 ; not enough measurements to fit *
* -> IERR = 103 ; invalid value in MES array *
* -> IERR = 104 ; invalid value of JSTART, JSTOP or JLAST *
* IERR = 8 ; covariance n.p.d. in FKLWM (1st call) *
* IERR = 9 ; covariance n.p.d. in FKLWM (2nd call) *
* IERR = 20 + n ; 2 < n < 10 iterations in point rejection *
* -> IERR = 130 ; 10 iterations in point rejection *
* IERR = 30 + n ; 2 < n < 10 iterations over fit sections *
* -> IERR = 140 ; 10 iterations over fit sections *
* IERR = 40 + n ; 1 < n < 10 restarts *
* -> IERR = 150 ; 10 restarts *
* -> IERR = 200 + ee ; fatal error ee from FKLPRO *
* -> IERR = 300 + ee ; fatal error ee from FKLFIL *
* -> IERR = 400 + ee ; fatal error ee from FKLSMO *
* -> IERR = 500 + ee ; fatal error ee from FKLPRS *
* -> IERR = 600 + ee ; fatal error ee from FKLPAS *
* *
* -> Fatal errors *
* *
**********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (IUTIL=0,IROUT=10)
*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.
*KEEP,FKNPL.
CHARACTER*5 CKDBG
PARAMETER (CKDBG='FKDBG')
PARAMETER (NPL=72)
LOGICAL LTRUE,LFIRST,LTRPL,LTRPLD
DOUBLE PRECISION TRUE,RTRUE,CHITRU,SPRO,CPRO,SFIL,CFIL
&, SSMT,CSMT,SSMTR,CSMTR,DPRO,CBPRO,QPRO,QGAIN
&, RPRO,CRPRO,RFIL,CRFIL,RSMT,CRSMT,CHIFIL,CHISMT
*
* Per-track values can go in H1WORK; note that LTRUE and LFIRST must
* be set at least per event.
*
* This is about 36k words long; the remaining common blocks are
* about 3.6k in total. Some of this could be in /H1WORK/, but the
* blocks would have to be reorganised.
*
COMMON /H1WORK/
* /FKPROJ/
& SPRO(5,NPL),CPRO(5,5,NPL)
* /FKFILT/
&, SFIL(5,NPL),CFIL(5,5,NPL)
* /FKSMTH/
&, SSMT(5,NPL),CSMT(5,5,NPL)
&, SSMTR(5,NPL),CSMTR(5,5,NPL)
* /FKINT/
&, DPRO(5,5,NPL),CBPRO(5,5,NPL),QPRO(5,5,NPL)
&, QGAIN(5,5,NPL),IAPROX,LFIRST
* /FKRSID/
&, RPRO(2,NPL),CRPRO(2,2,NPL),RFIL(2,NPL)
&, CRFIL(2,2,NPL),RSMT(2,NPL),CRSMT(2,2,NPL)
&, CHIFIL(NPL),CHISMT(NPL)
* /FKTRUE/
&, TRUE(5,NPL),RTRUE(5,NPL),CHITRU(NPL),LTRUE
* /FKDBG/
&, LTRPL(NPL),LTRPLD(NPL)
*KEEP,FKCNTL.
COMMON /FKCNTL/ LUN,IPR,ITR,IPL,JSTART,JSTOP,JLAST,JSTEP
*KEEP,FKFLAG.
LOGICAL LPRO,LFIL,LSMT,LMES,LRAD,LRPRO,LRFIL,LRSMT,LPOINT,LBLOCK
COMMON /FKFLAG/ LPRO(NPL),LFIL(NPL),LSMT(NPL),LMES(NPL)
&, LRAD(NPL),LRPRO,LRFIL,LRSMT,LPOINT,LBLOCK
*KEEP,FKCONS.
DOUBLE PRECISION ZPL,DZPL,RADL
COMMON /FKCONS/ ZPL(NPL),DZPL(NPL),RADL(NPL)
*KEEP,FKPROJ.
*KEEP,FKFILT.
*KEEP,FKSMTH.
*KEEP,FKRJCT.
DOUBLE PRECISION CHITOT,X2PCUT,X2CUTB,X2CUTA,X2CUTN
&, X2PCTI,X2CTBI,X2CTAI,X2CTNI
LOGICAL LWIRE,LPRINI
COMMON /FKRJCT/ X2PCUT,X2CUTB,X2CUTA,X2CUTN
&, X2PCTI,X2CTBI,X2CTAI,X2CTNI
&, CHITOT(NPL),NDF(NPL)
&, NBLOCK(NPL),NBADP(NPL),NBADB(NPL)
&, NFAILP(NPL),NFAILB(NPL),NNEWP(NPL)
&, NUNRJP(NPL),NUNRJB(NPL),NRERJP(NPL)
&, NCPRS,NBPRS,NCPAS,NPASS,IRJCT(NPL)
&, LWIRE(NPL),LPRINI
*KEEP,FKTRUE.
*KEEP,FKINT.
*KEND.
**********************************************************************
LOGICAL LFIXUP
DIMENSION STEMP(5),CTEMP(5,5)
**********************************************************************
IERR = 0
* Remember values of parameters which will be changed
KSTART = JSTART
KSTOP = JSTOP
KLAST = JLAST
CALL UCOPY(SPRO(1,JSTART),STEMP,10)
CALL FKCOPY(CPRO(1,1,JSTART),CTEMP)
* Initialise flags
LFIRST = .TRUE.
LFIXUP = .FALSE.
NTRY = 1
* Re-entry point for re-starts
1000 CONTINUE
* These need to be re-set for each pass
NSEC = 1
JEND = JLAST
JBEGIN = JSTOP
* Set the stop point to the start point
JSTOP = JSTART
* Switch off residuals and point rejection for the first pass
CALL FKPRSV(1)
CALL FKLFIT(IFAIL)
CALL FKPRSV(-1)
* Check errors
IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IFAIL/100,IFAIL,IERR)
IF (IFAIL.GT.100) GOTO 7000
* If this is getting too complicated, start again from scratch!
IF (LFIXUP .AND. NPASS.GT.1) GOTO 9000
* Now loop over the two sections until no points are altered
ISEC = 1
2000 CONTINUE
ISEC = -ISEC
* Remove the initial vector (to create a FILTERED vector)
CALL FKLWM(-1,SSMT(1,JSTART),CSMT(1,1,JSTART),SPRO(1,JSTART),
& CPRO(1,1,JSTART),SFIL(1,JSTART),CFIL(1,1,JSTART),IFAIL)
* If the vector can't be removed we have a serious problem
IF (IFAIL.GT.100) THEN
CALL FKERR(IUTIL,IROUT,IWARN,IOCV+(ISEC-5)/2,IERR)
* If we've rejected some points, start again from the beginning
IF (NSEC.GT.2) GOTO 9000
* First part is OK, so do half of the filter the other way round
IF (ISEC.EQ.1) GOTO 8000
* If it happens at the start, leave the initial vector in and carry on
CALL UCOPY(SSMT(1,JSTART),SFIL(1,JSTART),10)
CALL FKCOPY(CSMT(1,1,JSTART),CFIL(1,1,JSTART))
ENDIF
* Change direction
JSTEP = -JSTEP
* We've filtered at this point, so project into the other section
LPRO(JSTART+JSTEP) = .FALSE.
DZPL(JSTART) = ZPL(JSTART+JSTEP) - ZPL(JSTART)
CALL FKLPRO(JSTART,IFAIL)
IF (IFAIL.GT.100) THEN
CALL FKERR(IUTIL,IROUT,IFPRO,IFAIL,IERR)
GOTO 7000
ENDIF
* Flip the start and end points, and filter the other section
IF (ISEC.EQ.1) THEN
JLAST = JEND
ELSE
JLAST = JBEGIN
ENDIF
JSTART = JSTART + JSTEP
JSTOP = JSTART
CALL FKLFIT(IFAIL)
IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IFAIL/100,IFAIL,IERR)
IF (IFAIL.GT.100) GOTO 7000
* We've been once over every plane, so don't zero IRJCT again
LFIRST = .FALSE.
* If NPASS is > 1 we rejected some points, so quit and start again
IF (LFIXUP .AND. NPASS.GT.1) GOTO 9000
* If no point was rejected, and this isn't the first time, that's all
IF (LFIXUP) THEN
* Don't change error code
IF (NTRY.GT.1) CALL FKERR(IUTIL,IROUT,IWARN,IFREE2+NTRY,IER)
* Restore the saved vectors
CALL FKSAVE(-1,JBEGIN,JSTART-JSTEP)
ELSEIF (NPASS.EQ.1 .AND. NSEC.GT.1) THEN
* Pass back code 11 if it occurred
IER = 0
IF (NSEC.GT.2) CALL FKERR(IUTIL,IROUT,IWARN,IFREE1+NSEC,IER)
IF (NTRY.GT.1) CALL FKERR(IUTIL,IROUT,IWARN,IFREE2+NTRY,IER)
IF (IERR.NE.IOCV .AND. IER.GT.0) IERR = IER
ELSE
NSEC = NSEC + 1
IF (NSEC.LT.10) GOTO 2000
IF (NTRY.GT.1) CALL FKERR(IUTIL,IROUT,IWARN,IFREE2+NTRY,IERR)
CALL FKERR(IUTIL,IROUT,IFATAL,IFREE2,IERR)
ENDIF
7000 CONTINUE
* Reset the end-point flags so we don't screw up the calling routine
JSTART = KSTART
JSTOP = KSTOP
JLAST = KLAST
RETURN
8000 CONTINUE
*
* We get here if FKLWM fails on the second removal
* The fixup is to do half of the filter the other way round
*
LFIXUP = .TRUE.
* Save the smoothed vectors/residuals between JSTART and JLAST
CALL FKSAVE(1,JSTART,JLAST)
* Reset the starting values, and start again
JSTART = KSTART
JSTOP = KLAST
JLAST = KSTOP
GOTO 1000
* This resets everything and starts again
9000 CONTINUE
NTRY = NTRY + 1
IF (NTRY.GE.10) THEN
CALL FKERR(IUTIL,IROUT,IFATAL,IFREE3,IERR)
GOTO 7000
ENDIF
CALL UCOPY(STEMP,SPRO(1,JSTART),10)
CALL FKCOPY(CTEMP,CPRO(1,1,JSTART))
LFIXUP = .FALSE.
JSTART = KSTART
JSTOP = KSTOP
JLAST = KLAST
GOTO 1000
END
*