MODULE Space
   USE f90_kind
   PRIVATE 
   
! 2D coordinate
! point = Coordinate(X,Y)
! X = Pitch(point)
! Y = Length(point)
   PUBLIC :: Coordinate,   Pitch,    Length

! straight line
! track = Line(X,Phi)
! track = Line(point,phi)
! track = Line(pointA,pointB)
! track = Line(XYtrack,Y)
! getLine(track,X,Phi)
! X = getLineX(track,Y)
! Y = getLineY(track,X)
   PUBLIC :: Line,       getLine, getLineX, getLineY

! curve with radius
! XYtrack = XYTrack(track,invRad)
! getXYTrack(myTrack,X,Phi,invRad)
   PUBLIC :: XYTrack,    getXYTrack
   
! pad-types, pad can be rectangle or grid
! point = getCentre(pad)
   PUBLIC :: getCentre

! distance between two objects (1st argument has higher dimension):
! coordinate, line, rectangle, grid
! dist = Distance(A,B)
   PUBLIC :: Distance

! FillLoc(cX,cY,sX,sY,rect)
! FillLoc(cX,cY,sX,sY,mapNumber,grid)   
! AllocBitMap(NumMaps)   
! FillBitMap(mapNumber,bits)
   PUBLIC :: FillLoc, AllocBitMap, FillBitMap
! printGrid(grid)
   PUBLIC :: printGrid

! integral of Line over a pad (=rectangle or grid)   
! integral = intTrackPad(track,width,pad)
   PUBLIC :: intTrackPad
   
! just for fun
! == (.EQ.) defined for coordinate and line   
   PUBLIC :: OPERATOR(+), OPERATOR(-), OPERATOR(*), OPERATOR(==), ASSIGNMENT(=)
   
   PRIVATE :: assigna2i,assigna2r,assigna2d
   PRIVATE :: bitmap,recty
   PRIVATE :: FillRect,FillGrid
   PRIVATE :: getCentreRect,getCentreGrid,getCoordinate
   PRIVATE :: CoordinateD,CoordinateR,CoordinateRD,CoordinateDR
   PRIVATE :: assignCoord,assignCoordR,addCoord,subCoord,scaleCoordR,scaleCoordD
   PRIVATE :: prodCoord,sameCoord
   PRIVATE :: LineAngleD,LineAngleR,LinePktR,LinePktD,Line2Pkt
   PRIVATE :: LineXYTrackD,LineXYTrackR,XYTrackD,XYTrackR,sameLine
   PRIVATE :: getLineXD,getLineXR,getLineYD,getLineYR,sameXYTrack
   PRIVATE :: intTrackRect,eta,Corners,intTrackGrid
   PRIVATE :: DistXRect2Point,DistYRect2Point,DistXRect2Rect,DistYRect2Rect
   PRIVATE :: DistPoint2Point,SignDistLine2Point,DistLine2Point
   PRIVATE :: DistRect2Point,DistRect2Line,DistRect2Rect
   PRIVATE :: DistGrid2Point,DistGrid2Line,DistGrid2Grid
   
! if f90_kind is not available, select depending on compiler option
! kind parameter    sequential     byte
! integer: int8        1             1
!          int16       2             2
!          int32       3             4
!          int64       4             8
! real:    single      1             4
!          double      2             8
!          quad        3            16
! logical: byte        1             1
!          twobyte     2             2
!          word        3             4
! size of bitmap, can be adjusted
   INTEGER,PARAMETER,PRIVATE :: mapPrec=int32

   INTEGER,PARAMETER,PUBLIC :: spPrec=double
   REAL(KIND=spPrec),PRIVATE :: sqrt2=1.414213562 , sqrt2pi=2.506628275

! ###################### !
! all measurements in mm !
! ###################### !

   TYPE,PUBLIC :: Coordinate_type
      PRIVATE
      REAL(KIND=spPrec) :: X , Y
   END TYPE Coordinate_type   
   TYPE,PUBLIC :: Line_type
      PRIVATE
      REAL(KIND=spPrec) :: X0 , Phi
   END TYPE Line_type
   TYPE,PUBLIC :: XYTrack_type
      PRIVATE
! Slope: track at Y=0
! invRadius > 0 : left part of circle
! invRadius < 0 : right part of circle
! |invRadius| < 1E-6 : straight Track
! -> sagitta of 125um for 1m long track
      TYPE(Line_type)   :: Slope  
      REAL(KIND=spPrec) :: invRadius
   END TYPE XYTrack_type
 
   TYPE,PUBLIC :: Rectangle_type
      PRIVATE
      TYPE(Coordinate_type) :: centre , size
   END TYPE Rectangle_type
 
   TYPE,PUBLIC :: Grid_type
      PRIVATE
      TYPE(Coordinate_type) :: centre , size
      INTEGER :: mapID
   END TYPE Grid_type
   TYPE,PRIVATE :: bitmap_type
      PRIVATE
      INTEGER :: nx,ny
      INTEGER(KIND=mapPrec),ALLOCATABLE,DIMENSION(:) :: bmap
   END TYPE bitmap_type

   TYPE(bitmap_type),PRIVATE,SAVE,ALLOCATABLE,DIMENSION(:) :: gridMap
  
   INTERFACE ASSIGNMENT (=)
      MODULE PROCEDURE assignCoord,assignCoordR &
                      ,assigna2i,assigna2r,assigna2d
   END INTERFACE 
   INTERFACE OPERATOR (+)
      MODULE PROCEDURE addCoord
   END INTERFACE 
   INTERFACE OPERATOR (-)
      MODULE PROCEDURE subCoord
   END INTERFACE 
   INTERFACE OPERATOR (*)
      MODULE PROCEDURE prodCoord, scaleCoordR, scaleCoordD
   END INTERFACE 
   INTERFACE OPERATOR (==)
      MODULE PROCEDURE sameCoord, sameLine, sameXYtrack
   END INTERFACE 
   
   INTERFACE Coordinate
      MODULE PROCEDURE CoordinateD,CoordinateR,CoordinateDR,CoordinateRD
   END INTERFACE 
   
   INTERFACE Line
      MODULE PROCEDURE LineAngleD,LineAngleR,LinePktR,LinePktD,&
             Line2Pkt,LineXYTrackR,LineXYTrackD
   END INTERFACE 
      
   INTERFACE getLineX
      MODULE PROCEDURE getLineXR,getLineXD
   END INTERFACE 
      
   INTERFACE getLineY
      MODULE PROCEDURE getLineYR,getLineYD
   END INTERFACE 
      
   INTERFACE XYTrack
      MODULE PROCEDURE XYTrackR,XYTrackD
   END INTERFACE 
      
   INTERFACE Distance
      MODULE PROCEDURE DistPoint2Point,DistLine2Point,&
             DistRect2Point,DistRect2Line,DistRect2Rect,&
             DistGrid2Point,DistGrid2Line,DistGrid2Grid
   END INTERFACE 

! need soubroutines since the return type of a function 
! does not belong to the signature
   INTERFACE FillLoc
      MODULE PROCEDURE FillRect ,FillGrid
   END INTERFACE 
 
   INTERFACE getCentre
      MODULE PROCEDURE getCentreRect,getCentreGrid
   END INTERFACE 

   INTERFACE intTrackPad
      MODULE PROCEDURE intTrackRect,intTrackGrid
   END INTERFACE 
   
   CONTAINS

! *******************************************************************
   PURE SUBROUTINE assigna2i (i,a) 
      INTEGER,INTENT(OUT) :: i
      INTEGER,DIMENSION(1),INTENT(IN) :: a
      
      i = a(1)
   END SUBROUTINE assigna2i
   PURE SUBROUTINE assigna2r (r,a) 
      REAL,INTENT(OUT) :: r
      REAL,DIMENSION(1),INTENT(IN) :: a
      
      r = a(1)
   END SUBROUTINE assigna2r
   PURE SUBROUTINE assigna2d (d,a) 
      REAL(KIND=double),INTENT(OUT) :: d
      REAL(KIND=double),DIMENSION(1),INTENT(IN) :: a
      
      d = a(1)
   END SUBROUTINE assigna2d

! *******************************************************************

   PURE FUNCTION CoordinateD (X,Y) RESULT(point)
      TYPE(Coordinate_type) :: point
      REAL(KIND=spPrec),INTENT(IN) :: X,Y
      
      point%X = X
      point%Y = Y
   END FUNCTION CoordinateD
   PURE FUNCTION CoordinateR (X,Y) RESULT(point)
      TYPE(Coordinate_type) :: point
      REAL,INTENT(IN) :: X,Y
      
      point = CoordinateD(REAL(X,spPrec),REAL(Y,spPrec))
   END FUNCTION CoordinateR
   PURE FUNCTION CoordinateRD (X,Y) RESULT(point)
      TYPE(Coordinate_type) :: point
      REAL,INTENT(IN) :: X
      REAL(KIND=spPrec),INTENT(IN) :: Y
      
      point = CoordinateD(REAL(X,spPrec),Y)
   END FUNCTION CoordinateRD
   PURE FUNCTION CoordinateDR (X,Y) RESULT(point)
      TYPE(Coordinate_type) :: point
      REAL(KIND=spPrec),INTENT(IN) :: X
      REAL,INTENT(IN) :: Y
      
      point = CoordinateD(X,REAL(Y,spPrec))
   END FUNCTION CoordinateDR
   PURE SUBROUTINE assignCoord (point,X) 
      TYPE(Coordinate_type),INTENT(OUT) :: point
      REAL(KIND=spPrec),DIMENSION(2),INTENT(IN) :: X
      
      point%X = X(1)
      point%Y = X(2)
   END SUBROUTINE assignCoord
   PURE SUBROUTINE assignCoordR (point,X) 
      TYPE(Coordinate_type),INTENT(OUT) :: point
      REAL,DIMENSION(2),INTENT(IN) :: X
      
      point%X = X(1)
      point%Y = X(2)
   END SUBROUTINE assignCoordR
   
   PURE SUBROUTINE getCoordinate (point,X,Y)
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec),INTENT(OUT) :: X,Y
      
      X = point%X
      Y = point%Y
   END SUBROUTINE getCoordinate
   
   PURE FUNCTION Pitch(point) RESULT(w)
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec) :: w
      
      w = point%x
   END FUNCTION Pitch
   
   PURE FUNCTION Length(point) RESULT(h)
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec) :: h
      
      h = point%y
   END FUNCTION Length
   
   PURE FUNCTION addCoord (A,B) RESULT (AB)
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      TYPE(Coordinate_type) :: AB
      
      AB%X = A%X + B%X
      AB%Y = A%Y + B%Y
   END FUNCTION addCoord
  PURE  FUNCTION subCoord (A,B) RESULT (AB)
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      TYPE(Coordinate_type) :: AB
      
      AB%X = A%X - B%X
      AB%Y = A%Y - B%Y
   END FUNCTION subCoord
   PURE FUNCTION scaleCoordD (A,R) RESULT (AR)
      TYPE(Coordinate_type),INTENT(IN) :: A
      REAL(KIND=spPrec),INTENT(IN) :: R
      TYPE(Coordinate_type) :: AR
      
      AR%X = A%X * R
      AR%Y = A%Y * R
      
   END FUNCTION scaleCoordD
   PURE FUNCTION scaleCoordR (A,R) RESULT (AR)
      TYPE(Coordinate_type),INTENT(IN) :: A
      REAL,INTENT(IN) :: R
      TYPE(Coordinate_type) :: AR
      
      AR = scaleCoordD (A,REAL(R,kind=spPrec))
      
   END FUNCTION scaleCoordR
   PURE FUNCTION prodCoord (A,B) RESULT (AB)
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      REAL(KIND=spPrec) :: AB,dX,dY
      
      dX = A%X * B%X
      dY = A%Y * B%Y
      
      AB = dX + dY
   END FUNCTION prodCoord
   PURE FUNCTION sameCoord (A,B) RESULT (same)
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      LOGICAL :: same
      
      IF( A%X == B%X .and. A%Y == B%Y )THEN
         same = .TRUE.
      ELSE
         same = .FALSE.
      ENDIF
   END FUNCTION sameCoord
   
! *****************************************************************
   
   PURE FUNCTION LineAngleD(X,P) RESULT(myLine)
      TYPE(Line_type) :: myLine
      REAL(KIND=spPrec),INTENT(IN) :: X,P
      REAL(KIND=spPrec) :: PP,PI,PI2
      PARAMETER (PI=3.1415927_spPrec,PI2=PI/2.0_spPrec)
      
      PP = P
      IF( PP<-PI2 ) PP = -PI2
      IF( PP> PI2 ) PP =  PI2
      myLine%X0  = X
      myLine%Phi = PP        
   END FUNCTION LineAngleD  
   PURE FUNCTION LineAngleR(X,P) RESULT(myLine)
      TYPE(Line_type) :: myLine
      REAL,INTENT(IN) :: X,P

      myLine = LineAngleD(REAL(X,spPrec),REAL(P,spPrec))
   END FUNCTION LineAngleR
   
   PURE FUNCTION LinePktR(Point,P) RESULT(myLine)
      TYPE(Line_type) :: myLine
      TYPE(Coordinate_type),INTENT(IN) :: Point
      REAL,INTENT(IN) :: P
      
      myLine = LinePktD(Point,REAL(P,spPrec))
   END FUNCTION LinePktR
   PURE FUNCTION LinePktD(Point,P) RESULT(myLine)
      TYPE(Line_type) :: myLine
      TYPE(Coordinate_type),INTENT(IN) :: Point
      REAL(KIND=spPrec),INTENT(IN) :: P
      REAL(KIND=spPrec) :: X0,PP,PI,PI2
      PARAMETER (PI=3.1415927_spPrec,PI2=PI/2.0_spPrec)
         
      PP = P
      IF( PP<-PI2 ) PP = -PI2
      IF( PP> PI2 ) PP =  PI2
      X0 = Point%X - TAN(PP)*Point%Y
      myLine = LineAngleD(X0,PP)
   END FUNCTION LinePktD
   
   PURE FUNCTION Line2Pkt(A,B) RESULT(myLine)
      TYPE(Line_type) :: myLine
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      REAL(KIND=spPrec) :: Phi,dx,dy,PI,PI2
      PARAMETER (PI=3.1415927,PI2=PI/2.0)
         
      dx = A%X - B%X
      dy = A%Y - B%Y
      
      Phi = PI2-ATAN2(dy,dx) 
      IF( Phi > PI2 ) Phi = Phi - PI
      
      myLine = LinePktD(A,Phi)
   END FUNCTION Line2Pkt

   PURE FUNCTION LineXYTrackD(myTrack,Y) RESULT(myLine)
      TYPE(XYTrack_type),INTENT(IN) :: myTrack
      REAL(KIND=spPrec),INTENT(IN) :: Y
      TYPE(Line_type) :: myLine
      REAL(KIND=spPrec) :: Radius,PhiNew,X,RR,PhiTrk,X0Trk
      
      IF( ABS(myTrack%invRadius)<1.0E-16_spPrec )THEN 
         myLine = myTrack%Slope
      ELSE
         Radius  = 1.0_spPrec/myTrack%invRadius
         X0Trk   = myTrack%Slope%X0
         PhiTrk  = myTrack%Slope%Phi
! center of circle
!         CenterX = X0Trk + cos(PhiTrk) * Radius
!         CenterY = - sin(PhiTrk) * Radius
!         PhiNew = asin((Y-CenterY)/Radius)
!         X = CenterX - cos(PhiNew)*Radius
!  or     X = X0Trk + (cos(PhiTrk)-cos(PhiNew))*Radius
! numerically more stable:
         RR = Y/Radius + sin(PhiTrk)
         IF( abs(RR)<1.0 )THEN
            PhiNew = asin( Y/Radius + sin(PhiTrk) )
            X = X0Trk + Y * tan((PhiTrk+PhiNew)*0.5_spPrec)
            myLine = Line(Coordinate(X,Y),PhiNew)
         ELSE
! can't deal with curlers
            myLine = myTrack%Slope
         ENDIF
      ENDIF
   END FUNCTION LineXYTrackD
   PURE FUNCTION LineXYTrackR(myTrack,Y) RESULT(myLine)
      TYPE(XYTrack_type),INTENT(IN) :: myTrack
      REAL,INTENT(IN) :: Y
      TYPE(Line_type) :: myLine
      
      myLine = LineXYTrackD(myTrack,REAL(Y,spPrec))
   END FUNCTION LineXYTrackR

   PURE SUBROUTINE getLine(myLine,X,P)
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL(KIND=spPrec),INTENT(OUT) :: X,P
         
      X = myLine%X0
      P = myLine%Phi
   END SUBROUTINE getLine
    
   PURE FUNCTION getLineXD(myLine,Y) RESULT(X)
! line X at a given Y
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL(KIND=spPrec),INTENT(IN) :: Y
      REAL(KIND=spPrec) :: X
      
      X = myLine%X0 + TAN(myLine%Phi)*Y
   END FUNCTION getLineXD
   PURE FUNCTION getLineXR(myLine,Y) RESULT(X)
! line X at a given Y
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL,INTENT(IN) :: Y
      REAL(KIND=spPrec) :: X
      
      X = getLineXD(myLine,REAL(Y,spPrec))
   END FUNCTION getLineXR

   PURE FUNCTION getLineYD(myLine,X) RESULT(Y)
! line Y at a given X
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL(KIND=spPrec),INTENT(IN) :: X
      REAL(KIND=spPrec) :: Y
      
      Y = (X - myLine%X0)/TAN(myLine%Phi)
   END FUNCTION getLineYD
   PURE FUNCTION getLineYR(myLine,X) RESULT(Y)
! line Y at a given X
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL,INTENT(IN) :: X
      REAL(KIND=spPrec) :: Y
      
      Y = getLineYD(myLine,REAL(X,spPrec))
   END FUNCTION getLineYR
   
   FUNCTION sameLine (A,B) RESULT (same)
      TYPE(Line_type),INTENT(IN) :: A,B
      LOGICAL :: same
      
      IF( A%X0 == B%X0 .and. A%PHI == B%PHI )THEN
         same = .TRUE.
      ELSE
         same = .FALSE.
      ENDIF
   END FUNCTION sameLine

! *****************************************************************

   PURE FUNCTION XYTrackD(myLine,invRad) RESULT(myTrack)
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL(KIND=spPrec),INTENT(IN) :: invRad
      TYPE(XYTrack_type) :: myTrack
         
      myTrack%Slope = myLine
      myTrack%invRadius = invRad
   END FUNCTION XYTrackD
   PURE FUNCTION XYTrackR(myLine,invRad) RESULT(myTrack)
      TYPE(Line_type),INTENT(IN) :: myLine
      REAL,INTENT(IN) :: invRad
      TYPE(XYTrack_type) :: myTrack
         
      myTrack = XYTrackD(myLine,REAL(invRad,spPrec))
   END FUNCTION XYTrackR
   PURE SUBROUTINE getXYTrack(myTrack,X,P,iR)
      TYPE(XYTrack_type),INTENT(IN) :: myTrack
      REAL(KIND=spPrec),INTENT(OUT) :: X,P,iR
         
      CALL getLine(myTrack%Slope,X,P)
      iR = myTrack%invRadius
   END SUBROUTINE getXYTrack
   FUNCTION sameXYTrack (A,B) RESULT (same)
      TYPE(XYTrack_type),INTENT(IN) :: A,B
      LOGICAL :: same
      
      IF( A%slope == B%slope .and. A%invRadius == B%invRadius )THEN
         same = .TRUE.
      ELSE
         same = .FALSE.
      ENDIF
   END FUNCTION sameXYTrack

! *****************************************************************

   PURE SUBROUTINE FillRect(cX,cY,sX,sY,rect)
      REAL(KIND=spPrec),INTENT(IN) :: cX,cY,sX,sY
      TYPE(Rectangle_type),INTENT(OUT) :: rect
      
      rect%centre%X = cX
      rect%centre%Y = cY
      rect%size%X = sX
      rect%size%Y = sY
   END SUBROUTINE FillRect
   
   PURE FUNCTION getCentreRect(rect) RESULT(point)
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type) :: point
               
      point = rect%centre
   END FUNCTION getCentreRect
         
    
! *****************************************************************

   PURE SUBROUTINE FillGrid(cX,cY,sX,sY,mapNumber,grid)
      REAL(KIND=spPrec),INTENT(IN) :: cX,cY,sX,sY
      INTEGER ,INTENT(IN):: mapNumber
      TYPE(Grid_type),INTENT(OUT) :: grid
      
      grid%centre%X = cX
      grid%centre%Y = cY
      grid%size%X = sX
      grid%size%Y = sY
      grid%mapID    = mapNumber
   END SUBROUTINE FillGrid
      
   PURE FUNCTION getCentreGrid(grid) RESULT(point)
      TYPE(Grid_type),INTENT(IN) :: grid
      TYPE(Coordinate_type) :: point
      INTEGER :: ix,iy
      REAL(KIND=spPrec) :: nRect
      TYPE(bitmap_type) :: thismap
      INTEGER(KIND=mapPrec) :: mapline

! very short version         
!      X = grid%centre%X
!      Y = grid%centre%Y 
           
      thismap = gridMap(grid%mapID)
      
      Point = Coordinate(0.0,0.0)
      nRect = 0.0
      DO ix = 1,thismap%nx
         mapline = thismap%bmap(ix)
         DO iy = 1,thismap%ny
            IF( btest(mapline,iy) ) THEN
               Point = Point + getCentreRect(recty(grid,ix,iy))
               nRect = nRect + 1.0
            ENDIF
         ENDDO
      ENDDO

      nRect = 1.0 / nRect
      Point = Point * nRect
   END FUNCTION getCentreGrid

   
   SUBROUTINE printGrid(grid)
      TYPE(Grid_type),INTENT(IN) :: grid
      INTEGER :: ix,iy
      REAL(KIND=spPrec) :: X,Y
      TYPE(bitmap_type) :: thismap
      INTEGER(KIND=mapPrec) :: mapline
      CHARACTER(LEN=30) :: CHLINE
      CHLINE="                              "

      CALL getCoordinate(getCentreGrid(grid),X,Y)
      thismap = gridMap(grid%mapID)

      print*,thismap%nx," * ",thismap%ny," center:",x,y
      DO ix = 1,thismap%nx
         mapline = thismap%bmap(ix)
         DO iy = 1,thismap%ny
            IF( btest(mapline,iy) ) THEN
               CHLINE(iy:iy)="X"
            ELSE
               CHLINE(iy:iy)="."
            ENDIF
         ENDDO
         print*,CHLINE
      ENDDO
      print*,"----------------------------------------"
   END SUBROUTINE printGrid

   PURE FUNCTION bitmap(bits) RESULT(mymap)
      TYPE(bitmap_type) :: mymap
      CHARACTER(LEN=1),DIMENSION(:,:),INTENT(IN) :: bits
      INTEGER :: nx,ny,ix,iy
      INTEGER(KIND=mapPrec) :: mapline
      
      nx = size(bits,1)
      ny = size(bits,2)
      IF( bit_size(mapline) < ny)THEN
         mymap%nx = -1
         mymap%ny = -1
         RETURN
      ENDIF
      
      ALLOCATE( mymap%bmap(nx) )
      
      mymap%nx = nx
      mymap%ny = ny
      
      DO ix = 1,nx
         mapline = 0
         DO iy = 1,ny
            IF( bits(ix,iy)/=" " .AND.  bits(ix,iy)/=".") &
                mapline = ibset( mapline,iy )
         ENDDO
         mymap%bmap(ix) = mapline
      ENDDO
   END FUNCTION bitmap
      
   SUBROUTINE AllocBitMap(NumMaps)
      INTEGER,INTENT(IN) :: NumMaps

      ALLOCATE(gridMap(NumMaps))
   END SUBROUTINE AllocBitMap
   
   SUBROUTINE FillBitMap(mapNumber,bits)
      INTEGER,INTENT(IN) :: mapNumber
      CHARACTER(LEN=1),DIMENSION(:,:),INTENT(IN) :: bits
      
      gridMap(mapNumber) = bitmap(bits)
   END SUBROUTINE FillBitMap
   
   PURE FUNCTION recty(grid,ix,iy) RESULT(rect)
      TYPE(Grid_type),INTENT(IN) :: grid
      INTEGER,INTENT(IN) :: ix,iy
      TYPE(Rectangle_type) :: rect
      REAL(KIND=spPrec)  :: nx,ny,cX,cY,sX,sY
      
      nx = gridMap(grid%mapID)%nx
      ny = gridMap(grid%mapID)%ny
      sX = grid%size%X / nx
      sY = grid%size%Y / ny
      cX = grid%centre%X + (REAL(ix)-0.5-nx/2.0)*sX
      cY = grid%centre%Y + (REAL(iy)-0.5-ny/2.0)*sY
      CALL FillLoc(cX,cY,sX,sY,rect)
      
   END FUNCTION recty

! *** help routines for distance ***********************************  
  
   PURE SUBROUTINE Corners(rect,point)
!     returns the four corners of a rectangle
!     corner index:   3 4
!                     1 2
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type),DIMENSION(4),INTENT(OUT) :: point
      REAL(KIND=spPrec) :: dx , dy 
      
      dx = rect%size%X / 2.0
      dy = rect%size%Y / 2.0
      
      point(1)%x = rect%centre%X - dx
      point(1)%y = rect%centre%Y - dy
      
      point(2)%x = rect%centre%X + dx
      point(2)%y = rect%centre%Y - dy
      
      point(3)%x = rect%centre%X - dx
      point(3)%y = rect%centre%Y + dy
      
      point(4)%x = rect%centre%X + dx
      point(4)%y = rect%centre%Y + dy      
   END SUBROUTINE Corners
 
   FUNCTION DistXRect2Point(rect,point) RESULT(dist)
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type),INTENT(IN) :: point
      TYPE(Coordinate_type),DIMENSION(4) :: rectCorner
      REAL(KIND=spPrec) :: dist, distRight, distLeft
      
      CALL Corners(rect,rectCorner) 

      distRight =   point%X - rectCorner(2)%X
      distLeft  = - point%X + rectCorner(1)%X
      
      dist = max( distRight,distLeft )
   END FUNCTION DistXRect2Point
   FUNCTION DistYRect2Point(rect,point) RESULT(dist)
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type),INTENT(IN) :: point
      TYPE(Coordinate_type),DIMENSION(4) :: rectCorner
      REAL(KIND=spPrec) :: dist, distUp, distLow
      
      CALL Corners(rect,rectCorner) 

      distUp  =   point%Y - rectCorner(3)%Y
      distLow = - point%Y + rectCorner(1)%Y
      
      dist = max( distUp,distLow )
   END FUNCTION DistYRect2Point
   FUNCTION DistXRect2Rect(A,B) RESULT(mindist)
! distance  LC(RR) - RC(LR)  
!    (= left corner (right rectangle) - right corner (left rectangle))
!           LC(LR) - RC(RR) is always (more) negative
      TYPE(Rectangle_type),INTENT(IN) :: A,B
      TYPE(Coordinate_type),DIMENSION(4) :: PointA,PointB
      REAL(KIND=spPrec) :: mindist,dist1,dist2
      
      CALL Corners(A,PointA) ! get the 4 corners of A
      CALL Corners(B,PointB) ! get the 4 corners of B
      
      dist1 = PointA(1)%X - PointB(2)%X
      dist2 = PointB(1)%X - PointA(2)%X
      mindist = max(dist1,dist2)      
   END FUNCTION DistXRect2Rect
   FUNCTION DistYRect2Rect(A,B) RESULT(mindist)
! distance  LC(UR) - UC(LR)   
!    (= lower corner (upper rectangle) - upper corner (lower rectangle))
!           LC(LR) - UC(UR) is always (more) negative
      TYPE(Rectangle_type),INTENT(IN) :: A,B
      TYPE(Coordinate_type),DIMENSION(4) :: PointA,PointB
      REAL(KIND=spPrec) :: mindist,dist1,dist2
      
      CALL Corners(A,PointA) ! get the 4 corners of A
      CALL Corners(B,PointB) ! get the 4 corners of B
      
      dist1 = PointA(1)%Y - PointB(3)%Y
      dist2 = PointB(1)%Y - PointA(3)%Y
      mindist = max(dist1,dist2)      
   END FUNCTION DistYRect2Rect
   
   FUNCTION SignDistLine2Point(myline,point) RESULT(dist)
      TYPE(Line_type),INTENT(IN) :: myline
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec) :: dist,dx
! signed distance: line passes to the right >0 ; left <0
         
      dx = myline%X0 + point%Y*TAN(myline%Phi) - point%X
            
      dist = dx*COS(myline%Phi)         
   END FUNCTION SignDistLine2Point
   
! *** routines for Distance: ***************************************
   
   FUNCTION DistPoint2Point(A,B) RESULT(dist)
      TYPE(Coordinate_type),INTENT(IN) :: A,B
      TYPE(Coordinate_type) :: AsubB
      REAL(KIND=spPrec) :: dist
      
      AsubB = A - B
      dist = SQRT( AsubB*AsubB )
   END FUNCTION DistPoint2Point
   
   FUNCTION DistLine2Point(myline,point) RESULT(dist)
      TYPE(Line_type),INTENT(IN) :: myline
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec) :: dist
   
      dist = ABS( SignDistLine2Point(myline,point) )
      
   END FUNCTION DistLine2Point
   
   FUNCTION DistRect2Point(rect,point) RESULT(dist)
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec)  :: dist,dX,dY
   
      dX = DistXRect2Point(rect,point)
      dY = DistYRect2Point(rect,point)
      IF( dX<=0 .AND. dY<=0 )THEN
! point is inside rectangle, distance to border (negative)
         dist = max(dX,dY)
      ELSE      
         dX = max(dX,0.0_spPrec)  
         dY = max(dY,0.0_spPrec)
         dist = SQRT( dX**2 + dY**2 )
      ENDIF
   END FUNCTION DistRect2Point      
   
   FUNCTION DistRect2Line(rect,myline) RESULT(dist)
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Line_type),INTENT(IN) :: myline
      TYPE(Coordinate_type),DIMENSION(4) :: reCorner
      REAL(KIND=spPrec)  :: dist,mindist,cd,distA,distB,dx,dy,fac,PI,PI4
      REAL(KIND=spPrec),DIMENSION(4) :: cdist
      PARAMETER (PI=3.1415927,PI4=PI/4.0)
      INTEGER :: npos,nneg,I
      
      CALL Corners(rect,reCorner)
      mindist = huge(mindist)
      npos = 0
      nneg = 0
      DO I=1,4
         cd = SignDistLine2Point(myline,reCorner(I))
         cdist(I) = cd
         IF( cd > 1.0E-5 )THEN
            npos = npos+1
         ELSEIF( cd < -1.0E-5 )THEN
            nneg = nneg+1
         ENDIF
         mindist = MIN( mindist,ABS(cd) )
      ENDDO
            
      IF( npos == 4 .OR. nneg==4 )THEN
! all corners pointing to same direction
         dist = mindist
      ELSE
! line passes through rectangle ; it's supposed to be 
! the largest circle with centre on the line
! that fits into the rectangle   (negative)
         fac = COS( ABS(myline%Phi)-PI4 )*SQRT(2.0)
         IF( myline%Phi > 0.0 )THEN
            distA = -cdist(2)/fac
            distB =  cdist(3)/fac
         ELSE
            distA =  cdist(1)/fac
            distB = -cdist(4)/fac
         ENDIF

         IF( distA < 0.0 .OR. distB < 0.0 )THEN
            IF( distA < -1.0E-5 ) print*,"WARNING: DistRect2Line> distA should be >0",distA
            IF( distB < -1.0E-5 ) print*,"WARNING: DistRect2Line> distB should be >0",distB  
            dist = 0.0
         ELSE
            dx = rect%size%X / 2.0
            dy = rect%size%Y / 2.0
            dist = -min(distA,distB,dx,dy)
         ENDIF
      ENDIF
   END FUNCTION DistRect2Line      
   
   FUNCTION DistRect2Rect(A,B) RESULT(dist)
      TYPE(Rectangle_type),INTENT(IN) :: A,B
      REAL(KIND=spPrec) :: dist,dX,dY
   
      dX = DistXRect2Rect(A,B)
      dY = DistYRect2Rect(A,B)
      IF( dX<=0.0 .AND. dY<=0.0 )THEN
! rectangles overlap by   (negative)
         dist = max(dX,dY)
      ELSE      
         dX = max(dX,0.0_spPrec)  
         dY = max(dY,0.0_spPrec)
         dist = SQRT( dX**2 + dY**2 )
      ENDIF
   END FUNCTION DistRect2Rect

   FUNCTION DistGrid2Point(grid,point) RESULT(dist)
      TYPE(Grid_type),INTENT(IN) :: grid
      TYPE(Coordinate_type),INTENT(IN) :: point
      REAL(KIND=spPrec)  :: dist
      INTEGER :: ix,iy
      TYPE(bitmap_type) :: thismap
      INTEGER(KIND=mapPrec) :: mapline
      
      thismap = gridMap(grid%mapID)
      
      dist = huge(dist)
      DO ix = 1,thismap%nx
         mapline = thismap%bmap(ix)
         DO iy = 1,thismap%ny
            IF( btest(mapline,iy) ) &
            dist = min(dist,Distance( recty(grid,ix,iy),point ))
         ENDDO
      ENDDO
      
   END FUNCTION DistGrid2Point

   FUNCTION DistGrid2Line(grid,myline) RESULT(dist)
      TYPE(Grid_type),INTENT(IN) :: grid
      TYPE(Line_type),INTENT(IN) :: myline
      REAL(KIND=spPrec)  :: dist
      INTEGER :: ix,iy
      TYPE(bitmap_type) :: thismap
      INTEGER(KIND=mapPrec) :: mapline
      
      thismap = gridMap(grid%mapID)
      
      dist = huge(dist)
      DO ix = 1,thismap%nx
         mapline = thismap%bmap(ix)
         DO iy = 1,thismap%ny
            IF( btest(mapline,iy) ) &
            dist = min(dist,Distance( recty(grid,ix,iy),myline ))
         ENDDO
      ENDDO
      
   END FUNCTION DistGrid2Line

   FUNCTION DistGrid2Grid(A,B) RESULT(dist)
      TYPE(Grid_type),INTENT(IN) :: A,B
      REAL(KIND=spPrec)  :: dist
      INTEGER :: ixA,iyA,ixB,iyB
      TYPE(bitmap_type) :: Amap,Bmap
      TYPE(Rectangle_type) :: Arect,Brect
      INTEGER(KIND=mapPrec) :: Aline,Bline
      
      Amap = gridMap(A%mapID)
      Bmap = gridMap(B%mapID)
      
      dist = huge(dist)
      DO ixA = 1,Amap%nx
         Aline = Amap%bmap(ixA)
         DO iyA = 1,Amap%ny
            IF( btest(Aline,iyA) ) THEN
               Arect = recty(A,ixA,iyA)
            
               DO ixB = 1,Bmap%nx
                  Bline = Bmap%bmap(ixB)
                  DO iyB = 1,Bmap%ny
                     IF( btest(Bline,iyB) ) THEN
                        Brect = recty(B,ixB,iyB)
            
                        dist = min(dist,distance(Arect,Brect))
                     ENDIF
                  ENDDO
               ENDDO
            
            ENDIF
         ENDDO
      ENDDO
      
   END FUNCTION DistGrid2Grid

! *******************************************************************
!  integral of a track over a pad
   FUNCTION intTrackRect(track,width,rect) RESULT(integral)     
      TYPE(Line_type),INTENT(IN) :: track
      REAL(KIND=spPrec),INTENT(IN) :: width
      TYPE(Rectangle_type),INTENT(IN) :: rect
      TYPE(Coordinate_type),DIMENSION(4) :: point
      REAL(KIND=spPrec),DIMENSION(4) :: dPkt
      REAL(KIND=spPrec) :: tpcchr,integral,X,P,sinP,cosP,X1,X2,dfreq,mindist
      INTEGER :: I
      
      CALL Corners(rect,point)
      DO I=1,4
         dPkt(I) = SignDistLine2Point(track,point(I))
      ENDDO
      CALL getLine(track,X,P)
      
      sinP = sin(P)
      cosP = cos(P)
      
      mindist = minval(abs(dPkt))
      IF( mindist > 5.0 * width )THEN 
! track more than 5 sigma away from pad
         tpcchr = 0.0
      ELSEIF( mindist > 3.0*width .AND. mindist> 5.0*rect%size%x)THEN 
! track more than 3 sigma away from pad and more than 5 pads between
         tpcchr = 0.0
      ELSE
         IF( abs(sinP)<1.0E-5_spPrec )THEN
! track parallel to pad
            X1 = (dPkt(1)+dPkt(3)) / (2.0*width)
            X2 = (dPkt(2)+dPkt(4)) / (2.0*width)
            tpcchr = ( dfreq(X1)-dfreq(X2) )*rect%size%Y
         ELSE
! Deans formula
               tpcchr = eta(dPkt(3),width)-eta(dPkt(1),width) &
                       +eta(dPkt(2),width)-eta(dPkt(4),width)
               tpcchr = tpcchr / sinP/cosP
         ENDIF
      ENDIF
     
      IF( tpcchr<0.0 )THEN 
         print*,"WARNING: intTrackRect> numerical problems",P,tpcchr
         tpcchr = 0.0
      ENDIF
      
      integral = tpcchr
   END FUNCTION intTrackRect
   FUNCTION ETA(u,s) RESULT(x)
      REAL(KIND=spPrec),INTENT(IN) :: u,s
      REAL(KIND=spPrec) :: expo,x,derf

      expo = max(-100.0_spPrec,-u*u/s/s/2.0_spPrec)
      x = u/2.0_spPrec*derf(u/sqrt2/s) + s/sqrt2pi*exp(expo)
   END FUNCTION ETA
   
   FUNCTION intTrackGrid(track,width,grid) RESULT(integral)     
      TYPE(Line_type),INTENT(IN) :: track
      REAL(KIND=spPrec),INTENT(IN) :: width
      TYPE(Grid_type),INTENT(IN) :: grid
      TYPE(Rectangle_type) :: rect
      TYPE(bitmap_type) :: thismap
      INTEGER(KIND=mapPrec) :: mapline
      REAL(KIND=spPrec) :: dist,integral,X,P,nRect
      INTEGER :: ix,iy
      
      CALL getLine(track,X,P)      
      thismap = gridMap(grid%mapID)
      
      integral = 0.0
      nRect = 0.0
      DO ix = 1,thismap%nx
         mapline = thismap%bmap(ix)
         DO iy = 1,thismap%ny
            IF( btest(mapline,iy) ) THEN
               rect = recty(grid,ix,iy)
               dist = Distance(track,rect%centre)
               integral = integral + &
               exp(- (dist/width)**2/2.0_spPrec ) / width / sqrt2pi
               nRect = nRect + 1.0
            ENDIF
         ENDDO
      ENDDO
      integral = integral/nRect

   END FUNCTION intTrackGrid
END MODULE Space