*-- Author : S.Burke / J.V. Morris SUBROUTINE FKLSMO(JPL,IERR) ********************************************************************** * * * Kalman Smoothing from plane JPL+1 to plane JPL. * * * * ERROR CONDITIONS; * * IERR = 0 ; normal termination * * -> IERR = 101 ; no smoothed data at plane JPL+1 * * -> IERR = 102 ; no projected data at plane JPL+1 * * IERR = 5 ; smoothing already done ... but continue * * IERR = 11 ; smoothed covariance matrix n.p.d. * * IERR = 12 ; covariance of smoothed residuals n.p.d. * * IERR = 116 ; theta > pi/2: reset to pi/4 * * IERR = 17 ; theta > 1 (warning) * * * * -> Fatal errors * * * * NB Error 12 is not considered fatal, but the chi-sq will be zero * * * ********************************************************************** *KEEP,FKECODE. *KEND. * * Common block definitions * *KEEP,FKNPL. * * 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. * * /FKPROJ/ * /FKFILT/ * /FKSMTH/ * /FKINT/ * /FKRSID/ * /FKTRUE/ * /FKDBG/ *KEEP,FKCNTL. *KEEP,FKCONS. *KEEP,FKTRUE. *KEEP,FKFLAG. *KEEP,FKPROJ. *KEEP,FKFILT. *KEEP,FKSMTH. *KEEP,FKRSID. *KEEP,FKINT. *KEND. ********************************************************************** * * Local arrays etc ... * ********************************************************************** * * Initialisation and checks ... * * If JPL=JLAST then Smoothed is same as Filtered ... * Does filtered data exist at plane JPL ? CALL FKERR(IUTIL,IROUT,IFATAL,IINF1,IERR) CALL UCOPY(SFIL(1,JPL),SSMT(1,JPL),10) CALL FKCOPY(CFIL(1,1,JPL),CSMT(1,1,JPL)) * Number of next plane * Does smoothed data exist at previous plane? CALL FKERR(IUTIL,IROUT,IFATAL,IINF1,IERR) * RADL and LRAD are defined in an asymmetrical way * Does projected data exist at previous plane? (only needed with MS) CALL FKERR(IUTIL,IROUT,IFATAL,IINF2,IERR) * Has smoothing to JPL already been done? Continue anyway ..... !? IF (LSMT(JPL)) CALL FKERR(IUTIL,IROUT,IWARN,IDONE,IERR) ********************************************************************** * * Update the smoothed state vector at plane JPL+1 by removing the * estimated multiple scattering between planes JPL and JPL+1; * I think that this is probably the most `correct' way to do it. * If there is no multiple scattering QGAIN is zero, so this can be * skipped. Then transform the adjusted state vector from JPL+1 to JPL. * CALL FKADJ(SSMT(1,JPLN),SPRO(1,JPLN),QGAIN(1,1,JPL),SADJ) CALL FKNORM(SADJ,IFAIL) IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IFAIL/100,IFAIL,IERR) CALL FKTRAN(-DZPL(JPL),ZPL(JPLN),SADJ,SSMT(1,JPL),DINV) * * Compute the smoothed covariance. This may be numerically unstable if * the multiple scattering is large (probably not the case). * CALL FKDMQD(DINV,QGAIN(1,1,JPL),AGAIN) CALL FKMXM(CSMT(1,1,JPLN),AGAIN,CSMT(1,1,JPL)) CALL FKDQA(DINV,QPRO(1,1,JPL),AGAIN,CSMT(1,1,JPL)) * If there is no MS, smoothing is just back-extrapolation * Transform the state vector and covariance from JPL to JPL+1 ... CALL FKTRAN(-DZPL(JPL),ZPL(JPLN),SSMT(1,JPLN), CALL FKMUL(CSMT(1,1,JPLN),DINV,CSMT(1,1,JPL)) CALL FKNORM(SSMT(1,JPL),IFAIL) IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IFAIL/100,IFAIL,IERR) CALL FKRST(JPL,IFAIL) IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IWARN,IFAIL,IERR) CALL FKLRSD(JPL,SSMT(1,JPL),CSMT(1,1,JPL),-3, IF (IFAIL.NE.0) CALL FKERR(IUTIL,IROUT,IWARN,IFAIL,IERR) CALL VZERO(RSMT(1,JPL),4) CALL VZERO(CRSMT(1,1,JPL),8) * Set flag to show smoothing is done *