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