*-- Author : S.Burke / J.V. Morris
      SUBROUTINE FKSCAT(DZ,S,X0,D,Q)
**********************************************************************                                        
*                                                                    *                                        
* Compute the Multiple Scattering Matrix, Q, for transformation of   *                                        
* State Vector, S, from Z to Z+DZ.                                   *                                        
*                                                                    *                                        
**********************************************************************                                        
                                                                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               
                                                                        
      DIMENSION S(5),Q(5,5),D(5,5)                                      
                                                                        
* TT2MIN sets the lower limit allowed for tan**2(theta) (arbitrary)                                           
      PARAMETER (TT2MIN=1.D-20)                                         
      PARAMETER (G=0.5D0,G2=1.0D0/3.0D0,PSCALE=0.0141)                  
                                                                        
**********************************************************************                                        
                                                                        
*                                                                                                             
* Test for finite step length and non-zero radiation length ....                                              
*                                                                                                             
      IF (DZ.EQ.0.D0 .OR. X0.LE.0.D0) THEN                              
         CALL VZERO(Q,50)                                                                              
         RETURN                                                         
      ENDIF                                                             
                                                                        
* Fix for apparent singularity ... this should be checked further                                             
      TT2 = S(4)*S(4)                                                   
      IF (TT2.LT.TT2MIN) TT2 = TT2MIN                                   
      CT2I = 1.D0 + TT2                                                 
      ST2I = 1.D0/TT2 + 1.D0                                            
                                                                        
* Compute path length in medium = X ....                                                                      
      X = DABS(DZ*DSQRT(CT2I))                                          
                                                                        
*                                                                                                             
* compute the mean square scattering angle (SQMA) assuming Beta=1 ...                                         
* (this is the projection onto plane)                                                                         
*                                                                                                             
      SQMA   = (X/X0)*(PSCALE*S(3))**2                                  
                                                                        
      G2SQMA = G2*SQMA                                                  
      G2SQ54 = G2SQMA*D(5,4)                                            
      GSQMA  = G*SQMA                                                   
      GSQC2I = GSQMA*CT2I                                               
      GSQS2I = GSQMA*ST2I                                               
                                                                        
*                                                                                                             
* First compute the diagonal terms ...                                                                        
*                                                                                                             
      Q(1,1) = G2SQMA*(D(1,4)*D(1,4) + D(1,5)*D(1,5)*ST2I)              
      Q(2,2) = G2SQMA*(D(2,4)*D(2,4) + D(2,5)*D(2,5)*ST2I)              
      Q(4,4) = SQMA*CT2I*CT2I                                           
      Q(5,5) = G2SQ54*D(5,4) + SQMA*ST2I                                
                                                                        
*                                                                                                             
* Now the off-diagonal terms ...                                                                              
*                                                                                                             
      Q(2,1) = G2SQMA*(D(1,4)*D(2,4) + D(1,5)*D(2,5)*ST2I)              
      Q(4,1) = GSQC2I*D(1,4)                                            
      Q(5,1) = G2SQ54*D(1,4) + GSQS2I*D(1,5)                            
                                                                        
      Q(4,2) = GSQC2I*D(2,4)                                            
      Q(5,2) = G2SQ54*D(2,4) + GSQS2I*D(2,5)                            
                                                                        
      Q(5,4) = GSQC2I*D(5,4)                                            
                                                                        
      RETURN                                                            
      END                                                               
*