************************************************************************* * * Subroutine inter_B_inner.f (from tracker software analysis) * * 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 "gpfield.inc" 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