flight/src/inter_B.for

*************************************************************************
*     
*     Subroutine inter_B.f
*     
*     it computes the magnetic field in a chosen point x,y,z inside or 
*     outside the magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     if the point falls outside the interpolation volume, set the field to 0
*     
*     needs:
*     - ../common/common_B_inner.f common file for the inner magnetic field map
*     - ./inter_B_inner.f common file for the inner magnetic field map
*     - ../common/common_B_outer.f common file for the outer magnetic field map
*     - ./inter_B_outer.f common file for the outer magnetic field map
*     
*     to be called after ./read_B.f (magnetic field map reading subroutine)
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

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

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


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

      real*8 x,y,z              !point of interest
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)

      real*8 zl,zu
      real*8 resu(3),resl(3)

c------------------------------------------------------------------------
c     
c     set the field outside the interpolation volume to be 0
c     
c------------------------------------------------------------------------

      do ip=1,3
        res(ip)=0.
      enddo


c------------------------------------------------------------------------
c     
c     check if the point falls inside the interpolation volumes
c     
c------------------------------------------------------------------------

*     -----------------------
*     INNER MAP
*     -----------------------
      if(
     $          (x.ge.edgexmin).and.(x.le.edgexmax)
     $     .and.(y.ge.edgeymin).and.(y.le.edgeymax)
     $     .and.(z.ge.edgezmin).and.(z.le.edgezmax)
     $     ) then

        call inter_B_inner(x,y,z,res)
c        print*,'INNER - ',z,res

      endif

*     -----------------------
*     OUTER MAP
*     -----------------------
      if(
     $         ((x.ge.edgeuxmin).and.(x.le.edgeuxmax)
     $     .and.(y.ge.edgeuymin).and.(y.le.edgeuymax)
     $     .and.(z.ge.edgeuzmin).and.(z.le.edgeuzmax))
     $     .or.
     $         ((x.ge.edgelxmin).and.(x.le.edgelxmax)
     $     .and.(y.ge.edgelymin).and.(y.le.edgelymax)
     $     .and.(z.ge.edgelzmin).and.(z.le.edgelzmax))
     $     ) then
         
        call inter_B_outer(x,y,z,res)
c        res(2)=res(2)*10
c        print*,'OUTER - ',z,res

      endif

*     --------------------------------
*     GAP between INNER and OUTER MAPS
*     --------------------------------
      if(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgezmax).and.(z.lt.edgeuzmin)
     $     )then
         
         zu = edgeuzmin
         zl = edgezmax
         call inter_B_inner(x,y,zl,resu)
         call inter_B_outer(x,y,zu,resl)

         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP U - ',z,res

      elseif(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgelzmax).and.(z.lt.edgezmin)
     $     ) then
         
        
         zu = edgezmin
         zl = edgelzmax
         call inter_B_inner(x,y,zu,resu)
         call inter_B_outer(x,y,zl,resl)
         
         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP D - ',z,res

      endif

      return
      end


      include './inter_B_inner.for'
      include './inter_B_outer.for'
1	*************************************************************************
2	*
3	* Subroutine inter_B.f
4	*
5	* it computes the magnetic field in a chosen point x,y,z inside or
6	* outside the magnetic cavity, using a trilinear interpolation of
7	* B field measurements (read before by means of ./read_B.f)
8	* if the point falls outside the interpolation volume, set the field to 0
9	*
10	* needs:
11	* - ../common/common_B_inner.f common file for the inner magnetic field map
12	* - ./inter_B_inner.f common file for the inner magnetic field map
13	* - ../common/common_B_outer.f common file for the outer magnetic field map
14	* - ./inter_B_outer.f common file for the outer magnetic field map
15	*
16	* to be called after ./read_B.f (magnetic field map reading subroutine)
17	*
18	* input: coordinates in m
19	* output: magnetic field in T
20	*
21	*************************************************************************
22
23	subroutine inter_B(x,y,z,res) !coordinates in m, magnetic field in T
24
25	implicit double precision (a-h,o-z)
26	include './common_B_inner.for'
27	include './common_B_outer.for'
28
29
30	c------------------------------------------------------------------------
31	c
32	c local variables
33	c
34	c------------------------------------------------------------------------
35
36	real*8 x,y,z !point of interest
37	real*8 res(3) !interpolated B components: res = (Bx, By, Bz)
38
39	real*8 zl,zu
40	real*8 resu(3),resl(3)
41
42	c------------------------------------------------------------------------
43	c
44	c set the field outside the interpolation volume to be 0
45	c
46	c------------------------------------------------------------------------
47
48	do ip=1,3
49	res(ip)=0.
50	enddo
51
52
53	c------------------------------------------------------------------------
54	c
55	c check if the point falls inside the interpolation volumes
56	c
57	c------------------------------------------------------------------------
58
59	* -----------------------
60	* INNER MAP
61	* -----------------------
62	if(
63	$ (x.ge.edgexmin).and.(x.le.edgexmax)
64	$ .and.(y.ge.edgeymin).and.(y.le.edgeymax)
65	$ .and.(z.ge.edgezmin).and.(z.le.edgezmax)
66	$ ) then
67
68	call inter_B_inner(x,y,z,res)
69	c print*,'INNER - ',z,res
70
71	endif
72
73	* -----------------------
74	* OUTER MAP
75	* -----------------------
76	if(
77	$ ((x.ge.edgeuxmin).and.(x.le.edgeuxmax)
78	$ .and.(y.ge.edgeuymin).and.(y.le.edgeuymax)
79	$ .and.(z.ge.edgeuzmin).and.(z.le.edgeuzmax))
80	$ .or.
81	$ ((x.ge.edgelxmin).and.(x.le.edgelxmax)
82	$ .and.(y.ge.edgelymin).and.(y.le.edgelymax)
83	$ .and.(z.ge.edgelzmin).and.(z.le.edgelzmax))
84	$ ) then
85
86	call inter_B_outer(x,y,z,res)
87	c res(2)=res(2)*10
88	c print*,'OUTER - ',z,res
89
90	endif
91
92	* --------------------------------
93	* GAP between INNER and OUTER MAPS
94	* --------------------------------
95	if(
96	$ (x.gt.edgexmin).and.(x.lt.edgexmax)
97	$ .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
98	$ .and.(z.gt.edgezmax).and.(z.lt.edgeuzmin)
99	$ )then
100
101	zu = edgeuzmin
102	zl = edgezmax
103	call inter_B_inner(x,y,zl,resu)
104	call inter_B_outer(x,y,zu,resl)
105
106	do i=1,3
107	res(i) = z * ((resu(i)-resl(i))/(zu-zl))
108	$ + resu(i)
109	$ - zu * ((resu(i)-resl(i))/(zu-zl))
110	enddo
111	c print*,'GAP U - ',z,res
112
113	elseif(
114	$ (x.gt.edgexmin).and.(x.lt.edgexmax)
115	$ .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
116	$ .and.(z.gt.edgelzmax).and.(z.lt.edgezmin)
117	$ ) then
118
119
120	zu = edgezmin
121	zl = edgelzmax
122	call inter_B_inner(x,y,zu,resu)
123	call inter_B_outer(x,y,zl,resl)
124
125	do i=1,3
126	res(i) = z * ((resu(i)-resl(i))/(zu-zl))
127	$ + resu(i)
128	$ - zu * ((resu(i)-resl(i))/(zu-zl))
129	enddo
130	c print*,'GAP D - ',z,res
131
132	endif
133
134	return
135	end
136
137
138	include './inter_B_inner.for'
139	include './inter_B_outer.for'