ground/src/inter_B_inner.for

*************************************************************************
*     
*     Subroutine inter_B_inner.f
*     
*     it computes the magnetic field in a chosen point x,y,z inside the
*     magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     the value is computed for two different inner maps and then averaged
*     
*     needs:
*     - ../common/common_B_inner.f
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B_inner(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include './common_B_inner.for'


c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      real*8 x,y,z              !point of interpolation
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)
      real*8 res1(3),res2(3)    !interpolated B components for the two maps

      integer ic                !index for B components:
                                ! ic=1 ---> Bx
                                ! ic=2 ---> By
                                ! ic=3 ---> Bz

      integer cube(3)           !vector of indexes identifying the cube
                                ! containing the point of interpolation
                                ! (see later...)
      
      real*8 xl,xh,yl,yh,zl,zh  !cube vertexes coordinates

      real*8 xr,yr,zr           !reduced variables (coordinates of the 
                                ! point of interpolation inside the cube)

      real*8 Bp(8)              !vector of values of B component 
                                ! being computed, on the eight cube vertexes


c------------------------------------------------------------------------
c     
c     *** FIRST MAP ***
c     
c------------------------------------------------------------------------

      do ic=1,3                 !loops on the three B components

c------------------------------------------------------------------------
c     
c     chooses the coordinates interval containing the input point
c     
c------------------------------------------------------------------------
c     e.g.:
c     
c     x1    x2    x3    x4    x5...
c     |-----|-+---|-----|-----|----
c     ~~~~~~~~x
c     
c     in this case the right interval is identified by indexes 2-3, so the 
c     value assigned to cube variable is 2

        cube(1)=INT((nx-1)*(x-px1min(ic))/(px1max(ic)-px1min(ic))) +1
        cube(2)=INT((ny-1)*(y-py1min(ic))/(py1max(ic)-py1min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz1min(ic))/(pz1max(ic)-pz1min(ic))) +1
        
c------------------------------------------------------------------------
c     
c     if the point falls beyond the extremes of the grid...
c     
c------------------------------------------------------------------------
c     
c     ~~~~~~~~~~x1    x2    x3...
c     - - + - - |-----|-----|----
c     ~~~~x
c     
c     in the case cube = 1
c     
c     
c     ...nx-2  nx-1  nx
c     ----|-----|-----| - - - + - -
c     ~~~~~~~~~~~~~~~~~~~~~~~~x
c     
c     in this case cube = nx-1

        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1


c------------------------------------------------------------------------
c     
c     temporary variables definition for field computation
c     
c------------------------------------------------------------------------

        xl = px1(cube(1),ic)    !X coordinates of cube vertexes
        xh = px1(cube(1)+1,ic)

        yl = py1(cube(2),ic)    !Y coordinates of cube vertexes
        yh = py1(cube(2)+1,ic)

        zl = pz1(cube(3),ic)    !Z coordinates of cube vertexes
        zh = pz1(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)   !reduced variables
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b1(cube(1)  ,cube(2)  ,cube(3)  ,ic) !ic-th component of B
        Bp(2) = b1(cube(1)+1,cube(2)  ,cube(3)  ,ic) ! on the eight cube
        Bp(3) = b1(cube(1)  ,cube(2)+1,cube(3)  ,ic) ! vertexes
        Bp(4) = b1(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b1(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b1(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b1(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b1(cube(1)+1,cube(2)+1,cube(3)+1,ic)

c------------------------------------------------------------------------
c     
c     computes interpolated ic-th component of B in (x,y,z)
c     
c------------------------------------------------------------------------

        res1(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr


      enddo
      
c------------------------------------------------------------------------
c     
c     *** SECOND MAP ***
c     
c------------------------------------------------------------------------

c     second map is rotated by 180 degree along the Z axis. so change sign
c     of x and y coordinates and then change sign to Bx and By components
c     to obtain the correct result

      x=-x
      y=-y

      do ic=1,3

        cube(1)=INT((nx-1)*(x-px2min(ic))/(px2max(ic)-px2min(ic))) +1
        cube(2)=INT((ny-1)*(y-py2min(ic))/(py2max(ic)-py2min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz2min(ic))/(pz2max(ic)-pz2min(ic))) +1
        
        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1

        xl = px2(cube(1),ic)
        xh = px2(cube(1)+1,ic)

        yl = py2(cube(2),ic)
        yh = py2(cube(2)+1,ic)

        zl = pz2(cube(3),ic)
        zh = pz2(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b2(cube(1)  ,cube(2)  ,cube(3)  ,ic)
        Bp(2) = b2(cube(1)+1,cube(2)  ,cube(3)  ,ic)
        Bp(3) = b2(cube(1)  ,cube(2)+1,cube(3)  ,ic)
        Bp(4) = b2(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b2(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b2(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b2(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b2(cube(1)+1,cube(2)+1,cube(3)+1,ic)

        res2(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr

      enddo

c     change Bx and By components sign
      res2(1)=-res2(1)
      res2(2)=-res2(2)

c     change back the x and y coordinate signs
      x=-x
      y=-y


c------------------------------------------------------------------------
c     
c     average the two maps results
c     
c------------------------------------------------------------------------

      do ic=1,3
        res(ic)=(res1(ic)+res2(ic))/2
      enddo


      return
      end
1	mocchiut	1.1	*************************************************************************
2			*
3			* Subroutine inter_B_inner.f
4			*
5			* it computes the magnetic field in a chosen point x,y,z inside the
6			* magnetic cavity, using a trilinear interpolation of
7			* B field measurements (read before by means of ./read_B.f)
8			* the value is computed for two different inner maps and then averaged
9			*
10			* needs:
11			* - ../common/common_B_inner.f
12			*
13			* input: coordinates in m
14			* output: magnetic field in T
15			*
16			*************************************************************************
17
18			subroutine inter_B_inner(x,y,z,res) !coordinates in m, magnetic field in T
19
20			implicit double precision (a-h,o-z)
21			include './common_B_inner.for'
22
23
24			c------------------------------------------------------------------------
25			c
26			c local variables
27			c
28			c------------------------------------------------------------------------
29
30			real*8 x,y,z !point of interpolation
31			real*8 res(3) !interpolated B components: res = (Bx, By, Bz)
32			real*8 res1(3),res2(3) !interpolated B components for the two maps
33
34			integer ic !index for B components:
35			! ic=1 ---> Bx
36			! ic=2 ---> By
37			! ic=3 ---> Bz
38
39			integer cube(3) !vector of indexes identifying the cube
40			! containing the point of interpolation
41			! (see later...)
42
43			real*8 xl,xh,yl,yh,zl,zh !cube vertexes coordinates
44
45			real*8 xr,yr,zr !reduced variables (coordinates of the
46			! point of interpolation inside the cube)
47
48			real*8 Bp(8) !vector of values of B component
49			! being computed, on the eight cube vertexes
50
51
52			c------------------------------------------------------------------------
53			c
54			c * FIRST MAP *
55			c
56			c------------------------------------------------------------------------
57
58			do ic=1,3 !loops on the three B components
59
60			c------------------------------------------------------------------------
61			c
62			c chooses the coordinates interval containing the input point
63			c
64			c------------------------------------------------------------------------
65			c e.g.:
66			c
67			c x1 x2 x3 x4 x5...
68			c \|-----\|-+---\|-----\|-----\|----
69			c ~~~~~~~~x
70			c
71			c in this case the right interval is identified by indexes 2-3, so the
72			c value assigned to cube variable is 2
73
74			cube(1)=INT((nx-1)*(x-px1min(ic))/(px1max(ic)-px1min(ic))) +1
75			cube(2)=INT((ny-1)*(y-py1min(ic))/(py1max(ic)-py1min(ic))) +1
76			cube(3)=INT((nz-1)*(z-pz1min(ic))/(pz1max(ic)-pz1min(ic))) +1
77
78			c------------------------------------------------------------------------
79			c
80			c if the point falls beyond the extremes of the grid...
81			c
82			c------------------------------------------------------------------------
83			c
84			c ~~~~~~~~~~x1 x2 x3...
85			c - - + - - \|-----\|-----\|----
86			c ~~~~x
87			c
88			c in the case cube = 1
89			c
90			c
91			c ...nx-2 nx-1 nx
92			c ----\|-----\|-----\| - - - + - -
93			c ~~~~~~~~~~~~~~~~~~~~~~~~x
94			c
95			c in this case cube = nx-1
96
97			if (cube(1).le.0) cube(1) = 1
98			if (cube(2).le.0) cube(2) = 1
99			if (cube(3).le.0) cube(3) = 1
100			if (cube(1).ge.nx) cube(1) = nx-1
101			if (cube(2).ge.ny) cube(2) = ny-1
102			if (cube(3).ge.nz) cube(3) = nz-1
103
104
105			c------------------------------------------------------------------------
106			c
107			c temporary variables definition for field computation
108			c
109			c------------------------------------------------------------------------
110
111			xl = px1(cube(1),ic) !X coordinates of cube vertexes
112			xh = px1(cube(1)+1,ic)
113
114			yl = py1(cube(2),ic) !Y coordinates of cube vertexes
115			yh = py1(cube(2)+1,ic)
116
117			zl = pz1(cube(3),ic) !Z coordinates of cube vertexes
118			zh = pz1(cube(3)+1,ic)
119
120			xr = (x-xl) / (xh-xl) !reduced variables
121			yr = (y-yl) / (yh-yl)
122			zr = (z-zl) / (zh-zl)
123
124			Bp(1) = b1(cube(1) ,cube(2) ,cube(3) ,ic) !ic-th component of B
125			Bp(2) = b1(cube(1)+1,cube(2) ,cube(3) ,ic) ! on the eight cube
126			Bp(3) = b1(cube(1) ,cube(2)+1,cube(3) ,ic) ! vertexes
127			Bp(4) = b1(cube(1)+1,cube(2)+1,cube(3) ,ic)
128			Bp(5) = b1(cube(1) ,cube(2) ,cube(3)+1,ic)
129			Bp(6) = b1(cube(1)+1,cube(2) ,cube(3)+1,ic)
130			Bp(7) = b1(cube(1) ,cube(2)+1,cube(3)+1,ic)
131			Bp(8) = b1(cube(1)+1,cube(2)+1,cube(3)+1,ic)
132
133			c------------------------------------------------------------------------
134			c
135			c computes interpolated ic-th component of B in (x,y,z)
136			c
137			c------------------------------------------------------------------------
138
139			res1(ic) =
140			+ Bp(1)(1-xr)(1-yr)*(1-zr) +
141			+ Bp(2)xr(1-yr)*(1-zr) +
142			+ Bp(3)(1-xr)yr*(1-zr) +
143			+ Bp(4)xryr*(1-zr) +
144			+ Bp(5)(1-xr)(1-yr)*zr +
145			+ Bp(6)xr(1-yr)*zr +
146			+ Bp(7)(1-xr)yr*zr +
147			+ Bp(8)xryr*zr
148
149
150			enddo
151
152			c------------------------------------------------------------------------
153			c
154			c * SECOND MAP *
155			c
156			c------------------------------------------------------------------------
157
158			c second map is rotated by 180 degree along the Z axis. so change sign
159			c of x and y coordinates and then change sign to Bx and By components
160			c to obtain the correct result
161
162			x=-x
163			y=-y
164
165			do ic=1,3
166
167			cube(1)=INT((nx-1)*(x-px2min(ic))/(px2max(ic)-px2min(ic))) +1
168			cube(2)=INT((ny-1)*(y-py2min(ic))/(py2max(ic)-py2min(ic))) +1
169			cube(3)=INT((nz-1)*(z-pz2min(ic))/(pz2max(ic)-pz2min(ic))) +1
170
171			if (cube(1).le.0) cube(1) = 1
172			if (cube(2).le.0) cube(2) = 1
173			if (cube(3).le.0) cube(3) = 1
174			if (cube(1).ge.nx) cube(1) = nx-1
175			if (cube(2).ge.ny) cube(2) = ny-1
176			if (cube(3).ge.nz) cube(3) = nz-1
177
178			xl = px2(cube(1),ic)
179			xh = px2(cube(1)+1,ic)
180
181			yl = py2(cube(2),ic)
182			yh = py2(cube(2)+1,ic)
183
184			zl = pz2(cube(3),ic)
185			zh = pz2(cube(3)+1,ic)
186
187			xr = (x-xl) / (xh-xl)
188			yr = (y-yl) / (yh-yl)
189			zr = (z-zl) / (zh-zl)
190
191			Bp(1) = b2(cube(1) ,cube(2) ,cube(3) ,ic)
192			Bp(2) = b2(cube(1)+1,cube(2) ,cube(3) ,ic)
193			Bp(3) = b2(cube(1) ,cube(2)+1,cube(3) ,ic)
194			Bp(4) = b2(cube(1)+1,cube(2)+1,cube(3) ,ic)
195			Bp(5) = b2(cube(1) ,cube(2) ,cube(3)+1,ic)
196			Bp(6) = b2(cube(1)+1,cube(2) ,cube(3)+1,ic)
197			Bp(7) = b2(cube(1) ,cube(2)+1,cube(3)+1,ic)
198			Bp(8) = b2(cube(1)+1,cube(2)+1,cube(3)+1,ic)
199
200			res2(ic) =
201			+ Bp(1)(1-xr)(1-yr)*(1-zr) +
202			+ Bp(2)xr(1-yr)*(1-zr) +
203			+ Bp(3)(1-xr)yr*(1-zr) +
204			+ Bp(4)xryr*(1-zr) +
205			+ Bp(5)(1-xr)(1-yr)*zr +
206			+ Bp(6)xr(1-yr)*zr +
207			+ Bp(7)(1-xr)yr*zr +
208			+ Bp(8)xryr*zr
209
210			enddo
211
212			c change Bx and By components sign
213			res2(1)=-res2(1)
214			res2(2)=-res2(2)
215
216			c change back the x and y coordinate signs
217			x=-x
218			y=-y
219
220
221			c------------------------------------------------------------------------
222			c
223			c average the two maps results
224			c
225			c------------------------------------------------------------------------
226
227			do ic=1,3
228			res(ic)=(res1(ic)+res2(ic))/2
229			enddo
230
231
232			return
233			end