GEANT321/erpremc/trs1s2.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/03/06 15:37:36  mclareni
   6 * Add geane321 source directories
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.49  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE TRS1S2 (PD1,RD1,PD2,RD2,H,CH,IERR,SP1,SP2
  13      1,                  DJ1,DK1,DJ2,DK2)
  14 C
  15 C *** TRANSFORMS ERROR MATRIX
  16 C     FROM        VARIABLES (1/P,V1',W1',V1,W1)
  17 C      TO         VARIABLES (1/P,V2',W2',V2,W2)
  18 C
  19 C     Authors: A. Haas and W. Wittek
  20 C
  21 C
  22 C *** PD1(3)    1/P,V1',W1'                             INPUT
  23 C     PD2(3)    1/P,V2',W2'                            OUTPUT
  24 C     H(3)      MAGNETIC FIELD                          INPUT
  25 C     RD1(15)   ERROR MATRIX IN 1/P,V1',W1',V1,W1       INPUT      (TRIANGLE)
  26 C     RD2(15)   ERROR MATRIX IN 1/P,V2',W2',V2,W2      OUTPUT      (TRIANGLE)
  27 C     CH        CHARGE OF PARTICLE                      INPUT
  28 C               CHARGE AND MAGNETIC FIELD ARE NEEDED
  29 C               FOR CORRELATION TERMS (V2',V1),(V2',W1),(W2',V1),(W2',W1)
  30 C               THESE CORRELATION TERMS APPEAR BECAUSE RD1 IS ASSUMED
  31 C               TO BE THE ERROR MATRIX FOR FIXED U1
  32 C               AND RD2 FOR FIXED U2
  33 C     SP1       SIGN OF U1-COMPONENT OF PARTICLE MOMENTUM     INPUT
  34 C     SP2       SIGN OF U2-COMPONENT OF PARTICLE MOMENTUM    OUTPUT
  35 C     DJ1(3)    UNIT VECTOR IN V1-DIRECTION
  36 C     DK1(3)    UNIT VECTOR IN W1-DIRECTION    OF SYSTEM 1
  37 C     DJ2(3)    UNIT VECTOR IN V2-DIRECTION
  38 C     DK2(3)    UNIT VECTOR IN W2-DIRECTION    OF SYSTEM 2
  39 C
  40 C     IERR      = 0    TRANSFORMATION OK
  41 C               = 1    MOMENTUM PERPENDICULAR TO U2-DIRECTION (V2',W2' NOT DEFIN
  42 C               = 2    MOMENTUM PERPENDICULAR TO X-AXIS
  43 C
  44 C
  45 #if !defined(CERNLIB_SINGLE)
  46       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  47       REAL   PD1,PD2,RD1,RD2,H,CH,SP1,SP2,DJ1,DK1,DJ2,DK2
  48 #endif
  49 #include "geant321/trcom3.inc"
  50       DIMENSION PD1(3),PD2(3),RD1(15),RD2(15),H(3),DJ1(3),DK1(3)
  51      +,DJ2(3),DK2(3),UN(3),VN(3),DI1(3),DI2(3),TVW1(3),TVW2(3)
  52 C
  53 #if defined(CERNLIB_SINGLE)
  54       DATA CFACT8 / 2.997925 E-4 /
  55 #endif
  56 #if !defined(CERNLIB_SINGLE)
  57       DATA CFACT8 / 2.997925 D-4 /
  58 #endif
  59 C
  60       IERR=0
  61       PM=PD1(1)
  62       TVW1(1)=1./SQRT(1.+PD1(2)**2+PD1(3)**2)
  63       IF(SP1.LT.0.) TVW1(1)=-TVW1(1)
  64       TVW1(2)=PD1(2)*TVW1(1)
  65       TVW1(3)=PD1(3)*TVW1(1)
  66 C
  67       DI1(1)=DJ1(2)*DK1(3)-DJ1(3)*DK1(2)
  68       DI1(2)=DJ1(3)*DK1(1)-DJ1(1)*DK1(3)
  69       DI1(3)=DJ1(1)*DK1(2)-DJ1(2)*DK1(1)
  70 C
  71       DO 5 I=1,3
  72          TN(I)=TVW1(1)*DI1(I)+TVW1(2)*DJ1(I)+TVW1(3)*DK1(I)
  73     5 CONTINUE
  74 C
  75       DI2(1)=DJ2(2)*DK2(3)-DJ2(3)*DK2(2)
  76       DI2(2)=DJ2(3)*DK2(1)-DJ2(1)*DK2(3)
  77       DI2(3)=DJ2(1)*DK2(2)-DJ2(2)*DK2(1)
  78 C
  79       TVW2(1)=TN(1)*DI2(1)+TN(2)*DI2(2)+TN(3)*DI2(3)
  80       TVW2(2)=TN(1)*DJ2(1)+TN(2)*DJ2(2)+TN(3)*DJ2(3)
  81       TVW2(3)=TN(1)*DK2(1)+TN(2)*DK2(2)+TN(3)*DK2(3)
  82 C
  83       IF(TVW2(1).EQ.0.) GO TO 901
  84       TR=1./TVW2(1)
  85       SP2=1.
  86       IF(TVW2(1).LT.0.) SP2=-1.
  87       PD2(1)=PD1(1)
  88       PD2(2)=TVW2(2)*TR
  89       PD2(3)=TVW2(3)*TR
  90 C
  91       COSL=SQRT(ABS(1.-TN(3)**2))
  92       IF(COSL.EQ.0.) GO TO 902
  93       COSL1=1./COSL
  94       UN(1)=-TN(2)*COSL1
  95       UN(2)=TN(1)*COSL1
  96       UN(3)=0.
  97 C
  98       VN(1)=-TN(3)*UN(2)
  99       VN(2)=TN(3)*UN(1)
 100       VN(3)=COSL
 101 C
 102       UJ1=UN(1)*DJ1(1)+UN(2)*DJ1(2)+UN(3)*DJ1(3)
 103       UK1=UN(1)*DK1(1)+UN(2)*DK1(2)+UN(3)*DK1(3)
 104       VJ1=VN(1)*DJ1(1)+VN(2)*DJ1(2)+VN(3)*DJ1(3)
 105       VK1=VN(1)*DK1(1)+VN(2)*DK1(2)+VN(3)*DK1(3)
 106 C
 107       UJ2=UN(1)*DJ2(1)+UN(2)*DJ2(2)+UN(3)*DJ2(3)
 108       UK2=UN(1)*DK2(1)+UN(2)*DK2(2)+UN(3)*DK2(3)
 109       VJ2=VN(1)*DJ2(1)+VN(2)*DJ2(2)+VN(3)*DJ2(3)
 110       VK2=VN(1)*DK2(1)+VN(2)*DK2(2)+VN(3)*DK2(3)
 111 C
 112       J=0
 113       DO 10 I=1,5
 114          DO 4 K=I,5
 115             J=J+1
 116             A(I,K)=0.
 117             A(K,I)=0.
 118             S(J)=RD1(J)
 119     4    CONTINUE
 120    10 CONTINUE
 121 C
 122       IF(CH.EQ.0.) GO TO 6
 123       HA=SQRT(H(1)**2+H(2)**2+H(3)**2)
 124       HAM=HA*PM
 125       IF(HAM.EQ.0.) GO TO 6
 126       HM=CH/HA
 127 C
 128       Q=-HAM*CFACT8
 129       TT=-Q*TR**3
 130       SJ1I2=DJ1(1)*DI2(1)+DJ1(2)*DI2(2)+DJ1(3)*DI2(3)
 131       SK1I2=DK1(1)*DI2(1)+DK1(2)*DI2(2)+DK1(3)*DI2(3)
 132       SK2U=DK2(1)*UN(1)+DK2(2)*UN(2)+DK2(3)*UN(3)
 133       SK2V=DK2(1)*VN(1)+DK2(2)*VN(2)+DK2(3)*VN(3)
 134       SJ2U=DJ2(1)*UN(1)+DJ2(2)*UN(2)+DJ2(3)*UN(3)
 135       SJ2V=DJ2(1)*VN(1)+DJ2(2)*VN(2)+DJ2(3)*VN(3)
 136 C
 137       SINZ=-(H(1)*UN(1)+H(2)*UN(2)+H(3)*UN(3))*HM
 138       COSZ= (H(1)*VN(1)+H(2)*VN(2)+H(3)*VN(3))*HM
 139       A(2,4)=-TT*SJ1I2*(SK2U*SINZ-SK2V*COSZ)
 140       A(2,5)=-TT*SK1I2*(SK2U*SINZ-SK2V*COSZ)
 141       A(3,4)= TT*SJ1I2*(SJ2U*SINZ-SJ2V*COSZ)
 142       A(3,5)= TT*SK1I2*(SJ2U*SINZ-SJ2V*COSZ)
 143 C
 144     6 A(1,1)=1.
 145       A(4,4)=TR*(UJ1*VK2-VJ1*UK2)
 146       A(4,5)=TR*(UK1*VK2-VK1*UK2)
 147       A(5,4)=TR*(VJ1*UJ2-UJ1*VJ2)
 148       A(5,5)=TR*(VK1*UJ2-UK1*VJ2)
 149 C
 150       TS=TR*TVW1(1)
 151       A(2,2)=A(4,4)*TS
 152       A(2,3)=A(4,5)*TS
 153       A(3,2)=A(5,4)*TS
 154       A(3,3)=A(5,5)*TS
 155 C
 156       CALL SSMT5T(A,S,S)
 157 C
 158       J=0
 159       DO 25 I=1,5
 160          DO 20 K=I,5
 161             J=J+1
 162             RD2(J)=S(J)
 163    20    CONTINUE
 164    25 CONTINUE
 165 C
 166       RETURN
 167 C
 168 C *** ERROR EXITS
 169 C
 170   901 IERR=1
 171       GO TO 910
 172   902 IERR=2
 173   910 RETURN
 174       END
 175 *