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