/[PAMELA software]/DarthVader/TrackerLevel2/src/F77/interB.f
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/F77/interB.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Fri May 19 13:15:55 2006 UTC (18 years, 7 months ago) by mocchiut
Branch: DarthVader, MAIN
CVS Tags: v0r01, v0r02, v5r00, v4r00, v2r01, v10RED, v1r01beta, start, v9r00, v9r01, v1r00, v1r01, v10REDr01, v3r04, v3r05, v3r06, v3r00, v3r03, v6r01, v6r00, v2r00BETA, HEAD
Changes since 1.1: +0 -0 lines
Imported sources

1 mocchiut 1.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_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_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.f'
27    
28    
29     c------------------------------------------------------------------------
30     c
31     c local variables
32     c
33     c------------------------------------------------------------------------
34    
35     real*8 x,y,z !point of interest
36     real*8 res(3) !interpolated B components: res = (Bx, By, Bz)
37    
38     real*8 zl,zu
39     real*8 resu(3),resl(3)
40    
41     c------------------------------------------------------------------------
42     c
43     c set the field outside the interpolation volume to be 0
44     c
45     c------------------------------------------------------------------------
46    
47     do ip=1,3
48     res(ip)=0.
49     enddo
50    
51    
52     c------------------------------------------------------------------------
53     c
54     c check if the point falls inside the interpolation volumes
55     c
56     c------------------------------------------------------------------------
57    
58     * -----------------------
59     * INNER MAP
60     * -----------------------
61     if(
62     $ (x.ge.edgexmin).and.(x.le.edgexmax)
63     $ .and.(y.ge.edgeymin).and.(y.le.edgeymax)
64     $ .and.(z.ge.edgezmin).and.(z.le.edgezmax)
65     $ ) then
66    
67     call inter_B_inner(x,y,z,res)
68     c print*,'INNER - ',z,res
69    
70     endif
71    
72     * -----------------------
73     * OUTER MAP
74     * -----------------------
75     if(
76     $ ((x.ge.edgeuxmin).and.(x.le.edgeuxmax)
77     $ .and.(y.ge.edgeuymin).and.(y.le.edgeuymax)
78     $ .and.(z.ge.edgeuzmin).and.(z.le.edgeuzmax))
79     $ .or.
80     $ ((x.ge.edgelxmin).and.(x.le.edgelxmax)
81     $ .and.(y.ge.edgelymin).and.(y.le.edgelymax)
82     $ .and.(z.ge.edgelzmin).and.(z.le.edgelzmax))
83     $ ) then
84    
85     call inter_B_outer(x,y,z,res)
86     c res(2)=res(2)*10
87     c print*,'OUTER - ',z,res
88    
89     endif
90    
91     * --------------------------------
92     * GAP between INNER and OUTER MAPS
93     * --------------------------------
94     if(
95     $ (x.gt.edgexmin).and.(x.lt.edgexmax)
96     $ .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
97     $ .and.(z.gt.edgezmax).and.(z.lt.edgeuzmin)
98     $ )then
99    
100     zu = edgeuzmin
101     zl = edgezmax
102     call inter_B_inner(x,y,zl,resu)
103     call inter_B_outer(x,y,zu,resl)
104    
105     do i=1,3
106     res(i) = z * ((resu(i)-resl(i))/(zu-zl))
107     $ + resu(i)
108     $ - zu * ((resu(i)-resl(i))/(zu-zl))
109     enddo
110     c print*,'GAP U - ',z,res
111    
112     elseif(
113     $ (x.gt.edgexmin).and.(x.lt.edgexmax)
114     $ .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
115     $ .and.(z.gt.edgelzmax).and.(z.lt.edgezmin)
116     $ ) then
117    
118    
119     zu = edgezmin
120     zl = edgelzmax
121     call inter_B_inner(x,y,zu,resu)
122     call inter_B_outer(x,y,zl,resl)
123    
124     do i=1,3
125     res(i) = z * ((resu(i)-resl(i))/(zu-zl))
126     $ + resu(i)
127     $ - zu * ((resu(i)-resl(i))/(zu-zl))
128     enddo
129     c print*,'GAP D - ',z,res
130    
131     endif
132    
133     return
134     end
135    
136    
137     *************************************************************************
138     *
139     * Subroutine inter_B_inner.f
140     *
141     * it computes the magnetic field in a chosen point x,y,z inside the
142     * magnetic cavity, using a trilinear interpolation of
143     * B field measurements (read before by means of ./read_B.f)
144     * the value is computed for two different inner maps and then averaged
145     *
146     * needs:
147     * - ../common/common_B_inner.f
148     *
149     * input: coordinates in m
150     * output: magnetic field in T
151     *
152     *************************************************************************
153    
154     subroutine inter_B_inner(x,y,z,res) !coordinates in m, magnetic field in T
155    
156     implicit double precision (a-h,o-z)
157     include 'common_B.f'
158    
159    
160     c------------------------------------------------------------------------
161     c
162     c local variables
163     c
164     c------------------------------------------------------------------------
165    
166     real*8 x,y,z !point of interpolation
167     real*8 res(3) !interpolated B components: res = (Bx, By, Bz)
168     real*8 res1(3),res2(3) !interpolated B components for the two maps
169    
170     integer ic !index for B components:
171     ! ic=1 ---> Bx
172     ! ic=2 ---> By
173     ! ic=3 ---> Bz
174    
175     integer cube(3) !vector of indexes identifying the cube
176     ! containing the point of interpolation
177     ! (see later...)
178    
179     real*8 xl,xh,yl,yh,zl,zh !cube vertexes coordinates
180    
181     real*8 xr,yr,zr !reduced variables (coordinates of the
182     ! point of interpolation inside the cube)
183    
184     real*8 Bp(8) !vector of values of B component
185     ! being computed, on the eight cube vertexes
186    
187    
188     c------------------------------------------------------------------------
189     c
190     c *** FIRST MAP ***
191     c
192     c------------------------------------------------------------------------
193    
194     do ic=1,3 !loops on the three B components
195    
196     c------------------------------------------------------------------------
197     c
198     c chooses the coordinates interval containing the input point
199     c
200     c------------------------------------------------------------------------
201     c e.g.:
202     c
203     c x1 x2 x3 x4 x5...
204     c |-----|-+---|-----|-----|----
205     c ~~~~~~~~x
206     c
207     c in this case the right interval is identified by indexes 2-3, so the
208     c value assigned to cube variable is 2
209    
210     cube(1)=INT((nx-1)*(x-px1min(ic))/(px1max(ic)-px1min(ic))) +1
211     cube(2)=INT((ny-1)*(y-py1min(ic))/(py1max(ic)-py1min(ic))) +1
212     cube(3)=INT((nz-1)*(z-pz1min(ic))/(pz1max(ic)-pz1min(ic))) +1
213    
214     c------------------------------------------------------------------------
215     c
216     c if the point falls beyond the extremes of the grid...
217     c
218     c------------------------------------------------------------------------
219     c
220     c ~~~~~~~~~~x1 x2 x3...
221     c - - + - - |-----|-----|----
222     c ~~~~x
223     c
224     c in the case cube = 1
225     c
226     c
227     c ...nx-2 nx-1 nx
228     c ----|-----|-----| - - - + - -
229     c ~~~~~~~~~~~~~~~~~~~~~~~~x
230     c
231     c in this case cube = nx-1
232    
233     if (cube(1).le.0) cube(1) = 1
234     if (cube(2).le.0) cube(2) = 1
235     if (cube(3).le.0) cube(3) = 1
236     if (cube(1).ge.nx) cube(1) = nx-1
237     if (cube(2).ge.ny) cube(2) = ny-1
238     if (cube(3).ge.nz) cube(3) = nz-1
239    
240    
241     c------------------------------------------------------------------------
242     c
243     c temporary variables definition for field computation
244     c
245     c------------------------------------------------------------------------
246    
247     xl = px1(cube(1),ic) !X coordinates of cube vertexes
248     xh = px1(cube(1)+1,ic)
249    
250     yl = py1(cube(2),ic) !Y coordinates of cube vertexes
251     yh = py1(cube(2)+1,ic)
252    
253     zl = pz1(cube(3),ic) !Z coordinates of cube vertexes
254     zh = pz1(cube(3)+1,ic)
255    
256     xr = (x-xl) / (xh-xl) !reduced variables
257     yr = (y-yl) / (yh-yl)
258     zr = (z-zl) / (zh-zl)
259    
260     Bp(1) = b1(cube(1) ,cube(2) ,cube(3) ,ic) !ic-th component of B
261     Bp(2) = b1(cube(1)+1,cube(2) ,cube(3) ,ic) ! on the eight cube
262     Bp(3) = b1(cube(1) ,cube(2)+1,cube(3) ,ic) ! vertexes
263     Bp(4) = b1(cube(1)+1,cube(2)+1,cube(3) ,ic)
264     Bp(5) = b1(cube(1) ,cube(2) ,cube(3)+1,ic)
265     Bp(6) = b1(cube(1)+1,cube(2) ,cube(3)+1,ic)
266     Bp(7) = b1(cube(1) ,cube(2)+1,cube(3)+1,ic)
267     Bp(8) = b1(cube(1)+1,cube(2)+1,cube(3)+1,ic)
268    
269     c------------------------------------------------------------------------
270     c
271     c computes interpolated ic-th component of B in (x,y,z)
272     c
273     c------------------------------------------------------------------------
274    
275     res1(ic) =
276     + Bp(1)*(1-xr)*(1-yr)*(1-zr) +
277     + Bp(2)*xr*(1-yr)*(1-zr) +
278     + Bp(3)*(1-xr)*yr*(1-zr) +
279     + Bp(4)*xr*yr*(1-zr) +
280     + Bp(5)*(1-xr)*(1-yr)*zr +
281     + Bp(6)*xr*(1-yr)*zr +
282     + Bp(7)*(1-xr)*yr*zr +
283     + Bp(8)*xr*yr*zr
284    
285    
286     enddo
287    
288     c------------------------------------------------------------------------
289     c
290     c *** SECOND MAP ***
291     c
292     c------------------------------------------------------------------------
293    
294     c second map is rotated by 180 degree along the Z axis. so change sign
295     c of x and y coordinates and then change sign to Bx and By components
296     c to obtain the correct result
297    
298     x=-x
299     y=-y
300    
301     do ic=1,3
302    
303     cube(1)=INT((nx-1)*(x-px2min(ic))/(px2max(ic)-px2min(ic))) +1
304     cube(2)=INT((ny-1)*(y-py2min(ic))/(py2max(ic)-py2min(ic))) +1
305     cube(3)=INT((nz-1)*(z-pz2min(ic))/(pz2max(ic)-pz2min(ic))) +1
306    
307     if (cube(1).le.0) cube(1) = 1
308     if (cube(2).le.0) cube(2) = 1
309     if (cube(3).le.0) cube(3) = 1
310     if (cube(1).ge.nx) cube(1) = nx-1
311     if (cube(2).ge.ny) cube(2) = ny-1
312     if (cube(3).ge.nz) cube(3) = nz-1
313    
314     xl = px2(cube(1),ic)
315     xh = px2(cube(1)+1,ic)
316    
317     yl = py2(cube(2),ic)
318     yh = py2(cube(2)+1,ic)
319    
320     zl = pz2(cube(3),ic)
321     zh = pz2(cube(3)+1,ic)
322    
323     xr = (x-xl) / (xh-xl)
324     yr = (y-yl) / (yh-yl)
325     zr = (z-zl) / (zh-zl)
326    
327     Bp(1) = b2(cube(1) ,cube(2) ,cube(3) ,ic)
328     Bp(2) = b2(cube(1)+1,cube(2) ,cube(3) ,ic)
329     Bp(3) = b2(cube(1) ,cube(2)+1,cube(3) ,ic)
330     Bp(4) = b2(cube(1)+1,cube(2)+1,cube(3) ,ic)
331     Bp(5) = b2(cube(1) ,cube(2) ,cube(3)+1,ic)
332     Bp(6) = b2(cube(1)+1,cube(2) ,cube(3)+1,ic)
333     Bp(7) = b2(cube(1) ,cube(2)+1,cube(3)+1,ic)
334     Bp(8) = b2(cube(1)+1,cube(2)+1,cube(3)+1,ic)
335    
336     res2(ic) =
337     + Bp(1)*(1-xr)*(1-yr)*(1-zr) +
338     + Bp(2)*xr*(1-yr)*(1-zr) +
339     + Bp(3)*(1-xr)*yr*(1-zr) +
340     + Bp(4)*xr*yr*(1-zr) +
341     + Bp(5)*(1-xr)*(1-yr)*zr +
342     + Bp(6)*xr*(1-yr)*zr +
343     + Bp(7)*(1-xr)*yr*zr +
344     + Bp(8)*xr*yr*zr
345    
346     enddo
347    
348     c change Bx and By components sign
349     res2(1)=-res2(1)
350     res2(2)=-res2(2)
351    
352     c change back the x and y coordinate signs
353     x=-x
354     y=-y
355    
356    
357     c------------------------------------------------------------------------
358     c
359     c average the two maps results
360     c
361     c------------------------------------------------------------------------
362    
363     do ic=1,3
364     res(ic)=(res1(ic)+res2(ic))/2
365     enddo
366    
367    
368     return
369     end
370     *************************************************************************
371     *
372     * Subroutine inter_B_outer.f
373     *
374     * it computes the magnetic field in a chosen point x,y,z OUTSIDE the
375     * magnetic cavity, using a trilinear interpolation of
376     * B field measurements (read before by means of ./read_B.f)
377     * the value is computed for the outer map
378     *
379     * needs:
380     * - ../common/common_B_outer.f
381     *
382     * input: coordinates in m
383     * output: magnetic field in T
384     *
385     *************************************************************************
386    
387     subroutine inter_B_outer(x,y,z,res) !coordinates in m, magnetic field in T
388    
389     implicit double precision (a-h,o-z)
390     include 'common_B.f'
391    
392    
393     c------------------------------------------------------------------------
394     c
395     c local variables
396     c
397     c------------------------------------------------------------------------
398    
399     real*8 x,y,z !point of interpolation
400     real*8 res(3) !interpolated B components: res = (Bx, By, Bz)
401     real*8 zin
402    
403     integer ic
404     c !index for B components:
405     c ! ic=1 ---> Bx
406     c ! ic=2 ---> By
407     c ! ic=3 ---> Bz
408    
409     integer cube(3)
410     c !vector of indexes identifying the cube
411     c ! containing the point of interpolation
412     c ! (see later...)
413    
414     real*8 xl,xh,yl,yh,zl,zh !cube vertexes coordinates
415    
416     real*8 xr,yr,zr
417     c !reduced variables (coordinates of the
418     c ! point of interpolation inside the cube)
419    
420     real*8 Bp(8)
421     c !vector of values of B component
422     c ! being computed, on the eight cube vertexes
423    
424    
425     c LOWER MAP
426     c ---> up/down simmetry
427     zin=z
428     if(zin.le.edgelzmax)z=-z
429    
430     c------------------------------------------------------------------------
431     c
432     c *** MAP ***
433     c
434     c------------------------------------------------------------------------
435    
436     do ic=1,3 !loops on the three B components
437    
438     c------------------------------------------------------------------------
439     c
440     c chooses the coordinates interval containing the input point
441     c
442     c------------------------------------------------------------------------
443     c e.g.:
444     c
445     c x1 x2 x3 x4 x5... xN
446     c |-----|-+---|-----|-----|---- ... ----|-----|
447     c ~~~~~~~~x
448     c
449     c in this case the right interval is identified by indexes 2-3, so the
450     c value assigned to cube variable is 2
451    
452     cube(1)=INT((nox-1)*(x-poxmin(ic))/(poxmax(ic)-poxmin(ic))) +1
453     cube(2)=INT((noy-1)*(y-poymin(ic))/(poymax(ic)-poymin(ic))) +1
454     cube(3)=INT((noz-1)*(z-pozmin(ic))/(pozmax(ic)-pozmin(ic))) +1
455    
456     c------------------------------------------------------------------------
457     c
458     c if the point falls beyond the extremes of the grid...
459     c
460     c------------------------------------------------------------------------
461     c
462     c ~~~~~~~~~~x1 x2 x3...
463     c - - + - - |-----|-----|----
464     c ~~~~x
465     c
466     c in the case cube = 1
467     c
468     c
469     c ...nx-2 nx-1 nx
470     c ----|-----|-----| - - - + - -
471     c ~~~~~~~~~~~~~~~~~~~~~~~~x
472     c
473     c in this case cube = nx-1
474    
475     if (cube(1).le.0) cube(1) = 1
476     if (cube(2).le.0) cube(2) = 1
477     if (cube(3).le.0) cube(3) = 1
478     if (cube(1).ge.nox) cube(1) = nox-1
479     if (cube(2).ge.noy) cube(2) = noy-1
480     if (cube(3).ge.noz) cube(3) = noz-1
481    
482    
483     c------------------------------------------------------------------------
484     c
485     c temporary variables definition for field computation
486     c
487     c------------------------------------------------------------------------
488    
489     xl = pox(cube(1),ic) !X coordinates of cube vertexes
490     xh = pox(cube(1)+1,ic)
491    
492     yl = poy(cube(2),ic) !Y coordinates of cube vertexes
493     yh = poy(cube(2)+1,ic)
494    
495     zl = poz(cube(3),ic) !Z coordinates of cube vertexes
496     zh = poz(cube(3)+1,ic)
497    
498     xr = (x-xl) / (xh-xl) !reduced variables
499     yr = (y-yl) / (yh-yl)
500     zr = (z-zl) / (zh-zl)
501    
502     Bp(1) = bo(cube(1) ,cube(2) ,cube(3) ,ic) !ic-th component of B
503     Bp(2) = bo(cube(1)+1,cube(2) ,cube(3) ,ic) ! on the eight cube
504     Bp(3) = bo(cube(1) ,cube(2)+1,cube(3) ,ic) ! vertexes
505     Bp(4) = bo(cube(1)+1,cube(2)+1,cube(3) ,ic)
506     Bp(5) = bo(cube(1) ,cube(2) ,cube(3)+1,ic)
507     Bp(6) = bo(cube(1)+1,cube(2) ,cube(3)+1,ic)
508     Bp(7) = bo(cube(1) ,cube(2)+1,cube(3)+1,ic)
509     Bp(8) = bo(cube(1)+1,cube(2)+1,cube(3)+1,ic)
510    
511     c------------------------------------------------------------------------
512     c
513     c computes interpolated ic-th component of B in (x,y,z)
514     c
515     c------------------------------------------------------------------------
516    
517     res(ic) =
518     + Bp(1)*(1-xr)*(1-yr)*(1-zr) +
519     + Bp(2)*xr*(1-yr)*(1-zr) +
520     + Bp(3)*(1-xr)*yr*(1-zr) +
521     + Bp(4)*xr*yr*(1-zr) +
522     + Bp(5)*(1-xr)*(1-yr)*zr +
523     + Bp(6)*xr*(1-yr)*zr +
524     + Bp(7)*(1-xr)*yr*zr +
525     + Bp(8)*xr*yr*zr
526    
527    
528     enddo
529    
530     c LOWER MAP
531     c ---> up/down simmetry
532     if(zin.le.edgelzmax)then
533     z=-z !back to initial ccoordinate
534     res(3)=-res(3) !invert BZ component
535     endif
536    
537     return
538     end

  ViewVC Help
Powered by ViewVC 1.1.23