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