*-- 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 *
* *
**********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
REAL R(3),R2(3),B2(3),DB(3),BR(3),BZEND
DIMENSION S1(5),S2(5)
*KEEP,FKPIDP.
DOUBLE PRECISION PI,TWOPI,PIBY2
PARAMETER (PI=3.141592653589793238)
PARAMETER (TWOPI=PI*2.0D0,PIBY2=PI/2.0D0)
*KEND.
PARAMETER (BFACT=2.997886E-4)
**********************************************************************
*
* Thresholds for various approximations (*** UNOPTIMISED ***)
*
* Max z step over which field is assumed constant
PARAMETER (DZMAX=1.0D1)
* If Dphi CALL GUFLD(R,BR)
IF (ABS(BR(3)).LT.0.01) BR(3) = SIGN(0.01,BR(3))
ZR = Z1
ENDIF
* Work out Delta phi
DPHI = -DZSECT*BR(3)*BFACT*S1(3)
S2(5) = S2(5) + DPHI
* Choose the order of approximation
IF (DABS(DZ).GT.DZMAX) THEN
* Must consider change in field over DZ
IAPROX = 5
ELSEIF (DABS(DZ*DPHI).GT.DXYMIN) THEN
* Must take helix (x and y affected by dphi)
IAPROX = 4
ELSE
* Ignore change in phi for DX and DY
S2(1) = S2(1) + DZ*CPHI1*S1(4)
S2(2) = S2(2) + DZ*SPHI1*S1(4)
GOTO 9100
ENDIF
IF (DABS(DPHI).LT.DPHMAX) THEN
* Small angle approximations for cos and sin (to order DPHI**2)
DPHI2 = 0.5D0*DPHI
CDPHI = 1.D0 - DPHI2*DPHI
SPHI2 = SPHI1*CDPHI + CPHI1*DPHI
CPHI2 = CPHI1*CDPHI - SPHI1*DPHI
DXOTT = DZ*(CPHI1 - SPHI1*DPHI2)
DYOTT = DZ*(SPHI1 + CPHI1*DPHI2)
ELSE
* Full expressions
CPHI2 = DCOS(S2(5))
SPHI2 = DSIN(S2(5))
FAC = DZ/DPHI
DXOTT = FAC*(SPHI2 - SPHI1)
DYOTT = FAC*(CPHI1 - CPHI2)
* Remember values
PHIR = S2(5)
CPHIR = CPHI2
SPHIR = SPHI2
ENDIF
* Change in x and y
S2(1) = S2(1) + DXOTT*S1(4)
S2(2) = S2(2) + DYOTT*S1(4)
*
* 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 ......
*
DPHIOB = DZSECT*S1(3)*BFACT
CON1 = DPHIOB*DZ*0.5D0
CON2 = DPHIOB*OPTSQ
IF (S1(4).NE.0.D0) THEN
CON3 = DPHIOB/S1(4)
ELSE
* If theta = 0, phi is undefined, so CON3 can be set to any value
CON3 = 0.D0
ENDIF
S2(1) = S2(1) - CON1*BR(2)
S2(2) = S2(2) + CON1*BR(1)
S2(4) = S2(4) + CON2*(BR(1)*SPHI1 - BR(2)*CPHI1)
S2(5) = S2(5) + CON3*(BR(1)*CPHI1 + BR(2)*SPHI1)
* Are corrections for field variation required?
IF (IAPROX.LT.5) GOTO 9000
* Get field at end point ...
B2(1) = BR(1)
B2(2) = BR(2)
B2(3) = BR(3)
R2(1) = S2(1)
R2(2) = S2(2)
R2(3) = Z1 + DZ
ZR = Z1 + DZ
CALL GUFLD(R2,BR)
IF (ABS(BR(3)).LT.0.01) BR(3) = SIGN(0.01,BR(3))
DB(1) = 0.5*(BR(1) - B2(1))
DB(2) = 0.5*(BR(2) - B2(2))
DB(3) = 0.5*(BR(3) - B2(3))
*
* 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.
*
DBZTT = DB(3)*S1(4)
S2(1) = S2(1) + CON1*(DBZTT*SPHI1 - DB(2))
S2(2) = S2(2) + CON1*(DB(1) - DBZTT*CPHI1)
S2(4) = S2(4) + CON2*(DB(1)*SPHI1 - DB(2)*CPHI1)
S2(5) = S2(5) + CON3*(DB(1)*CPHI1 + DB(2)*SPHI1) - DPHIOB*DB(3)
9000 CONTINUE
IF (S2(4).LT.0.D0) THEN
S2(4) = DABS(S2(4))
S2(5) = S2(5) + PI
ENDIF
9100 CONTINUE
IF (DABS(S2(5)).GT.TWOPI) S2(5)=DMOD(S2(5),TWOPI)
IF (S2(5).LT.0.D0) S2(5) = S2(5) + TWOPI
* Stop values getting too large (i.e. stop overflows)
IF (ABS(S2(1)).GT.1.D4) S2(1) = SIGN(1.D4,S2(1))
IF (ABS(S2(2)).GT.1.D4) S2(2) = SIGN(1.D4,S2(2))
IF (S2(4).GT.1.D6) S2(4) = 1.D6
BZEND = BR(3)
RETURN
END