! From Aurelio.Bay@ch.unil.ipn Mon Aug 14 14:28:42 1995
! Subject: vectors and 4-vectors
!----------------------------------------------------------------
!   Some vector and four-vector algebra.
!   The possibility of defining vector entities with specific algebra
!   makes the code more readable, compared to f77.
!   A couple of examples:
! 
!   scalar product of two vectors:
!   F90:
!        s = v1*v2
!   F77:
!        s = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)
!
!   vector product of two vectors:
!   F90:
!        V = A .X. B
!   F77:
!      v(1) = A(2)*B(3) - B(2)*A(3)
!      V(2) = A(3)*B(1) - B(3)*A(1)
!      V(3) = A(1)*B(2) - B(1)*A(2)
! 
!
! this example is of interests for physics :
!   mass of particle with four momentum p
!   F90:
!        m = sqrt(p*p)
!   F77:
!        s = p(1)**2 + p(2)**2 + p(3)**2
!        m = sqrt( p(4)**2 - s )
!
!===========================================================================
!===========================================================================
!  MODULE P3V_HANDLING
!   type defined:    P3V : vector with components (x, y, z)
!  MODULE P4V_HANDLING
!   type defined:    P4V : 4-vector with components (E, V)=(E, V%x, V%y, V%z)
!                          metrix=diag(1,-1,-1,-1)
!
!
!   First version:  19-AUG-1992   A. Bay          BAYA@VXCERN.CERN.CH
!                                 T. Monteiro
!   Last update:    14-Aug-1995
!===========================================================================
!  operations defined on vectors (v, v1, v2):
!     v1 + v2    add two vectors
!     v1 - v2    subtract
!     v1*v2      scalar product
!     v1 .X. v2  vector product
!     v*R        vector*real
!     v = R      3 components of the vector set to same real value R
!     v = R(3)   vector equated to real(3)
!     R(3) = v   real(3) equated to vector
!
! functions: 
!    v2 = norm_3v(v1)       vector v2 is v1 normalized to 1.
!    ang = find_ang_3v(v1,v2)  angle in [rad] between v1 and v2
!
!................................................................ 
!  operations on 4-vectors :
!     v1 + v2
!     v1 - v2 
!     v1*v2     scalar product  
!     v*R       4vector*real
!     v = R     4 components of the 4vector set to same real value R
!     v = R(4)  vector equated to real(4)
!     R(4) = v  real(4) equated to vector
!
! functions: 
!    find_ang(v1,v2) angle in rad between v1%V and v2%V
!
!===========================================================================

MODULE P3V_HANDLING

      TYPE P3V
        REAL X
        REAL Y
        REAL Z
      END TYPE P3V

      INTERFACE ASSIGNMENT(=)
        MODULE PROCEDURE P3V_EQU_P3V, P3V_EQU_R, P3V_EQU_V3, V3_EQU_P3V
      END INTERFACE
      INTERFACE OPERATOR(+)
        MODULE PROCEDURE ADD_P3V
      END INTERFACE
      INTERFACE OPERATOR(-)
        MODULE PROCEDURE SUB_P3V
      END INTERFACE
      INTERFACE OPERATOR(*)
        MODULE PROCEDURE MULT_P3V,MULT_P3V_R,MULT_R_P3V
      END INTERFACE
      INTERFACE OPERATOR(.X.)
        MODULE PROCEDURE VEC_P3V
      END INTERFACE


CONTAINS
  SUBROUTINE P3V_EQU_P3V(A,B)
      TYPE(P3V), INTENT(OUT) :: A
      TYPE(P3V), INTENT(IN) :: B
      A%X = B%X;      A%Y = B%Y;      A%Z = B%Z
  END SUBROUTINE P3V_EQU_P3V
  SUBROUTINE P3V_EQU_R(A,B)
      TYPE(P3V), INTENT(OUT) :: A
      REAL     , INTENT(IN)  :: B
      A%X = B;      A%Y = B;      A%Z = B
  END SUBROUTINE P3V_EQU_R
  SUBROUTINE P3V_EQU_V3(A,R)
      TYPE(P3V), INTENT(OUT):: A
      REAL, DIMENSION(3), INTENT(IN) :: R
      A%X  = R(1); A%Y = R(2); A%Z = R(3)
  END SUBROUTINE P3V_EQU_V3
  SUBROUTINE V3_EQU_P3V(R,A)
      TYPE(P3V), INTENT(IN):: A
      REAL, DIMENSION(3), INTENT(OUT) :: R
      R(1) = A%X; R(2) = A%Y; R(3) = A%Z
  END SUBROUTINE V3_EQU_P3V


  FUNCTION ADD_P3V(A,B)
      TYPE(P3V) ADD_P3V
      TYPE(P3V) , INTENT(IN):: A,B
      ADD_P3V%X = A%X + B%X;  ADD_P3V%Y = A%Y + B%Y;  ADD_P3V%Z = A%Z + B%Z
  END FUNCTION ADD_P3V
  FUNCTION SUB_P3V(A,B)
      TYPE(P3V) SUB_P3V
      TYPE(P3V) , INTENT(IN):: A,B
      SUB_P3V%X = A%X - B%X;  SUB_P3V%Y = A%Y - B%Y;  SUB_P3V%Z = A%Z - B%Z
  END FUNCTION SUB_P3V
  FUNCTION MULT_P3V(A,B)
      REAL      MULT_P3V
      TYPE(P3V) , INTENT(IN):: A,B
      MULT_P3V = A%X*B%X + A%Y*B%Y + A%Z*B%Z
  END FUNCTION MULT_P3V
  FUNCTION MULT_P3V_R(A,R)
      TYPE(P3V) MULT_P3V_R
      TYPE(P3V) , INTENT(IN) :: A
      REAL      , INTENT(IN) :: R
      MULT_P3V_R%X = A%X * R; MULT_P3V_R%Y = A%Y * R; MULT_P3V_R%Z = A%Z * R
  END FUNCTION MULT_P3V_R
  FUNCTION MULT_R_P3V(R,A)
      TYPE(P3V) MULT_R_P3V
      TYPE(P3V) , INTENT(IN) :: A
      REAL      , INTENT(IN) :: R
      MULT_R_P3V%X = A%X * R;  MULT_R_P3V%Y = A%Y * R; MULT_R_P3V%Z = A%Z * R
  END FUNCTION MULT_R_P3V
  FUNCTION VEC_P3V(A,B)
      TYPE(P3V) VEC_P3V
      TYPE(P3V) , INTENT(IN) :: A,B
      VEC_P3V%X=A%Y*B%Z-B%Y*A%Z
      VEC_P3V%Y=A%Z*B%X-B%Z*A%X
      VEC_P3V%Z=A%X*B%Y-B%X*A%Y
  END FUNCTION VEC_P3V

  FUNCTION NORM_3V(V)
   TYPE (P3V)              ::NORM_3V
   TYPE (P3V),INTENT(IN)   ::V
   REAL :: N
    N  = SQRT(V*V)
    IF (N>0) THEN
     NORM_3V%X = V%X/N
     NORM_3V%Y = V%Y/N
     NORM_3V%Z = V%Z/N
    ENDIF
  END  FUNCTION NORM_3V

!-------------------------------------- Angles
  FUNCTION FIND_ANG_3V(A,B)
    REAL                    ::FIND_ANG_3V
    TYPE (P3V), INTENT(IN)  ::A,B
    REAL                    ::AB,D,X
    AB = A*B;   D = (A*A)*(B*B)
    IF (D==0.) THEN
      FIND_ANG_3V = -99999.
    ELSE
      X = (AB/SQRT(D))
      IF (ABS(X)>1.) THEN
        WRITE(*,*)'Error in Find_ang', X
        WRITE(*,*)'A=',A
        WRITE(*,*)'B=',B
        IF (X<-1.) X = -1.
        IF (X>1.)  X = 1.
      ENDIF
    FIND_ANG_3V = ACOS(X)
    ENDIF
  END FUNCTION FIND_ANG_3V

END MODULE P3V_HANDLING

!=============================================================

MODULE P4V_HANDLING
USE P3V_HANDLING

      TYPE P4V
        TYPE(P3V) V
        REAL E
      END TYPE P4V

      INTERFACE ASSIGNMENT(=)
        MODULE PROCEDURE P4V_EQU_P4V, P4V_EQU_R,P4V_EQU_V4,V4_EQU_P4V
      END INTERFACE
      INTERFACE OPERATOR(+)
        MODULE PROCEDURE ADD_P4V
      END INTERFACE
      INTERFACE OPERATOR(-)
        MODULE PROCEDURE SUB_P4V
      END INTERFACE
      INTERFACE OPERATOR(*)
        MODULE PROCEDURE MULT_P4V,MULT_P4V_R,MULT_R_P4V
      END INTERFACE

      INTERFACE FIND_ANG
        MODULE PROCEDURE FIND_ANG_3V,FIND_ANG_4V
      END INTERFACE


CONTAINS
  SUBROUTINE P4V_EQU_P4V(A,B)
      TYPE(P4V) , INTENT(IN) :: B
      TYPE(P4V) , INTENT(OUT)  :: A
      A%V = B%V;  A%E = B%E
  END SUBROUTINE P4V_EQU_P4V
  SUBROUTINE P4V_EQU_R(A,B)
      TYPE(P4V) , INTENT(OUT) :: A
      REAL      , INTENT(IN) :: B
      A%V%X = B;      A%V%Y = B;      A%V%Z = B
      A%E = B
  END SUBROUTINE P4V_EQU_R
  SUBROUTINE P4V_EQU_V4(A,R)
      TYPE(P4V), INTENT(OUT):: A
      REAL, DIMENSION(4), INTENT(IN) :: R(4)
      A%V%X  = R(1);      A%V%Y  = R(2);      A%V%Z  = R(3)
      A%E    = R(4)
  END SUBROUTINE P4V_EQU_V4
  SUBROUTINE V4_EQU_P4V(R,A)
      TYPE(P4V), INTENT(IN) :: A
      REAL, DIMENSION(4), INTENT(OUT) :: R
      R(1) = A%V%X;      R(2) = A%V%Y;      R(3) = A%V%Z
      R(4) = A%E
  END SUBROUTINE V4_EQU_P4V


  FUNCTION ADD_P4V(A,B)
      TYPE(P4V) ADD_P4V
      TYPE(P4V) , INTENT(IN) :: A,B
      ADD_P4V%V = A%V + B%V
      ADD_P4V%E = A%E + B%E
  END FUNCTION ADD_P4V
  FUNCTION SUB_P4V(A,B)
      TYPE(P4V) SUB_P4V
      TYPE(P4V) , INTENT(IN) :: A,B
      SUB_P4V%V = A%V - B%V
      SUB_P4V%E = A%E - B%E
  END FUNCTION SUB_P4V
  FUNCTION MULT_P4V(A,B)
      REAL      MULT_P4V
      TYPE(P4V), INTENT(IN) :: A,B
      MULT_P4V = - (A%V*B%V) + A%E*B%E
  END FUNCTION MULT_P4V
  FUNCTION MULT_P4V_R(A,R)
      TYPE(P4V) MULT_P4V_R
      TYPE(P4V) , INTENT(IN) :: A
      REAL      , INTENT(IN) :: R
      MULT_P4V_R%V = A%V * R; MULT_P4V_R%E = A%E * R
  END FUNCTION MULT_P4V_R
  FUNCTION MULT_R_P4V(R,A)
      TYPE(P4V) MULT_R_P4V
      TYPE(P4V) , INTENT(IN) :: A
      REAL      , INTENT(IN) :: R
      MULT_R_P4V%V = A%V * R; MULT_R_P4V%E = A%E * R
  END FUNCTION MULT_R_P4V

! angles ----------------------------------------

  FUNCTION FIND_ANG_4V(A,B)
   REAL                    ::FIND_ANG_4V
   TYPE (P4V), INTENT(IN)  ::A,B
   FIND_ANG_4V = FIND_ANG_3V(A%V,B%V)
  END  FUNCTION FIND_ANG_4V

END MODULE P4V_HANDLING