*
* $Id: gufld.F,v 3.1.1.1 2002/07/11 16:02:01 cafagna Exp $
*
* $Log: gufld.F,v $
* Revision 3.1.1.1  2002/07/11 16:02:01  cafagna
* First GPAMELA release on CVS
*
*
*CMZ :  2.01/00 06/03/2000  13.07.03  by  Francesco Cafagna
*CMZ :  2.00/00 03/03/2000  15.39.05  by  Francesco Cafagna
*CMZU:  1.01/00 26/04/96  15.12.30  by  Paolo (The Magic) Papini
*-- Author :
      SUBROUTINE GUFLD(V,F)
************************************************************************
*                                                                      *
* To map Pamela magnetic field                                         *
* A bit of housekeeping: delete print statement etc. etc., by F. Caf.  *
*                                                                      *
* Variables definition:                                                *
* IN:                                                                  *
*  V , vector with coordinates in MARS                                 *
* OUT:                                                                 *
*  F , Magnetic field components along X, Y and Z                      *
*                                                                      *
* Called by: GHELIX, GRKUTA                                            *
* Author: Paolo Papini   16/02/96                                      *
*                                                                      *
************************************************************************
#include "gpfield.inc"
      REAL*8 VVINT(3),FFINT(3)
      REAL V(3),F(3)
      REAL*8 CM_TO_M , TESLA_TO_KGAUSS
      PARAMETER(CM_TO_M=1.D-2 , TESLA_TO_KGAUSS = 1.D1)
    

C*
C      INTEGER II,III
C      REAL DISM,F0X,F0Y,F0Z,F1X,F1Y,F1Z,F2X,F2Y,F2Z,
C     +     F3X,F3Y,F3Z
C      REAL V(3),F(3),AV(3)
C*
C* Transform coordinates to Spectrometer frame
C*
C      CALL GPMASPE(V)
C*
C* Take just the absolute value for the coordinates
C*
C      DO I=1,3
C         AV(I) = ABS( V(I) )
C      ENDDO
C      F(1)=0.
C      F(2)=0.
C      F(3)=0.
C*
C* Check if we are outside the map
C*
C      IF( (AV(1).GE.20).OR.(AV(2).GE.20).OR.(AV(3).GE.60.) )
C     +    GOTO 10
C      IV(1)=INT(AV(1)*2.)+1
C      IV(2)=INT(AV(2)*2.)+1
C      IV(3)=INT(AV(3)/2.)+1
C      DO I1=0,1
C         DO I2=0,1
C            DO I3=0,1
C               II=I1*4+I2*2+I3+1
C               VV(II,1)=FLOAT(IV(1)+I1-1)*0.5
C               VV(II,2)=FLOAT(IV(2)+I2-1)*0.5
C               VV(II,3)=FLOAT(IV(3)+I3-1)*2.
C               IVV(II,1)=IV(1)+I1
C               IVV(II,2)=IV(2)+I2
C               IVV(II,3)=IV(3)+I3
C               DD(II)=(VV(II,1)-AV(1))**2 + (VV(II,2)-AV(2))**2 +
C     +         (VV(II,3)-AV(3))**2
C            ENDDO
C         ENDDO
C      ENDDO
C* --- v0
C      DISM=1.E9
C      II=0
C      DO I=1,8
C         IF(DD(I).LT.DISM) THEN
C            DISM=DD(I)
C            II=I
C         END IF
C      END DO
C      DO I=1,3
C         V0(I)=VV(II,I)
C      END DO
C      F0X=FX(IVV(II,1),IVV(II,2),IVV(II,3))
C      F0Y=FY(IVV(II,1),IVV(II,2),IVV(II,3))
C      F0Z=FZ(IVV(II,1),IVV(II,2),IVV(II,3))
C* --- v1
C      V1(2)=V0(2)
C      V1(3)=V0(3)
C      IF(AV(1).GE.V0(1)) THEN
C         III=IVV(II,1)+1
C         V1(1)=V0(1)+0.5
C      ELSE
C         III=IVV(II,1)-1
C         V1(1)=V0(1)-0.5
C      END IF
C      F1X=FX(III,IVV(II,2),IVV(II,3))
C      F1Y=FY(III,IVV(II,2),IVV(II,3))
C      F1Z=FZ(III,IVV(II,2),IVV(II,3))
C* --- v2
C      V2(1)=V0(1)
C      V2(3)=V0(3)
C      IF(AV(2).GE.V0(2)) THEN
C         III=IVV(II,2)+1
C         V2(2)=V0(2)+0.5
C      ELSE
C         III=IVV(II,2)-1
C         V2(2)=V0(2)-0.5
C      END IF
C      F2X=FX(IVV(II,1),III,IVV(II,3))
C      F2Y=FY(IVV(II,1),III,IVV(II,3))
C      F2Z=FZ(IVV(II,1),III,IVV(II,3))
C* --- v3
C      V3(1)=V0(1)
C      V3(2)=V0(2)
C      IF(AV(3).GE.V0(3)) THEN
C         III=IVV(II,3)+1
C         V3(3)=V0(3)+2.
C      ELSE
C         III=IVV(II,3)-1
C         V3(3)=V0(3)-2.
C      END IF
C      F3X=FX(IVV(II,1),IVV(II,2),III)
C      F3Y=FY(IVV(II,1),IVV(II,2),III)
C      F3Z=FZ(IVV(II,1),IVV(II,2),III)
C* --- linear interpolation, magnetic field calculation
C      CALL FLIN3(V0,V1,V2,V3,F0X,F1X,F2X,F3X,AV,F(1))
C      CALL FLIN3(V0,V1,V2,V3,F0Y,F1Y,F2Y,F3Y,AV,F(2))
C      CALL FLIN3(V0,V1,V2,V3,F0Z,F1Z,F2Z,F3Z,AV,F(3))
C* --- mirroing
C      IF(V(2).LT.0.) THEN
C         F(1)=-1.*F(1)
C         F(3)=-1.*F(3)
C      END IF
C      IF(V(1).LT.0.) F(1)=-1.*F(1)
C      IF(V(3).LT.0.) F(3)=-1.*F(3)
C*
C* Transform coordinates back to MARS
C*
C   10 CALL GPSPEMA(V)
C      RETURN
C      END

*

   
*
* Transform coordinates to Spectrometer frame
*
      CALL GPMASPE(V)
*
* INTERFACE TO TRACKER FIELD ROUTINES
*      
      DO I=1,3
         VVINT(I) = DBLE(V(I)) * CM_TO_M
      ENDDO
      
      CALL inter_B(VVINT(1),VVINT(2),VVINT(3),FFINT) !coordinates in m, Field in Tesla

      DO I=1,3
         F(I) = REAL( FFINT(I) * TESLA_TO_KGAUSS )
      ENDDO
      
*
* Transform coordinates back to MARS
*
   10 CALL GPSPEMA(V)
      RETURN
      END
