src/F77/bdll.f


***   
*     subroutine to compute the Integral(B dl) and 
*     L along the trajectory
***   


***   
*     
*     created by paolo papini, 11 oct 2005
*     
*     c$$$ max:  modified by massimo bongi, 11 oct 2005 (later...)
*     
***


      subroutine CalcBdL(Istep,BdL,L,IFAIL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      integer Istep,IFAIL
      real*8 BdL,L,Lstep
      real vvv(3),fff(3)
      real*8 f(3)
      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure
*     
*     set parameters for GRKUTA from ZINI to the first plane
*     
c$$$  max    IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
c$$$  max    IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
c$$$      IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
      IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
      IF(AL(5).EQ.0) CHARGE=1.
      VOUT(1)=AL(1)
      VOUT(2)=AL(2)
      VOUT(3)=ZINI              !  DBLE(Z0)-DBLE(ZSPEC)
      VOUT(4)=AL(3)*DCOS(AL(4))
      VOUT(5)=AL(3)*DSIN(AL(4))
      VOUT(6)=-1.*DSQRT(1.-AL(3)**2)
      IF(AL(5).NE.0.) VOUT(7)=DABS(1./AL(5))
      IF(AL(5).EQ.0.) VOUT(7)=1.E8
*     
*     first integrate the track from reference plane to first trk plane
*
*      zout=zm(1)     
      zout=zv(1)
      step=vout(3)-zout
 10   DO J=1,7
         VECT(J)=VOUT(J)
         VECTINI(J)=VOUT(J)
      ENDDO
 11   continue
      CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
      IF(VOUT(3).GT.VECT(3)) THEN
         IFAIL=1
         PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
         print*,'charge',charge
         print*,'vect',vect
         print*,'vout',vout
         print*,'step',step
         RETURN
      ENDIF
      Z=VOUT(3)
      IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 100
      IF(Z.GT.Zout+TOLL) GOTO 10      
      IF(Z.LE.Zout-TOLL) THEN
         STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
         DO J=1,7
            VECT(J)=VECTINI(J)
         ENDDO
         GOTO 11
      ENDIF
*
*     then integrate B o dL along the track, from the first to the last plane
*
 100  continue
      BdL=0.
      L=0.
*      dz=(zm(1)-zm(6))/float(Istep)
      dz = (zv(1)-zv(6))/float(Istep)
      zout = zv(6)
      
      do i=1,10*Istep
         BdLstep=0.
         Lstep=0.
         vvv(1)=sngl(vout(1))
         vvv(2)=sngl(vout(2))
         vvv(3)=sngl(vout(3))
         CALL GUFLD(VVV,FFF)
         do j=1,3
            f(j)=dble(fff(j))
         enddo
c         zout=vout(3)-dz
         step = dz
c 20      DO J=1,7
         DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 21      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
c$$$         BdLstep=BdLstep+
c$$$     $        (vout(1)-vect(1))*f(1)+
c$$$     $        (vout(2)-vect(2))*f(2)+
c$$$     $        (vout(3)-vect(3))*f(3)
         BdLstep=BdLstep+
     $        dsqrt(
     $        ((vout(2)-vect(2))*f(3)-(vout(3)-vect(3))*f(2))**2+
     $        ((vout(3)-vect(3))*f(1)-(vout(1)-vect(1))*f(3))**2+
     $        ((vout(1)-vect(1))*f(2)-(vout(2)-vect(2))*f(1))**2)
         Lstep=Lstep+
     $        dsqrt(
     $        (vout(1)-vect(1))**2+
     $        (vout(2)-vect(2))**2+
     $        (vout(3)-vect(3))**2)
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
c 200     continue
         continue
         BdL = BdL + (BdLstep*1.D-3) !Tesla * meters
         L = L + (Lstep*1.D-2) !meters
c$$$         print*,i,' #### ',BdL,L

         Z=VOUT(3)
         IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 2000 !end
         IF(Z.LE.Zout-TOLL) THEN
            STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            BdL=BdL-(BdLstep*1.D-3) !Tesla * meters
            L=L-(Lstep*1.D-2)    !meters
            BdLstep=0.
            Lstep=0.
            GOTO 21
         ENDIF
      enddo
 2000 continue
c$$$      print*,' #### ',BdL
c$$$      print*,' #### ',L
c$$$      print*,' #### ',xgood
c$$$      print*,' #### ',ygood
c$$$      print*,'========================================='
      return
      end


1
2	***
3	* subroutine to compute the Integral(B dl) and
4	* L along the trajectory
5	***
6
7
8	***
9	*
10	* created by paolo papini, 11 oct 2005
11	*
12	* c$$$ max: modified by massimo bongi, 11 oct 2005 (later...)
13	*
14	***
15
16
17	subroutine CalcBdL(Istep,BdL,L,IFAIL)
18	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
19	integer Istep,IFAIL
20	real*8 BdL,L,Lstep
21	real vvv(3),fff(3)
22	real*8 f(3)
23	include 'commontracker.f' !tracker general common
24	include 'common_mini_2.f' !common for the tracking procedure
25	*
26	* set parameters for GRKUTA from ZINI to the first plane
27	*
28	c$$$ max IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
29	c$$$ max IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
30	c$$$ IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
31	IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
32	IF(AL(5).EQ.0) CHARGE=1.
33	VOUT(1)=AL(1)
34	VOUT(2)=AL(2)
35	VOUT(3)=ZINI ! DBLE(Z0)-DBLE(ZSPEC)
36	VOUT(4)=AL(3)*DCOS(AL(4))
37	VOUT(5)=AL(3)*DSIN(AL(4))
38	VOUT(6)=-1.DSQRT(1.-AL(3)*2)
39	IF(AL(5).NE.0.) VOUT(7)=DABS(1./AL(5))
40	IF(AL(5).EQ.0.) VOUT(7)=1.E8
41	*
42	* first integrate the track from reference plane to first trk plane
43	*
44	* zout=zm(1)
45	zout=zv(1)
46	step=vout(3)-zout
47	10 DO J=1,7
48	VECT(J)=VOUT(J)
49	VECTINI(J)=VOUT(J)
50	ENDDO
51	11 continue
52	CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
53	IF(VOUT(3).GT.VECT(3)) THEN
54	IFAIL=1
55	PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
56	print*,'charge',charge
57	print*,'vect',vect
58	print*,'vout',vout
59	print*,'step',step
60	RETURN
61	ENDIF
62	Z=VOUT(3)
63	IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 100
64	IF(Z.GT.Zout+TOLL) GOTO 10
65	IF(Z.LE.Zout-TOLL) THEN
66	STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
67	DO J=1,7
68	VECT(J)=VECTINI(J)
69	ENDDO
70	GOTO 11
71	ENDIF
72	*
73	* then integrate B o dL along the track, from the first to the last plane
74	*
75	100 continue
76	BdL=0.
77	L=0.
78	* dz=(zm(1)-zm(6))/float(Istep)
79	dz = (zv(1)-zv(6))/float(Istep)
80	zout = zv(6)
81
82	do i=1,10*Istep
83	BdLstep=0.
84	Lstep=0.
85	vvv(1)=sngl(vout(1))
86	vvv(2)=sngl(vout(2))
87	vvv(3)=sngl(vout(3))
88	CALL GUFLD(VVV,FFF)
89	do j=1,3
90	f(j)=dble(fff(j))
91	enddo
92	c zout=vout(3)-dz
93	step = dz
94	c 20 DO J=1,7
95	DO J=1,7
96	VECT(J)=VOUT(J)
97	VECTINI(J)=VOUT(J)
98	ENDDO
99	21 continue
100	CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
101	c$$$ BdLstep=BdLstep+
102	c$$$ $ (vout(1)-vect(1))*f(1)+
103	c$$$ $ (vout(2)-vect(2))*f(2)+
104	c$$$ $ (vout(3)-vect(3))*f(3)
105	BdLstep=BdLstep+
106	$ dsqrt(
107	$ ((vout(2)-vect(2))f(3)-(vout(3)-vect(3))f(2))**2+
108	$ ((vout(3)-vect(3))f(1)-(vout(1)-vect(1))f(3))**2+
109	$ ((vout(1)-vect(1))f(2)-(vout(2)-vect(2))f(1))**2)
110	Lstep=Lstep+
111	$ dsqrt(
112	$ (vout(1)-vect(1))**2+
113	$ (vout(2)-vect(2))**2+
114	$ (vout(3)-vect(3))**2)
115	IF(VOUT(3).GT.VECT(3)) THEN
116	IFAIL=1
117	PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
118	print*,'charge',charge
119	print*,'vect',vect
120	print*,'vout',vout
121	print*,'step',step
122	RETURN
123	ENDIF
124	c 200 continue
125	continue
126	BdL = BdL + (BdLstep1.D-3) !Tesla meters
127	L = L + (Lstep*1.D-2) !meters
128	c$$$ print*,i,' #### ',BdL,L
129
130	Z=VOUT(3)
131	IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 2000 !end
132	IF(Z.LE.Zout-TOLL) THEN
133	STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
134	DO J=1,7
135	VECT(J)=VECTINI(J)
136	ENDDO
137	BdL=BdL-(BdLstep1.D-3) !Tesla meters
138	L=L-(Lstep*1.D-2) !meters
139	BdLstep=0.
140	Lstep=0.
141	GOTO 21
142	ENDIF
143	enddo
144	2000 continue
145	c$$$ print*,' #### ',BdL
146	c$$$ print*,' #### ',L
147	c$$$ print*,' #### ',xgood
148	c$$$ print*,' #### ',ygood
149	c$$$ print*,'========================================='
150	return
151	end
152
153