*
* $Id$
*
* $Log$
* Revision 1.1  1996/03/06 15:37:34  mclareni
* Initial revision
*
*
*      MODIFIED FOR MAGFIELD GRADIENT 02-APR-82
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.49  by  S.Giani
*-- Author :
C
      SUBROUTINE TRPROP(X1,P1,H1,X2,P2,H2,CH,XL,R,MVAR,IFLAG,ITRAN,IERR)
C
C *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
C     ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
C     AND (HX,HY,HZ) ARE RATHER CONSTANT. DELTA(PHI) MUST NOT BE TOO LARGE
C
C     Authors: A. Haas and W. Wittek
C
C *** IFLAG  =  -1   INITIALIZATION, TRANSFORMATION OF ERROR MATRIX FROM
C                    EXTERNAL TO SC VARIABLES
C            =   0   ERROR PROPAGATION FROM X1 TO X2
C            =   1   TRANSFORMATION OF ERROR MATRIX FROM SC TO
C                    EXTERNAL VARIABLES
C
C     ITRAN          USED FOR IFLAG = 0 OR 1 ONLY
C            =   0   TRANSFORMATION MATRIX IS UPDATED ,BUT ERROR MATRIX IS NOT
C                    TRANSFORMED
C           =    1   TRANSF. MATRIX IS UPDATED  AND ERROR MATRIX IS TRANSFORMED
C
C     MVAR           SPECIFIES TYPE OF EXTERNAL VARIABLES
C            =   0   ( 1/P,LAMBDA,PHI,YT, ZT ;   SC   )
C            =   1   ( 1/P,  Y',  Z',  Y,  Z ; SPLINE )
C
C *** X1, P1, H1     X,Y,Z COMPONENTS OF POSITION, MOMENTUM AND MAGNETIC   INPUT
C                    FIELD VECTOR/GRADIENT AT STARTING POINT OF INTERVAL
C     X2, P2, H2     ......  AT END POINT OF INTERVAL                      INPUT
C     CH             CHARGE OF PARTICLE                                    INPUT
C     XL             PATHLENGTH FROM X1 TO X2   ( NEGATIVE IF OPPOSITE
C                    TO ACTUAL MOVEMENT OF PARTICLE )                      INPUT
C     R              ERROR MATRIX  (TRIANGLE)                       INPUT/OUTPUT
C     B              5 * 5 TRANSFORMATION MATRIX FOR ERRORS IN
C                    SC VARIABLES                                         OUTPUT
C
C *** IERR   =  1    ILLEGAL VALUE OF MVAR                                OUTPUT
C               2    MOMENTUM IS ZERO
C               3    H*ALFA/P AT X1 AND X2 DIFFER TOO MUCH
C                    OR DELTA PHI IS TOO LARGE
C               4    PARTICLE MOVES IN Z - DIRECTION
C
#if !defined(CERNLIB_SINGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL  X1,P1,H1,X2,P2,H2,R,CH,PS,PC,XL,SPX
#endif
#include "geant321/trcom3.inc"
#include "geant321/gcunit.inc"
      DIMENSION X1(3),P1(3),H1(9),X2(3),P2(3),H2(9)
      DIMENSION R(15),PS(3),PC(3),HN(9)
C
      SAVE INIT,DELHP6,CFACT8
*
      DATA INIT/0/
#if defined(CERNLIB_SINGLE)
      DATA DELHP6/300./,DELFI6/0.1/
      DATA CFACT8 / 2.997925 E-4 /
#endif
#if !defined(CERNLIB_SINGLE)
      DATA DELHP6/300.D0/,DELFI6/0.1D0/
      DATA CFACT8 / 2.997925 D-4 /
#endif
C
      IERR=0
      IF(IFLAG) 10, 20, 80
C
C *** TRANSFORM ERROR MATRIX FROM EXTERNAL TO INTERNAL VARIABLES;
C
   10 NEW=1
      IF(MVAR.NE.1) GO TO 11
      PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
      IF(PA1.EQ.0.) GO TO 902
      PS(1)=1./PA1
      IF(P1(1).EQ.0.) GO TO 904
      PS(2)=P1(2)/P1(1)
      PS(3)=P1(3)/P1(1)
      SPX=1.
      IF(P1(1).LT.0.) SPX=-1.
      CALL TRSPSC(PS,R,PC,R,H1,CH,IERR,SPX)
      GO TO 19
C
   11 IF(MVAR.NE.0) GO TO 901
   19 GO TO 900
C
C *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
 
C
   20 PA1=SQRT(P1(1)**2+P1(2)**2+P1(3)**2)
      PA2=SQRT(P2(1)**2+P2(2)**2+P2(3)**2)
      IF(PA1*PA2.EQ.0.) GO TO 902
      PM1=1./PA1
      PM2=1./PA2
C
      TN(1)=P1(1)+P2(1)
      TN(2)=P1(2)+P2(2)
      TN(3)=P1(3)+P2(3)
      PM12=1./SQRT(TN(1)**2+TN(2)**2+TN(3)**2)
      TN(1)=TN(1)*PM12
      TN(2)=TN(2)*PM12
      TN(3)=TN(3)*PM12
C
      SINL=TN(3)
      COSL=SQRT(ABS(1.-SINL**2))
      IF(COSL.EQ.0.) GO TO 904
      COSL1=1./COSL
      SINP=TN(2)*COSL1
      COSP=TN(1)*COSL1
C
C *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
C *** NEUTRAL PARTICLE OR FIELDFREE REGION
C
      DO 26 I=1,5
         DO 15 K=1,5
            A(I,K)=0.
   15    CONTINUE
         A(I,I)=1.
   26 CONTINUE
      A(4,3)=XL*COSL
      A(5,2)=XL
C
      IF(CH.EQ.0.) GO TO 45
      HA1=SQRT(H1(1)**2+H1(2)**2+H1(3)**2)
      HA2=SQRT(H2(1)**2+H2(2)**2+H2(3)**2)
      HAM1=HA1*PM1
      HAM2=HA2*PM2
      HAMX = MAX(HAM1,HAM2)
      IF(HAMX.EQ.0.) GO TO 45
C
C *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
C
      PM12=(PM1+PM2)*0.5
      P12=1./(2.*PM12)
      HN(1)=(H1(1)*PM1+H2(1)*PM2)*P12*CH*CFACT8
      HN(2)=(H1(2)*PM1+H2(2)*PM2)*P12*CH*CFACT8
      HN(3)=(H1(3)*PM1+H2(3)*PM2)*P12*CH*CFACT8
      HN(4)=(H1(4)*PM1+H2(4)*PM2)*P12*CH*CFACT8
      HN(5)=(H1(5)*PM1+H2(5)*PM2)*P12*CH*CFACT8
      HN(6)=(H1(6)*PM1+H2(6)*PM2)*P12*CH*CFACT8
      HN(7)=(H1(7)*PM1+H2(7)*PM2)*P12*CH*CFACT8
      HN(8)=(H1(8)*PM1+H2(8)*PM2)*P12*CH*CFACT8
      HN(9)=(H1(9)*PM1+H2(9)*PM2)*P12*CH*CFACT8
C
      B0=HN(1)*COSP+HN(2)*SINP
      B2=-HN(1)*SINP+HN(2)*COSP
      B3=-B0*SINL+HN(3)*COSL
      TGL=SINL*COSL1
C
C
C *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
C     AND WHETHER CHANGE OF TRACK DIRECTION DUE TO MAG.FIELD IS TOO LARGE
C
C
      IF(HA2.EQ.0.) GO TO 29
 
      GAM=(H2(1)*TN(1)+H2(2)*TN(2)+H2(3)*TN(3))/HA2
      GO TO 28
   29 GAM=(H1(1)*TN(1)+H1(2)*TN(2)+H1(3)*TN(3))/HA1
   28 CONTINUE
      ALFA=SQRT(ABS(1.-GAM**2))
C
      DH2=(H1(1)*PM1-H2(1)*PM2)**2+
     1    (H1(2)*PM1-H2(2)*PM2)**2+
     1    (H1(3)*PM1-H2(3)*PM2)**2
      IF(DH2*ALFA**2.GT.DELHP6**2) GO TO 903
      ALFAQ=-ALFA*CFACT8*(HAM1+HAM2)*0.5
      DFI=ABS(XL*ALFAQ)
      IF(DFI.GT.DELFI6) GO TO 903
C
C *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
C *** TAKING INTO ACCOUNT  FIELD GRADIENT PERPENDICULAR TO TRACK
C
   30 COSP2=COSP*COSP
      SINP2=SINP*SINP
      COSIP=COSP*SINP
C
      G22=SINP2*HN(9)+COSP2*HN(8)-2.0*COSIP*HN(7)
      G33=SINL*SINL*(COSP2*HN(9)+SINP2*HN(8)+2.0*COSIP*HN(7))
     ++COSL*(COSL*HN(6)-2.0*SINL*(COSP*HN(4)+SINP*HN(5)))
      G23=SINL*(COSIP*(HN(9)-HN(8))+(SINP2-COSP2)*HN(7))
     ++COSL*(COSP*HN(5)-SINP*HN(4))
C
      A(2,1)=XL*B2
      A(2,3)=-B0*XL*PM12
      A(2,4)=(B2*B3*PM12+G22)*XL*PM12
      A(2,5)=(-B2*B2*PM12+G23)*XL*PM12
C
      A(3,1)=-XL*B3*COSL1
      A(3,2)=B0*XL*PM12*COSL1**2
      A(3,3)=1.+TGL*B2*XL*PM12
      A(3,4)=(-B3*B3*PM12-G23)*XL*PM12*COSL1
      A(3,5)=(B3*B2*PM12-G33)*XL*PM12*COSL1
C
      A(4,5)=-B3*TGL*XL*PM12
      A(5,4)=B3*TGL*XL*PM12
C
   45 CONTINUE
C
C *** NEW = 0  TRANSFORMATION MATRIX IS UPDATED
C           1  TRANSFORMATION MATRIX IS INITIALIZED
C
      IF(NEW.EQ.0) GO TO 23
      NEW=0
      DO 25 I=1,5
         DO 24 K=1,5
            B(I,K)=A(I,K)
   24    CONTINUE
   25 CONTINUE
      GO TO 27
   23 CONTINUE
C
      CALL XMM55(A,B,B)
C
   27 CONTINUE
   80 IF(ITRAN.EQ.0) GO TO 90
C
C
      J=0
      DO 22 I=1,5
         DO 21 K=I,5
            J=J+1
            S(J)=R(J)
   21    CONTINUE
   22 CONTINUE
C
C
C *** TRANSFORM ERROR MATRIX
C
      CALL SSMT5T(B,S,S)
C
      NEW=1
 
      J=0
      DO 41 I=1,5
         DO 40 K=I,5
            J=J+1
            R(J)=S(J)
   40    CONTINUE
   41 CONTINUE
C
   90 IF(IFLAG.LE.0) GO TO 900
C
C
C *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
C
C
      NEW=1
      IF(MVAR.NE.1) GO TO 91
      PC(1)=PM2
      PC(2)=ASIN(P2(3)*PC(1))
      IF (ABS (P2(1)) .LT. 1.E-30) P2(1) = 1.E-30
      PC(3) = ATAN2 (P2(2),P2(1))
      CALL TRSCSP(PC,R,PS,R,H2,CH,IERR,SPX)
      GO TO 900
C
   91 IF(MVAR.NE.0) GO TO 901
      GO TO 900
C
C *** ERROR EXITS
C
  901 IERR=1
      GO TO 999
  902 IERR=2
      GO TO 999
  903 IERR=3
      IF(INIT.NE.0) GO TO 30
C     WRITE (LOUT, 998) DH2,DFI,ALFA,XL
C 998 FORMAT(1H0,48H *** S/R TRPROP   DELTA(H*ALFA/P)  OR DELTA(PHI),5X
C    1,22HEXCEEDS TOLERANCE     /1H0,4E12.5//16H **********    ,
C    251HATTENTION !   NO FURTHER WARNINGS WILL BE GIVEN    ///)
      INIT=1
      GO TO 30
  904 IERR=4
  999 WRITE (LOUT, 1000) IERR
 1000 FORMAT(1H ,' *** S/R ERPROP   IERR =',I5)
C
  900 RETURN
      END
