C THREEKIN. THREE PARTICLE KINENATICS BASED ON FORMULAE OF C G.G. OHLSEN NIM 37(1965)240, PROGRAMMED ON ROCHESTER C PDP 6 BY RAINER OST, TRANSFERRED TO 4070 AND MODIFIED TO C RUN BY BRIAN FULTON.........IT'LL NEVER RUN!!! COMMON E1LMAX,E1LMIN,E2LMAX,E2LMIN,TH1MAX,TH2MAX,E0, FXMP,XMT,XM1,XM2,XM3,Q,EPL,DELE1L,TH1L,TH2L,PHI2L, FLMAX,LMIN,ICODE COMMON/COMON1/ E1L(250),E2L(250,2),E23(250),E31(250,2), FE12(250,2),TH1C(250),TH2C(250,2),TH3C(250,2),THCP1(250,2) COMMON/COMON2/THCP2(250,2),THCP3(250,2),RHO1(250,2), FRHO2(250,2),TH3L(250,2),PHI3L(250,2), FIM(5),IZ(5),EX(5),ER(5) C DIMENSION EEXX(10),EPLX(10),TH1X(10) C FIRST SET DEFAULT VALUES CHARACTER*20 FILEOU, OVERW*1 LOGICAL EXI RAD=3.141592654/180. EPL=0. E0=0. DELE1L=0. TH1L=0. TH2L=0. PHI2L=0. PHI1L=180. IPAP=0 C C INPUT /OUTPUT STREAMS C ISIP=5 ISOP=6 IDOP=6 c c Ask for output filename. Carriage return implies output to standard out. 51 WRITE(ISOP,*)'Please enter name of output file (20 chars MAX).' WRITE(ISOP,*)'Carriage return for output to screen.' READ(ISIP,'(A)') FILEOU IF(FILEOU.EQ.'') THEN GOTO 59 ELSE IDOP=14 INQUIRE(FILE=FILEOU,IOSTAT=IOS,ERR=53,EXIST=EXI) IF(EXI.EQ..TRUE.)THEN WRITE(ISOP,*)'File exists. Overwrite ?' READ(ISIP,'(A)') OVERW IF(OVERW.EQ.'Y'.OR.OVERW.EQ.'y')GOTO 52 GOTO 51 ENDIF 52 OPEN(IDOP,IOSTAT=IOS,ERR=54,FILE=FILEOU) GOTO 59 ENDIF 53 WRITE(ISOP,*)'Error in INQUIREing about file. IOSTAT = ',IOS GOTO 51 54 WRITE(ISOP,*)'Error in OPENing file. IOSTAT = ',IOS GOTO 51 C C NOW READ IN INPUT 59 WRITE(ISOP,200) 200 FORMAT(' (1) FOR WIDE PAPER'/) READ(ISIP,*) IPAP 6 CONTINUE WRITE(ISOP,202) 202 FORMAT(' PARTICLE MASSES'/' PROJECTILE,TARGET,1ST OUT,2ND OUT') READ(ISIP,*) (IM(I),I=1,4) WRITE(ISOP,204) 204 FORMAT(' PARTICLE CHARGE'/' PROJECTILE,TARGET,1ST OUT,2ND OUT') READ(ISIP,*) (IZ(I),I=1,4) IM(5)=IM(1)+IM(2)-(IM(3)+IM(4)) IZ(5)=IZ(1)+IZ(2)-(IZ(3)+IZ(4)) XMP=FLOAT(IM(1)) XMT=FLOAT(IM(2)) XM1=FLOAT(IM(3)) XM2=FLOAT(IM(4)) XM3=FLOAT(IM(5)) WRITE(ISOP,206) 206 FORMAT(' BEAM ENERGY,QVALUE,EXCITATION ENERGY (MEV)'/) READ(ISIP,*) EPL,QU,EEX Q=QU-EEX 207 WRITE(ISOP,208) 208 FORMAT(' LAB ANGLE 1, LAB ANGLE 2'/) READ(ISIP,*) TH1L,TH2L WRITE(ISOP,210) 210 FORMAT(' OUTPUT ENERGY START, OUTPUT ENERGY STEP (MEV)'/) READ(ISIP,*) E0,DELE1L IBRA=0 WRITE(ISOP,212) 212 FORMAT(' OUTPUT FOR UPPER(0), LOWER(1), BOTH(-1) BRANCHES'/) READ(ISIP,*) IBRA 9 CONTINUE C CALL KDRE C IF(ICODE.NE.1) GO TO 6 WRITE(IDOP,101)IM(1),IM(2),IM(3),IM(4),IM(5),EPL,QU,EEX WRITE(IDOP,104)IZ(1),IZ(2),IZ(3),IZ(4),IZ(5) WRITE(IDOP,107) WRITE(IDOP,105)TH1L,TH2L WRITE(IDOP,106) PHI1L,PHI2L WRITE(IDOP,107) WRITE(IDOP,102) TH1MAX,TH2MAX WRITE(IDOP,103) E1LMIN,E1LMAX,E2LMIN,E2LMAX IA=1 IB=1 IF(IBRA)3,2,1 3 IB=2 GO TO 2 1 IA=2 IB=2 2 CONTINUE C VTOT=SQRT(2.*XMP*EPL)/(XMP+XMT) A1=SQRT(0.5*XM1*VTOT*VTOT) A2=SQRT(0.5*XM2*VTOT*VTOT) A3=SQRT(0.5*XM3*VTOT*VTOT) CO1=COS(TH1L*RAD) CO2=COS(TH2L*RAD) ETOT=Q+(XMT*EPL)/(XMP+XMT) WRITE(IDOP,107) WRITE(IDOP,113)VTOT C DO 4 II=IA,IB WRITE(IDOP,107) WRITE(IDOP,111) II WRITE(IDOP,107) 221 FORMAT(' LAB ENERGY RELATIVE ENERGY',7X,'PARTICLE ANGLES' 1,4X,'DENSITY STATES'/15X,'(EXCITATION-BREAKUP)',10X,'(CM)', 112X,'E1L AXIS'/) 220 FORMAT(' LAB ENERGY RELATIVE ENERGY',7X, 1'PARTICLE ANGLES',4X,'ANGLE BETWEEN COORD',5X,'ANGLE BETWEEN ', 16X,'DENSITY OF STATES'/ 115X,'(EXCITATION-BREAKUP)',10X,'(CM)',13X,'SYSTEMS (CM)',8X, 1'PARTICLES (CM)',5X,'E1L AXIS E AXIS'/58X, 1'(SEQUENTIAL DECAY)'/) IF(IPAP.EQ.0) GO TO 11 WRITE(IDOP,220) WRITE(IDOP,112) GO TO 12 11 WRITE(IDOP,221) WRITE(IDOP,118) 12 DO 5 L=LMIN,LMAX IF(E2L(L,II).LE.0.) GO TO 5 IF(PHI3L(L,II).LT.90.) TH3C(L,II)=-TH3C(L,II) IF(PHI2L.LT.90.) TH2C(L,II)=-TH2C(L,II) T1=180.+TH1C(L) T2=180.+TH2C(L,II) T3=180.+TH3C(L,II) DIF1=ABS(T1-T2) DIF2=ABS(T1-T3) DIF3=ABS(T2-T3) IF(DIF1.GT.180.) DIF1=360.-DIF1 IF(DIF2.GT.180.) DIF2=360.-DIF2 IF(DIF3.GT.180.) DIF3=360.-DIF3 E1C=E1L(L)-2.*A1*SQRT(E1L(L))*CO1+A1*A1 E2C=E2L(L,II)-2.*A2*SQRT(E2L(L,II))*CO2+A2*A2 E3C=ETOT-E1C-E2C P1=SQRT((2.*E1C)/XM1)*XM1 P2=SQRT((2.*E2C)/XM2)*XM2 P3=SQRT((2.*E3C)/XM3)*XM3 IF(IPAP.EQ.0) GO TO 10 WRITE(IDOP,110) E1L(L),E2L(L,II),E12(L,II),E31(L,II),E23(L), F TH1C(L),TH2C(L,II),TH3C(L,II),THCP1(L,II),THCP2(L,II) F ,THCP3(L,II),DIF1,DIF2,DIF3,RHO1(L,II),RHO2(L,II) GO TO 5 10 CONTINUE WRITE(IDOP,119)E1L(L),E2L(L,II),E12(L,II),E31(L,II),E23(L) F ,TH1C(L),TH2C(L,II),TH3C(L,II),RHO1(L,II) 5 CONTINUE 4 CONTINUE WRITE(IDOP,107) WRITE(IDOP,107) C C WRITE(ISOP,100) 100 FORMAT(' (1) NEW REACTION (2) NEW ANGLES (0) TO STOP'/) READ(ISIP,*) I IF(I.EQ.1) GOTO 6 IF(I.EQ.2) GOTO 207 STOP 101 FORMAT(1H ,I3,3H + ,I3,5H -- ,I3,3H + ,I3,3H + ,I3,5X, C F7.2,11HMEV Q = ,F6.2,4X,5HEX = ,F6.2) 102 FORMAT(25H PARTICLE 1 MAXIMUM ANGLE ,F9.3, F 25H PARTICLE 2 MAXIMUM ANGLE ,F9.3) 103 FORMAT(7H E1LMIN ,F8.3,7H E1LMAX ,F8.3,7H E2LMIN,F8.3, F 7H E2LMAX ,F8.3) 104 FORMAT(1H ,I3,3X,I3,5X,I3,3X,I3,3X,I3) 105 FORMAT(14H ANGLES POL ,2F6.2) 106 FORMAT(14H ANGLES AZM ,2F6.2) 107 FORMAT(1H//) 110 FORMAT(1H ,2F6.2,2X,3F6.1,2X,9F7.1,3X,E8.2,2X,E8.2) 111 FORMAT(1H ,I1,8H.BRANCH ) 112 FORMAT(2X,36H E1L E2L E12 E31 E32 , F 53HTHC1 THC2 THC3 T1C23 T2C13 T3C12 TC12 TC13 , F 8H TC23 ,3X,16H RHO1 RHO2/) 113 FORMAT(14H CM VELOCITY ,F7.2) 118 FORMAT(32H E1L E2L E12 E31 E23 , F 4X,26HTHC1 THC2 THC3 RHO1 /) 119 FORMAT(1H ,2F6.2,2X,3F6.2,1X,3F7.2,2X,E8.2) END C C SUBROUTINE KDRE C C G.G.OHLSEN,NUCL.INSTR.AND METH. 37(1965)240 C COMMON E1LMAX,E1LMIN,E2LMAX,E2LMIN,TH1MAX,TH2MAX,E0, F XMP,XMT,XM1,XM2,XM3,Q,EPL,DELE1L,TH1L,TH2L,PHI2L,LMAX, F LMIN,ICODE COMMON/COMON1/ E1L(250),E2L(250,2),E23(250),E31(250,2),E12(250,2), F TH1C(250),TH2C(250,2),TH3C(250,2),THCP1(250,2) COMMON/COMON2/THCP2(250,2), F THCP3(250,2),RHO1(250,2),RHO2(250,2),TH3L(250,2),PHI3L(250,2) F ,IM(5),IZ(5),EX(5),ER(5) DIMENSION PCX(3),PCY(3),PCZ(3),PC(3),THC(3),THCP(3), F VRELX(3),VRELY(3),VRELZ(3),VREL(3) ICODE=1 DO 1 I=1,250 E2L(I,1)=0.0 1 CONTINUE RAD=3.141592654/180. TH1LR=TH1L*RAD TH2LR=TH2L*RAD PHI2LR=PHI2L*RAD XM=XM1+XM2+XM3 XM31=XM1+XM3 XM12=XM1+XM2 XM23=XM2+XM3 PP=SQRT(2.0*XMP*EPL) VV=PP/(XMP+XMT) A1=VV*SQRT(XM1/2.0) A2=VV*SQRT(XM2/2.0) ETOT=Q+XMT*EPL/(XMP+XMT) IF(ETOT) 2,2,3 2 WRITE(ISOP,101) ICODE=2 RETURN 3 C1L=COS(TH1LR) S1L=SIN(TH1LR) C2L=COS(TH2LR) S2L=SIN(TH2LR) CPHI2L=COS(PHI2LR) SPHI2L=SIN(PHI2LR) C C COMPUTE MAXIMUM POSSIBLE ENERGIES AND ANGLES C E1CMAX=XM23*ETOT/XM E2CMAX=XM31*ETOT/XM IF(A1**2-E1CMAX) 4,4,5 4 TH1MAX=180. GO TO 6 5 TH1MAX=ARCTAN(SQRT(E1CMAX),SQRT(A1**2-E1CMAX)) 6 IF(A2**2-E2CMAX)7,7,8 7 TH2MAX=180.0 GO TO 9 8 TH2MAX=ARCTAN(SQRT(E2CMAX),SQRT(A2**2-E2CMAX)) 9 IF(TH1MAX-TH1L)11,11,10 10 IF(TH2MAX-TH2L)11,11,12 11 WRITE(ISOP,100) ICODE=2 RETURN 12 T1=SQRT(E1CMAX-(A1*S1L)**2) T2=SQRT(E2CMAX-(A2*S2L)**2) E1LMAX=(A1*C1L+T1)**2 IF(A1*C1L-T1)13,13,14 13 E1LMIN=0.0 GO TO 15 14 E1LMIN=(A1*C1L-T1)**2 15 E2LMAX=(A2*C2L+T2)**2 IF(A2*C2L-T2)16,16,17 16 E2LMIN=0.0 GO TO 18 17 E2LMIN=(A2*C2L-T2)**2 C C MAIN LOOP C 18 DO 30 L=1,250 LEL=L-1 E1L(L)=DELE1L*FLOAT(LEL)+E0 IF(E1L(L).LE.0.) GO TO 30 IF(E1L(L)-E1LMAX) 19,19,31 19 E1C=E1L(L)-2.0*SQRT(E1L(L))*A1*C1L+A1**2 C1C=(SQRT(E1L(L))*C1L-A1)/SQRT(E1C) E123=E1C*XM/XM23 E23(L)=ETOT-E123 IF(E23(L))30,30,22 22 DVSQ=2.0*XM3*XM*(A1*C1L/SQRT(E1L(L))-1.0)/(XM2*XM23**2) V23C=SQRT(2.0*XM1*E1C)/XM23 UX=V23C*SQRT(1.0-C1C*C1C) UZ=VV-V23C*C1C ALPHA=UX*S2L*CPHI2L+UZ*C2L BETA=2.0*XM3*E23(L)/(XM2*XM23)-UX*UX-UZ*UZ UXP=S1L*SQRT(XM1/(2.0*E1L(L)))/XM23 UZP=-C1L*SQRT(XM1/(2.0*E1L(L)))/XM23 ALPHAP=UXP*S2L*CPHI2L+UZP*C2L BETAP=DVSQ-2.0*(UX*UXP+UZ*UZP) P1=SQRT(2.0*XM1*E1L(L)) PCX(1)=-P1*S1L PCY(1)=0.0 PCZ(1)=P1*C1L-A1*SQRT(2.0*XM1) PC(1)=SQRT(PCX(1)**2+PCZ(1)**2) IF(ALPHA**2+BETA)30,23,23 23 IF(BETA)24,25,25 24 N=2 GO TO 26 25 N=1 E2L(L,2)=-0.0 E31(L,2)=-0.0 E12(L,2)=-0.0 TH2C(L,2)=-0.0 TH3C(L,2)=-0.0 RHO1(L,2)=-0.0 RHO2(L,2)=-0.0 THCP1(L,2)=-0.0 THCP2(L,2)=-0.0 THCP3(L,2)=-0.0 26 DO 29 K=1,N RTAB=SQRT(ALPHA**2+BETA) SGN=-(-1.0)**K E2L(L,K)=0.5*XM2*(ALPHA+SGN*RTAB)**2 P2=SQRT(2.0*XM2*E2L(L,K)) POP2=PP*P2*C2L P1OP2=P1*P2*(C1L*C2L-S1L*S2L*CPHI2L) DENOM=SQRT((XM23-(POP2-P1OP2)*XM2/P2**2)**2) RHO1(L,K)=XM1*XM2*XM3*P1*P2/DENOM DE2=XM2*(ALPHA+SGN*RTAB)*(ALPHAP+SGN*(2.0*ALPHA*ALPHAP+ F BETAP)/(2.0*RTAB)) RHO2(L,K)=RHO1(L,K)/SQRT(1.0+DE2*DE2) TH1C(L)=ARCTAN(SQRT(1.0-C1C*C1C),C1C) PCX(2)=P2*S2L*CPHI2L PCY(2)=P2*S2L*SPHI2L PCZ(2)=P2*C2L-A2*SQRT(2.0*XM2) PC(2)=SQRT(PCX(2)**2+PCY(2)**2+PCZ(2)**2) PCX(3)=-PCX(1)-PCX(2) PCY(3)=-PCY(1)-PCY(2) PCZ(3)=-PCZ(1)-PCZ(2) PC(3)=SQRT(PCX(3)**2+PCY(3)**2+PCZ(3)**2) E31(L,K)=ETOT-XM*PC(2)**2/(2.0*XM2*XM31) E12(L,K)=ETOT-XM*PC(3)**2/(2.0*XM3*XM12) VRELX(1)=PCX(2)/XM2-PCX(3)/XM3 VRELY(1)=PCY(2)/XM2-PCY(3)/XM3 VRELZ(1)=PCZ(2)/XM2-PCZ(3)/XM3 VRELX(2)=PCX(3)/XM3-PCX(1)/XM1 VRELY(2)=PCY(3)/XM3-PCY(1)/XM1 VRELZ(2)=PCZ(3)/XM3-PCZ(1)/XM1 VRELX(3)=PCX(1)/XM1-PCX(2)/XM2 VRELY(3)=PCY(1)/XM1-PCY(2)/XM2 VRELZ(3)=PCZ(1)/XM1-PCZ(2)/XM2 VREL(1)=SQRT(VRELX(1)**2+VRELY(1)**2+VRELZ(1)**2) VREL(2)=SQRT(VRELX(2)**2+VRELY(2)**2+VRELZ(2)**2) VREL(3)=SQRT(VRELX(3)**2+VRELY(3)**2+VRELZ(3)**2) DO 27 J=1,3 CTH=(PCX(J)*VRELX(J)+PCY(J)*VRELY(J)+PCZ(J)*VRELZ(J))/(PC(J) F *VREL(J)) THCP(J)=ARCTAN(SQRT(1.0-CTH*CTH),CTH) 27 CONTINUE THCP1(L,K)=THCP(1) THCP2(L,K)=THCP(2) THCP3(L,K)=THCP(3) DO 28 I=2,3 CIC=PCZ(I)/PC(I) THC(I)=ARCTAN(SQRT(1.0-CIC*CIC),CIC) 28 CONTINUE TH2C(L,K)=THC(2) TH3C(L,K)=THC(3) TH3L(L,K)=ARCTAN(SQRT(PCX(3)**2+PCY(3)**2),(PCZ(3)+XM3*VV)) PHI3L(L,K)=ARCTAN(PCY(3),PCX(3)) 29 CONTINUE 30 CONTINUE C C DETERMINE LMAX AND LMIN C 31 DO 34 I=1,250 IF(E2L(I,1))34,34,35 34 CONTINUE 35 LMIN=I DO 38 I=LMIN,250 IF(E2L(I,1)) 39,39,38 38 CONTINUE 39 LMAX=I-1 RETURN 100 FORMAT(20H UNACCEPTABLE ANGLES ) 101 FORMAT(16H BELOW THRESHOLD ) END C C FUNCTION ARCTAN(Y,X) PI=3.141592654 IF(X)1,4,8 1 IF(Y)2,3,3 2 ARCTAN=ATAN(Y/X)-PI GO TO 9 3 ARCTAN=ATAN(Y/X)+PI GO TO 9 4 IF(Y)5,6,7 5 ARCTAN=-PI/2.0 GO TO 9 6 ARCTAN=0. GO TO 9 7 ARCTAN=PI/2.0 GO TO 9 8 ARCTAN=ATAN(Y/X) 9 ARCTAN=ARCTAN*180./PI RETURN END