C RELKIN. RELATIVISTIC KINEMATICS ON THE GEC4070. FORMULA C FROM MINESOTA PROGRAM. SYMBOLIC ENTRY OF REACTION WITH C MASS AND QVALUE RETURNED FROM COPY OF THE NUCLEAR DATA C TABLES (1983) TABULATION OBTAINED FROM MANCHESTER. C Sept. 84 Code converted to Fortran-77. C Input/output completely rewritten. C Subroutine REACT rewritten. C v/c, TOF and sig(LAB)/sig(CM) info. added. C J.Lukasiak, Univ. of Manchester C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL FIRST,NEWREA,NEWELB CHARACTER CHAR CHARACTER *24 LINEIN, TODAY DIMENSION RM(4),EM(4),IA(4),IZ(4) REAL ANGLE(3),EXCEN(15) EQUIVALENCE (A1,RM(1)),(A2,RM(2)),(A3,RM(3)),(A4,RM(4)) PARAMETER ( RADIAN = 0.01745329D0 , SPEED = 2.99792458D-1, $ PI = 3.1415927D0 ) DATA FIRST /.TRUE./ C C ENTRY OF SYMBOLIC REACTION 1010 CALL NOCRLF(4) WRITE(4,*) $ 'Reaction e.g. 28Si(C12,16O) : ' READ(3,4) LINEIN 4 FORMAT( A ) CALL REACT(RM,QGS,QERROR,EM,IA,IZ,IFLAG,LINEIN) IF (IFLAG.LE.1) GOTO 1019 C C Pathological cases C IF (IFLAG .EQ. 2) THEN WRITE (4,*) 'Symbol not an element name or string faulty.' , $ ' Retype reaction.' GO TO 1010 ELSE IF (IFLAG .EQ. 3) THEN DO 15 I=1,4 IF (RM(I) .EQ. 0.0) THEN WRITE (4,800) IZ(I),IA(I) 800 FORMAT('No mass found for Z,A = ', 2I4) 16 WRITE (4,810) 810 FORMAT ('Type / to retype reaction,' / $ 'type M to specify mass,' / $ 'type 0 to quit' ) READ (3,4) CHAR IF (CHAR .EQ. '/') GO TO 1010 IF (CHAR .EQ. '0') STOP IF (CHAR .EQ. 'M' ) THEN RM(I) = ACCPTR('Mass [amu] = ') QGS = ACCPTR('Q(gs) [MeV] = ') GO TO 15 END IF C Ask again for option GO TO 16 END IF 15 CONTINUE ELSE WRITE (4,*) 'IZ(4) or IA(4) negative. Retype reaction.' GO TO 1010 END IF 1019 CONTINUE NEWREA = .TRUE. IF (.NOT. FIRST) GO TO 1000 1020 CALL NOCRLF(4) WRITE (4,*) 'Level energies [MeV] : ' CALL GETRL(EXCEN, NEXCEN) IF (.NOT. FIRST) GO TO 1000 1030 ELAB = ACCPTR('Beam energy [MeV] = ') NEWELB = .TRUE. IF (.NOT. FIRST) GO TO 1000 1040 CALL NOCRLF(4) WRITE (4,*) 'Lab angles: START [,STOP [,STEP]] = ' CALL GETRL(ANGLE,NANGLE) ANG1 = ANGLE(1) IF (NANGLE .GE. 2) THEN ANG2 = ANGLE(2) ELSE ANG2 = ANG1 END IF IF (NANGLE .EQ. 3) THEN ANG3 = ANGLE(3) ELSE ANG3 = 1. END IF FIRST = .FALSE. C C End of input, ask for OPTION C GO TO 1000 C C ------------------------------------- C Do C 1050 CONTINUE ECM = ELAB*A2/(A1+A2) ETHRGS = -QGS*(A1+A2)/A2 IF (NEWREA .OR. NEWELB) THEN CALL INDATE (TODAY,24) WRITE (6,2000) LINEIN,TODAY, A2,A1,A3,A4,QGS,QERROR, $ ETHRGS,ELAB,ECM 2000 FORMAT (1H1 / 1X , A , 60X, A // $ 1X, 'masses :' , 4F11.6, 5X, $ 'Qgs =', F7.3, ' +/-' , F6.3, ' MeV' , 5X, $ 'threshold energy (g.s) = ', F7.3 // $ ' Elab =', F6.2, ' MeV Ecm =', F7.3, ' MeV' //) IF(IFLAG.EQ.2) THEN WRITE(4,65) WRITE(6,65) 65 FORMAT $ (' **** WARNING: you are using interpolated masses'/) END IF WRITE(6,2070) 2070 FORMAT (/ 1 68X, ' ASSOCIATED PARTICLE'/ $ ' LAB level Q ', $ ' CM E3 LAB/CM beta3 TOF3' , $ 5X , 'LAB', 6X, 'LAB LAB/CM beta4 TOF4' / $ ' angle' , 20X, 'angle' ,5X, 'LAB ' , 2 5X , 'ratio %c ns/m', 4X, 'angle', 3X, 'energy', 3 ' ratio %c ns/m' /) NEWREA = .FALSE. NEWELB = .FALSE. END IF WRITE(4,70) 70 FORMAT (/ $' LAB ANGLE LEVEL CM ANGLE LAB ENERGY LAB/CM RATIO' 1 ,' ASSOCIATED PARTICLE'/ 2 58X, 'LAB ANGLE LAB ENERGY'/) C C CALCULATION 100 ANGL=ANG1 C C LOOP THROUGH ANGLES 110 ANGLR=ANGL*RADIAN C Loop through levels DO 200 NE=1,NEXCEN Q=QGS-EXCEN(NE) IF (ECM+Q .LE. 0.0) THEN ETHRES=-Q*(A1+A2)/A2 WRITE(4,90) ANGL,ETHRES 90 FORMAT( 3X, F6.2, ' REACTION BELOW THRESHOLD (',F6.3,' MeV)'/) WRITE(6,2090) ANGL,EXCEN(NE),ETHRES 2090 FORMAT( 3X, F6.2, 4X, F6.3, 5X, 10('*'), 5X, $ ' REACTION BELOW THRESHOLD (',F6.3,' MeV)' ) GO TO 200 END IF CALL LABTCM(A1,A2,A3,A4,ELAB,Q,ANGLR,ANGCM1,T3L1,G1, 1ANGCM2,T3L2,G2,ANG4L1,ANG4L2,T4L1,T4L2,NS) IF(ANGCM1.EQ.0.AND.ANGL.NE.0.) THEN WRITE (4,160) ANGL,EXCEN(NE) WRITE (6,2160) ANGL,EXCEN(NE) 160 FORMAT (3X, F6.2, F7.3, 5X, 10('*'), $ ' Maximum scattering angle exceeded') 2160 FORMAT (3X, F6.2, F10.3, 5X, 10('*'), $ ' Maximum scattering angle exceeded') GO TO 200 END IF WRITE(4,120) $ ANGL,EXCEN(NE),ANGCM1/RADIAN,T3L1,G1,ANG4L1/RADIAN,T4L1 120 FORMAT(3X,F6.2,3X,F6.3,4X,F6.2,2X,F10.3,5X,F6.3,6X,F6.2,5X,F6.2) C beta, TOF & G4 (solid anlgle ratio) calculation put in by JL G4 = SIN(ANG4L1)**2 / SIN(PI-ANGCM1)**2 * COS(PI-ANGCM1-ANG4L1) G4 = ABS(G4) BETA3 = BETA(T3L1,A3) BETA4 = BETA(T4L1,A4) TOF3 = 1./(BETA3 * SPEED) TOF4 = 1./(BETA4 * SPEED) WRITE(6,2120) ANGL,EXCEN(NE),Q, $ ANGCM1/RADIAN,T3L1,G1,100*BETA3,TOF3,ANG4L1/RADIAN,T4L1, $ G4,100*BETA4,TOF4 2120 FORMAT $(F7.2,F8.3,F9.3,2X,F6.2,2X,F8.3,2X,F8.3,2F7.1,F9.2,F9.3, $ 2X, F6.3, 1X, 2F6.1) IF(NS.EQ.2) THEN G4 = SIN(ANG4L2)**2 / SIN(PI-ANGCM2)**2 * COS(PI-ANGCM2-ANG4L2) G4 = ABS(G4) BETA3 = BETA(T3L2,A3) BETA4 = BETA(T4L2,A4) TOF3 = 1./(BETA3 * SPEED) TOF4 = 1./(BETA4 * SPEED) WRITE(4,130) ANGCM2/RADIAN,T3L2,G2,ANG4L2/RADIAN,T4L2 WRITE(6,2130) $ ANGCM2/RADIAN,T3L2,G2,100*BETA3,TOF3,ANG4L2/RADIAN,T4L2, $ G4,100*BETA4,TOF4 130 FORMAT(22X,F6.2,4X,F8.3,2X,F9.3,6X,F6.2,5X,F6.2) 2130 FORMAT $( 24X, 2X,F6.2,2X,F8.3,2X,F8.3,2F7.1,F9.2,F9.3, $ 2X, F6.3, 1X, 2F6.1) END IF 200 CONTINUE ANGL=ANGL+ANG3 IF(ANGL.LE.ANG2) GOTO 110 1000 CONTINUE WRITE (4,2010) 2010 FORMAT (//'Options available: R(eaction)' / $ 19X, 'A(ngle)' / $ 19X, 'L(evels)' / $ 19X, 'B(eam energy)' / $ 19X, 'S(how)' / $ 19X, 'D(o)' / $ 19X, 'X(it)' ) CALL NOCRLF(4) WRITE (4,*) 'Option : ' READ (3,4) CHAR IF (CHAR .EQ. 'X') S T O P IF (CHAR .EQ. 'R') GO TO 1010 IF (CHAR .EQ. 'A') GO TO 1040 IF (CHAR .EQ. 'L') GO TO 1020 IF (CHAR .EQ. 'B') GO TO 1030 IF (CHAR .EQ. 'D') GO TO 1050 IF (CHAR .EQ. 'S') THEN WRITE (4,2020) LINEIN,QGS,-QGS*(A1+A2)/A2, $ ELAB,ELAB*A2/(A1+A2),ANG1,ANG2,ANG3, (EXCEN(I), I=1,NEXCEN) 2020 FORMAT $ (/A,2X,'Qgs =', F7.3, ' threshold (g.s.) = ', F7.3, ' MeV' / $ 'Elab = ', F7.3, ' MeV' ,' E(cm) = ', F7.3, ' MeV' / $ 'Lab. angle: from' , F7.2, ' to', F7.2, $ ' in steps of', F6.2, ' degrees' / $ 'Levels:', 5F7.3, (/7X, 5F7.3)) GO TO 1000 END IF GO TO 1000 END C C SUBROUTINE LABTCM TO RETURN CM ANGLE, ENERGY AND JACOBIAN C FOR EJECTILE IN GIVEN REACTION. NS=2 IF DOUBLE SOLUTION. C MASSES ARE IN AMU, ENERGIES IN MEV C THIS ROUTINE FROM MINNESOTA AND PROVIDED BY J.S.LILLEY SUBROUTINE LABTCM(AA1,AA2,AA3,AA4,TIL,Q,ANGL, 1ANGCM1,T3L1,G1,ANGCM2,T3L2,G2,ANG4L1,ANG4L2,T4L1,T4L2,NS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA AMU,PI/931.5016D0,3.1415926D0/ A1=AA1*AMU A2=AA2*AMU A3=AA3*AMU A4=AA4*AMU Y=A1+A2 BETAC=SQRT(TIL*(TIL+2.*A1))/(Y+TIL) GAMCI2=1.-BETAC*BETAC ECMI=SQRT(Y*Y+2.*TIL*A2) ECMF=ECMI+Q-Y+A3+A4 E3CM=(ECMF*ECMF+A3*A3-A4*A4)/(2.*ECMF) GAM3C2=(E3CM/A3)*(E3CM/A3) BETA3C=SQRT(1.-1./GAM3C2) Y=GAM3C2*GAMCI2 COSAGL=COS(ANGL) B=-BETAC*COSAGL A=Y+B*B C=1.-Y NS=0 ANGCM1=0. T3L1=0. G1=0. ANGCM2=0. T3L2=0. G2=0. ANG4L1=0. ANG4L2=0. D=B*B-A*C IF(D) 10,10,1 1 D=SQRT(D) DO 9 KJ=1,3,2 J=KJ-2 FJ=-J BETA3L=(-B+SIGN(D,FJ))/A IF(BETA3L-0.0000001) 10,10,2 2 GAMM3L=1./SQRT(1.-BETA3L*BETA3L) TL=(GAMM3L-1.)*A3 Y=BETA3L*COSAGL Y=(Y-BETAC)/((1.-BETAC*Y)*BETA3C) IF(ABS(Y)-1.) 5,5,3 3 IF(ABS(Y)-1.00001) 4,9,9 4 Y= SIGN(1D0,Y) 5 AC=DACOS(Y) G=BETA3C*(1.+B/BETA3L)/(1.+BETAC*BETA3C*Y)/BETA3L IF(NS) 7,6,7 6 ANGCM1=AC T3L1=TL G1=G GOTO 8 7 ANGCM2=AC T3L2=TL G2=-G 8 NS=NS+1 C ASSOCIATED PARTICLE ANGLE CALCULATION C THIS CODE PUT IN BY B.R.FULTON 30/7/81 T4=TIL+Q-TL P32C2=TL*(TL+2.*A3) P42C2=T4*(T4+2.*A4) ANGL4=ASIN(SQRT(P32C2/P42C2)*SIN(ANGL)) P12C2=TIL*(TIL+2.*A1) TEST=P32C2*COSAGL*COSAGL IF(TEST.GT.P12C2.AND.ANGL.LT.PI/2) ANGL4=PI-ANGL4 IF(NS.EQ.1.) ANG4L1=ANGL4 IF(NS.EQ.2.) ANG4L2=ANGL4 IF(NS.EQ.1.) T4L1=T4 IF(NS.EQ.2.) T4L2=T4 9 CONTINUE 10 CONTINUE RETURN END FUNCTION BETA(EKIN,A) C Calculates relativistic v/c for a nucleus with atomic C mass A and kinetic energy EKIN. C beta**2 = p2/(p2+m**2) C where p2 is momentum**2 and m is the rest mass IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA AMU /931.5016D0/ AMEV = A*AMU P2 = EKIN * (EKIN + 2*AMEV) BETA = SQRT( P2/(P2+AMEV**2) ) RETURN END C SUBROUTINE REACT TO INTERPRET A REACTION IN 28SI(12C,16O) C SYMBOLIC FORMAT AND RETURN REAL*8 VALUES OF THE MASSES, C QVALUE AND ERROR AS WELL AS THE RESIDUAL PARTICLE SYMBOL, C AND A FEW OTHER RELEVENT DETAILS. C INPUT:LINE REACTION (CHARACTER*24) C OUTPUT:RM(4) MASSES(AMU)................DOUBLE PRECISION C Q REACTION Q-VALUE (MEV)........DOUBLE PRECISION C EQ ERROR ON Q-VALUE (MEV)........DOUBLE PRECISION C EM(4) ERROR IN MASSES (KEV)......DOUBLE PRECISION C IA(4) MASSES C IZ(4) ATOMIC NUMBER C IFLAG = 0 O.K. ,all masses are experimental ones C 1 O.K., one or more masses are interpolated C 2 ONE OR MORE REQUIRED SYMBOLS NOT IN LIST C or string improperly constructed C 3 no isotope for a certain Z C 4 IZ(4) or IA(4) negative C C SUBROUTINE REACT(RM,Q,EQ,EM,IA,IZ,IFLAG,LINE) CHARACTER *(*) LINE CHARACTER*2 NUCNAM(110),SYM(3) CHARACTER*4 FORMAT(3) INTEGER*2 IAMIN,IAMAX,MSTART,NAME DOUBLE PRECISION RM(4),EM(4),Q,EQ,AM(4),X,AMU DIMENSION IA(4),IZ(4) C PARAMETER (MAXIZ = 109) PARAMETER (AMU = 931501.6D0) C DATA NUCNAM/' N','H ','HE', +'LI','BE','B ','C ','N ','O ','F ','NE', +'NA','MG','AL','SI','P ','S ','CL','AR', +'K ','CA','SC','TI','V ','CR','MN','FE','CO', +'NI','CU','ZN','GA','GE','AS','SE','BR','KR', +'RB','SR','Y ','ZR','NB','MO','TC','RU','RH', +'PD','AG','CD','IN','SN','SB','TE','I ','XE', +'CS','BA','LA','CE','PR','ND','PM','SM', +'EU','GD','TB','DY','HO','ER','TM','YB', +'LU','HF','TA','W ','RE','OS','IR','PT', +'AU','HG','TL','PB','BI','PO','AT','RN', +'FR','RA','AC','TH','PA','U ','NP','PU', +'AM','CM','BK','CF','ES','FM','MD','NO', +'LR','RF','HA','NH','NS','UO','UE'/ DATA FORMAT /'(I1)' , '(I2)' , '(I3)'/ C IF (IFIRST.NE.0) GOTO 10 CALL CONECT('LOCAL.USERPROG.DATA.MASS83 ',5) IFIRST=1 10 IFLAG=0 Q = 0.0 EQ = 0.0 DO 20 I=1,4 RM(I)=0.0 EM(I)=0.0 IA(I)=0 IZ(I)=0 20 CONTINUE C C *READING ALONG REACTION TO FIND NUCLEI* C C Find positions of ( , ) I1 = INDEX(LINE, '(' ) I2 = INDEX(LINE, ',' ) I3 = INDEX(LINE, ')' ) C Zero value of Ii means improper string construction IF (I1 * I2 * I3 .EQ. 0) THEN IFLAG = 2 RETURN END IF C Get masses and symbols from the input string CALL SYMDEC (LINE,1,I1-1,IA(1),SYM(1)) CALL SYMDEC (LINE,I1+1,I2-1,IA(2),SYM(2)) CALL SYMDEC (LINE,I2+1,I3-1,IA(3),SYM(3)) DO 40 I=1,3 IF (IA(I) .EQ. 0) THEN C Is it one of the special cases : A,P,N,D,T ? IF (SYM(I) .EQ. 'A ') THEN IZ(I) = 2 IA(I) = 4 ELSE IF (SYM(I) .EQ. 'N ') THEN IZ(I) = 0 IA(I) = 1 ELSE IF (SYM(I) .EQ. 'P ') THEN IZ(I) = 1 IA(I) = 1 ELSE IF (SYM(I) .EQ. 'D ') THEN IZ(I) = 1 IA(I) = 2 ELSE IF (SYM(I) .EQ. 'T ') THEN IZ(I) = 1 IA(I) = 3 ELSE C None of the special cases and no mass given, C therefore it must be an error; quit IFLAG = 2 RETURN END IF C End of mass=0 case ELSE C Mass given, so get Z DO 50 J=1,MAXIZ 50 IF (NUCNAM(J) .EQ. SYM(I)) GO TO 55 IFLAG = 2 RETURN 55 IZ(I) = J-1 END IF 40 CONTINUE C C *CALCULATING A AND Z OF 4TH NUCLEI* C 130 IZ(4)=IZ(1)+IZ(2)-IZ(3) IA(4)=IA(1)+IA(2)-IA(3) IF (IZ(4).LT.0.OR.IA(4).LT.0) THEN IFLAG = 4 RETURN END IF C C *OBTAINING DATA FROM MASS FILE* C DO 170 J=1,4 KEYPTR=IZ(J)+1 READ(5,REC=KEYPTR)IAMIN,IAMAX,MSTART,NAME IF (IA(J).LT.IAMIN.OR.IA(J).GT.IAMAX) IFLAG = 3 KEYPTR=IA(J)-IAMIN+MSTART READ(5,REC=KEYPTR)RMASS,ERRM C C *CALCULATION DATA PUT INTO ARRAYS* C IF (ERRM.EQ.1000.0) IFLAG=1 RM(J)=DBLE(RMASS)/AMU+DBLE (IA(J)) AM(J)=RMASS EM(J)=ERRM 170 CONTINUE C C *CONVERTING Z OF 4TH NUCLEI INTO NUCLEI SYMBOL* C * AND CHANGING INTEGER A INTO CHARACTER A * C C get number of significant digits NSD = LOG10 (REAL(IA(4))+.01) + 1 WRITE (LINE(I3+1:I3+NSD), FORMAT(NSD)) IA(4) C symbol LINE(I3+NSD+1:I3+NSD+2) = NUCNAM(IZ(4)+1) C C *CALCULATING QVALUE AND ERROR* C Q=(AM(1)+AM(2)-AM(3)-AM(4))/1000.0D0 EQ=SQRT(EM(1)**2+EM(2)**2+EM(3)**2+EM(4)**2)/1000.0D0 C C *INTERCHANGING TARGET AND PROJECTILE POSITIONS IN ARRAYS* C 260 X=RM(1) RM(1)=RM(2) RM(2)=X X=EM(1) EM(1)=EM(2) EM(2)=X IX=IA(1) IA(1)=IA(2) IA(2)=IX IX=IZ(1) IZ(1)=IZ(2) IZ(2)=IX RETURN END SUBROUTINE SYMDEC(STRING,FROM,TO,MASS,SYMBOL) C Searches STRING(FROM:TO) for mass and chemical symbol, C ignoring spaces. C Returns 0 and/or ' ' if unsuccessfull. C IMPLICIT INTEGER (A-Z) CHARACTER SPACE , CHAR CHARACTER *(*) STRING CHARACTER*2 SYMBOL, SYM PARAMETER (SPACE = ' ') C MASS = 0 SYM = ' ' DO 100 I=FROM,TO CHAR = STRING(I:I) IF (CHAR .EQ. SPACE) GO TO 100 C IF (CHAR .GE. '0' .AND. CHAR. LE. '9') THEN C Digit found READ (CHAR, '(I1)') DIGIT MASS = MASS*10 + DIGIT ELSE IF (SYM(1:1) .EQ. SPACE) THEN SYM(1:1) = CHAR ELSE SYM(2:2) = CHAR END IF END IF 100 CONTINUE SYMBOL = SYM RETURN END