source/analysis/BdL.f


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


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


      subroutine CalcBdL(Nstep,BdL,IFAIL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      real BdL,vvv(3),fff(3)
      real*8 f(3)
      include '../common/commontracker.f' !tracker general common
      include '../common/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
*
*      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
*
 100  continue
      BdL=0.
*      dz=(zm(1)-zm(6))/float(Nstep)
      dz = (zv(1)-zv(6))/float(Nstep)
      zout = zv(6)
      
      do i=1,10*Nstep
         BdLstep=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
 20      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+
     $        sqrt(
     $        ((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
     $        )
c$$$         if(mod(i,10).eq.0)print*,i,' **** ',vect(1),vect(2),vect(3)
c$$$     $        ,' ---- ',f(1),f(2),f(3)
         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$$$         Z=VOUT(3)
c$$$         IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 200
c$$$         IF(Z.GT.Zout+TOLL) GOTO 20      
c$$$         IF(Z.LE.Zout-TOLL) THEN
c$$$            STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
c$$$            DO J=1,7
c$$$               VECT(J)=VECTINI(J)
c$$$            ENDDO
c$$$            BdLstep=0.
c$$$            GOTO 21
c$$$         ENDIF
 200     continue
         BdL=BdL+sngl(BdLstep*1.D-3) !Tesla * meters
c         print*,i,' #### ',BdL

         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-sngl(BdLstep*1.D-3) !Tesla * meters            
            BdLstep=0.
            GOTO 21
         ENDIF
      enddo
 2000 continue
c$$$      print*,' #### ',BdL
c$$$      print*,' #### ',xgood
c$$$      print*,' #### ',ygood
c$$$      print*,'========================================='
      return
      end


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