SUBROUTINE MAMRKUT(XSTART,DDXST,SIGP,H,IENTR,XEND,DDXEND,DELS,IER) C ===================================================================== C C Runge Kutta solver of the equation of motions. Mods to standard C routine. C C Implicit None REAL XSTART(3),DDXST(3),XEND(3),DDXEND(3) REAL SIGP, H, DELS INTEGER IER, IENTR C C INPUT PARAMETERS C XSTART(3) COORDINATES OF START POINT C DDXST (3) DX/DS,DY/DS,DZ/DS AT START POINT C SINP MOMENTUM (UNCHANGED ALONG TRACK) C GIVE NEGATIVE SIGN IF TRACK TURNS C CLOCKWISE WHEN LOOKING ALONG B C H STEP SIZE IN R C IENTR =1,2,3 OR 4 IF EQUATIONS ARE C TO BE INTEGRATED ALONG X,Y,Z OR R C RESPECTIVELY C C OUTPUT PARAMETERS C XEND (3) COORDINATES OF END POINT C DDXEND(3) DX/DS,DY/DS,DZ/DS AT END POINT C (THESE VARIABLES MAY OVERLAP WITH START POINT VAR.) C DELS CURVED ARC LENGTH DIFFERENCE C ALONG TRACK FOR THIS STEP C IER AN ERROR FLAG SIGNALLING C THE (UN-) SUITABLILITY OF THE C INTEGRATION VARIABLE. C =1 IF DDXST(IENTR) IS LESS THAN C TESTX (WHEN IENTR=1,2 OR 3) OR C DS/DR GREATER THAN TEST C (WHEN IENTR=4) IN ABSOLUTE VALUE, C =2 IF DDXEND(IENTR) IS LESS THAN C TESTX (WHEN IENTR=1,2 OR 3) OR C DS/DR GREATER THAN TEST C (WHEN IENTR=4) IN ABSOLUTE VALUE, C =0 IF O.K. C C USER ROUTINE CALLED C CALL MAMAGF (FIELD(1...3),X,Y,Z) C gb: 7/03/95 C CALL MAMAGF(FIELD,X,Y,Z, ier* ) C C ASSUMPTION C SPEED OF LIGHT IS .029.., WANTING GEV/C, M AND KGAUSS C C WARNING C INTEGRATION ALONG R (IENTR=4) IS NOT RECOMMENDED FOR C SMALL R,CORRECT PROTECTION CAN NOT BE GUARANTED. C ALSO,THERE IS A MORE ROBUST ALGORITHM AVAILABLE (RKRAD), C WHICH IS NOT USING THE NYSTROM INTEGRATION AND IS FASTER, C UNLESS THE FIELD ROUTINE IS THE LIMITING FACTOR C C LOCAL VARIABLES C C-- DECLARATIONS ADDED A.DAKE 05.06.92 C REAL SECYX(4),SECZX(4),F(3),FF(3),XPT(3) REAL SECXR(4),SECYR(4),SECZR(4) REAL SECSX(4),SECSR(4) REAL XT, YT, ZT EQUIVALENCE (XPT(1),XT),(XPT(2),YT),(XPT(3),ZT) REAL PINV, H1, H2, H4, H6, H8, X, Y, Z, A, B, C REAL PH, E, D, DSDR, ROOT, AT, BT, DXDS, R, CT, S REAL PHI, SGPHI, PHI1, AA, BB, CC, DD, PSI, RT INTEGER IX, IY, IZ, JX, JY, JZ REAL TESTX, TESTR, THIRD real testp C CGB> REAL CS,DS,FMOD CGB< DATA THIRD/.333333333333333/ DATA TESTX/0.1/ DATA TESTR/10./ data testp /0.1/ C C INITIALISATION. C IER = 0 CGB> IF(SIGP.EQ.0) GOTO 1400 if(abs(sigp).lt.testp) goto 200 CGB< C PINV = 0.02997925/SIGP GO TO (10,20,30,50),IENTR C C INTEGRATION WITH RESPECT TO X,Y OR Z C 10 IX = 1 IY = 2 IZ = 3 JX = 1 JY = 2 JZ = 3 GO TO 40 20 IX = 2 IY = 3 IZ = 1 JX = 3 JY = 1 JZ = 2 GO TO 40 30 IX = 3 IY = 1 IZ = 2 JX = 2 JY = 3 JZ = 1 40 CONTINUE IF (ABS(DDXST(IX)).LE.TESTX) GO TO 160 C C PRESERVE START POINT VALUES C X = XSTART(IX) Y = XSTART(IY) Z = XSTART(IZ) A = DDXST(IY)/DDXST(IX) B = DDXST(IZ)/DDXST(IX) H2 = 0.5*H H4 = 0.25*H CALL MAMAGF(F,XSTART(1),XSTART(2),XSTART(3),ier) if(ier.gt.0) goto 210 CGB> FMOD=F(1)*F(1)+F(2)*F(2)+F(3)*F(3) IF(FMOD.EQ.0) GOTO 1400 CGB< PH = PINV*H2 IF (DDXST(IX).LT.0.) PH = -PH C C FIRST INTERMEDIATE POINT C C = A*A+1. E = B*B D = E+1. DSDR = SQRT(E+C) ROOT = DSDR*PH SECYX(1) = ((F(IY)*A+F(IX))*B-F(IZ)*C)*ROOT SECZX(1) = ((-F(IZ)*B-F(IX))*A+F(IY)*D)*ROOT SECSX(1) = (E+C)*(B*F(IY)-A*F(IZ))*PH c print *,' first int' c print *, ' XSTART ',XSTART c print *, ' DDXST ',DDXST c print *,' a,b ',a,b c print *,' sec yx/zx,sx ',secyx(1),seczx(1),secsx(1) c print *, 'P ',sigp XT = X+H2 YT = H2*A+H4*SECYX(1)+Y ZT = H2*B+H4*SECZX(1)+Z AT = A+SECYX(1) C = AT*AT+1.0 BT = B+SECZX(1) E = BT*BT D = E+1.0 ROOT = SQRT(E+C)*PH C C SECOND INTERMEDIATE POINT C CALL MAMAGF(F,XPT(JX),XPT(JY),XPT(JZ),ier) if(ier.gt.0) goto 210 SECYX(2) = ((F(IY)*AT+F(IX))*BT-F(IZ)*C)*ROOT SECZX(2) = ((-F(IZ)*BT-F(IX))*AT+F(IY)*D)*ROOT SECSX(2) = (E+C)*(BT*F(IY)-AT*F(IZ))*PH c print *,' second int' c print *,' at,bt ',at,bt c print *,' sec yx/zx,sx ',secyx(2),seczx(2),secsx(2) c print *, 'P ',sigp C C TENTATIVE STEP SECOND HALF C AT = A+SECYX(2) C = AT*AT+1.0 BT = B+SECZX(2) E = BT*BT D = E+1.0 ROOT = SQRT(E+C)*PH SECYX(3) = ((F(IY)*AT+F(IX))*BT-F(IZ)*C)*ROOT SECZX(3) = ((-F(IZ)*BT-F(IX))*AT+F(IY)*D)*ROOT SECSX(3) = (E+C)*(BT*F(IY)-AT*F(IZ))*PH c print *,' tentative 2nd half' c print *,' at,bt ',at,bt c print *,' sec yx/zx,sx ',secyx(3),seczx(3),secsx(3) c print *, 'p ',sigp C C REACH FINAL X-VALUE OF THIS STEP C XT = X+H YT = (SECYX(3)+A)*H+Y ZT = (SECZX(3)+B)*H+Z AT = 2.0*SECYX(3)+A C = AT*AT+1.0 BT = B+2.0*SECZX(3) E = BT*BT D = E+1.0 ROOT = SQRT(E+C)*PH CALL MAMAGF(F,XPT(JX),XPT(JY),XPT(JZ),ier) if(ier.gt.0) goto 210 SECYX(4) = ((F(IY)*AT+F(IX))*BT-F(IZ)*C)*ROOT SECZX(4) = ((-F(IZ)*BT-F(IX))*AT+F(IY)*D)*ROOT C C FINAL VARIABLES THIS STEP C XEND(IY) = ((SECYX(1)+SECYX(2)+SECYX(3))*THIRD+A)*H+Y XEND(IZ) = ((SECZX(1)+SECZX(2)+SECZX(3))*THIRD+B)*H+Z DELS = ((SECSX(1)+SECSX(2)+SECSX(3))*THIRD+DSDR)*H IF (DDXST(IX).LT.0.) DELS=-DELS A = ((SECYX(2)+SECYX(3))*2.0+SECYX(1)+SECYX(4))*THIRD+A B = ((SECZX(2)+SECZX(3))*2.0+SECZX(1)+SECZX(4))*THIRD+B XEND(IX) = XT DXDS = SIGN(1./SQRT(1.+A*A+B*B),DDXST(IX)) IF (ABS(DXDS).LE.TESTX) GO TO 170 DDXEND(IX) = DXDS DDXEND(IY) = DXDS*A DDXEND(IZ) = DXDS*B RETURN C C INTEGRATION WITH RESPECT TO R C 50 CONTINUE C C PRESERVE START POINT VALUES C X = XSTART(1) Y = XSTART(2) Z = XSTART(3) A = DDXST(1) B = DDXST(2) C = DDXST(3) R = SQRT(Y*Y+Z*Z) IF (R.EQ.0.) GO TO 60 PHI=Y*B+Z*C IF (PHI.EQ.0.) GOTO 160 PHI=R/PHI SGPHI = SIGN(1.,PHI) H1 = SGPHI*H H2 = 0.5*H1 H6 = H1/6. H8 = 0.125*H1*H1 PHI = ABS(PHI) GO TO 70 60 PHI=B*B+C*C IF (PHI.EQ.0.) GOTO 160 PHI=1./PHI 70 IF (PHI.GT.TESTR) GO TO 160 PHI1 = PHI A = A*PHI B = B*PHI C = C*PHI CALL MAMAGF(F,X,Y,Z,ier) if(ier.gt.0) goto 210 CGB> FMOD=F(1)*F(1)+F(2)*F(2)+F(3)*F(3) if(FMOD.EQ.0) GOTO 1400 CGB< C C FIRST INTERMEDIATE POINT C AA = PINV*(B*F(3)-C*F(2)) BB = PINV*(C*F(1)-A*F(3)) CC = PINV*(A*F(2)-B*F(1)) DD = B*B+C*C PHI = SQRT(A*A+DD) IF (R.EQ.0.) GO TO 80 PSI = (1.-DD-PHI*(Y*BB+Z*CC))/R GO TO 90 80 PSI = -PHI*(B*BB+C*CC) 90 PSI = SGPHI*PSI SECXR(1) = PSI*A+PHI*AA SECYR(1) = PSI*B+PHI*BB SECZR(1) = PSI*C+PHI*CC SECSR(1) = PHI*PSI RT = R+H2 XT = X+H2*A+H8*SECXR(1) YT = Y+H2*B+H8*SECYR(1) ZT = Z+H2*C+H8*SECZR(1) AT = A+H2*SECXR(1) BT = B+H2*SECYR(1) CT = C+H2*SECZR(1) C C SECOND INTERMEDIATE POINT C CALL MAMAGF(F,XT,YT,ZT,ier) if(ier.gt.0) goto 210 AA = PINV*(BT*F(3)-CT*F(2)) BB = PINV*(CT*F(1)-AT*F(3)) CC = PINV*(AT*F(2)-BT*F(1)) DD = BT*BT+CT*CT PHI = SQRT(AT*AT+DD) IF (RT.EQ.0.) GO TO 100 PSI = (1.-DD-PHI*(YT*BB+ZT*CC))/RT GO TO 110 100 PSI = -PHI*(BT*BB+CT*CC) 110 PSI = SGPHI*PSI SECXR(2) = PSI*AT+PHI*AA SECYR(2) = PSI*BT+PHI*BB SECZR(2) = PSI*CT+PHI*CC SECSR(2) = PHI*PSI C C TENTATIVE STEP SECOND HALF C AT = A+H2*SECXR(2) BT = B+H2*SECYR(2) CT = C+H2*SECZR(2) AA = PINV*(BT*F(3)-CT*F(2)) BB = PINV*(CT*F(1)-AT*F(3)) CC = PINV*(AT*F(2)-BT*F(1)) DD = BT*BT+CT*CT PHI = SQRT(AT*AT+DD) IF (RT.EQ.0.) GO TO 120 PSI = (1.-DD-PHI*(YT*BB+ZT*CC))/RT GO TO 130 120 PSI = -PHI*(BT*BB+CT*CC) 130 PSI = SGPHI*PSI SECXR(3) = PSI*AT+PHI*AA SECYR(3) = PSI*BT+PHI*BB SECZR(3) = PSI*CT+PHI*CC SECSR(3) = PHI*PSI C C REACH FINAL R-VALUE OF THIS STEP C RT = R+H1 XT = X+H1*(A+H1*SECXR(3)) YT = Y+H1*(B+H1*SECYR(3)) ZT = Z+H1*(C+H1*SECZR(3)) AT = A+H1*SECXR(3) BT = B+H1*SECYR(3) CT = C+H1*SECZR(3) CALL MAMAGF(F,XT,YT,ZT,ier) if(ier.gt.0) goto 210 AA = PINV*(BT*F(3)-CT*F(2)) BB = PINV*(CT*F(1)-AT*F(3)) CC = PINV*(AT*F(2)-BT*F(1)) DD = BT*BT+CT*CT PHI = SQRT(AT*AT+DD) IF (PHI.GT.TESTR) GO TO 170 IF (RT.EQ.0.) GO TO 140 PSI = (1.-DD-PHI*(YT*BB+ZT*CC))/RT GO TO 150 140 PSI = -PHI*(BT*BB+CT*CC) 150 PSI = SGPHI*PSI SECXR(4) = PSI*AT+PHI*AA SECYR(4) = PSI*BT+PHI*BB SECZR(4) = PSI*CT+PHI*CC C C FINAL VARIABLES THIS STEP C R = RT X = X+H1*(A+H6*(SECXR(1)+SECXR(2)+SECXR(3))) Y = Y+H1*(B+H6*(SECYR(1)+SECYR(2)+SECYR(3))) Z = Z+H1*(C+H6*(SECZR(1)+SECZR(2)+SECZR(3))) DELS = H1*(PHI1+H6*(SECSR(1)+SECSR(2)+SECSR(3))) A = A+H6*(SECXR(1)+SECXR(4)+2.*(SECXR(2)+SECXR(3))) B = B+H6*(SECYR(1)+SECYR(4)+2.*(SECYR(2)+SECYR(3))) C = C+H6*(SECZR(1)+SECZR(4)+2.*(SECZR(2)+SECZR(3))) XEND(1) = X XEND(2) = Y XEND(3) = Z S = SQRT(A*A+B*B+C*C) DDXEND(1) = A/S DDXEND(2) = B/S DDXEND(3) = C/S RETURN C C CGB> 1400 goto(1401,1402,1403,1404) ientr 1401 if(abs(ddxst(1)).lt.testx) goto 180 XEND(1)=XSTART(1)+H H1=DDXST(2)/DDXST(1) H2=DDXST(3)/DDXST(1) XEND(2)=XSTART(2)+H*H1 XEND(3)=XSTART(3)+H*H2 DDXEND(1)=DDXST(1) DDXEND(2)=DDXST(2) DDXEND(3)=DDXST(3) DELS=H*SQRT(1.+H1*H1+H2*H2) RETURN 1402 if(abs(ddxst(2)).lt.testx) goto 180 XEND(2)=XSTART(2)+H H1=DDXST(1)/DDXST(2) H2=DDXST(3)/DDXST(2) XEND(1)=XSTART(1)+H*H1 XEND(2)=XSTART(3)+H*H2 DDXEND(1)=DDXST(1) DDXEND(2)=DDXST(2) DDXEND(3)=DDXST(3) DELS=H*SQRT(1.+H1*H1+H2*H2) RETURN 1403 if(abs(ddxst(3)).lt.testx) goto 180 XEND(3)=XSTART(3)+H H1=DDXST(1)/DDXST(3) H2=DDXST(2)/DDXST(3) XEND(1)=XSTART(1)+H*H1 XEND(2)=XSTART(2)+H*H2 DDXEND(1)=DDXST(1) DDXEND(2)=DDXST(2) DDXEND(3)=DDXST(3) DELS=H*SQRT(1.+H1*H1+H2*H2) RETURN 1404 R = XSTART(2)*XSTART(2)+XSTART(3)*XSTART(3) D = DDXST(2)*DDXST(2)+DDXST(3)*DDXST(3) C = XSTART(2)*DDXST(2)+XSTART(3)*DDXST(3) CS = C*C+D*(H*H+2.*H*SQRT(R)) IF(CS.LT.0.) GOTO 190 R = SQRT(CS) CS=(R-C)/D DS=(-R-C)/D IF(ABS(DS).LT.ABS(CS)) CS=DS XEND(1)=XSTART(1)+DDXST(1)*CS XEND(2)=XSTART(2)+DDXST(2)*CS XEND(3)=XSTART(3)+DDXST(3)*CS DDXEND(1)=DDXST(1) DDXEND(2)=DDXST(2) DDXEND(3)=DDXST(3) DELS=CS*SQRT(DDXST(1)**2+DDXST(2)**2+DDXST(3)**2) RETURN CGB< C C C ERROR EXITS, NO USABLE RESULT C 160 IER = 1 write(6,*) 'RKUT. IX,DDXST ',IX,DDXST(IX) RETURN 170 IER = 2 write(6,*) 'RKUT. DXDS ',DXDS RETURN CGB> 180 ier = 3 write(6,*) 'RKUT. LINEAR IER=3' return 190 ier =4 return 200 ier=5 return 210 ier=6 CGB< END C