/[PAMELA software]/gpamela/gpfield/gufld.F
ViewVC logotype

Diff of /gpamela/gpfield/gufld.F

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3.1.1.1 by cafagna, Thu Jul 11 16:02:01 2002 UTC revision 3.2 by pam-ba, Mon Dec 5 12:15:20 2005 UTC
# Line 1  Line 1 
1  *  *
2  * $Id$  * $Id: gufld.F,v 3.1.1.1 2002/07/11 16:02:01 cafagna Exp $
3    *
4    * $Log: gufld.F,v $
5    * Revision 3.1.1.1  2002/07/11 16:02:01  cafagna
6    * First GPAMELA release on CVS
7  *  *
 * $Log$  
8  *  *
9  *CMZ :  2.01/00 06/03/2000  13.07.03  by  Francesco Cafagna  *CMZ :  2.01/00 06/03/2000  13.07.03  by  Francesco Cafagna
10  *CMZ :  2.00/00 03/03/2000  15.39.05  by  Francesco Cafagna  *CMZ :  2.00/00 03/03/2000  15.39.05  by  Francesco Cafagna
# Line 24  Line 27 
27  *                                                                      *  *                                                                      *
28  ************************************************************************  ************************************************************************
29  #include "gpfield.inc"  #include "gpfield.inc"
30          REAL*8 VVINT(3),FFINT(3)
31          REAL V(3),F(3)
32          REAL*8 CM_TO_M , TESLA_TO_KGAUSS
33          PARAMETER(CM_TO_M=1.D-2 , TESLA_TO_KGAUSS = 1.D1)
34        
35    
36    C*
37    C      INTEGER II,III
38    C      REAL DISM,F0X,F0Y,F0Z,F1X,F1Y,F1Z,F2X,F2Y,F2Z,
39    C     +     F3X,F3Y,F3Z
40    C      REAL V(3),F(3),AV(3)
41    C*
42    C* Transform coordinates to Spectrometer frame
43    C*
44    C      CALL GPMASPE(V)
45    C*
46    C* Take just the absolute value for the coordinates
47    C*
48    C      DO I=1,3
49    C         AV(I) = ABS( V(I) )
50    C      ENDDO
51    C      F(1)=0.
52    C      F(2)=0.
53    C      F(3)=0.
54    C*
55    C* Check if we are outside the map
56    C*
57    C      IF( (AV(1).GE.20).OR.(AV(2).GE.20).OR.(AV(3).GE.60.) )
58    C     +    GOTO 10
59    C      IV(1)=INT(AV(1)*2.)+1
60    C      IV(2)=INT(AV(2)*2.)+1
61    C      IV(3)=INT(AV(3)/2.)+1
62    C      DO I1=0,1
63    C         DO I2=0,1
64    C            DO I3=0,1
65    C               II=I1*4+I2*2+I3+1
66    C               VV(II,1)=FLOAT(IV(1)+I1-1)*0.5
67    C               VV(II,2)=FLOAT(IV(2)+I2-1)*0.5
68    C               VV(II,3)=FLOAT(IV(3)+I3-1)*2.
69    C               IVV(II,1)=IV(1)+I1
70    C               IVV(II,2)=IV(2)+I2
71    C               IVV(II,3)=IV(3)+I3
72    C               DD(II)=(VV(II,1)-AV(1))**2 + (VV(II,2)-AV(2))**2 +
73    C     +         (VV(II,3)-AV(3))**2
74    C            ENDDO
75    C         ENDDO
76    C      ENDDO
77    C* --- v0
78    C      DISM=1.E9
79    C      II=0
80    C      DO I=1,8
81    C         IF(DD(I).LT.DISM) THEN
82    C            DISM=DD(I)
83    C            II=I
84    C         END IF
85    C      END DO
86    C      DO I=1,3
87    C         V0(I)=VV(II,I)
88    C      END DO
89    C      F0X=FX(IVV(II,1),IVV(II,2),IVV(II,3))
90    C      F0Y=FY(IVV(II,1),IVV(II,2),IVV(II,3))
91    C      F0Z=FZ(IVV(II,1),IVV(II,2),IVV(II,3))
92    C* --- v1
93    C      V1(2)=V0(2)
94    C      V1(3)=V0(3)
95    C      IF(AV(1).GE.V0(1)) THEN
96    C         III=IVV(II,1)+1
97    C         V1(1)=V0(1)+0.5
98    C      ELSE
99    C         III=IVV(II,1)-1
100    C         V1(1)=V0(1)-0.5
101    C      END IF
102    C      F1X=FX(III,IVV(II,2),IVV(II,3))
103    C      F1Y=FY(III,IVV(II,2),IVV(II,3))
104    C      F1Z=FZ(III,IVV(II,2),IVV(II,3))
105    C* --- v2
106    C      V2(1)=V0(1)
107    C      V2(3)=V0(3)
108    C      IF(AV(2).GE.V0(2)) THEN
109    C         III=IVV(II,2)+1
110    C         V2(2)=V0(2)+0.5
111    C      ELSE
112    C         III=IVV(II,2)-1
113    C         V2(2)=V0(2)-0.5
114    C      END IF
115    C      F2X=FX(IVV(II,1),III,IVV(II,3))
116    C      F2Y=FY(IVV(II,1),III,IVV(II,3))
117    C      F2Z=FZ(IVV(II,1),III,IVV(II,3))
118    C* --- v3
119    C      V3(1)=V0(1)
120    C      V3(2)=V0(2)
121    C      IF(AV(3).GE.V0(3)) THEN
122    C         III=IVV(II,3)+1
123    C         V3(3)=V0(3)+2.
124    C      ELSE
125    C         III=IVV(II,3)-1
126    C         V3(3)=V0(3)-2.
127    C      END IF
128    C      F3X=FX(IVV(II,1),IVV(II,2),III)
129    C      F3Y=FY(IVV(II,1),IVV(II,2),III)
130    C      F3Z=FZ(IVV(II,1),IVV(II,2),III)
131    C* --- linear interpolation, magnetic field calculation
132    C      CALL FLIN3(V0,V1,V2,V3,F0X,F1X,F2X,F3X,AV,F(1))
133    C      CALL FLIN3(V0,V1,V2,V3,F0Y,F1Y,F2Y,F3Y,AV,F(2))
134    C      CALL FLIN3(V0,V1,V2,V3,F0Z,F1Z,F2Z,F3Z,AV,F(3))
135    C* --- mirroing
136    C      IF(V(2).LT.0.) THEN
137    C         F(1)=-1.*F(1)
138    C         F(3)=-1.*F(3)
139    C      END IF
140    C      IF(V(1).LT.0.) F(1)=-1.*F(1)
141    C      IF(V(3).LT.0.) F(3)=-1.*F(3)
142    C*
143    C* Transform coordinates back to MARS
144    C*
145    C   10 CALL GPSPEMA(V)
146    C      RETURN
147    C      END
148    
149  *  *
150        INTEGER II,III  
151        REAL DISM,F0X,F0Y,F0Z,F1X,F1Y,F1Z,F2X,F2Y,F2Z,    
      +     F3X,F3Y,F3Z  
       REAL V(3),F(3),AV(3)  
152  *  *
153  * Transform coordinates to Spectrometer frame  * Transform coordinates to Spectrometer frame
154  *  *
155        CALL GPMASPE(V)        CALL GPMASPE(V)
156  *  *
157  * Take just the absolute value for the coordinates  * INTERFACE TO TRACKER FIELD ROUTINES
158  *  *      
159        DO I=1,3        DO I=1,3
160           AV(I) = ABS( V(I) )           VVINT(I) = DBLE(V(I)) * CM_TO_M
161        ENDDO        ENDDO
162        F(1)=0.        
163        F(2)=0.        CALL inter_B(VVINT(1),VVINT(2),VVINT(3),FFINT) !coordinates in m, Field in Tesla
164        F(3)=0.  
 *  
 * Check if we are outside the map  
 *  
       IF( (AV(1).GE.20).OR.(AV(2).GE.20).OR.(AV(3).GE.60.) )  
      +    GOTO 10  
       IV(1)=INT(AV(1)*2.)+1  
       IV(2)=INT(AV(2)*2.)+1  
       IV(3)=INT(AV(3)/2.)+1  
       DO I1=0,1  
          DO I2=0,1  
             DO I3=0,1  
                II=I1*4+I2*2+I3+1  
                VV(II,1)=FLOAT(IV(1)+I1-1)*0.5  
                VV(II,2)=FLOAT(IV(2)+I2-1)*0.5  
                VV(II,3)=FLOAT(IV(3)+I3-1)*2.  
                IVV(II,1)=IV(1)+I1  
                IVV(II,2)=IV(2)+I2  
                IVV(II,3)=IV(3)+I3  
                DD(II)=(VV(II,1)-AV(1))**2 + (VV(II,2)-AV(2))**2 +  
      +         (VV(II,3)-AV(3))**2  
             ENDDO  
          ENDDO  
       ENDDO  
 * --- v0  
       DISM=1.E9  
       II=0  
       DO I=1,8  
          IF(DD(I).LT.DISM) THEN  
             DISM=DD(I)  
             II=I  
          END IF  
       END DO  
165        DO I=1,3        DO I=1,3
166           V0(I)=VV(II,I)           F(I) = REAL( FFINT(I) * TESLA_TO_KGAUSS )
167        END DO        ENDDO
168        F0X=FX(IVV(II,1),IVV(II,2),IVV(II,3))        
       F0Y=FY(IVV(II,1),IVV(II,2),IVV(II,3))  
       F0Z=FZ(IVV(II,1),IVV(II,2),IVV(II,3))  
 * --- v1  
       V1(2)=V0(2)  
       V1(3)=V0(3)  
       IF(AV(1).GE.V0(1)) THEN  
          III=IVV(II,1)+1  
          V1(1)=V0(1)+0.5  
       ELSE  
          III=IVV(II,1)-1  
          V1(1)=V0(1)-0.5  
       END IF  
       F1X=FX(III,IVV(II,2),IVV(II,3))  
       F1Y=FY(III,IVV(II,2),IVV(II,3))  
       F1Z=FZ(III,IVV(II,2),IVV(II,3))  
 * --- v2  
       V2(1)=V0(1)  
       V2(3)=V0(3)  
       IF(AV(2).GE.V0(2)) THEN  
          III=IVV(II,2)+1  
          V2(2)=V0(2)+0.5  
       ELSE  
          III=IVV(II,2)-1  
          V2(2)=V0(2)-0.5  
       END IF  
       F2X=FX(IVV(II,1),III,IVV(II,3))  
       F2Y=FY(IVV(II,1),III,IVV(II,3))  
       F2Z=FZ(IVV(II,1),III,IVV(II,3))  
 * --- v3  
       V3(1)=V0(1)  
       V3(2)=V0(2)  
       IF(AV(3).GE.V0(3)) THEN  
          III=IVV(II,3)+1  
          V3(3)=V0(3)+2.  
       ELSE  
          III=IVV(II,3)-1  
          V3(3)=V0(3)-2.  
       END IF  
       F3X=FX(IVV(II,1),IVV(II,2),III)  
       F3Y=FY(IVV(II,1),IVV(II,2),III)  
       F3Z=FZ(IVV(II,1),IVV(II,2),III)  
 * --- linear interpolation, magnetic field calculation  
       CALL FLIN3(V0,V1,V2,V3,F0X,F1X,F2X,F3X,AV,F(1))  
       CALL FLIN3(V0,V1,V2,V3,F0Y,F1Y,F2Y,F3Y,AV,F(2))  
       CALL FLIN3(V0,V1,V2,V3,F0Z,F1Z,F2Z,F3Z,AV,F(3))  
 * --- mirroing  
       IF(V(2).LT.0.) THEN  
          F(1)=-1.*F(1)  
          F(3)=-1.*F(3)  
       END IF  
       IF(V(1).LT.0.) F(1)=-1.*F(1)  
       IF(V(3).LT.0.) F(3)=-1.*F(3)  
169  *  *
170  * Transform coordinates back to MARS  * Transform coordinates back to MARS
171  *  *
Legend:
Removed from v.3.1.1.1  
changed lines
  Added in v.3.2
  ViewVC Help
Powered by ViewVC 1.1.23