FFTRF COMMENTS
*-- Author :    Stephen Burke   26/11/96
      SUBROUTINE FFTRF(DZ,Z1,S1,S2,BZEND)
**********************************************************************                                        
*                                                                    *                                        
*  Transform track vector (S1) through DZ to (S2) in magnetic field. *                                        
*                                                                    *                                        
*  The field in kG must be supplied through call to GUFLD, and       *                                        
*  is assumed to satisfy   |Bz| >> SQRT( Bx**2 + By**2 )  and        *                                        
*                          |B2-B1| << |B1|                           *                                        
*  ie that the field is constant and aligned with the +/- Z axis.    *                                        
*  A helix track model is used to compute S2. First order            *                                        
*  corrections are applied if necessary:                             *                                        
*                                                                    *                                        
*       IAPROX = 1 -> Straight line extrapolation (unused)           *                                        
*       IAPROX = 2 -> As 1 for x and y, but Delta phi varying        *                                        
*       IAPROX = 3 -> Helix model. Bz=constant, Bx = By = 0 (unused) *                                        
*       IAPROX = 4 -> 1st order corrections assuming Bx and By are   *                                        
*                     small but constant across DZ.                  *                                        
*       IAPROX = 5 -> as 4 plus 1st order corrections for change     *                                        
*                     in field along trajectory.                     *                                        
*                                                                    *                                        
* INPUT;   DZ    = step length in cm (can be negative)               *                                        
*          Z1    = Z coordinate of starting point in cm              *                                        
*          S1    = (X,Y,q/P,Tan(Theta),Phi) = state vector at Z1     *                                        
*                                                                    *                                        
* OUTPUT;  S2    = state vector at Z1 + DZ                           *                                        
*          BZEND = Bz at track end                                   *                                        
*                                                                    *                                        
**********************************************************************                                        
*KEEP,FKPIDP.                                                                                                 
*KEND.                                                                                                        
**********************************************************************                                        
*                                                                                                             
* Thresholds for various approximations (*** UNOPTIMISED ***)                                                 
*                                                                                                             
* Max z step over which field is assumed constant                                                             
* If Dphi CALL GUFLD(R,BR)                                                                              
* Work out Delta phi                                                                                          
* Choose the order of approximation                                                                           
* Must consider change in field over DZ                                                                       
* Must take helix (x and y affected by dphi)                                                                  
* Ignore change in phi for DX and DY                                                                          
* Small angle approximations for cos and sin (to order DPHI**2)                                               
* Full expressions                                                                                            
* Remember values                                                                                             
* Change in x and y                                                                                           
*                                                                                                             
* Apply 1st order corrections assuming that Bx and By are                                                     
* small but non-zero ........ (Bx,By,Bz) is still assumed                                                     
* to be constant across DZ ......                                                                             
*                                                                                                             
* If theta = 0, phi is undefined, so CON3 can be set to any value                                             
* Are corrections for field variation required?                                                               
* Get field at end point ...                                                                                  
      CALL GUFLD(R2,BR)                                                                                
*                                                                                                             
* DB is the mean difference from the assumed field along the trajectory,                                      
* assuming the field components vary linearly with Z over DZ.                                                 
* ie B and dB/dZ are non-zero but  d2B/dZ2  etc are assumed zero.                                             
*                                                                                                             
* Stop values getting too large (i.e. stop overflows)