C++
C**********  P  D  O  T  B  ********** 
C
      LOGICAL FUNCTION PdotB(PART,GEOM,ELEMENT,Pup,Pdown,Direction)
C
C  Compute the projection of the polarization of the particle PART on the
C    direction of the magnetic field at its location, using geometry data 
C    from GEOM and ELEMENT.  The ortho and para projections are returned as 
C    Pup and Pdown (sum = 1), and Direction is a unit vector in the field 
C    direction.  If the field is 0, then Pup=Pdown=1/2, Direction=0, and 
C    the function is set to .FALSE.
C
C  P. A. Seeger, 19 August, 2002.
C    06 Nov 2009: changed calling sequence to BFIELD to include P [PAS]
C    14 Nov 2009: converted from STRUCTUREs to TYPEs, "." to "%" [PAS]
C
C  Definitions of TYPEs:
      IMPLICIT NONE
      INCLUDE 'mc_geom.inc'
C
C  Variables in calling sequence:
C    PART = record containing particle coordinates (input)
C    GEOM = complete geometry of the simulation (input)
C    ELEMENT = all element parameters (input)
C    Pup,Pdown = probabilities of projections of P on B (output)
C    Direction = unit vector in B direction (or 0 if |B|=0) (output)
      TYPE(PARTICLE)   PART
      TYPE(MC_GEOM)    GEOM
      TYPE(MC_ELEMENT) ELEMENT
      REAL*8 Pup, Pdown, Direction(3)
C
C  Externals
C    BFIELD    BREGION
C--
C  Local variables
      TYPE(PARTICLE) LocalPART
      REAL*8   B(3), BSQ, P(3), F, PdotN
      INTEGER  MREG, J
C
C     Advance neutron 10 ns to avoid roundoff errors
      LocalPART     = PART
      LocalPART%X   = PART%X   + 0.01*PART%VX
      LocalPART%Y   = PART%Y   + 0.01*PART%VY
      LocalPART%Z   = PART%Z   + 0.01*PART%VZ
      LocalPART%TOF = PART%TOF + 0.01
C
C     Find out if in a field region, and if there is a field
      BSQ = 0.d0
      CALL BREGION(LocalPART, GEOM, ELEMENT, MREG, J, F)
      IF (MREG .NE. 0) THEN
C        In a field region, find the field
         F = 0.d0
         P(1) = PART%PX
         P(2) = PART%PY
         P(3) = PART%PZ
         CALL BFIELD(LocalPART, P, F, GEOM, MREG, ELEMENT%NAME(MREG),   &
     &               ELEMENT%PARAM(ELEMENT%INDEX(MREG)), B, J, F)
         BSQ = B(1)**2 + B(2)**2 + B(3)**2
      END IF
C
      IF (BSQ .GT. 0.d0) THEN
C        Use the field direction to find projections of P
         F     = DSQRT(BSQ)
         PdotN = (P(1)*B(1) + P(2)*B(2) + P(3)*B(3))/F
C        Relative probabilities of spin up or spin down
         Pup   = (1.d0 + PdotN)/2.d0
         Pdown = (1.d0 - PdotN)/2.d0
         Direction(1) = B(1)/F
         Direction(2) = B(2)/F
         Direction(3) = B(3)/F
         PdotB = .TRUE.
      ELSE
C        No field, null result
         Pup   = 0.5d0
         Pdown = 0.5d0
         Direction(1) = 0.d0
         Direction(2) = 0.d0
         Direction(3) = 0.d0
         PdotB = .FALSE.
      END IF
C
      RETURN
      END
