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

Diff of /DarthVader/TrackerLevel2/src/F77/grkuta.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Fri May 19 13:15:55 2006 UTC revision 1.6 by pam-fi, Tue Nov 27 11:43:51 2007 UTC
# Line 28  C.    *                                 Line 28  C.    *                                
28  C.    ******************************************************************  C.    ******************************************************************
29  C.  C.
30        IMPLICIT DOUBLE PRECISION(A-H,O-Z)        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
31          COMMON/DELTAB/DELTA0,DELTA1
32  *  *
33        REAL VVV(3),FFF(3)        REAL VVV(3),FFF(3)
34        REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4)        REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4)
# Line 41  C. Line 42  C.
42        PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)        PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
43        PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)        PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
44        PARAMETER (PISQUA=.986960440109D+01)        PARAMETER (PISQUA=.986960440109D+01)
45        PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)        PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
46    
47          REAL*8 DELTAB(3)
48    
49  *.  *.
50  *.    ------------------------------------------------------------------  *.    ------------------------------------------------------------------
51  *.  *.
52  *             This constant is for units CM,GEV/C and KGAUSS  *             This constant is for units CM,GEV/C and KGAUSS
53  *  *
54    
55        ITER = 0        ITER = 0
56        NCUT = 0        NCUT = 0
57        DO 10 J=1,7        DO 10 J=1,7
# Line 69  C. Line 73  C.
73        DO I=1,3        DO I=1,3
74         F(I)=DBLE(FFF(I))         F(I)=DBLE(FFF(I))
75        ENDDO        ENDDO
76          DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2))
77          F(2) = F(2)+DELTAB(2)
78    cPP   -----------------
79  *  *
80  *             Start of integration  *             Start of integration
81  *  *
# Line 106  C. Line 113  C.
113        CALL GUFLD(VVV,FFF)        CALL GUFLD(VVV,FFF)
114        DO I=1,3        DO I=1,3
115         F(I)=DBLE(FFF(I))         F(I)=DBLE(FFF(I))
116        ENDDO            ENDDO  
117          DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2))
118          F(2) = F(2)+DELTAB(2)
119    cPP   -----------------
120  C      CALL GUFLD(XYZT,F)  C      CALL GUFLD(XYZT,F)
121        AT     = A + SECXS(1)        AT     = A + SECXS(1)
122        BT     = B + SECYS(1)        BT     = B + SECYS(1)
# Line 141  C      CALL GUFLD(XYZT,F) Line 151  C      CALL GUFLD(XYZT,F)
151        DO I=1,3        DO I=1,3
152         F(I)=DBLE(FFF(I))         F(I)=DBLE(FFF(I))
153        ENDDO        ENDDO
154          DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2))
155          F(2) = F(2)+DELTAB(2)
156    cPP   -----------------
157  C      CALL GUFLD(XYZT,F)  C      CALL GUFLD(XYZT,F)
158  *  *
159        Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H        Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
# Line 237  C      CALL GUFLD(XYZT,F) Line 250  C      CALL GUFLD(XYZT,F)
250  *      *    
251  *      *    
252    
253    **********************************************************************
254    *
255    *
256    *     routine per tracciare la particella di uno STEP
257    *     *** extended version ***
258    *     it return also the track-length
259    *
260          SUBROUTINE GRKUTA2 (CHARGE,STEP,VECT,VOUT)
261    C.
262    C.    ******************************************************************
263    C.    *                                                                *
264    C.    *  Runge-Kutta method for tracking a particle through a magnetic *
265    C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
266    C.    *  Standards, procedure 25.5.20)                                 *
267    C.    *                                                                *
268    C.    *  Input parameters                                              *
269    C.    *       CHARGE    Particle charge                                *
270    C.    *       STEP      Step size                                      *
271    C.    *       VECT      Initial co-ords,direction cosines,momentum     *
272    C.    *  Output parameters                                             *
273    C.    *       VOUT      Output co-ords,direction cosines,momentum      *
274    C.    *  User routine called                                           *
275    C.    *       CALL GUFLD(X,F)                                          *
276    C.    *                                                                *
277    C.    *    ==>Called by : <USER>, GUSWIM                               *
278    C.    *       Authors    R.Brun, M.Hansroul  *********                 *
279    C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
280    C.    *                                                                *
281    C.    *                                                                *
282    C.    ******************************************************************
283    C.
284          IMPLICIT DOUBLE PRECISION(A-H,O-Z)
285    *
286          REAL VVV(3),FFF(3)
287          REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4)
288          REAL*8 XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
289          DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
290          EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
291         +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
292    *
293          PARAMETER (MAXIT = 1992, MAXCUT = 11)
294          PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
295          PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
296          PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
297          PARAMETER (PISQUA=.986960440109D+01)
298          PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
299    
300    *     track length
301          REAL*8 DL
302    
303    *.
304    *.    ------------------------------------------------------------------
305    *.
306    *             This constant is for units CM,GEV/C and KGAUSS
307    *
308          ITER = 0
309          NCUT = 0
310          DO 10 J=1,8
311             VOUT(J)=VECT(J)
312       10 CONTINUE
313          PINV   = EC * CHARGE / VECT(7)
314          TL = 0.
315          H      = STEP
316          
317    c      print*,'===================== START GRKUTA2'
318          
319    *
320    *
321       20 REST  = STEP-TL
322          IF (DABS(H).GT.DABS(REST)) H = REST
323          DO I=1,3
324           VVV(I)=SNGL(VOUT(I))
325          ENDDO
326          
327          CALL GUFLD(VVV,FFF)
328    *      print*,'GRKUTA Bx,By,Bz: ',(FFF(i),i=1,3)
329          DO I=1,3
330           F(I)=DBLE(FFF(I))
331          ENDDO
332    *
333    *             Start of integration
334    *
335          X      = VOUT(1)
336          Y      = VOUT(2)
337          Z      = VOUT(3)
338          A      = VOUT(4)
339          B      = VOUT(5)
340          C      = VOUT(6)
341    
342          DL     = VOUT(8)
343    
344    *
345          H2     = HALF * H
346          H4     = HALF * H2
347          PH     = PINV * H
348          PH2    = HALF * PH
349          SECXS(1) = (B * F(3) - C * F(2)) * PH2
350          SECYS(1) = (C * F(1) - A * F(3)) * PH2
351          SECZS(1) = (A * F(2) - B * F(1)) * PH2
352          ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
353          IF (ANG2.GT.PISQUA) GO TO 40
354          DXT    = H2 * A + H4 * SECXS(1)
355          DYT    = H2 * B + H4 * SECYS(1)
356          DZT    = H2 * C + H4 * SECZS(1)
357          XT     = X + DXT
358          YT     = Y + DYT
359          ZT     = Z + DZT
360    *
361    *              Second intermediate point
362    *
363          EST = DABS(DXT)+DABS(DYT)+DABS(DZT)
364          IF (EST.GT.H) GO TO 30
365    
366          DO I=1,3
367           VVV(I)=SNGL(XYZT(I))
368          ENDDO
369          CALL GUFLD(VVV,FFF)
370          DO I=1,3
371           F(I)=DBLE(FFF(I))
372          ENDDO    
373    C      CALL GUFLD(XYZT,F)
374          AT     = A + SECXS(1)
375          BT     = B + SECYS(1)
376          CT     = C + SECZS(1)
377    *
378          SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
379          SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
380          SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
381          AT     = A + SECXS(2)
382          BT     = B + SECYS(2)
383          CT     = C + SECZS(2)
384          SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
385          SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
386          SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
387          DXT    = H * (A + SECXS(3))
388          DYT    = H * (B + SECYS(3))
389          DZT    = H * (C + SECZS(3))
390          XT     = X + DXT
391          YT     = Y + DYT
392          ZT     = Z + DZT
393          AT     = A + TWO*SECXS(3)
394          BT     = B + TWO*SECYS(3)
395          CT     = C + TWO*SECZS(3)
396    *
397          EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
398          IF (EST.GT.2.*ABS(H)) GO TO 30
399    
400          DO I=1,3
401           VVV(I)=SNGL(XYZT(I))
402          ENDDO
403          CALL GUFLD(VVV,FFF)
404          DO I=1,3
405           F(I)=DBLE(FFF(I))
406          ENDDO
407    C      CALL GUFLD(XYZT,F)
408    *      
409          Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
410          Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
411          X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
412    *
413          SECXS(4) = (BT*F(3) - CT*F(2))* PH2
414          SECYS(4) = (CT*F(1) - AT*F(3))* PH2
415          SECZS(4) = (AT*F(2) - BT*F(1))* PH2
416          A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
417          B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
418          C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
419    *
420          EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
421         ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
422         ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
423    *
424          IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
425    
426          ITER = ITER + 1
427          NCUT = 0
428    *               If too many iterations, go to HELIX
429          IF (ITER.GT.MAXIT) GO TO 40
430    *
431          DL     = VOUT(8) +
432         $     DSQRT( 0
433         $     + (X-VOUT(1))**2
434         $     + (Y-VOUT(2))**2
435         $     + (Z-VOUT(3))**2
436         $     )
437    c      print*,'- ',VOUT(3),z,VOUT(1),x,VOUT(2),y,DL
438    *
439          TL = TL + H
440          IF (EST.LT.(DLT32)) THEN
441             H = H*TWO
442          ENDIF
443          CBA    = ONE/ SQRT(A*A + B*B + C*C)
444          VOUT(1) = X
445          VOUT(2) = Y
446          VOUT(3) = Z
447          VOUT(4) = CBA*A
448          VOUT(5) = CBA*B
449          VOUT(6) = CBA*C
450          VOUT(8) = DL
451          REST = STEP - TL
452          IF (STEP.LT.0.) REST = -REST
453          IF (REST .GT. 1.E-5*DABS(STEP)) GO TO 20
454    *
455          GO TO 999
456    *
457    **              CUT STEP
458       30 NCUT = NCUT + 1
459    *               If too many cuts , go to HELIX
460          IF (NCUT.GT.MAXCUT)       GO TO 40
461          H = H*HALF
462          GO TO 20
463    *
464    **              ANGLE TOO BIG, USE HELIX
465       40 F1  = F(1)
466          F2  = F(2)
467          F3  = F(3)
468          F4  = DSQRT(F1**2+F2**2+F3**2)
469          RHO = -F4*PINV
470          TET = RHO * STEP
471          IF(TET.NE.0.) THEN
472             HNORM = ONE/F4
473             F1 = F1*HNORM
474             F2 = F2*HNORM
475             F3 = F3*HNORM
476    *
477             HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
478             HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
479             HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
480    
481             HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
482    *
483             RHO1 = ONE/RHO
484             SINT = DSIN(TET)
485             COST = TWO*DSIN(HALF*TET)**2
486    *
487             G1 = SINT*RHO1
488             G2 = COST*RHO1
489             G3 = (TET-SINT) * HP*RHO1
490             G4 = -COST
491             G5 = SINT
492             G6 = COST * HP
493    
494             VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
495             VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
496             VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
497    
498             VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
499             VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
500             VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
501    *
502          ELSE
503             VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
504             VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
505             VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
506    *
507          ENDIF
508    *     TEMP !!! TEMP !!! TEMP !!! TEMP !!! TEMP !!! TEMP !!!
509    *     devo mettere la lunghezza dell'elica!!!!!!!!!!!!!!
510    *     ma non mi riesce :-(
511          VOUT(8) = DSQRT( 0
512         $     +(VOUT(IX)-VECT(IX))**2
513         $     +(VOUT(IY)-VECT(IY))**2
514         $     +(VOUT(IZ)-VECT(IZ))**2
515         $     )
516    c      print*,'WARNING: GRKUTA2 --> '
517    c     $     ,'helix :-( ... length evaluated with straight line'
518    
519    *
520      999 END
521    *    
522    *    
523    
524  **********************************************************************  **********************************************************************
525  *      *    
# Line 252  C      CALL GUFLD(XYZT,F) Line 534  C      CALL GUFLD(XYZT,F)
534        real*8 vv(3),ff(3)        !inter_B.f works in double precision        real*8 vv(3),ff(3)        !inter_B.f works in double precision
535    
536    
537          do i=1,3        do i=1,3
538            vv(i)=v(i)/100.       !inter_B.f works in meters           vv(i)=v(i)/100.        !inter_B.f works in meters
539          enddo        enddo
540  c       inter_B: coordinates in m, B field in Tesla  c     inter_B: coordinates in m, B field in Tesla
541          call inter_B(vv(1),vv(2),vv(3),ff)  c$$$      print*,'GUFLD: v ',v
542          do i=1,3                !change back the field in kGauss        call inter_B(vv(1),vv(2),vv(3),ff)
543            f(i)=ff(i)*10.        do i=1,3                  !change back the field in kGauss
544          enddo           f(i)=ff(i)*10.
545          enddo
546    c$$$      print*,'GUFLD: b ',f
547          
548        return        return
549        end        end
550    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23