From 9a99de0fb1c0c4546683caccecf39f64da102b98 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Thu, 21 Feb 2013 10:01:02 +0000 Subject: [PATCH] Cleanup code in progress. --- trimsp/src/trimsp7l.F | 2539 +++++++++++++++++++---------------------- 1 file changed, 1201 insertions(+), 1338 deletions(-) diff --git a/trimsp/src/trimsp7l.F b/trimsp/src/trimsp7l.F index 2cce871..a8eb93e 100644 --- a/trimsp/src/trimsp7l.F +++ b/trimsp/src/trimsp7l.F @@ -1153,8 +1153,7 @@ C GET MAX AND MIN INDEX OF TEST FAILURES IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 4106 GO TO 4104 - 4106 DO 4108 IV=1,IH1 -C IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 4106 ' + 4106 DO IV=1,IH1 ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) SQE=DSQRT(EPS(IV)) CC=(.009611D0+SQE)/(.005175D0+SQE) @@ -1163,33 +1162,34 @@ C IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 4106 ' DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0) C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0) C2(IV)=DMIN1(1.0D0,C*C) - 4108 S2(IV)=1.D0-(1.D0*C2(IV)) + S2(IV)=1.D0-(1.D0*C2(IV)) + ENDDO GO TO 4103 -C + 4102 IF(IPOT.NE.3) GO TO 4103 C ZBL POTENTIAL C CALL MAGICZBL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) - 5104 DO 5105 IV=IVMIN,IVMAX + 5104 DO IV=IVMIN,IVMAX IF(R(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 5104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - R(IV)=0.00001D0 - ENDIF - EX1(IV)=DEXP(-.20162D0*R(IV)) - EX2(IV)=DEXP(-.4029D0*R(IV)) - EX3(IV)=DEXP(-.94229D0*R(IV)) - EX4(IV)=DEXP(-3.1998D0*R(IV)) -C IF(R(IV).EQ.0.D0)WRITE(99,'(A50)')' R nach Label 5104 ' - RR1=1.D0/R(IV) - V(IV)=(.02817D0*EX1(IV)+.28022D0*EX2(IV)+.50986D0*EX3(IV)+ - 1 .18175D0*EX4(IV))*RR1 - FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) - V1(IV)=-(V(IV)+.0056796354D0*EX1(IV)+.112900638D0*EX2(IV)+ - 1 .4804359794D0*EX3(IV)+.58156365D0*EX4(IV))*RR1 - FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 - Q=FR/FR1 - R(IV)=R(IV)-Q - TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 - 5105 CONTINUE + WRITE(99,'(A70)' + & )'in DO 5104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + R(IV)=0.00001D0 + ENDIF + EX1(IV)=DEXP(-.20162D0*R(IV)) + EX2(IV)=DEXP(-.4029D0*R(IV)) + EX3(IV)=DEXP(-.94229D0*R(IV)) + EX4(IV)=DEXP(-3.1998D0*R(IV)) + RR1=1.D0/R(IV) + V(IV)=(.02817D0*EX1(IV)+.28022D0*EX2(IV)+.50986D0*EX3(IV)+ + & .18175D0*EX4(IV))*RR1 + FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) + V1(IV)=-(V(IV)+.0056796354D0*EX1(IV)+.112900638D0*EX2(IV)+ + & .4804359794D0*EX3(IV)+.58156365D0*EX4(IV))*RR1 + FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 + Q=FR/FR1 + R(IV)=R(IV)-Q + TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) @@ -1197,31 +1197,30 @@ C GET MAX AND MIN INDEX OF TEST FAILURES IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 5106 GO TO 5104 - 5106 DO 5108 IV=1,IH1 - IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 5106 ' - ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) - SQE=DSQRT(EPS(IV)) - CC=(.011615D0+SQE)/(.0071222D0+SQE) - AA=2.D0*EPS(IV)*(1.D0+(0.99229D0/SQE))*B(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPS(IV))/(14.813D0+EPS(IV))) - DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0) - C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0) - C2(IV)=DMIN1(1.0D0,C*C) - 5108 S2(IV)=1.D0-(1.D0*C2(IV)) + 5106 DO IV=1,IH1 + IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 5106 ' + ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) + SQE=DSQRT(EPS(IV)) + CC=(.011615D0+SQE)/(.0071222D0+SQE) + AA=2.D0*EPS(IV)*(1.D0+(0.99229D0/SQE))*B(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPS(IV))/(14.813D0+EPS(IV) + & )) + DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0) + C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0) + C2(IV)=DMIN1(1.0D0,C*C) + S2(IV)=1.D0-(1.D0*C2(IV)) + ENDDO 4103 CONTINUE C C END OF MAGIC C - DO 65 IV=1,IH1 + DO IV=1,IH1 DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV) -C TAU(IV)=CVMGT(P(IV)*DSQRT(S2(IV)/C2(IV)),0.,KK.EQ.4) IF(C2(IV).LT.1.D-10) THEN -c WRITE(*,*)C2(IV),S2(IV) WRITE(99,'(A50)')' C2 < 10^-10, C2,S2,DEN resettet ' C2(IV)=1.D-10 S2(IV)=1.D0-(1.D0*C2(IV)) DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV) -c WRITE(*,*)C2(IV),S2(IV) ENDIF TAU(IV)=CVMGT(P(IV)*DSQRT(DABS(S2(IV)/C2(IV))),0.D0,KK.EQ.0) TAU(IV)=DMIN1(TAU(IV),LM(LLL(IV))) @@ -1234,59 +1233,59 @@ c WRITE(*,*)C2(IV),S2(IV) CPSI(IV)=CVMGT(TA2,-TA2,CU.GT.0.D0) SPSI(IV)=DABS(TA)*TA2 DEEOR=CVMGT(KOR1(JJJ(IV))*DSQRT(DABS(E(IV)))*EX1(IV),0.D0, - # KDEE1.EQ.2.OR.KDEE1.EQ.3) + & KDEE1.EQ.2.OR.KDEE1.EQ.3) DENS(IV)=DENS(IV)+DEN(IV) DEES(IV)=DEES(IV)+DEEOR - 65 CONTINUE + ENDDO C C DETERMINATION OF NEW FLIGHT DIRECTIONS C CALL DIRCOS(COSX(1),COSY(1),COSZ(1),SINE(1),CPSI(1),SPSI(1) - * ,CPHI(1),SPHI(1),IH1) + & ,CPHI(1),SPHI(1),IH1) 245 CONTINUE C C END OF COLLISION LOOP C C INELASTIC ENERGY LOSS( 5 POSSIBILITIES) C - DO 14 IV=1,IH1 + DO IV=1,IH1 ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV)) TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV)) - 14 CONTINUE + ENDDO GO TO(15,16,17,18,19),KDEE1 - 15 DO 26 IV=1,IH1 - DEE(IV)=CVMGT(0.D0,KLM1(LLL(IV))*ASIGT(IV)*DSQRT(E(IV)), - # X(IV).LT.HLM.OR.X(IV).GT.HLMT) - 26 CONTINUE + 15 DO IV=1,IH1 + DEE(IV)=CVMGT(0.D0,KLM1(LLL(IV))*ASIGT(IV)*DSQRT(E(IV)),X(IV) + & .LT.HLM.OR.X(IV).GT.HLMT) + ENDDO GO TO 40 - 16 DO 21 IV=1,IH1 + 16 DO IV=1,IH1 DEE(IV)=DEES(IV) - 21 CONTINUE + ENDDO GO TO 40 - 17 DO 22 IV=1,IH1 - DEE(IV)=CVMGT(DEES(IV),0.5D0*(KLM1(LLL(IV))*ASIGT(IV)* - # DSQRT(E(IV))+DEES(IV)),X(IV).LT.HLM.OR.X(IV).GT.HLMT) - 22 CONTINUE + 17 DO IV=1,IH1 + DEE(IV)=CVMGT(DEES(IV),0.5D0*(KLM1(LLL(IV))*ASIGT(IV) + & * DSQRT(E(IV))+DEES(IV)),X(IV).LT.HLM.OR.X(IV).GT.HLMT) + ENDDO GO TO 40 - 18 DO 23 IV=1,IH1 + 18 DO IV=1,IH1 SM(IV)=0.D0 EM(IV)=E(IV)*0.001D0/M1 - 23 CONTINUE - DO 66 IV=1,IH1 - DO 66 J=1,NJ(LLL(IV)) - SH(IV,J)=CVMGT(CH1(LLL(IV),J)*DSQRT(EM(IV)) - # ,CH2(LLL(IV),J)*EM(IV)**0.45D0*(CH3(LLL(IV),J)/EM(IV)) - # *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/ - # EM(IV))+CH5(LLL(IV),J)*EM(IV))) - # /(CH2(LLL(IV),J)*EM(IV)**0.45D0+(CH3(LLL(IV),J)/EM(IV)) - # *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/ - # EM(IV))+CH5(LLL(IV),J)*EM(IV)))) - # ,EM(IV).LT.10.D0) - 66 CONTINUE - DO 73 IV=1,IH1 - DO 73 J=1,NJ(LLL(IV)) + ENDDO + DO IV=1,IH1 + DO J=1,NJ(LLL(IV)) + SH(IV,J)=CVMGT(CH1(LLL(IV),J)*DSQRT(EM(IV)) ,CH2(LLL(IV),J) + & *EM(IV)**0.45D0*(CH3(LLL(IV),J)/EM(IV)) *DLOG(DABS(1.D0 + & +(CH4(LLL(IV),J)/ EM(IV))+CH5(LLL(IV),J)*EM(IV))) + & /(CH2(LLL(IV),J)*EM(IV)**0.45D0+(CH3(LLL(IV),J)/EM(IV)) + & *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/ EM(IV))+CH5(LLL(IV),J) + & *EM(IV)))) ,EM(IV).LT.10.D0) + ENDDO + ENDDO + DO IV=1,IH1 + DO J=1,NJ(LLL(IV)) SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) - 73 CONTINUE + ENDDO + ENDDO DO IV=1,IH1 DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10 & .D0) @@ -1428,7 +1427,7 @@ C KO(NREC1,JJR(NREC1,1),1)=1 INOUT(NREC1,1)=SIGN(1.D0,CX(IV)) NPA=NPA+1 - 6 CONTINUE + 6 CONTINUE C IF(NREC1.LT.NUM) GO TO 27 C @@ -1497,7 +1496,7 @@ C ENDDO KK2=KK0R - DO 235 KKR=KK2,0,-1 + DO KKR=KK2,0,-1 C C CHOICE OF COLLISION PARTNERS C @@ -1507,294 +1506,280 @@ C & ,1,DBLE(random))+JT(LRR(IREC1,2)) ENDDO - DO IREC1=1,NREC2 - call ranlux(ran2, 2) - PHIPR=PI2*DBLE(ran2(1)) - CPHIR(IREC1,2)=DCOS(PHIPR) - SPHIR(IREC1,2)=DSIN(PHIPR) - PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(ran2(2))+KKR) - X2(IREC1)=XR(IREC1,2)-PR(IREC1)*CPHIR(IREC1,2)*SXR(IREC1) - PR(IREC1)=CVMGT(1.D10,PR(IREC1),X2(IREC1).LT.0.0D0 .OR.X2(IREC1 - & ).GT.TT) - BR(IREC1)=PR(IREC1)/A(JJR(IREC1,2),JJR(IREC1,1)) - EPSR(IREC1)=F(JJR(IREC1,2),JJR(IREC1,1))*ER(IREC1,2) - ENDDO + DO IREC1=1,NREC2 + call ranlux(ran2, 2) + PHIPR=PI2*DBLE(ran2(1)) + CPHIR(IREC1,2)=DCOS(PHIPR) + SPHIR(IREC1,2)=DSIN(PHIPR) + PR(IREC1)=PDMAX(LRR(IREC1,2))*DSQRT(DBLE(ran2(2))+KKR) + X2(IREC1)=XR(IREC1,2)-PR(IREC1)*CPHIR(IREC1,2)*SXR(IREC1) + PR(IREC1)=CVMGT(1.D10,PR(IREC1),X2(IREC1).LT.0.0D0 .OR + & .X2(IREC1).GT.TT) + BR(IREC1)=PR(IREC1)/A(JJR(IREC1,2),JJR(IREC1,1)) + EPSR(IREC1)=F(JJR(IREC1,2),JJR(IREC1,1))*ER(IREC1,2) + ENDDO - CALL SCOPY(NREC2,BR,1,RR,1) - IVMIN=1 - IVMAX=NREC2 + CALL SCOPY(NREC2,BR,1,RR,1) + IVMIN=1 + IVMAX=NREC2 C C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.) C - IF(IPOTR.NE.1) GO TO 4201 + IF(IPOTR.NE.1) GO TO 4201 C KR-C POTENTIAL C CALL MAGICKRC(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) - 205 DO IV=IVMIN,IVMAX - IF(RR(IV).LT.1.D-20)THEN - WRITE(99,'(A70)' - & )'in DO 205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - RR(IV)=0.00001D0 - ENDIF - EX1R(IV)=DEXP(-.278544D0*RR(IV)) - EX2R(IV)=DEXP(-.637174D0*RR(IV)) - EX3R(IV)=DEXP(-1.919249D0*RR(IV)) - RRR1=1./RR(IV) - VR(IV)=(.190945D0*EX1R(IV)+.473674D0*EX2R(IV)+.335381D0 - & *EX3R(IV))*RRR1 - FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) - V1R(IV)=-(VR(IV)+.053186584080D0*EX1R(IV)+.301812757276D0 - & *EX2R(IV)+.643679648869D0*EX3R(IV))*RRR1 - FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/EPSR(IV) - & -1.D0 - Q=FR/FR1 - RR(IV)=RR(IV)-Q - TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 - ENDDO + 205 DO IV=IVMIN,IVMAX + IF(RR(IV).LT.1.D-20)THEN + WRITE(99,'(A70)' + & )'in DO 205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + RR(IV)=0.00001D0 + ENDIF + EX1R(IV)=DEXP(-.278544D0*RR(IV)) + EX2R(IV)=DEXP(-.637174D0*RR(IV)) + EX3R(IV)=DEXP(-1.919249D0*RR(IV)) + RRR1=1./RR(IV) + VR(IV)=(.190945D0*EX1R(IV)+.473674D0*EX2R(IV)+.335381D0 + & *EX3R(IV))*RRR1 + FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) + V1R(IV)=-(VR(IV)+.053186584080D0*EX1R(IV)+.301812757276D0 + & *EX2R(IV)+.643679648869D0*EX3R(IV))*RRR1 + FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/EPSR(IV + & )-1.D0 + Q=FR/FR1 + RR(IV)=RR(IV)-Q + TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES - IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) - IF(IVMIN.GT.IVMAX) GO TO 207 - IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) - IF(IVMIN.GT.IVMAX) GO TO 207 - GO TO 205 - 207 DO IV=1,NREC2 - ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) - SQE=DSQRT(EPSR(IV)) - CC=(.235800D0+SQE)/(.126000D0+SQE) - AA=2.D0*EPSR(IV)*(1.D0+(1.0144D0/SQE))*BR(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((69350.D0+EPSR(IV)) /(83550.D0 - & +EPSR(IV))) - DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) - C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) - C2R(IV)=DMIN1(1.0D0,C*C) - S2R(IV)=1.D0-C2R(IV) - ENDDO - GO TO 4203 + IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) + IF(IVMIN.GT.IVMAX) GO TO 207 + IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) + IF(IVMIN.GT.IVMAX) GO TO 207 + GO TO 205 + 207 DO IV=1,NREC2 + ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) + SQE=DSQRT(EPSR(IV)) + CC=(.235800D0+SQE)/(.126000D0+SQE) + AA=2.D0*EPSR(IV)*(1.D0+(1.0144D0/SQE))*BR(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((69350.D0+EPSR(IV)) /(83550.D0 + & +EPSR(IV))) + DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) + C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) + C2R(IV)=DMIN1(1.0D0,C*C) + S2R(IV)=1.D0-C2R(IV) + ENDDO + GO TO 4203 - 4201 IF(IPOTR.NE.2) GO TO 4202 + 4201 IF(IPOTR.NE.2) GO TO 4202 C MOLIERE POTENTIAL C CALL MAGICMOL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) - 4205 DO IV=IVMIN,IVMAX - IF(RR(IV).LT.1.D-20)THEN - WRITE(99,'(A70)' - & )'in DO 4205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - RR(IV)=0.00001D0 - ENDIF - EX1R(IV)=DEXP(-.3D0*RR(IV)) - EX2R(IV)=DEXP(-1.2D0*RR(IV)) - EX3R(IV)=DEXP(-6.0D0*RR(IV)) - RRR1=1.D0/RR(IV) - VR(IV)=(.35D0*EX1R(IV)+.55D0*EX2R(IV)+.10D0*EX3R(IV))*RRR1 - FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) - V1R(IV)=-(VR(IV)+.105D0*EX1R(IV)+ .66D0*EX2R(IV)+.6D0*EX3R(IV)) - & *RRR1 - FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ EPSR(IV) - & -1.D0 - Q=FR/FR1 - RR(IV)=RR(IV)-Q - TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 - ENDDO + 4205 DO IV=IVMIN,IVMAX + IF(RR(IV).LT.1.D-20)THEN + WRITE(99,'(A70)' + & )'in DO 4205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + RR(IV)=0.00001D0 + ENDIF + EX1R(IV)=DEXP(-.3D0*RR(IV)) + EX2R(IV)=DEXP(-1.2D0*RR(IV)) + EX3R(IV)=DEXP(-6.0D0*RR(IV)) + RRR1=1.D0/RR(IV) + VR(IV)=(.35D0*EX1R(IV)+.55D0*EX2R(IV)+.10D0*EX3R(IV))*RRR1 + FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) + V1R(IV)=-(VR(IV)+.105D0*EX1R(IV)+ .66D0*EX2R(IV)+.6D0 + & *EX3R(IV))*RRR1 + FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ + & EPSR(IV)-1.D0 + Q=FR/FR1 + RR(IV)=RR(IV)-Q + TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES - IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) - IF(IVMIN.GT.IVMAX) GO TO 4207 - IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) - IF(IVMIN.GT.IVMAX) GO TO 4207 - GO TO 4205 - 4207 DO 4208 IV=1,NREC2 - ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) - SQE=DSQRT(EPSR(IV)) - CC=(.009611D0+SQE)/(.005175D0+SQE) - AA=2.D0*EPSR(IV)*(1.D0+(0.6743D0/SQE))*BR(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPSR(IV))/(10.+EPSR(IV))) - DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) - C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) - C2R(IV)=DMIN1(1.0D0,C*C) - 4208 S2R(IV)=1.D0-C2R(IV) - GO TO 4203 -C - 4202 IF(IPOTR.NE.3) GO TO 4203 + IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) + IF(IVMIN.GT.IVMAX) GO TO 4207 + IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) + IF(IVMIN.GT.IVMAX) GO TO 4207 + GO TO 4205 + 4207 DO IV=1,NREC2 + ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) + SQE=DSQRT(EPSR(IV)) + CC=(.009611D0+SQE)/(.005175D0+SQE) + AA=2.D0*EPSR(IV)*(1.D0+(0.6743D0/SQE))*BR(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPSR(IV))/(10.+EPSR(IV)) + & ) + DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) + C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) + C2R(IV)=DMIN1(1.0D0,C*C) + S2R(IV)=1.D0-C2R(IV) + ENDDO + GO TO 4203 + + 4202 IF(IPOTR.NE.3) GO TO 4203 C ZBL POTENTIAL C CALL MAGICZBL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) - 5205 DO 5206 IV=IVMIN,IVMAX - IF(RR(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 5205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - RR(IV)=0.00001D0 - ENDIF - EX1R(IV)=DEXP(-.20162D0*RR(IV)) - EX2R(IV)=DEXP(-.40290D0*RR(IV)) - EX3R(IV)=DEXP(-.94229D0*RR(IV)) - EX4R(IV)=DEXP(-3.1998D0*RR(IV)) - RRR1=1./RR(IV) - VR(IV)=(.02817D0*EX1R(IV)+.28022D0*EX2R(IV)+.50986D0*EX3R(IV)+ - 1 .18175D0*EX4R(IV))*RRR1 - FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) - V1R(IV)=-(VR(IV)+.0056796354D0*EX1R(IV)+.112900638D0*EX2R(IV)+ - 1 .4804359794D0*EX3R(IV)+.581563650D0*EX4R(IV))*RRR1 - FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ - 1 EPSR(IV)-1.D0 - Q=FR/FR1 - RR(IV)=RR(IV)-Q - TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 - 5206 CONTINUE + 5205 DO 5206 IV=IVMIN,IVMAX + IF(RR(IV).LT.1.D-20)THEN + WRITE(99,'(A70)' + & )'in DO 5205 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + RR(IV)=0.00001D0 + ENDIF + EX1R(IV)=DEXP(-.20162D0*RR(IV)) + EX2R(IV)=DEXP(-.40290D0*RR(IV)) + EX3R(IV)=DEXP(-.94229D0*RR(IV)) + EX4R(IV)=DEXP(-3.1998D0*RR(IV)) + RRR1=1./RR(IV) + VR(IV)=(.02817D0*EX1R(IV)+.28022D0*EX2R(IV)+.50986D0*EX3R(IV + & )+.18175D0*EX4R(IV))*RRR1 + FR=BR(IV)*BR(IV)*RRR1+VR(IV)*RR(IV)/EPSR(IV)-RR(IV) + V1R(IV)=-(VR(IV)+.0056796354D0*EX1R(IV)+.112900638D0*EX2R(IV + & )+.4804359794D0*EX3R(IV)+.581563650D0*EX4R(IV))*RRR1 + FR1=-BR(IV)*BR(IV)*RRR1*RRR1+(VR(IV)+V1R(IV)*RR(IV))/ + & EPSR(IV)-1.D0 + Q=FR/FR1 + RR(IV)=RR(IV)-Q + TEST1(IV)=DABS(Q/RR(IV)).GT.0.001D0 + 5206 CONTINUE C GET MAX AND MIN INDEX OF TEST FAILURES - IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) - IF(IVMIN.GT.IVMAX) GO TO 5207 - IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) - IF(IVMIN.GT.IVMAX) GO TO 5207 - GO TO 5205 - 5207 DO 5208 IV=1,NREC2 - ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) - SQE=DSQRT(EPSR(IV)) - CC=(.011615D0+SQE)/(.0071222D0+SQE) - AA=2.*EPSR(IV)*(1.D0+(0.99229D0/SQE))*BR(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPSR(IV)) - # /(14.813D0+EPSR(IV))) - DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) - C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) - C2R(IV)=DMIN1(1.0D0,C*C) - 5208 S2R(IV)=1.D0-C2R(IV) - 4203 CONTINUE + IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),1) + IF(IVMIN.GT.IVMAX) GO TO 5207 + IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST1(IVMIN),-1) + IF(IVMIN.GT.IVMAX) GO TO 5207 + GO TO 5205 + 5207 DO IV=1,NREC2 + ROCINV=-.5D0*V1R(IV)/(EPSR(IV)-VR(IV)) + SQE=DSQRT(EPSR(IV)) + CC=(.011615D0+SQE)/(.0071222D0+SQE) + AA=2.*EPSR(IV)*(1.D0+(0.99229D0/SQE))*BR(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((9.3066D0+EPSR(IV)) /(14.813D0 + & +EPSR(IV))) + DELTA=(RR(IV)-BR(IV))*AA*FF/(FF+1.D0) + C=(ROCINV*(BR(IV)+DELTA)+1.D0)/(ROCINV*RR(IV)+1.D0) + C2R(IV)=DMIN1(1.0D0,C*C) + S2R(IV)=1.D0-C2R(IV) + ENDDO + 4203 CONTINUE C - DO 237 IREC1=1,NREC2 - T(IREC1)=ER(IREC1,2)*S2R(IREC1)*EC(JJR(IREC1,2),JJR(IREC1,1)) - TS(IREC1)=TS(IREC1)+T(IREC1) - T1=CVMGT(T(IREC1),0.D0,KKR.EQ.3) - TR1=TR1+T1 - DEEORR=CVMGT(0.D0,KOR(JJR(IREC1,2),JJR(IREC1,1))* - # DSQRT(ER(IREC1,2))*EX1R(IREC1),KDEE2.EQ.1) - DEERS(IREC1)=DEERS(IREC1)+DEEORR - TAUR(IREC1)=CVMGT(PR(IREC1)*DSQRT(S2R(IREC1)/C2R(IREC1)),0.D0, - 1 KKR.EQ.0) - TAUR(IREC1)=DMIN1(TAUR(IREC1),LM(LRR(IREC1,2))) - CTR(IREC1)=C2R(IREC1)+C2R(IREC1)-1.D0 - STR(IREC1)=DSQRT(1.D0-CTR(IREC1)*CTR(IREC1)) - CUR(IREC1) = CTR(IREC1)+MU(JJR(IREC1,2),JJR(IREC1,1)) - CUR(IREC1) = CVMGT(CUR(IREC1),1.0D-8,DABS(CUR(IREC1)).GE.1.0D-8) - TAR=STR(IREC1)/CUR(IREC1) - TAR2=1./DSQRT(1.D0+TAR*TAR) - CPSIR(IREC1,2)=TAR2 - SPSIR(IREC1,2)=TAR*TAR2 - 237 CONTINUE -C - CALL DIRCOS(CSXR(1,2),CSYR(1,2),CSZR(1,2),SNXR(1,2), - 1 CPSIR(1,2),SPSIR(1,2),CPHIR(1,2),SPHIR(1,2),NREC2) -C - 235 CONTINUE + DO IREC1=1,NREC2 + T(IREC1)=ER(IREC1,2)*S2R(IREC1)*EC(JJR(IREC1,2),JJR(IREC1,1) + & ) + TS(IREC1)=TS(IREC1)+T(IREC1) + T1=CVMGT(T(IREC1),0.D0,KKR.EQ.3) + TR1=TR1+T1 + DEEORR=CVMGT(0.D0,KOR(JJR(IREC1,2),JJR(IREC1,1)) + & * DSQRT(ER(IREC1,2))*EX1R(IREC1),KDEE2.EQ.1) + DEERS(IREC1)=DEERS(IREC1)+DEEORR + TAUR(IREC1)=CVMGT(PR(IREC1)*DSQRT(S2R(IREC1)/C2R(IREC1)),0 + & .D0,KKR.EQ.0) + TAUR(IREC1)=DMIN1(TAUR(IREC1),LM(LRR(IREC1,2))) + CTR(IREC1)=C2R(IREC1)+C2R(IREC1)-1.D0 + STR(IREC1)=DSQRT(1.D0-CTR(IREC1)*CTR(IREC1)) + CUR(IREC1) = CTR(IREC1)+MU(JJR(IREC1,2),JJR(IREC1,1)) + CUR(IREC1) = CVMGT(CUR(IREC1),1.0D-8,DABS(CUR(IREC1)).GE.1 + & .0D-8) + TAR=STR(IREC1)/CUR(IREC1) + TAR2=1./DSQRT(1.D0+TAR*TAR) + CPSIR(IREC1,2)=TAR2 + SPSIR(IREC1,2)=TAR*TAR2 + ENDDO + CALL DIRCOS(CSXR(1,2),CSYR(1,2),CSZR(1,2),SNXR(1,2), CPSIR(1,2) + & ,SPSIR(1,2),CPHIR(1,2),SPHIR(1,2),NREC2) + ENDDO C C CREATE SECONDARY KNOCK-ON ATOMS C DO 246 IREC1=1,NREC2 -cc IF(T(IREC1).LE.SFE) GO TO 246 - IF(T(IREC1).LE.ERC) GO TO 246 - IF(X2(IREC1).GT.RD.AND.X2(IREC1).LT.RT) GO TO 246 -C -C CALL NEWREC(NREC1,T(IREC1),XR(IREC1,2),YR(IREC1,2),ZR(IREC1,2), -C 1 CXR(IREC1),CYR(IREC1),CZR(IREC1),SXR(IREC1), -C 2 CTR(IREC1),STR(IREC1),PHIR(IREC1,2),PR(IREC1), -C 3 ER(I,1),XR(1,1),YR(1,1),ZR(1,1),PHIR(1,1),PSIR(1,1) -C 4 ,CSXR(1,1),CSYR(1,1),CSZR(1,1),SNXR(1,1),L(1,1) - NREC1=NREC1+1 - ER(NREC1,1)=T(IREC1)-EP(JJR(IREC1,1)) - XR(NREC1,1)=X2(IREC1) - YR(NREC1,1)=YR(IREC1,2)-PR(IREC1)*(SPHIR(IREC1,2)*CZR(IREC1)- - 1 CPHIR(IREC1,2)*CYR(IREC1)*CXR(IREC1))/SXR(IREC1) - ZR(NREC1,1)=ZR(IREC1,2)+PR(IREC1)*(SPHIR(IREC1,2)*CYR(IREC1)+ - 1 CPHIR(IREC1,2)*CXR(IREC1)*CZR(IREC1))/SXR(IREC1) - CSXR(NREC1,1)=CXR(IREC1) - CSYR(NREC1,1)=CYR(IREC1) - CSZR(NREC1,1)=CZR(IREC1) - SNXR(NREC1,1)=SXR(IREC1) - CPHIR(NREC1,1)=-CPHIR(IREC1,2) - SPHIR(NREC1,1)=-SPHIR(IREC1,2) - CTR(NREC1)=DMIN1(CTR(IREC1),.99999999D0) - TAR=STR(IREC1)/(1.D0-CTR(NREC1)) - TAR2=1./DSQRT(1.D0+TAR*TAR) - CPSIR(NREC1,1)=TAR2 - SPSIR(NREC1,1)=TAR*TAR2 - TAUPSR(NREC1,1)=0.D0 - KO(NREC1,JJR(IREC1,1),1)=0 - INOUT(NREC1,1)=INOUT(IREC1,2) - JJR(NREC1,1)=JJR(IREC1,1) - NSA=NSA+1 - 246 CONTINUE + IF(T(IREC1).LE.ERC) GO TO 246 + IF(X2(IREC1).GT.RD.AND.X2(IREC1).LT.RT) GO TO 246 + NREC1=NREC1+1 + ER(NREC1,1)=T(IREC1)-EP(JJR(IREC1,1)) + XR(NREC1,1)=X2(IREC1) + YR(NREC1,1)=YR(IREC1,2)-PR(IREC1)*(SPHIR(IREC1,2)*CZR(IREC1)- + & CPHIR(IREC1,2)*CYR(IREC1)*CXR(IREC1))/SXR(IREC1) + ZR(NREC1,1)=ZR(IREC1,2)+PR(IREC1)*(SPHIR(IREC1,2)*CYR(IREC1)+ + & CPHIR(IREC1,2)*CXR(IREC1)*CZR(IREC1))/SXR(IREC1) + CSXR(NREC1,1)=CXR(IREC1) + CSYR(NREC1,1)=CYR(IREC1) + CSZR(NREC1,1)=CZR(IREC1) + SNXR(NREC1,1)=SXR(IREC1) + CPHIR(NREC1,1)=-CPHIR(IREC1,2) + SPHIR(NREC1,1)=-SPHIR(IREC1,2) + CTR(NREC1)=DMIN1(CTR(IREC1),.99999999D0) + TAR=STR(IREC1)/(1.D0-CTR(NREC1)) + TAR2=1./DSQRT(1.D0+TAR*TAR) + CPSIR(NREC1,1)=TAR2 + SPSIR(NREC1,1)=TAR*TAR2 + TAUPSR(NREC1,1)=0.D0 + KO(NREC1,JJR(IREC1,1),1)=0 + INOUT(NREC1,1)=INOUT(IREC1,2) + JJR(NREC1,1)=JJR(IREC1,1) + NSA=NSA+1 + 246 CONTINUE C C INELASTIC ENERGY LOSS C - DO 238 IREC1=1,NREC2 - ASIGTR(IREC1)=(LM(LRR(IREC1,2))-TAUR(IREC1)+ - # TAUPSR(IREC1,2))*ARHO(LRR(IREC1,2)) - TAUPSR(IREC1,2)=TAUR(IREC1)*DABS(CPSIR(IREC1,2)) - 238 CONTINUE + DO IREC1=1,NREC2 + ASIGTR(IREC1)=(LM(LRR(IREC1,2))-TAUR(IREC1)+ TAUPSR(IREC1,2)) + & *ARHO(LRR(IREC1,2)) + TAUPSR(IREC1,2)=TAUR(IREC1)*DABS(CPSIR(IREC1,2)) + ENDDO GO TO(115,116,117),KDEE2 - 115 DO 241 IREC1=1,NREC2 - DEER(IREC1)=CVMGT(0.D0,KLM(LRR(IREC1,2), - & JJR(IREC1,2))*ASIGTR(IREC1)*DSQRT(ER(IREC1,2)), - # XR(IREC1,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) -cTR DEER(IREC1)=CVMGT(0.,KLM(LRR(IREC1,2),JJR(IREC1,2))*ASIGTR(IREC1)* -cTR # DSQRT(ER(IREC1,2)),XR(IREC1,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) - 241 CONTINUE + 115 DO IREC1=1,NREC2 + DEER(IREC1)=CVMGT(0.D0,KLM(LRR(IREC1,2), JJR(IREC1,2)) + & *ASIGTR(IREC1)*DSQRT(ER(IREC1,2)), XR(IREC1,2).LT.HLM.OR + & .XR(IREC1,2).GT.HLMT) + ENDDO GO TO 242 - 116 DO 243 IREC1=1,NREC2 - DEER(IREC1)=DEERS(IREC1) - 243 CONTINUE + 116 DO IREC1=1,NREC2 + DEER(IREC1)=DEERS(IREC1) + ENDDO GO TO 242 - 117 DO 244 IREC1=1,NREC2 - DEER(IREC1)=CVMGT(DEERS(IREC1),.5*(KLM(LRR(IREC1,2),JJR(IREC1,2))* - 1 ASIGTR(IREC1)*DSQRT(ER(IREC1,2))+DEERS(IREC1)), - 2 XR(IREC1,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) - 244 CONTINUE - 242 CONTINUE -C - DO 344 IREC1=1,NREC2 - DELR=DMAX1(1.0D-20,TS(IREC1)+DEER(IREC1)) - TS(IREC1)=CVMGT(ER(IREC1,2)*TS(IREC1)/DELR,TS(IREC1) - 1 ,DELR.GT.ER(IREC1,2)) - DEER(IREC1)=CVMGT(ER(IREC1,2)*DEER(IREC1)/DELR,DEER(IREC1) - 1 ,DELR.GT.ER(IREC1,2)) - 344 CONTINUE + 117 DO IREC1=1,NREC2 + DEER(IREC1)=CVMGT(DEERS(IREC1),.5*(KLM(LRR(IREC1,2),JJR(IREC1,2 + & ))*ASIGTR(IREC1)*DSQRT(ER(IREC1,2))+DEERS(IREC1)),XR(IREC1 + & ,2).LT.HLM.OR.XR(IREC1,2).GT.HLMT) + ENDDO + 242 CONTINUE C + DO IREC1=1,NREC2 + DELR=DMAX1(1.0D-20,TS(IREC1)+DEER(IREC1)) + TS(IREC1)=CVMGT(ER(IREC1,2)*TS(IREC1)/DELR,TS(IREC1) ,DELR.GT + & .ER(IREC1,2)) + DEER(IREC1)=CVMGT(ER(IREC1,2)*DEER(IREC1)/DELR,DEER(IREC1) + & ,DELR.GT.ER(IREC1,2)) + ENDDO + DO 252 IREC1=1,NREC2 -C IF(XR(IREC1,2).LT.0.0.OR.XR(IREC1,2).GT.TT) GO TO 252 - I=MAX0(MIN0(IDINT(X2(IREC1)/CW+1.D0),100),1) -C I=DMAX(MIN(INT(XR(IREC1,2)/CW+1.),100),1) - DENTR(I)=DENTR(I)+TS(IREC1) - DMGNR(I)=DMGNR(I)+T(IREC1) - IONR(I)=IONR(I)+DEER(IREC1) -C ELER(I,JJR(IREC1,2))=ELER(I,JJR(IREC1,2))+T(IREC1) -C ELIR(I,JJR(IREC1,2))=ELIR(I,JJR(IREC1,2))+DEER(IREC1) - IF(T(IREC1).LE.DI(JJR(IREC1,1))) GO TO 84 -CC EPS(IV)=F1(JJJ(IV))*DEN(IV) -CC G=EPS(IV)+.40244*EPS(IV)**.75+3.4008*EPS(IV)**.16667 -CC MOT=DEN(IV)/(1.+K2(LLL(IV))*G) -CC CASMOT(I)=CASMOT(I)+MOT - ELGDR(I)=ELGDR(I)+T(IREC1) -C ELDR(I,JJR(IREC1,2))=ELDR(I,JJR(IREC1,2))+T(IREC1) - ICDR(I,JJR(IREC1,2))=ICDR(I,JJR(IREC1,2))+1 - ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))= - 1 ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))+1 - GO TO 252 - 84 PHONR(I)=PHONR(I)+T(IREC1) -C ELPR(I,JJJ(IREC1,2))=ELPR(I,JJJ(IREC1,2))+T(IREC1) - 252 CONTINUE - DO 82 IREC1=1,NREC2 - ICDIR=ICDIR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1)))) - ICSBR=ICSBR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.SB(1))) - ICSUMR=ICSUMR+IDINT(CVMGT(1.D0,0.D0,TS(IREC1).GT.0.D0)) - ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= - 1 ICDIRJ(JJR(IREC1,2),JJR(IREC1,1)) - 2 +IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1)))) - 82 CONTINUE - DO 253 IREC1=1,NREC2 - TEL=TEL+TS(IREC1) - TR=TR+T(IREC1) - EI=EI+DEER(IREC1) - ER(IREC1,2)=ER(IREC1,2)-DEER(IREC1)-TS(IREC1)+1.0D-10 - XR(IREC1,2)=XR(IREC1,2)-TAUR(IREC1)*CXR(IREC1) - YR(IREC1,2)=YR(IREC1,2)-TAUR(IREC1)*CYR(IREC1) - ZR(IREC1,2)=ZR(IREC1,2)-TAUR(IREC1)*CZR(IREC1) - TESTR(IREC1)=ER(IREC1,2).LE.SFE.OR.XR(IREC1,2).LT.(-SU) - 1 .OR.XR(IREC1,2).GT.SUT - 2 .OR.(XR(IREC1,2).GT.RD.AND.XR(IREC1,2).LT.RT) - 253 CONTINUE + I=MAX0(MIN0(IDINT(X2(IREC1)/CW+1.D0),100),1) + DENTR(I)=DENTR(I)+TS(IREC1) + DMGNR(I)=DMGNR(I)+T(IREC1) + IONR(I)=IONR(I)+DEER(IREC1) + IF(T(IREC1).LE.DI(JJR(IREC1,1))) GO TO 84 + ELGDR(I)=ELGDR(I)+T(IREC1) + ICDR(I,JJR(IREC1,2))=ICDR(I,JJR(IREC1,2))+1 + ICDIRI(I,JJR(IREC1,2),JJR(IREC1,1))= ICDIRI(I,JJR(IREC1,2) + & ,JJR(IREC1,1))+1 + GO TO 252 + 84 PHONR(I)=PHONR(I)+T(IREC1) + 252 CONTINUE + DO IREC1=1,NREC2 + ICDIR=ICDIR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.DI(JJR(IREC1,1))) + & ) + ICSBR=ICSBR+IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT.SB(1))) + ICSUMR=ICSUMR+IDINT(CVMGT(1.D0,0.D0,TS(IREC1).GT.0.D0)) + ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= ICDIRJ(JJR(IREC1,2) + & ,JJR(IREC1,1)) +IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT + & .DI(JJR(IREC1,1)))) + ENDDO + DO IREC1=1,NREC2 + TEL=TEL+TS(IREC1) + TR=TR+T(IREC1) + EI=EI+DEER(IREC1) + ER(IREC1,2)=ER(IREC1,2)-DEER(IREC1)-TS(IREC1)+1.0D-10 + XR(IREC1,2)=XR(IREC1,2)-TAUR(IREC1)*CXR(IREC1) + YR(IREC1,2)=YR(IREC1,2)-TAUR(IREC1)*CYR(IREC1) + ZR(IREC1,2)=ZR(IREC1,2)-TAUR(IREC1)*CZR(IREC1) + TESTR(IREC1)=ER(IREC1,2).LE.SFE.OR.XR(IREC1,2).LT.(-SU) .OR + & .XR(IREC1,2).GT.SUT .OR.(XR(IREC1,2).GT.RD.AND.XR(IREC1,2) + & .LT.RT) + ENDDO C C CHECK TO SEE IF ANY RECOIL ATOM IS SPUTTERED OR C IF THE ENERGY OF ANY RECOIL ATOM IS TOO LOW @@ -1804,894 +1789,818 @@ C IVMAX=NREC2-ILLZ(NREC2,TESTR,-1) C DO 248 IREC1=IVMIN,IVMAX - 249 IF(IREC1.GT.NREC2) GO TO 247 - IF(XR(IREC1,2).LT.(-SU)) GO TO 251 - IF(XR(IREC1,2).GT.SUT) GO TO 250 -cc IF(ER(IREC1,2).LE.SFE) GO TO 255 - IF(ER(IREC1,2).LE.ERC) GO TO 255 - IF(XR(IREC1,2).GT.RD.AND.XR(IREC1,2).LT.RT) GO TO 255 - GO TO 248 - 251 ENOR=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1) - IF(ENOR.GT.SB(1)) GO TO 254 + 249 IF(IREC1.GT.NREC2) GO TO 247 + IF(XR(IREC1,2).LT.(-SU)) GO TO 251 + IF(XR(IREC1,2).GT.SUT) GO TO 250 + IF(ER(IREC1,2).LE.ERC) GO TO 255 + IF(XR(IREC1,2).GT.RD.AND.XR(IREC1,2).LT.RT) GO TO 255 + GO TO 248 + 251 ENOR=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1) + IF(ENOR.GT.SB(1)) GO TO 254 C C RECOIL ATOM IS REFLECTED BACK INTO THE SOLID BY THE C POTENTIAL BARRIER C - XR(IREC1,2)=-1.D0*SU - CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2) - KIS=KIS+1 - GO TO 248 + XR(IREC1,2)=-1.D0*SU + CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2) + KIS=KIS+1 + GO TO 248 C C RECOIL ATOM IS SPUTTERED (BACKWARD) C - 254 ESP=ER(IREC1,2)-SB(1) + 254 ESP=ER(IREC1,2)-SB(1) CC254 ESP=ER(IREC1,2)-SB(LRR(IREC1,2)) C C NUMBER, ENERGY AND MOMENTS OF ALL SPUTTERED PARTICLES C - IBSP(JJR(IREC1,2))=IBSP(JJR(IREC1,2))+1 - EBSP(JJR(IREC1,2))=EBSP(JJR(IREC1,2))+ESP - SPE2=ESP*ESP - SPE3=SPE2*ESP - SPE2S(JJR(IREC1,2))=SPE2S(JJR(IREC1,2))+SPE2 - SPE3S(JJR(IREC1,2))=SPE3S(JJR(IREC1,2))+SPE3 - SPE4S(JJR(IREC1,2))=SPE4S(JJR(IREC1,2))+SPE2*SPE2 - SPE5S(JJR(IREC1,2))=SPE5S(JJR(IREC1,2))+SPE3*SPE2 - SPE6S(JJR(IREC1,2))=SPE6S(JJR(IREC1,2))+SPE3*SPE3 - IF(ESP.LT.0.1) GO TO 256 - IBSPL(JJR(IREC1,2))=IBSPL(JJR(IREC1,2))+1 - SPE1SL(JJR(IREC1,2))=SPE1SL(JJR(IREC1,2))+DLOG10(DABS(ESP)) - SPE2L=(DLOG10(DABS(ESP)))**2.D0 - SPE3L=SPE2L*DLOG10(DABS(ESP)) - SPE2SL(JJR(IREC1,2))=SPE2SL(JJR(IREC1,2))+SPE2L - SPE3SL(JJR(IREC1,2))=SPE3SL(JJR(IREC1,2))+SPE3L - SPE4SL(JJR(IREC1,2))=SPE4SL(JJR(IREC1,2))+SPE2L*SPE2L - SPE5SL(JJR(IREC1,2))=SPE5SL(JJR(IREC1,2))+SPE3L*SPE2L - SPE6SL(JJR(IREC1,2))=SPE6SL(JJR(IREC1,2))+SPE3L*SPE3L - 256 CONTINUE + IBSP(JJR(IREC1,2))=IBSP(JJR(IREC1,2))+1 + EBSP(JJR(IREC1,2))=EBSP(JJR(IREC1,2))+ESP + SPE2=ESP*ESP + SPE3=SPE2*ESP + SPE2S(JJR(IREC1,2))=SPE2S(JJR(IREC1,2))+SPE2 + SPE3S(JJR(IREC1,2))=SPE3S(JJR(IREC1,2))+SPE3 + SPE4S(JJR(IREC1,2))=SPE4S(JJR(IREC1,2))+SPE2*SPE2 + SPE5S(JJR(IREC1,2))=SPE5S(JJR(IREC1,2))+SPE3*SPE2 + SPE6S(JJR(IREC1,2))=SPE6S(JJR(IREC1,2))+SPE3*SPE3 + IF(ESP.LT.0.1) GO TO 256 + IBSPL(JJR(IREC1,2))=IBSPL(JJR(IREC1,2))+1 + SPE1SL(JJR(IREC1,2))=SPE1SL(JJR(IREC1,2))+DLOG10(DABS(ESP)) + SPE2L=(DLOG10(DABS(ESP)))**2.D0 + SPE3L=SPE2L*DLOG10(DABS(ESP)) + SPE2SL(JJR(IREC1,2))=SPE2SL(JJR(IREC1,2))+SPE2L + SPE3SL(JJR(IREC1,2))=SPE3SL(JJR(IREC1,2))+SPE3L + SPE4SL(JJR(IREC1,2))=SPE4SL(JJR(IREC1,2))+SPE2L*SPE2L + SPE5SL(JJR(IREC1,2))=SPE5SL(JJR(IREC1,2))+SPE3L*SPE2L + SPE6SL(JJR(IREC1,2))=SPE6SL(JJR(IREC1,2))+SPE3L*SPE3L + 256 CONTINUE C C SURFACE REFRACTION C - EXIR=DSQRT((ENOR-SB(1))/ESP) - IF ( EXIR .GE. 1.D0 ) EXIR = .999999D0 + EXIR=DSQRT((ENOR-SB(1))/ESP) + IF ( EXIR .GE. 1.D0 ) EXIR = .999999D0 C C TOTAL ANGULAR DISTRIBUTIONS C - IAG=IDINT(EXIR*20.D0+1.D0) - KADS(IAG)=KADS(IAG)+1 - KADSJ(IAG,JJR(IREC1,2))=KADSJ(IAG,JJR(IREC1,2))+1 -C + IAG=IDINT(EXIR*20.D0+1.D0) + KADS(IAG)=KADS(IAG)+1 + KADSJ(IAG,JJR(IREC1,2))=KADSJ(IAG,JJR(IREC1,2))+1 +C C 4 GROUPS:ION IN , PKA ;ION IN , SKA ;ION OUT, PKA ;ION OUT, SKA C - KOI=KO(IREC1,JJR(IREC1,2),2) - IF(INOUT(IREC1,2).EQ.-1) GO TO 61 - IF(KOI.EQ.0) GO TO 62 - ISPIP(JJR(IREC1,2))=ISPIP(JJR(IREC1,2))+1 - ESPIP(JJR(IREC1,2))=ESPIP(JJR(IREC1,2))+ESP - KADRIP(IAG,JJR(IREC1,2))=KADRIP(IAG,JJR(IREC1,2))+1 - GO TO 164 - 62 ISPIS(JJR(IREC1,2))=ISPIS(JJR(IREC1,2))+1 - ESPIS(JJR(IREC1,2))=ESPIS(JJR(IREC1,2))+ESP - KADRIS(IAG,JJR(IREC1,2))=KADRIS(IAG,JJR(IREC1,2))+1 - GO TO 164 - 61 IF(KOI.EQ.0) GO TO 163 - ISPOP(JJR(IREC1,2))=ISPOP(JJR(IREC1,2))+1 - ESPOP(JJR(IREC1,2))=ESPOP(JJR(IREC1,2))+ESP - KADROP(IAG,JJR(IREC1,2))=KADROP(IAG,JJR(IREC1,2))+1 - GO TO 164 - 163 ISPOS(JJR(IREC1,2))=ISPOS(JJR(IREC1,2))+1 - ESPOS(JJR(IREC1,2))=ESPOS(JJR(IREC1,2))+ESP - KADROS(IAG,JJR(IREC1,2))=KADROS(IAG,JJR(IREC1,2))+1 - 164 CONTINUE + KOI=KO(IREC1,JJR(IREC1,2),2) + IF(INOUT(IREC1,2).EQ.-1) GO TO 61 + IF(KOI.EQ.0) GO TO 62 + ISPIP(JJR(IREC1,2))=ISPIP(JJR(IREC1,2))+1 + ESPIP(JJR(IREC1,2))=ESPIP(JJR(IREC1,2))+ESP + KADRIP(IAG,JJR(IREC1,2))=KADRIP(IAG,JJR(IREC1,2))+1 + GO TO 164 + 62 ISPIS(JJR(IREC1,2))=ISPIS(JJR(IREC1,2))+1 + ESPIS(JJR(IREC1,2))=ESPIS(JJR(IREC1,2))+ESP + KADRIS(IAG,JJR(IREC1,2))=KADRIS(IAG,JJR(IREC1,2))+1 + GO TO 164 + 61 IF(KOI.EQ.0) GO TO 163 + ISPOP(JJR(IREC1,2))=ISPOP(JJR(IREC1,2))+1 + ESPOP(JJR(IREC1,2))=ESPOP(JJR(IREC1,2))+ESP + KADROP(IAG,JJR(IREC1,2))=KADROP(IAG,JJR(IREC1,2))+1 + GO TO 164 + 163 ISPOS(JJR(IREC1,2))=ISPOS(JJR(IREC1,2))+1 + ESPOS(JJR(IREC1,2))=ESPOS(JJR(IREC1,2))+ESP + KADROS(IAG,JJR(IREC1,2))=KADROS(IAG,JJR(IREC1,2))+1 + 164 CONTINUE C C OUTPUT MATRICES OF BACKWARD SPUTTERED ATOMS C - IF(JJR(IREC1,2).GT.JT(3)) GO TO 255 - IESP = IDINT(ESP*E0DE+2.0D0) - IESP = MIN0(101,IESP) - IESLOG=2 - IF(ESP.LE.0.1D0) GO TO 75 - IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESP))+3.D0) - IESLOG=MIN0(IESLOG,75) - 75 CONTINUE - IF(EXIR.GT.1.0D0)WRITE(99,'(A50)')' EXIR nach Label 75' - IA=IDINT(DAW*DACOS(EXIR)+2.D0) - IAGS = IAG+1 - IG=2 -C IF(SXR(I).EQ.0.) GO TO 182 -C SXR(I)=DMAX1(SXR(I),1.0E-12) - SXR(IREC1)=DMAX1(SXR(IREC1),1.0D-12) - U=CYR(IREC1)/SXR(IREC1) - IF(DABS(U).GT.1.D0) U = SIGN(1.D0,U) - IF(U.GT.1.0D0)WRITE(99,'(A50)')' U vor Label 182 CON..' - ACS=DACOS(U) - IG=IDINT(DGW*ACS+2.D0) -C 182 CONTINUE - IGG=IDINT(DGIK*ACS+1.D0) - IGG=MAX0(MIN0(NGIK,IGG),1) - MAGS(IG,IAGS,JJR(IREC1,2)) = MAGS(IG,IAGS,JJR(IREC1,2))+1 - MAGS(IG,22,JJR(IREC1,2)) = MAGS(IG,22,JJR(IREC1,2))+1 - MAGS(NG,IAGS,JJR(IREC1,2)) = MAGS(NG,IAGS,JJR(IREC1,2))+1 - MAGS(NG,22,JJR(IREC1,2)) = MAGS(NG,22,JJR(IREC1,2))+1 - MEAS(IESP,IAGS,JJR(IREC1,2)) = MEAS(IESP,IAGS,JJR(IREC1,2))+1 - MEAS(102,IAGS,JJR(IREC1,2)) = MEAS(102,IAGS,JJR(IREC1,2))+1 - MEAS(IESP,22,JJR(IREC1,2)) = MEAS(IESP,22,JJR(IREC1,2))+1 - MEAS(102,22,JJR(IREC1,2)) = MEAS(102,22,JJR(IREC1,2))+1 - MEASL(IESLOG,IAG,JJR(IREC1,2)) = MEASL(IESLOG,IAG,JJR(IREC1,2))+1 - MEASL(IESLOG,21,JJR(IREC1,2)) = MEASL(IESLOG,21,JJR(IREC1,2))+1 - MEASL(75,IAG,JJR(IREC1,2)) = MEASL(75,IAG,JJR(IREC1,2))+1 - MEASL(75,21,JJR(IREC1,2)) = MEASL(75,21,JJR(IREC1,2))+1 - IF(ALPHA.LT.1.0) GO TO 181 - MEAGS(IESP,IGG,IAGS,JJR(IREC1,2)) = - * MEAGS(IESP,IGG,IAGS,JJR(IREC1,2))+1 - MEAGS(102,IGG,IAGS,JJR(IREC1,2)) = - * MEAGS(102,IGG,IAGS,JJR(IREC1,2))+1 - MEAGS(IESP,IGG,22,JJR(IREC1,2)) = - * MEAGS(IESP,IGG,22,JJR(IREC1,2))+1 - MEAGS(102,IGG,22,JJR(IREC1,2)) = - * MEAGS(102,IGG,22,JJR(IREC1,2))+1 -C MEAGSL(IESLOG,IGG,IAG)=MEAGSL(IESLOG,IGG,IAG)+1 -C MEAGSL(IESLOG,IGG,21)=MEAGSL(IESLOG,IGG,21)+1 -C MEAGSL(75,IGG,IAG)=MEAGSL(75,IGG,IAG)+1 - MAGSA(IG,IA,JJR(IREC1,2)) = MAGSA(IG,IA,JJR(IREC1,2))+1 - MAGSA(IG,32,JJR(IREC1,2)) = MAGSA(IG,32,JJR(IREC1,2))+1 - MAGSA(NG,IA,JJR(IREC1,2)) = MAGSA(NG,IA,JJR(IREC1,2))+1 - MAGSA(NG,32,JJR(IREC1,2)) = MAGSA(NG,32,JJR(IREC1,2))+1 - 181 CONTINUE - GO TO 255 + IF(JJR(IREC1,2).GT.JT(3)) GO TO 255 + IESP = IDINT(ESP*E0DE+2.0D0) + IESP = MIN0(101,IESP) + IESLOG=2 + IF(ESP.LE.0.1D0) GO TO 75 + IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESP))+3.D0) + IESLOG=MIN0(IESLOG,75) + 75 CONTINUE + IF(EXIR.GT.1.0D0)WRITE(99,'(A50)')' EXIR nach Label 75' + IA=IDINT(DAW*DACOS(EXIR)+2.D0) + IAGS = IAG+1 + IG=2 + SXR(IREC1)=DMAX1(SXR(IREC1),1.0D-12) + U=CYR(IREC1)/SXR(IREC1) + IF(DABS(U).GT.1.D0) U = SIGN(1.D0,U) + IF(U.GT.1.0D0)WRITE(99,'(A50)')' U vor Label 182 CON..' + ACS=DACOS(U) + IG=IDINT(DGW*ACS+2.D0) + IGG=IDINT(DGIK*ACS+1.D0) + IGG=MAX0(MIN0(NGIK,IGG),1) + MAGS(IG,IAGS,JJR(IREC1,2))=MAGS(IG,IAGS,JJR(IREC1,2))+1 + MAGS(IG,22,JJR(IREC1,2))=MAGS(IG,22,JJR(IREC1,2))+1 + MAGS(NG,IAGS,JJR(IREC1,2))=MAGS(NG,IAGS,JJR(IREC1,2))+1 + MAGS(NG,22,JJR(IREC1,2))=MAGS(NG,22,JJR(IREC1,2))+1 + MEAS(IESP,IAGS,JJR(IREC1,2))=MEAS(IESP,IAGS,JJR(IREC1,2))+1 + MEAS(102,IAGS,JJR(IREC1,2))=MEAS(102,IAGS,JJR(IREC1,2))+1 + MEAS(IESP,22,JJR(IREC1,2))=MEAS(IESP,22,JJR(IREC1,2))+1 + MEAS(102,22,JJR(IREC1,2))=MEAS(102,22,JJR(IREC1,2))+1 + MEASL(IESLOG,IAG,JJR(IREC1,2))=MEASL(IESLOG,IAG,JJR(IREC1,2))+1 + MEASL(IESLOG,21,JJR(IREC1,2))=MEASL(IESLOG,21,JJR(IREC1,2))+1 + MEASL(75,IAG,JJR(IREC1,2))=MEASL(75,IAG,JJR(IREC1,2))+1 + MEASL(75,21,JJR(IREC1,2))=MEASL(75,21,JJR(IREC1,2))+1 + IF(ALPHA.LT.1.0) GO TO 181 + MEAGS(IESP,IGG,IAGS,JJR(IREC1,2)) = MEAGS(IESP,IGG,IAGS + & ,JJR(IREC1,2))+1 + MEAGS(102,IGG,IAGS,JJR(IREC1,2)) = MEAGS(102,IGG,IAGS + & ,JJR(IREC1,2))+1 + MEAGS(IESP,IGG,22,JJR(IREC1,2)) = MEAGS(IESP,IGG,22,JJR(IREC1 + & ,2))+1 + MEAGS(102,IGG,22,JJR(IREC1,2)) = MEAGS(102,IGG,22,JJR(IREC1 + & ,2))+1 + MAGSA(IG,IA,JJR(IREC1,2)) = MAGSA(IG,IA,JJR(IREC1,2))+1 + MAGSA(IG,32,JJR(IREC1,2)) = MAGSA(IG,32,JJR(IREC1,2))+1 + MAGSA(NG,IA,JJR(IREC1,2)) = MAGSA(NG,IA,JJR(IREC1,2))+1 + MAGSA(NG,32,JJR(IREC1,2)) = MAGSA(NG,32,JJR(IREC1,2))+1 + 181 CONTINUE + GO TO 255 C - 250 ENORT=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1) - IF(ENORT.GT.SB(L)) GO TO 257 + 250 ENORT=ER(IREC1,2)*CXR(IREC1)*CXR(IREC1) + IF(ENORT.GT.SB(L)) GO TO 257 C C RECOIL ATOM IS REFLECTED BACK INTO THE SOLID BY THE C POTENTIAL BARRIER C - XR(IREC1,2)=SUT - CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2) - KIST=KIST+1 - GO TO 248 + XR(IREC1,2)=SUT + CSXR(IREC1,2)=-1.D0*CSXR(IREC1,2) + KIST=KIST+1 + GO TO 248 C C RECOIL ATOM IS SPUTTERED (TRANSMISSION) C - 257 ESPT=ER(IREC1,2)-SB(L) + 257 ESPT=ER(IREC1,2)-SB(L) C C NUMBER AND ENERGY OF ALL SPUTTERED PARTICLES C - ITSP(JJR(IREC1,2))=ITSP(JJR(IREC1,2))+1 - ETSP(JJR(IREC1,2))=ETSP(JJR(IREC1,2))+ESPT + ITSP(JJR(IREC1,2))=ITSP(JJR(IREC1,2))+1 + ETSP(JJR(IREC1,2))=ETSP(JJR(IREC1,2))+ESPT C C SURFACE REFRACTION C - EXIRT=DSQRT((ENORT-SB(L))/ESPT) - IF ( EXIRT .GE. 1.D0 ) EXIRT = .999999D0 + EXIRT=DSQRT((ENORT-SB(L))/ESPT) + IF ( EXIRT .GE. 1.D0 ) EXIRT = .999999D0 C C TOTAL ANGULAR DISTRIBUTIONS C - IAG=IDINT(EXIRT*20.D0+1.D0) - KADST(IAG)=KADST(IAG)+1 - KDSTJ(IAG,JJR(IREC1,2))=KDSTJ(IAG,JJR(IREC1,2))+1 + IAG=IDINT(EXIRT*20.D0+1.D0) + KADST(IAG)=KADST(IAG)+1 + KDSTJ(IAG,JJR(IREC1,2))=KDSTJ(IAG,JJR(IREC1,2))+1 C C 4 GROUPS:ION IN , PKA ;ION IN , SKA ;ION OUT, PKA ;ION OUT, SKA C - KOI=KO(IREC1,JJR(IREC1,2),2) - IF(INOUT(IREC1,2).EQ.-1) GO TO 85 - IF(KOI.EQ.0) GO TO 86 - ISPIPT(JJR(IREC1,2))=ISPIPT(JJR(IREC1,2))+1 - ESPIPT(JJR(IREC1,2))=ESPIPT(JJR(IREC1,2))+ESPT -C KADRIP(IAG,JJR(IREC1,2))=KADRIP(IAG,JJR(IREC1,2))+1 - GO TO 88 - 86 ISPIST(JJR(IREC1,2))=ISPIST(JJR(IREC1,2))+1 - ESPIST(JJR(IREC1,2))=ESPIST(JJR(IREC1,2))+ESPT -C KADRIS(IAG,JJR(IREC1,2))=KADRIS(IAG,JJR(IREC1,2))+1 - GO TO 88 - 85 IF(KOI.EQ.0) GO TO 87 - ISPOPT(JJR(IREC1,2))=ISPOPT(JJR(IREC1,2))+1 - ESPOPT(JJR(IREC1,2))=ESPOPT(JJR(IREC1,2))+ESPT -C KADROP(IAG,JJR(IREC1,2))=KADROP(IAG,JJR(IREC1,2))+1 - GO TO 88 - 87 ISPOST(JJR(IREC1,2))=ISPOST(JJR(IREC1,2))+1 - ESPOST(JJR(IREC1,2))=ESPOST(JJR(IREC1,2))+ESPT -C KADROS(IAG,JJR(IREC1,2))=KADROS(IAG,JJR(IREC1,2))+1 - 88 CONTINUE + KOI=KO(IREC1,JJR(IREC1,2),2) + IF(INOUT(IREC1,2).EQ.-1) GO TO 85 + IF(KOI.EQ.0) GO TO 86 + ISPIPT(JJR(IREC1,2))=ISPIPT(JJR(IREC1,2))+1 + ESPIPT(JJR(IREC1,2))=ESPIPT(JJR(IREC1,2))+ESPT + GO TO 88 + 86 ISPIST(JJR(IREC1,2))=ISPIST(JJR(IREC1,2))+1 + ESPIST(JJR(IREC1,2))=ESPIST(JJR(IREC1,2))+ESPT + GO TO 88 + 85 IF(KOI.EQ.0) GO TO 87 + ISPOPT(JJR(IREC1,2))=ISPOPT(JJR(IREC1,2))+1 + ESPOPT(JJR(IREC1,2))=ESPOPT(JJR(IREC1,2))+ESPT + GO TO 88 + 87 ISPOST(JJR(IREC1,2))=ISPOST(JJR(IREC1,2))+1 + ESPOST(JJR(IREC1,2))=ESPOST(JJR(IREC1,2))+ESPT + 88 CONTINUE C C OUTPUT MATRICES OF FORWARD SPUTTERED ATOMS C - JRT=JJR(IREC1,2) - IF(L.EQ.3) JRT=JJR(IREC1,2)-NJ(1) - IF(JRT.LT.1) GO TO 255 - IESPT = IDINT(ESPT*E0DE+2.0D0) - IESPT = MIN0(101,IESPT) - IESLOG=2 - IF(ESPT.LE.0.1D0) GO TO 76 - IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESPT))+3.D0) - IESLOG=MIN0(IESLOG,75) - 76 CONTINUE - IAGS = IAG+1 - MAGST(IG,IAGS,JRT) = MAGST(IG,IAGS,JRT)+1 - MAGST(IG,22,JRT) = MAGST(IG,22,JRT)+1 - MAGST(NG,IAGS,JRT) = MAGST(NG,IAGS,JRT)+1 - MAGST(NG,22,JRT) = MAGST(NG,22,JRT)+1 - MEAST(IESPT,IAGS,JRT) = MEAST(IESPT,IAGS,JRT)+1 - MEAST(102,IAGS,JRT) = MEAST(102,IAGS,JRT)+1 - MEAST(IESPT,22,JRT) = MEAST(IESPT,22,JRT)+1 - MEAST(102,22,JRT) = MEAST(102,22,JRT)+1 - MEASTL(IESLOG,IAG,JRT) = MEASTL(IESLOG,IAG,JRT)+1 - MEASTL(IESLOG,21,JRT) = MEASTL(IESLOG,21,JRT)+1 - MEASTL(75,IAG,JRT) = MEASTL(75,IAG,JRT)+1 - MEASTL(75,21,JRT) = MEASTL(75,21,JRT)+1 -C IF(ALPHA.LT.1.0) GO TO 181 -C MEAGS(IESPT,IGG,IAGS) = MEAGS(IESPT,IGG,IAGS)+1 -C MEAGS(102,IGG,IAGS) = MEAGS(102,IGG,IAGS)+1 -C MEAGS(IESPT,IGG,22) = MEAGS(IESPT,IGG,22)+1 -C MEAGSL(IESLOG,IGG,IAG)=MEAGSL(IESLOG,IGG,IAG)+1 -C MEAGSL(IESLOG,IGG,21)=MEAGSL(IESLOG,IGG,21)+1 -C MEAGSL(75,IGG,IAG)=MEAGSL(75,IGG,IAG)+1 -C 181 CONTINUE + JRT=JJR(IREC1,2) + IF(L.EQ.3) JRT=JJR(IREC1,2)-NJ(1) + IF(JRT.LT.1) GO TO 255 + IESPT = IDINT(ESPT*E0DE+2.0D0) + IESPT = MIN0(101,IESPT) + IESLOG=2 + IF(ESPT.LE.0.1D0) GO TO 76 + IESLOG=IDINT(12.D0*DLOG10(DABS(10.D0*ESPT))+3.D0) + IESLOG=MIN0(IESLOG,75) + 76 CONTINUE + IAGS = IAG+1 + MAGST(IG,IAGS,JRT) = MAGST(IG,IAGS,JRT)+1 + MAGST(IG,22,JRT) = MAGST(IG,22,JRT)+1 + MAGST(NG,IAGS,JRT) = MAGST(NG,IAGS,JRT)+1 + MAGST(NG,22,JRT) = MAGST(NG,22,JRT)+1 + MEAST(IESPT,IAGS,JRT) = MEAST(IESPT,IAGS,JRT)+1 + MEAST(102,IAGS,JRT) = MEAST(102,IAGS,JRT)+1 + MEAST(IESPT,22,JRT) = MEAST(IESPT,22,JRT)+1 + MEAST(102,22,JRT) = MEAST(102,22,JRT)+1 + MEASTL(IESLOG,IAG,JRT) = MEASTL(IESLOG,IAG,JRT)+1 + MEASTL(IESLOG,21,JRT) = MEASTL(IESLOG,21,JRT)+1 + MEASTL(75,IAG,JRT) = MEASTL(75,IAG,JRT)+1 + MEASTL(75,21,JRT) = MEASTL(75,21,JRT)+1 C C REARRANGEMENT OF PARTICLES IN LIST 2 C - 255 ER(IREC1,2)=ER(NREC2,2) - XR(IREC1,2)=XR(NREC2,2) - YR(IREC1,2)=YR(NREC2,2) - ZR(IREC1,2)=ZR(NREC2,2) - CSXR(IREC1,2)=CSXR(NREC2,2) - CSYR(IREC1,2)=CSYR(NREC2,2) - CSZR(IREC1,2)=CSZR(NREC2,2) - SNXR(IREC1,2)=SNXR(NREC2,2) - CPHIR(IREC1,2)=CPHIR(NREC2,2) - SPHIR(IREC1,2)=SPHIR(NREC2,2) - CPSIR(IREC1,2)=CPSIR(NREC2,2) - SPSIR(IREC1,2)=SPSIR(NREC2,2) - TAUPSR(IREC1,2)=TAUPSR(NREC2,2) - JJR(IREC1,2)=JJR(NREC2,2) - KO(IREC1,JJR(IREC1,2),2)=KO(NREC2,JJR(NREC2,2),2) - INOUT(IREC1,2)=INOUT(NREC2,2) - NREC2=NREC2-1 -C - IF(IREC1.EQ.NREC2+1) GO TO 247 + 255 ER(IREC1,2)=ER(NREC2,2) + XR(IREC1,2)=XR(NREC2,2) + YR(IREC1,2)=YR(NREC2,2) + ZR(IREC1,2)=ZR(NREC2,2) + CSXR(IREC1,2)=CSXR(NREC2,2) + CSYR(IREC1,2)=CSYR(NREC2,2) + CSZR(IREC1,2)=CSZR(NREC2,2) + SNXR(IREC1,2)=SNXR(NREC2,2) + CPHIR(IREC1,2)=CPHIR(NREC2,2) + SPHIR(IREC1,2)=SPHIR(NREC2,2) + CPSIR(IREC1,2)=CPSIR(NREC2,2) + SPSIR(IREC1,2)=SPSIR(NREC2,2) + TAUPSR(IREC1,2)=TAUPSR(NREC2,2) + JJR(IREC1,2)=JJR(NREC2,2) + KO(IREC1,JJR(IREC1,2),2)=KO(NREC2,JJR(NREC2,2),2) + INOUT(IREC1,2)=INOUT(NREC2,2) + NREC2=NREC2-1 +C + IF(IREC1.EQ.NREC2+1) GO TO 247 C THE NREC2 PARTICLE FAILS THE TEST - IF(NREC2+1.GT.IVMAX) GO TO 248 - GO TO 249 - 248 CONTINUE -C - 247 CONTINUE -C + IF(NREC2+1.GT.IVMAX) GO TO 248 + GO TO 249 + 248 CONTINUE + 247 CONTINUE IF(NREC1+NREC2.EQ.0) GO TO 27 IF(NREC2.GE.NUM.OR.IH1.EQ.0) GO TO 83 C C END OF RECOIL ATOM SECTION C - 27 CONTINUE -C + 27 CONTINUE IF(IH1.EQ.0.AND.IH.EQ.NH) GO TO 140 C C PROJECTILE CANDIDATE FOR REFLECTION C - DO 29 IV=1,IH1 - E(IV)=E(IV)-DEE(IV)-DENS(IV)+1.0D-10 - X(IV)=X(IV)-TAU(IV)*CX(IV) - Y(IV)=Y(IV)-TAU(IV)*CY(IV) - Z(IV)=Z(IV)-TAU(IV)*CZ(IV) - PL(IV)=PL(IV)-TAU(IV) - TEST(IV)=E(IV).LE.EF.OR.X(IV).LT.-1.D0*SU.OR.X(IV).GT.SUT - 29 CONTINUE + DO IV=1,IH1 + E(IV)=E(IV)-DEE(IV)-DENS(IV)+1.0D-10 + X(IV)=X(IV)-TAU(IV)*CX(IV) + Y(IV)=Y(IV)-TAU(IV)*CY(IV) + Z(IV)=Z(IV)-TAU(IV)*CZ(IV) + PL(IV)=PL(IV)-TAU(IV) + TEST(IV)=E(IV).LE.EF.OR.X(IV).LT.-1.D0*SU.OR.X(IV).GT.SUT + ENDDO IVMIN=1+ILLZ(IH1,TEST,1) IF(IVMIN.GT.IH1) GO TO 90 IVMAX=IH1-ILLZ(IH1,TEST,-1) DO 120 IV=IVMIN,IVMAX - 160 IF(IV.GT.IH1) GO TO 90 - IF(X(IV).LT.-SU) GO TO 8 - IF(X(IV).GT.SUT) GO TO 9 - IF(E(IV).GT.EF) GO TO 125 - IF(E(IV).GT.ESB.AND.X(IV).LT.0.D0) GO TO 125 - IF(E(IV).GT.ESB.AND.X(IV).GT.TT) GO TO 125 + 160 IF(IV.GT.IH1) GO TO 90 + IF(X(IV).LT.-SU) GO TO 8 + IF(X(IV).GT.SUT) GO TO 9 + IF(E(IV).GT.EF) GO TO 125 + IF(E(IV).GT.ESB.AND.X(IV).LT.0.D0) GO TO 125 + IF(E(IV).GT.ESB.AND.X(IV).GT.TT) GO TO 125 C C PROJECTILE HAS STOPPED (PATHLENGTH,RANGE,SPREAD AND MOMENTS) C -C IF(X(IV).LT.0..OR.X(IV).GT.TT) GO TO 110 - IP = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) - IPL(IP)=IPL(IP)+1 - I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), 101), 0) - IRP(I1)=IRP(I1)+1 -c -c Berechnung der gestoppten Teilchen im jeweiligen Layer -c - LowTiefe = 0.D0 - UpTiefe = DX(1) -c - DO laufzahl=1,l - IF(X(IV).GT.LowTiefe.AND.X(IV).LE.UpTiefe) THEN - Number_in_Layer(laufzahl)=Number_in_Layer(laufzahl)+1 - ENDIF - LowTiefe = UpTiefe - UpTiefe = UpTiefe+DX(laufzahl+1) - ENDDO -c - PL2=PL(IV)*PL(IV) - PL3=PL2*PL(IV) - PLSUM=PLSUM+PL(IV) - PL2SUM=PL2SUM+PL2 - PL3SUM=PL3SUM+PL3 - PL4SUM=PL4SUM+PL2*PL2 - PL5SUM=PL5SUM+PL3*PL2 - PL6SUM=PL6SUM+PL3*PL3 - IF(X(IV).LT.0.D0.OR.X(IV).GT.TT) GO TO 111 - XQ=X(IV)*X(IV) - XQ3=XQ*X(IV) - XSUM=XSUM+X(IV) - X2SUM=X2SUM+XQ - X3SUM=X3SUM+XQ3 - X4SUM=X4SUM+XQ*XQ - X5SUM=X5SUM+XQ3*XQ - X6SUM=X6SUM+XQ3*XQ3 - RQ=Y(IV)*Y(IV)+Z(IV)*Z(IV) - RQW=DSQRT(RQ) - RQ3=RQ*RQW - RSUM=RSUM+RQW - R2SUM=R2SUM+RQ - R3SUM=R3SUM+RQ3 - R4SUM=R4SUM+RQ*RQ - R5SUM=R5SUM+RQ3*RQ - R6SUM=R6SUM+RQ3*RQ3 - 111 CONTINUE - ENUCLI=ENUCLI+ENUCL(IV) - ENL2I=ENL2I+ENUCL(IV)*ENUCL(IV) - ENUCL(IV)=0.D0 - EINELI=EINELI+EINEL(IV) - EIL2I=EIL2I+EINEL(IV)*EINEL(IV) - EINEL(IV)=0.D0 - GO TO 110 - 8 ENO=E(IV)*CX(IV)*CX(IV) - IF(ENO.LE.ESB) GO TO 24 + IP = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) + IPL(IP)=IPL(IP)+1 + I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), 101), 0) + IRP(I1)=IRP(I1)+1 +C +C Berechnung der gestoppten Teilchen im jeweiligen Layer +C + LowTiefe = 0.D0 + UpTiefe = DX(1) + DO laufzahl=1,l + IF(X(IV).GT.LowTiefe.AND.X(IV).LE.UpTiefe) THEN + Number_in_Layer(laufzahl)=Number_in_Layer(laufzahl)+1 + ENDIF + LowTiefe = UpTiefe + UpTiefe = UpTiefe+DX(laufzahl+1) + ENDDO + + PL2=PL(IV)*PL(IV) + PL3=PL2*PL(IV) + PLSUM=PLSUM+PL(IV) + PL2SUM=PL2SUM+PL2 + PL3SUM=PL3SUM+PL3 + PL4SUM=PL4SUM+PL2*PL2 + PL5SUM=PL5SUM+PL3*PL2 + PL6SUM=PL6SUM+PL3*PL3 + IF(X(IV).LT.0.D0.OR.X(IV).GT.TT) GO TO 111 + XQ=X(IV)*X(IV) + XQ3=XQ*X(IV) + XSUM=XSUM+X(IV) + X2SUM=X2SUM+XQ + X3SUM=X3SUM+XQ3 + X4SUM=X4SUM+XQ*XQ + X5SUM=X5SUM+XQ3*XQ + X6SUM=X6SUM+XQ3*XQ3 + RQ=Y(IV)*Y(IV)+Z(IV)*Z(IV) + RQW=DSQRT(RQ) + RQ3=RQ*RQW + RSUM=RSUM+RQW + R2SUM=R2SUM+RQ + R3SUM=R3SUM+RQ3 + R4SUM=R4SUM+RQ*RQ + R5SUM=R5SUM+RQ3*RQ + R6SUM=R6SUM+RQ3*RQ3 + 111 CONTINUE + ENUCLI=ENUCLI+ENUCL(IV) + ENL2I=ENL2I+ENUCL(IV)*ENUCL(IV) + ENUCL(IV)=0.D0 + EINELI=EINELI+EINEL(IV) + EIL2I=EIL2I+EINEL(IV)*EINEL(IV) + EINEL(IV)=0.D0 + GO TO 110 + 8 ENO=E(IV)*CX(IV)*CX(IV) + IF(ENO.LE.ESB) GO TO 24 C C PROJECTILE IS BACKSCATTERED C - IB=IB+1 - ES=E(IV)-ESB -C IJKLMN=IJKLMN+1 -C ESVDL(IJKLMN)=ES - ESQ=ES*ES - ES3=ESQ*ES - EB=EB+ES - EB2SUM=EB2SUM+ESQ - EB3SUM=EB3SUM+ES3 - EB4SUM=EB4SUM+ESQ*ESQ - EB5SUM=EB5SUM+ES3*ESQ - EB6SUM=EB6SUM+ES3*ES3 - IF(ES.LT.0.1D0) GO TO 112 - IBL=IBL+1 - ESQL=(DLOG10(DABS(ES)))**2.D0 - ES3L=ESQL*DLOG10(DABS(ES)) - EB1SUL=EB1SUL+DLOG10(DABS(ES)) - EB2SUL=EB2SUL+ESQL - EB3SUL=EB3SUL+ES3L - EB4SUL=EB4SUL+ESQL*ESQL - EB5SUL=EB5SUL+ES3L*ESQL - EB6SUL=EB6SUL+ES3L*ES3L - 112 CONTINUE - IPB = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) - IPLB(IPB)=IPLB(IPB)+1 - PLQB=PL(IV)*PL(IV) - PL3B=PLQB*PL(IV) - PLSB=PLSB+PL(IV) - PL2SB=PL2SB+PLQB - PL3SB=PL3SB+PL3B - PL4SB=PL4SB+PLQB*PLQB - PL5SB=PL5SB+PL3B*PLQB - PL6SB=PL6SB+PL3B*PL3B - ENUCLB=ENUCLB+ENUCL(IV) - ENL2B=ENL2B+ENUCL(IV)*ENUCL(IV) - ENUCL(IV)=0.D0 - EINELB=EINELB+EINEL(IV) - EIL2B=EIL2B+EINEL(IV)*EINEL(IV) - EINEL(IV)=0.D0 + IB=IB+1 + ES=E(IV)-ESB + ESQ=ES*ES + ES3=ESQ*ES + EB=EB+ES + EB2SUM=EB2SUM+ESQ + EB3SUM=EB3SUM+ES3 + EB4SUM=EB4SUM+ESQ*ESQ + EB5SUM=EB5SUM+ES3*ESQ + EB6SUM=EB6SUM+ES3*ES3 + IF(ES.LT.0.1D0) GO TO 112 + IBL=IBL+1 + ESQL=(DLOG10(DABS(ES)))**2.D0 + ES3L=ESQL*DLOG10(DABS(ES)) + EB1SUL=EB1SUL+DLOG10(DABS(ES)) + EB2SUL=EB2SUL+ESQL + EB3SUL=EB3SUL+ES3L + EB4SUL=EB4SUL+ESQL*ESQL + EB5SUL=EB5SUL+ES3L*ESQL + EB6SUL=EB6SUL+ES3L*ES3L + 112 CONTINUE + IPB = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) + IPLB(IPB)=IPLB(IPB)+1 + PLQB=PL(IV)*PL(IV) + PL3B=PLQB*PL(IV) + PLSB=PLSB+PL(IV) + PL2SB=PL2SB+PLQB + PL3SB=PL3SB+PL3B + PL4SB=PL4SB+PLQB*PLQB + PL5SB=PL5SB+PL3B*PLQB + PL6SB=PL6SB+PL3B*PL3B + ENUCLB=ENUCLB+ENUCL(IV) + ENL2B=ENL2B+ENUCL(IV)*ENUCL(IV) + ENUCL(IV)=0.D0 + EINELB=EINELB+EINEL(IV) + EIL2B=EIL2B+EINEL(IV)*EINEL(IV) + EINEL(IV)=0.D0 C C SURFACE REFRACTION C - EXI=DSQRT((ENO-ESB)/ES) - exi1s=exi1s+exi - exiq=exi*exi - exic=exiq*exi - exi2s=exi2s+exiq - exi3s=exi3s+exic - exi4s=exi4s+exiq*exiq - exi5s=exi5s+exic*exiq - exi6s=exi6s+exic*exic + EXI=DSQRT((ENO-ESB)/ES) + exi1s=exi1s+exi + exiq=exi*exi + exic=exiq*exi + exi2s=exi2s+exiq + exi3s=exi3s+exic + exi4s=exi4s+exiq*exiq + exi5s=exi5s+exic*exiq + exi6s=exi6s+exic*exic C C DIVISIONS FOR VECTORS AND MATRICES C - IE = IDINT(E0DE*ES+2.D0) - IE = MAX0( MIN0( IE,NE1), 2) - IERLOG = 2 - IF(ES.LE.0.1D0) GO TO 4 - IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*ES))+3.D0) - IERLOG=MIN0(IERLOG,75) - 4 CONTINUE - IAG=IDINT(EXI*20.D0+1.D0) - IAG = MIN0( IAG, 20) - IAGB = IAG+1 - KADB(IAG)=KADB(IAG)+1 - IG=2 - COSSIN=CY(IV)/SX(IV) - COSSIN=DMIN1(COSSIN,1.D0) - COSSIN=DMAX1(COSSIN,-1.D0) - coss1s=coss1s+cossin - cossq=cossin*cossin - cosst=cossq*cossin - coss2s=coss2s+cossq - coss3s=coss3s+cosst - coss4s=coss4s+cossq*cossq - coss5s=coss5s+cosst*cossq - coss6s=coss6s+cosst*cosst - IF(COSSIN.GT.1.0D0)WRITE(99,'(A50)')' nach coss6s' - AC=DACOS(COSSIN) - IG=IDINT(DAW*AC+2.D0) - IGG=IDINT(DGIK*AC+1.D0) - IF(IGG.GT.NGIK) IGG=NGIK - IPB1=IPB+1 - MEABL(IERLOG,IAG) = MEABL(IERLOG,IAG)+1 - MEABL(IERLOG,21) = MEABL(IERLOG,21)+1 - MEABL(75,IAG) = MEABL(75,IAG)+1 - MAGB(IG,IAGB) = MAGB(IG,IAGB)+1 - MAGB(NG,IAGB) = MAGB(NG,IAGB)+1 - MAGB(IG,22) = MAGB(IG,22)+1 - MEAB(IE,IAGB) = MEAB(IE,IAGB)+1 - MEAB(NE,IAGB) = MEAB(NE,IAGB)+1 - MEAB(IE,22) = MEAB(IE,22)+1 -C IF(ALPHA.LT.1.0) GO TO 183 - MEAGB(IE,IGG,IAGB) = MEAGB(IE,IGG,IAGB)+1 - MEAGB(102,IGG,IAGB) = MEAGB(102,IGG,IAGB)+1 - MEAGB(IE,IGG,22) = MEAGB(IE,IGG,22)+1 - MEAGB(102,IGG,22) = MEAGB(102,IGG,22)+1 -C 183 CONTINUE - MEPB(IE,IPB1) = MEPB(IE,IPB1)+1 - MEPB(NE,IPB1) = MEPB(NE,IPB1)+1 - MEPB(IE,102) = MEPB(IE,102)+1 - EMA(IG,IAGB) = EMA(IG,IAGB)+ES - EMA(IG,22) = EMA(IG,22)+ES - EMA(NG,IAGB) = EMA(NG,IAGB)+ES - GO TO 110 + IE = IDINT(E0DE*ES+2.D0) + IE = MAX0( MIN0( IE,NE1), 2) + IERLOG = 2 + IF(ES.LE.0.1D0) GO TO 4 + IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*ES))+3.D0) + IERLOG=MIN0(IERLOG,75) + 4 CONTINUE + IAG=IDINT(EXI*20.D0+1.D0) + IAG = MIN0( IAG, 20) + IAGB = IAG+1 + KADB(IAG)=KADB(IAG)+1 + IG=2 + COSSIN=CY(IV)/SX(IV) + COSSIN=DMIN1(COSSIN,1.D0) + COSSIN=DMAX1(COSSIN,-1.D0) + coss1s=coss1s+cossin + cossq=cossin*cossin + cosst=cossq*cossin + coss2s=coss2s+cossq + coss3s=coss3s+cosst + coss4s=coss4s+cossq*cossq + coss5s=coss5s+cosst*cossq + coss6s=coss6s+cosst*cosst + IF(COSSIN.GT.1.0D0)WRITE(99,'(A50)')' nach coss6s' + AC=DACOS(COSSIN) + IG=IDINT(DAW*AC+2.D0) + IGG=IDINT(DGIK*AC+1.D0) + IF(IGG.GT.NGIK) IGG=NGIK + IPB1=IPB+1 + MEABL(IERLOG,IAG) = MEABL(IERLOG,IAG)+1 + MEABL(IERLOG,21) = MEABL(IERLOG,21)+1 + MEABL(75,IAG) = MEABL(75,IAG)+1 + MAGB(IG,IAGB) = MAGB(IG,IAGB)+1 + MAGB(NG,IAGB) = MAGB(NG,IAGB)+1 + MAGB(IG,22) = MAGB(IG,22)+1 + MEAB(IE,IAGB) = MEAB(IE,IAGB)+1 + MEAB(NE,IAGB) = MEAB(NE,IAGB)+1 + MEAB(IE,22) = MEAB(IE,22)+1 + MEAGB(IE,IGG,IAGB) = MEAGB(IE,IGG,IAGB)+1 + MEAGB(102,IGG,IAGB) = MEAGB(102,IGG,IAGB)+1 + MEAGB(IE,IGG,22) = MEAGB(IE,IGG,22)+1 + MEAGB(102,IGG,22) = MEAGB(102,IGG,22)+1 + MEPB(IE,IPB1) = MEPB(IE,IPB1)+1 + MEPB(NE,IPB1) = MEPB(NE,IPB1)+1 + MEPB(IE,102) = MEPB(IE,102)+1 + EMA(IG,IAGB) = EMA(IG,IAGB)+ES + EMA(IG,22) = EMA(IG,22)+ES + EMA(NG,IAGB) = EMA(NG,IAGB)+ES + GO TO 110 C C PROJECTILE IS REFLECTED BACK INTO THE TARGET BY THE SURF. BARRIER C - 24 X(IV)=-1.D0*SU - COSX(IV)=-1.D0*COSX(IV) - KIB=KIB+1 - GO TO 125 + 24 X(IV)=-1.D0*SU + COSX(IV)=-1.D0*COSX(IV) + KIB=KIB+1 + GO TO 125 C C PROJECTILE IS TRANSMITTED C - 9 ENOT=E(IV)*CX(IV)*CX(IV) - IF(ENOT.LE.ESB) GO TO 517 - IT=IT+1 - EST=E(IV)-ESB - ET=ET+EST - ETQ=EST*EST - ET3=ETQ*EST - ET2SUM=ET2SUM+ETQ - ET3SUM=ET3SUM+ET3 - ET4SUM=ET4SUM+ETQ*ETQ - ET5SUM=ET5SUM+ET3*ETQ - ET6SUM=ET6SUM+ET3*ET3 - IPT = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) - IPLT(IP)=IPLT(IP)+1 - PLQT=PL(IV)*PL(IV) - PL3T=PLQT*PL(IV) - PLST=PLST+PL(IV) - PL2ST=PL2ST+PLQT - PL3ST=PL3ST+PL3T - PL4ST=PL4ST+PLQT*PLQT - PL5ST=PL5ST+PL3T*PLQT - PL6ST=PL6ST+PL3T*PL3T - ENUCLT=ENUCLT+ENUCL(IV) - ENL2T=ENL2T+ENUCL(IV)*ENUCL(IV) - ENUCL(IV)=0.D0 - EINELT=EINELT+EINEL(IV) - EIL2T=EIL2T+EINEL(IV)*EINEL(IV) - EINEL(IV)=0.D0 -C + 9 ENOT=E(IV)*CX(IV)*CX(IV) + IF(ENOT.LE.ESB) GO TO 517 + IT=IT+1 + EST=E(IV)-ESB + ET=ET+EST + ETQ=EST*EST + ET3=ETQ*EST + ET2SUM=ET2SUM+ETQ + ET3SUM=ET3SUM+ET3 + ET4SUM=ET4SUM+ETQ*ETQ + ET5SUM=ET5SUM+ET3*ETQ + ET6SUM=ET6SUM+ET3*ET3 + IPT = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) + IPLT(IP)=IPLT(IP)+1 + PLQT=PL(IV)*PL(IV) + PL3T=PLQT*PL(IV) + PLST=PLST+PL(IV) + PL2ST=PL2ST+PLQT + PL3ST=PL3ST+PL3T + PL4ST=PL4ST+PLQT*PLQT + PL5ST=PL5ST+PL3T*PLQT + PL6ST=PL6ST+PL3T*PL3T + ENUCLT=ENUCLT+ENUCL(IV) + ENL2T=ENL2T+ENUCL(IV)*ENUCL(IV) + ENUCL(IV)=0.D0 + EINELT=EINELT+EINEL(IV) + EIL2T=EIL2T+EINEL(IV)*EINEL(IV) + EINEL(IV)=0.D0 +C C SURFACE REFRACTION C - EXI=DSQRT((ENOT-ESB)/EST) + EXI=DSQRT((ENOT-ESB)/EST) C C DIVISIONS FOR VECTORS AND MATRICES C - IE=IDINT(E0DE*EST+2.D0) - IERLOG = 2 - IF(EST.LE.0.1D0) GO TO 5 - IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*EST))+3.D0) - IERLOG=MIN0(IERLOG,75) - 5 CONTINUE - IAG=IDINT(EXI*20.D0+1.D0) - IAG = MIN0( IAG, 20) - IAGB = IAG+1 - KADT(IAG)=KADT(IAG)+1 - IG=2 - COSSIN=CY(IV)/SX(IV) - COSSIN=DMIN1(COSSIN,1.D0) - COSSIN=DMAX1(COSSIN,-1.D0) - IF(COSSIN.GT.1.0D0) WRITE(99,'(A50)')' nach COSSIN' - AC=DACOS(COSSIN) - IG=IDINT(DAW*AC+2.D0) - IGG=IDINT(DGIK*AC+1.D0) - IF(IGG.GT.NGIK) IGG=NGIK - MEATL(IERLOG,IAG) = MEATL(IERLOG,IAG)+1 - MEATL(IERLOG,21) = MEATL(IERLOG,21)+1 - MEATL(75,IAG) = MEATL(75,IAG)+1 - MAGT(IG,IAGB) = MAGT(IG,IAGB)+1 - MAGT(NG,IAGB) = MAGT(NG,IAGB)+1 - MAGT(IG,22) = MAGT(IG,22)+1 - MEAT(IE,IAGB) = MEAT(IE,IAGB)+1 - MEAT(NE,IAGB) = MEAT(NE,IAGB)+1 - MEAT(IE,22) = MEAT(IE,22)+1 -C IF(ALPHA.LT.1.0) GO TO 183 - MEAGT(IE,IGG,IAGB) = MEAGT(IE,IGG,IAGB)+1 - MEAGT(102,IGG,IAGB) = MEAGT(102,IGG,IAGB)+1 - MEAGT(IE,IGG,22) = MEAGT(IE,IGG,22)+1 - MEAGT(102,IGG,22) = MEAGT(102,IGG,22)+1 -C 183 CONTINUE - MEPT(IE,IPT) = MEPT(IE,IPT)+1 - MEPT(NE,IPT) = MEPT(NE,IPT)+1 - MEPT(IE,102) = MEPT(IE,102)+1 - EMAT(IG,IAGB) = EMAT(IG,IAGB)+ES - EMAT(IG,22) = EMAT(IG,22)+ES - EMAT(NG,IAGB) = EMAT(NG,IAGB)+ES - GO TO 110 + IE=IDINT(E0DE*EST+2.D0) + IERLOG = 2 + IF(EST.LE.0.1D0) GO TO 5 + IERLOG = IDINT(12.D0*DLOG10(DABS(10.D0*EST))+3.D0) + IERLOG=MIN0(IERLOG,75) + 5 CONTINUE + IAG=IDINT(EXI*20.D0+1.D0) + IAG = MIN0( IAG, 20) + IAGB = IAG+1 + KADT(IAG)=KADT(IAG)+1 + IG=2 + COSSIN=CY(IV)/SX(IV) + COSSIN=DMIN1(COSSIN,1.D0) + COSSIN=DMAX1(COSSIN,-1.D0) + IF(COSSIN.GT.1.0D0) WRITE(99,'(A50)')' nach COSSIN' + AC=DACOS(COSSIN) + IG=IDINT(DAW*AC+2.D0) + IGG=IDINT(DGIK*AC+1.D0) + IF(IGG.GT.NGIK) IGG=NGIK + MEATL(IERLOG,IAG) = MEATL(IERLOG,IAG)+1 + MEATL(IERLOG,21) = MEATL(IERLOG,21)+1 + MEATL(75,IAG) = MEATL(75,IAG)+1 + MAGT(IG,IAGB) = MAGT(IG,IAGB)+1 + MAGT(NG,IAGB) = MAGT(NG,IAGB)+1 + MAGT(IG,22) = MAGT(IG,22)+1 + MEAT(IE,IAGB) = MEAT(IE,IAGB)+1 + MEAT(NE,IAGB) = MEAT(NE,IAGB)+1 + MEAT(IE,22) = MEAT(IE,22)+1 + MEAGT(IE,IGG,IAGB) = MEAGT(IE,IGG,IAGB)+1 + MEAGT(102,IGG,IAGB) = MEAGT(102,IGG,IAGB)+1 + MEAGT(IE,IGG,22) = MEAGT(IE,IGG,22)+1 + MEAGT(102,IGG,22) = MEAGT(102,IGG,22)+1 + MEPT(IE,IPT) = MEPT(IE,IPT)+1 + MEPT(NE,IPT) = MEPT(NE,IPT)+1 + MEPT(IE,102) = MEPT(IE,102)+1 + EMAT(IG,IAGB) = EMAT(IG,IAGB)+ES + EMAT(IG,22) = EMAT(IG,22)+ES + EMAT(NG,IAGB) = EMAT(NG,IAGB)+ES + GO TO 110 C C PROJECTILE IS REFLECTED BACK INTO THE TARGET BY THE SURF. BARRIER C - 517 X(IV)=SUT - COSX(IV)=-1.D0*COSX(IV) - KIT=KIT+1 - GO TO 125 -C - 110 IF(IH.EQ.NH) GO TO 130 -C - IH=IH+1 - IF(E0.GE.0.D0) GO TO 702 - IF(ALPHA.LT.0.D0) GO TO 703 + 517 X(IV)=SUT + COSX(IV)=-1.D0*COSX(IV) + KIT=KIT+1 + GO TO 125 + 110 IF(IH.EQ.NH) GO TO 130 + IH=IH+1 + IF(E0.GE.0.D0) GO TO 702 + IF(ALPHA.LT.0.D0) GO TO 703 C C MAXWELLIAN VELOCITY DISTRIBUTION C - CALL VELOC(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV)) + CALL VELOC(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV)) EMX = EMX+E(IV) - ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) - me(ne) = me(ne)+1 - GO TO 707 + ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) + me(ne) = me(ne)+1 + GO TO 707 C C MAXWELLIAN ENERGY DISTRIBUTION C - 703 CALL ENERG(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV)) -CC703 CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,1) + 703 CALL ENERG(E(IV),COSX(IV),COSY(IV),COSZ(IV),SINE(IV)) EMX = EMX+E(IV) -CC WRITE(6,*) E(IV) - GO TO 707 -C - 702 IF (EQUAL(Esig,0.D0)) THEN + GO TO 707 + 702 IF (EQUAL(Esig,0.D0)) THEN C FIXED PROJECTILE ENERGY -C WRITE(*,*)' Da Esig=0 ist E=E0' - E(IV)=E0 + E(IV)=E0 C GAUSSIAN ENERGY DISTRIBUTION - ELSE -7020 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) - tryE = tryE+1 - IF (Epar.LE.0.D0) THEN - negE = negE+1 - GO TO 7020 - ENDIF - E(IV)=Epar -C WRITE(*,*)E(IV),Epar,E0 - ENDIF -C - TAUPSI(IV)=0.D0 -C -C IF(ALPHA.EQ.-2.) GO TO 705 -C IF(ALPHA.EQ.-1.) GO TO 706 - IF(EQUAL(ALPHA,-2.D0)) GO TO 705 - IF(EQUAL(ALPHA,-1.D0)) GO TO 706 -C - IF(EQUAL(ALPHASIG,0.D0))THEN + ELSE + 7020 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) + tryE = tryE+1 + IF (Epar.LE.0.D0) THEN + negE = negE+1 + GO TO 7020 + ENDIF + E(IV)=Epar + ENDIF + TAUPSI(IV)=0.D0 + IF(EQUAL(ALPHA,-2.D0)) GO TO 705 + IF(EQUAL(ALPHA,-1.D0)) GO TO 706 + IF(EQUAL(ALPHASIG,0.D0))THEN C FIXED PROJECTILE ANGLE -C WRITE(88,*)ALPHA,CALFA,SALFA - COSX(IV)=CALFA - COSY(IV)=SALFA - COSZ(IV)=0.D0 - SINE(IV)=COSY(IV) - ELSE + COSX(IV)=CALFA + COSY(IV)=SALFA + COSZ(IV)=0.D0 + SINE(IV)=COSY(IV) + ELSE C C 1D-GAUSSIAN DISTRIBUTION PROJECTILE ANGLE C - CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, - + CALFA,SALFA,BW) -C WRITE(88,'(5(F12.5))')ALPHA,ALPHASIG,ALPHApar,CALFA,SALFA - COSX(IV) = CALFA - COSY(IV) = SALFA - COSZ(IV) = 0.D0 - SINE(IV) = COSY(IV) - ENDIF -C - GO TO 707 + CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, CALFA + & ,SALFA,BW) + COSX(IV) = CALFA + COSY(IV) = SALFA + COSZ(IV) = 0.D0 + SINE(IV) = COSY(IV) + ENDIF + GO TO 707 C C COSINE ANGLE DISTRIBUTION C -CC705 RPHI=PI2*RANF() -CC705 RPHI=PI2*DRAND48() -CC705 RPHI=PI2*DBLE(RAN(ISEED)) -705 call ranlux(ran2, 2) - RPHI=PI2*DBLE(ran2(1)) -CC RTHETA=RANF() -CC RTHETA=DRAND48() -CC RTHETA=DBLE(RAN(ISEED)) - RTHETA=DBLE(ran2(1)) - COSX(IV)=DSQRT(RTHETA) - SINE(IV)=DSQRT(1.D0-RTHETA) - COSY(IV)=SINE(IV)*DCOS(RPHI) - COSZ(IV)=SINE(IV)*DSIN(RPHI) - GO TO 707 + 705 call ranlux(ran2, 2) + RPHI=PI2*DBLE(ran2(1)) + RTHETA=DBLE(ran2(1)) + COSX(IV)=DSQRT(RTHETA) + SINE(IV)=DSQRT(1.D0-RTHETA) + COSY(IV)=SINE(IV)*DCOS(RPHI) + COSZ(IV)=SINE(IV)*DSIN(RPHI) + GO TO 707 C C RANDOM DISTRIBUTION C - 706 IF(X0.GT.0.D0) GO TO 709 -C -CC RPHI=PI2*RANF() -CC RPHI=PI2*DRAND48() -CC RPHI=PI2*DBLE(RAN(ISEED)) - call ranlux(ran2, 2) - RPHI=PI2*DBLE(ran2(1)) -CC RTHETA=RANF() -CC RTHETA=DRAND48() -cc RTHETA=DBLE(RAN(ISEED)) - RTHETA=DBLE(ran2(2)) - COSX(IV)=1.D0-RTHETA - SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) - COSY(IV)=SINE(IV)*DSIN(RPHI) - COSZ(IV)=SINE(IV)*DCOS(RPHI) - GO TO 707 -C -CC709 RPHI=PI2*RANF() -CC709 RPHI=PI2*DRAND48() -CC709 RPHI=PI2*DBLE(RAN(ISEED)) - 709 call ranlux(ran2, 2) - RPHI=PI2*DBLE(ran2(1)) -CC RTHETA=RANF() -CC RTHETA=DRAND48() -CC RTHETA=DBLE(RAN(ISEED)) - RTHETA=DBLE(ran2(2)) - COSX(IV)=1.D0-2.D0*RTHETA - SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) - COSY(IV)=SINE(IV)*DSIN(RPHI) - COSZ(IV)=SINE(IV)*DCOS(RPHI) - GO TO 708 -C - 707 IF(X0.GT.0.D0) GO TO 708 + 706 IF(X0.GT.0.D0) GO TO 709 + call ranlux(ran2, 2) + RPHI=PI2*DBLE(ran2(1)) + RTHETA=DBLE(ran2(2)) + COSX(IV)=1.D0-RTHETA + SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) + COSY(IV)=SINE(IV)*DSIN(RPHI) + COSZ(IV)=SINE(IV)*DCOS(RPHI) + GO TO 707 + 709 call ranlux(ran2, 2) + RPHI=PI2*DBLE(ran2(1)) + RTHETA=DBLE(ran2(2)) + COSX(IV)=1.D0-2.D0*RTHETA + SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) + COSY(IV)=SINE(IV)*DSIN(RPHI) + COSZ(IV)=SINE(IV)*DCOS(RPHI) + GO TO 708 + 707 IF(X0.GT.0.D0) GO TO 708 C C EXTERNAL START C - SINA=SINE(IV) - COSX(IV)=DSQRT((E(IV)*COSX(IV)*COSX(IV)+ESB)/(E(IV)+ESB)) - SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) - COSY(IV)=COSY(IV)*SINE(IV)/SINA - COSZ(IV)=COSZ(IV)*SINE(IV)/SINA - E(IV)=E(IV)+ESB + SINA=SINE(IV) + COSX(IV)=DSQRT((E(IV)*COSX(IV)*COSX(IV)+ESB)/(E(IV)+ESB)) + SINE(IV)=DSQRT(1.D0-COSX(IV)*COSX(IV)) + COSY(IV)=COSY(IV)*SINE(IV)/SINA + COSZ(IV)=COSZ(IV)*SINE(IV)/SINA + E(IV)=E(IV)+ESB C C LOCUS OF FIRST COLLISION C - 708 LLL(IV)=ISRCHFGT(L,XX(1),1,X0) -CC RA1=CVMGT(RANF(),1.,X0.LE.0.) -CC RA1=CVMGT(DRAND48(),1.,X0.LE.0.) -CC RA1=CVMGT(DBLE(RAN(ISEED)),1.D0,X0.LE.0.D0) - call ranlux(random, 1) - RA1=CVMGT(DBLE(random),1.D0,X0.LE.0.D0) - X(IV)=XC+LM(LLL(IV))*RA1*COSX(IV) - Y(IV)=LM(LLL(IV))*RA1*COSY(IV) - Z(IV)=LM(LLL(IV))*RA1*COSZ(IV) - PL(IV)=CVMGT(0.D0,LM(LLL(IV))*RA1,X0.LE.0.D0) - GO TO 120 + 708 LLL(IV)=ISRCHFGT(L,XX(1),1,X0) + call ranlux(random, 1) + RA1=CVMGT(DBLE(random),1.D0,X0.LE.0.D0) + X(IV)=XC+LM(LLL(IV))*RA1*COSX(IV) + Y(IV)=LM(LLL(IV))*RA1*COSY(IV) + Z(IV)=LM(LLL(IV))*RA1*COSZ(IV) + PL(IV)=CVMGT(0.D0,LM(LLL(IV))*RA1,X0.LE.0.D0) + GO TO 120 C C COUNTING DOWN IH1 , ONLY LESS THAN (NH-IH) HAVE TO BE PROCESSED -C - 130 CONTINUE - E(IV)=E(IH1) - COSX(IV)=COSX(IH1) - COSY(IV)=COSY(IH1) - COSZ(IV)=COSZ(IH1) - SINE(IV)=SINE(IH1) - X(IV)=X(IH1) - Y(IV)=Y(IH1) - Z(IV)=Z(IH1) - PL(IV)=PL(IH1) - TAU(IV)=TAU(IH1) - TAUPSI(IV)=TAUPSI(IH1) - CPSI(IV)=CPSI(IH1) - ENUCL(IV)=ENUCL(IH1) - EINEL(IV)=EINEL(IH1) - IH1=IH1-1 - IF(IV.EQ.IH1+1) GO TO 90 - IF(IH1+1.GT.IVMAX) GO TO 125 - GO TO 160 - 125 CONTINUE - X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) - Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) - Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) - PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) - * ,X(IV).LE.0.D0) - LLL(IV)=MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) - 120 CONTINUE - 90 CONTINUE +C + 130 CONTINUE + E(IV)=E(IH1) + COSX(IV)=COSX(IH1) + COSY(IV)=COSY(IH1) + COSZ(IV)=COSZ(IH1) + SINE(IV)=SINE(IH1) + X(IV)=X(IH1) + Y(IV)=Y(IH1) + Z(IV)=Z(IH1) + PL(IV)=PL(IH1) + TAU(IV)=TAU(IH1) + TAUPSI(IV)=TAUPSI(IH1) + CPSI(IV)=CPSI(IH1) + ENUCL(IV)=ENUCL(IH1) + EINEL(IV)=EINEL(IH1) + IH1=IH1-1 + IF(IV.EQ.IH1+1) GO TO 90 + IF(IH1+1.GT.IVMAX) GO TO 125 + GO TO 160 + 125 CONTINUE + X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) + Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) + Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) + PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0 + & .D0) + LLL(IV)=MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) + 120 CONTINUE + 90 CONTINUE C C INCREMENT OF PROJECTILE ENERGY AND POSITION C OF PARTICLES NOT HANDLED IN LOOP 120 C -CC IF(IVMIN.LE.1) GO TO 134 DO 128 IV=1,IVMIN-1 - LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) - 128 CONTINUE - DO 129 IV=1,IVMIN-1 - X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) - Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) - Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) - PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) - * ,X(IV).LE.0.D0) - 129 CONTINUE - 134 CONTINUE + LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) + 128 CONTINUE + DO IV=1,IVMIN-1 + X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) + Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) + Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) + PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0 + & .D0) + ENDDO + 134 CONTINUE C IF(IVMAX.LT.IVMIN) GO TO 132 -CC IF(IVMAX.LT.IH1) GO TO 132 - DO 133 IV=IVMAX+1,IH1 - LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) - 133 CONTINUE - DO 131 IV=IVMAX+1,IH1 - X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) - Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) - Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) - PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) - * ,X(IV).LE.0.D0) - 131 CONTINUE - 132 CONTINUE -C + DO IV=IVMAX+1,IH1 + LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) + ENDDO + DO IV=IVMAX+1,IH1 + X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV) + Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV) + Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV) + PL(IV)=CVMGT(PL(IV),PL(IV)+LM(LLL(IV))+TAUPSI(IV) ,X(IV).LE.0 + & .D0) + ENDDO + 132 CONTINUE GO TO 1 -C - 140 IF(NREC1+NREC2.GT.0) GO TO 83 + 140 IF(NREC1+NREC2.GT.0) GO TO 83 C C C PRINTOUT C C do ima = 5000,1,-1 - if(me(ima).ne.0) goto 1010 + if(me(ima).ne.0) goto 1010 enddo ima = 1 1010 ima = MIN0(ima+2,5000) open(20,file='edist') do ne=1,ima - write(20,1020) ne,me(ne) + write(20,1020) ne,me(ne) enddo 1020 format(1x,2i6) close(20) -c -c Berechnung der part. reflec. coeff. nach Thomas et al. -c - E0keV=E0/1.D3 - EsigkeV=Esig/1.D3 -c - IF(ZT(1,2).LT.1.0D-3) THEN - epsilon = 32.55D0*(MT(1,1)/M1)/(1.D0+(MT(1,1)/M1))* - 1 1.D0/(Z1*ZT(1,1)*DSQRT(Z1**(2.D0/3.D0)+ZT(1,1)**(2.D0/3.D0)))* - 2 E0keV -cTR 1 1.D0/(Z1*ZT(1,1)*DSQRT(Z1**2.D0/3.D0+ZT(1,1)**2.D0/3.D0))* -cTR 2 E0keV - prcoeff = prc(1)*DLOG(prc(2)*epsilon+DEXP(1.D0))/ - 1 (1.D0+(PRC(3)*epsilon**PRC(4))+(PRC(5)*epsilon**PRC(6))) - ELSE - epsilon = 0.D0 - prcoeff = 0.D0 +C +C Berechnung der part. reflec. coeff. nach Thomas et al. +C + E0keV=E0/1.D3 + EsigkeV=Esig/1.D3 + + IF(ZT(1,2).LT.1.0D-3) THEN + epsilon = 32.55D0*(MT(1,1)/M1)/(1.D0+(MT(1,1)/M1))* 1.D0/(Z1 + & *ZT(1,1)*DSQRT(Z1**(2.D0/3.D0)+ZT(1,1)**(2.D0/3.D0))) + & * E0keV + prcoeff = prc(1)*DLOG(prc(2)*epsilon+DEXP(1.D0))/ (1.D0+(PRC(3) + & *epsilon**PRC(4))+(PRC(5)*epsilon**PRC(6))) + ELSE + epsilon = 0.D0 + prcoeff = 0.D0 ENDIF C 2nd CALL DATE_AND_TIME CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3), - 1 Date_Time) -C - IF(Date_Time(2).EQ.1) THEN - month_stop='Jan.' - days_total_stop=Date_Time(3) - ELSEIF(Date_Time(2).EQ.2) THEN - month_stop='Feb.' - days_stop_total=31+Date_Time(3) - ELSEIF(Date_Time(2).EQ.3) THEN - month_stop='Mar.' - days_stop_total=31+28+Date_Time(3) - ELSEIF(Date_Time(2).EQ.4) THEN - month_stop='Apr.' - days_stop_total=31+28+31+Date_Time(3) - ELSEIF(Date_Time(2).EQ.5) THEN - month_stop='May ' - days_stop_total=31+28+31+30+Date_Time(3) - ELSEIF(Date_Time(2).EQ.6) THEN - month_stop='Jun.' - days_stop_total=31+28+31+30+31+Date_Time(3) - ELSEIF(Date_Time(2).EQ.7) THEN - month_stop='Jul.' - days_stop_total=31+28+31+30+31+30+Date_Time(3) - ELSEIF(Date_Time(2).EQ.8) THEN - month_stop='Aug.' - days_stop_total=31+28+31+30+31+30+31+Date_Time(3) - ELSEIF(Date_Time(2).EQ.9) THEN - month_stop='Sep.' - days_stop_total=31+28+31+30+31+30+31+31+Date_Time(3) - ELSEIF(Date_Time(2).EQ.10) THEN - month_stop='Oct.' - days_stop_total=31+28+31+30+31+30+31+31+30+Date_Time(3) - ELSEIF(Date_Time(2).EQ.11) THEN - month_stop='Nov.' - days_stop_total=31+28+31+30+31+30+31+31+30+31+Date_Time(3) - ELSE - month_stop='Dec.' - days_stop_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3) - ENDIF -C - READ(Real_Clock(1)(1:4),'(A4)')year_stop - READ(Real_Clock(1)(7:8),'(A2)')day_stop - READ(Real_Clock(2)(1:2),'(A2)')hour_stop - READ(Real_Clock(2)(3:4),'(A2)')min_stop - READ(Real_Clock(2)(5:6),'(A2)')sec_stop + & Date_Time) + IF(Date_Time(2).EQ.1) THEN + month_stop='Jan.' + days_total_stop=Date_Time(3) + ELSEIF(Date_Time(2).EQ.2) THEN + month_stop='Feb.' + days_stop_total=31+Date_Time(3) + ELSEIF(Date_Time(2).EQ.3) THEN + month_stop='Mar.' + days_stop_total=31+28+Date_Time(3) + ELSEIF(Date_Time(2).EQ.4) THEN + month_stop='Apr.' + days_stop_total=31+28+31+Date_Time(3) + ELSEIF(Date_Time(2).EQ.5) THEN + month_stop='May ' + days_stop_total=31+28+31+30+Date_Time(3) + ELSEIF(Date_Time(2).EQ.6) THEN + month_stop='Jun.' + days_stop_total=31+28+31+30+31+Date_Time(3) + ELSEIF(Date_Time(2).EQ.7) THEN + month_stop='Jul.' + days_stop_total=31+28+31+30+31+30+Date_Time(3) + ELSEIF(Date_Time(2).EQ.8) THEN + month_stop='Aug.' + days_stop_total=31+28+31+30+31+30+31+Date_Time(3) + ELSEIF(Date_Time(2).EQ.9) THEN + month_stop='Sep.' + days_stop_total=31+28+31+30+31+30+31+31+Date_Time(3) + ELSEIF(Date_Time(2).EQ.10) THEN + month_stop='Oct.' + days_stop_total=31+28+31+30+31+30+31+31+30+Date_Time(3) + ELSEIF(Date_Time(2).EQ.11) THEN + month_stop='Nov.' + days_stop_total=31+28+31+30+31+30+31+31+30+31+Date_Time(3) + ELSE + month_stop='Dec.' + days_stop_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3) + ENDIF + + READ(Real_Clock(1)(1:4),'(A4)')year_stop + READ(Real_Clock(1)(7:8),'(A2)')day_stop + READ(Real_Clock(2)(1:2),'(A2)')hour_stop + READ(Real_Clock(2)(3:4),'(A2)')min_stop + READ(Real_Clock(2)(5:6),'(A2)')sec_stop C C how many seconds are needed for the simulation ?? C - seconds_stop_total=Date_Time(7)+(Date_Time(6)*60)+ - 1 (Date_Time(5)*3600)+(days_stop_total-1)*86400 -C - WRITE(21,*) - WRITE(21,10051)day_stop,month_stop,year_stop, - 1 hour_stop,min_stop,sec_stop -10051 FORMAT(1x,' TrimSp simulation ended at: ',A2,'.',A4,1x,A4, - 1 1x,A2,':',A2,':',A2) - WRITE(21,*) - WRITE(21,10052)nh,(seconds_stop_total-seconds_start_total) + seconds_stop_total=Date_Time(7)+(Date_Time(6)*60)+ (Date_Time(5) + & *3600)+(days_stop_total-1)*86400 + + WRITE(21,*) + WRITE(21,10051)day_stop,month_stop,year_stop, hour_stop,min_stop + & ,sec_stop +10051 FORMAT(1x,' TrimSp simulation ended at: ',A2,'.',A4,1x,A4, 1x,A2 + & ,':',A2,':',A2) + WRITE(21,*) + WRITE(21,10052)nh,(seconds_stop_total-seconds_start_total) 10052 FORMAT(1x,' Simulation needed for ',I7,' muons ',I7,' seconds') -C + WRITE(21,1402)innam 1402 FORMAT(//30X,'* INPUT DATA *',5X,A12) WRITE(21,1404) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC 1404 FORMAT(//,7X,2HZ1,8X,2HM1,10X,2HE0,6X,4HEsig,7X,5HALPHA,7X - 1 ,8HALPHASIG,7X,2HEF,7X - 2 ,3HESB,6X,6HSHEATH,5X,3HERC/2F10.2,1F13.2,7F10.2) + & ,8HALPHASIG,7X,2HEF,7X ,3HESB,6X,6HSHEATH,5X,3HERC/2F10.2 + & ,1F13.2,7F10.2) WRITE(21,1406) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT - 1 ,IPOTR,IRL + & ,IPOTR,IRL 1406 FORMAT(/7X,2HNH,8X,2HRI,5X,3HRI2,5X,3HRI3,11X,2HX0,8X,2HRD,8X,2HCW - 1 ,8X,2HCA - 1 ,7X,3HKK0,3X,4HKK0R,3X,5HKDEE1,2X,5HKDEE2,2X,4HIPOT,3X,5HIPOTR - 2 ,3X,3HIRL/I10,3F10.2,1F13.2,3F10.2,1X,7I7) + & ,8X,2HCA ,7X,3HKK0,3X,4HKK0R,3X,5HKDEE1,2X,5HKDEE2,2X,4HIPOT + & ,3X,5HIPOTR ,3X,3HIRL/I10,3F10.2,1F13.2,3F10.2,1X,7I7) WRITE(21,1408) 1408 FORMAT(//13X,2HDX,6X,3HRHO,4X,2HCK,2X - 1 ,5HZ(,1),1X,5HZ(,2),1X,5HZ(,3),1X,5HZ(,4),1X,5HZ(,5),2X - 2 ,5HM(,1),2X,5HM(,2),2X,5HM(,3),2X,5HM(,4),2X,5HM(,5),1X - 3 ,5HC(,1),1X,5HC(,2),1X,5HC(,3),1X,5HC(,4),1X,5HC(,5)) - DO 1410 I=1,L - WRITE(21,1412) I,DX(I),RHO(I),CK(I),(ZT(I,J),J=1,5) - 1 ,(MT(I,J),J=1,5),(CO(I,J),J=1,5) - 1412 FORMAT(/1X,I1,6H.LAYER,1X,1F8.2,2F7.2,5F6.0,5F7.2,5F6.3) - 1410 CONTINUE + & ,5HZ(,1),1X,5HZ(,2),1X,5HZ(,3),1X,5HZ(,4),1X,5HZ(,5),2X + & ,5HM(,1),2X,5HM(,2),2X,5HM(,3),2X,5HM(,4),2X,5HM(,5),1X + & ,5HC(,1),1X,5HC(,2),1X,5HC(,3),1X,5HC(,4),1X,5HC(,5)) + DO I=1,L + WRITE(21,1412) I,DX(I),RHO(I),CK(I),(ZT(I,J),J=1,5) ,(MT(I,J),J + & =1,5),(CO(I,J),J=1,5) + 1412 FORMAT(/1X,I1,6H.LAYER,1X,1F8.2,2F7.2,5F6.0,5F7.2,5F6.3) + ENDDO WRITE(21,1414) 1414 FORMAT(//27X,'***',2X,'SBE(LAYER,ELEMENT)',2X,'***',5X - 1 ,'***',5X,'ED(LAYER,ELEMENT)',5X,'***',5X - 2 ,'***',3X,'BE(LAYER,ELEMENT)',2X,'***') + & ,'***',5X,'ED(LAYER,ELEMENT)',5X,'***',5X + & ,'***',3X,'BE(LAYER,ELEMENT)',2X,'***') DO 1416 I=1,L - WRITE(21,1418) I,(SBE(I,J),J=1,5),(ED(I,J),J=1,5),(BE(I,J),J=1,5) - 1418 FORMAT(/1X,I1,6H.LAYER,17X,5F6.2,3X,5F7.2,3X,5F6.2) + WRITE(21,1418) I,(SBE(I,J),J=1,5),(ED(I,J),J=1,5),(BE(I,J),J=1 + & ,5) + 1418 FORMAT(/1X,I1,6H.LAYER,17X,5F6.2,3X,5F7.2,3X,5F6.2) 1416 CONTINUE IF(KDEE1.LT.4) GO TO 1421 WRITE(21,1419) 1419 FORMAT(//30X,'CH1',10X,'CH2',10X,'CH3',10X,'CH4',10X,'CH5') DO 1417 I=1,L - WRITE(21,1415) I,CH1(I,1),CH2(I,1),CH3(I,1),CH4(I,1),CH5(I,1) - IF(NJ(I).LT.2) GO TO 1417 - WRITE(21,1423) CH1(I,2),CH2(I,2),CH3(I,2),CH4(I,2),CH5(I,2) - IF(NJ(I).LT.3) GO TO 1417 - WRITE(21,1423) CH1(I,3),CH2(I,3),CH3(I,3),CH4(I,3),CH5(I,3) - IF(NJ(I).LT.4) GO TO 1417 - WRITE(21,1423) CH1(I,4),CH2(I,4),CH3(I,4),CH4(I,4),CH5(I,4) - IF(NJ(I).LT.5) GO TO 1417 - WRITE(21,1423) CH1(I,5),CH2(I,5),CH3(I,5),CH4(I,5),CH5(I,5) + WRITE(21,1415) I,CH1(I,1),CH2(I,1),CH3(I,1),CH4(I,1),CH5(I,1) + IF(NJ(I).LT.2) GO TO 1417 + WRITE(21,1423) CH1(I,2),CH2(I,2),CH3(I,2),CH4(I,2),CH5(I,2) + IF(NJ(I).LT.3) GO TO 1417 + WRITE(21,1423) CH1(I,3),CH2(I,3),CH3(I,3),CH4(I,3),CH5(I,3) + IF(NJ(I).LT.4) GO TO 1417 + WRITE(21,1423) CH1(I,4),CH2(I,4),CH3(I,4),CH4(I,4),CH5(I,4) + IF(NJ(I).LT.5) GO TO 1417 + WRITE(21,1423) CH1(I,5),CH2(I,5),CH3(I,5),CH4(I,5),CH5(I,5) 1417 CONTINUE 1415 FORMAT(/1X,I1,6H.LAYER,17X,5F13.6) 1423 FORMAT(/25X,5F13.6) @@ -2703,8 +2612,8 @@ C IF(IPOTR.EQ.2) DPOTR='MOLIERE POTENTIAL' IF(IPOTR.EQ.3) DPOTR='ZBL POTENTIAL' WRITE(21,1411) DPOT,DPOTR - 1411 FORMAT(//7X,'INTERACTION POTENTIAL : PROJECTILE-TARGET : ',A18,' - 1 TARGET-TARGET : ',A18) + 1411 FORMAT(//7X,'INTERACTION POTENTIAL : PROJECTILE-TARGET : ',A18 + & ,' TARGET-TARGET : ',A18) IF(KDEE1.EQ.1) DKDEE1='LINDHARD-SCHARFF' IF(KDEE1.EQ.2) DKDEE1='OEN-ROBINSON' IF(KDEE1.EQ.3) DKDEE1='50% LS 50% OR' @@ -2714,38 +2623,38 @@ C IF(KDEE2.EQ.2) DKDEE2='OEN-ROBINSON' IF(KDEE2.EQ.3) DKDEE2='50% LS 50% OR' WRITE(21,1413) DKDEE1,DKDEE2 - 1413 FORMAT(//7X,'INELASTIC LOSS MODEL : PROJECTILE-TARGET : ',A18,' - 1 TARGET-TARGET : ',A18) + 1413 FORMAT(//7X,'INELASTIC LOSS MODEL : PROJECTILE-TARGET : ',A18 + & ,' TARGET-TARGET : ',A18) IF(E0.GT.0.D0) GO TO 1420 IF(ALPHA.LT.0.D0) GO TO 1405 WRITE(21,1422) TI,ZARG,VELC,EMX - 1422 FORMAT(//6X,'MAXWELLIAN DISTRIBUTION',7X,2HTI,5X,4HZARG,5X - 1 ,4HVELC,8X,3HEMX/29X,1F10.2,2F9.4,1E14.6) + 1422 FORMAT(//6X,'MAXWELLIAN DISTRIBUTION',7X,2HTI,5X,4HZARG,5X ,4HVELC + & ,8X,3HEMX/29X,1F10.2,2F9.4,1E14.6) GO TO 1427 1405 ALPHAM=-ALPHA WRITE(21,1407) TI,SHEATH,ALPHAM,EMX 1407 FORMAT(//6X,'MAXWELLIAN DISTRIBUTION (ENERGY)',7X,'TI',5X - 1 ,'SHEATH',5X,'ALPHAM',8X,'EMX'/38X,3F10.2,2X,1E14.6) + & ,'SHEATH',5X,'ALPHAM',8X,'EMX'/38X,3F10.2,2X,1E14.6) GO TO 1427 1420 IF(ALPHA.EQ.-1.) WRITE(21,1424) 1424 FORMAT(//6X,'RANDOM DISTRIBUTION'/) IF(ALPHA.EQ.-2.) WRITE(21,1426) 1426 FORMAT(//6X,'COSINE DISTRIBUTION'/) 1427 CONTINUE - IF(EQUAL(Esig,0.D0)) THEN - WRITE(21,14271) -14271 FORMAT(//6X,'fixed PROJECTILE ENERGY'/) - ELSE - WRITE(21,14272) -14272 FORMAT(//6X,'PROJECTILE ENERGY has GAUSSIAN DISTRIBUTION '/) - ENDIF - IF(EQUAL(ALPHASIG,0.D0)) THEN - WRITE(21,14273) -14273 FORMAT(//6X,'fixed PROJECTILE ANGLE'/) - ELSE - WRITE(21,14274) -14274 FORMAT(//6X,'PROJECTILE ANGLE has 1D GAUSSIAN DISTRIBUTION '/) - ENDIF + IF(EQUAL(Esig,0.D0)) THEN + WRITE(21,14271) +14271 FORMAT(//6X,'fixed PROJECTILE ENERGY'/) + ELSE + WRITE(21,14272) +14272 FORMAT(//6X,'PROJECTILE ENERGY has GAUSSIAN DISTRIBUTION '/) + ENDIF + IF(EQUAL(ALPHASIG,0.D0)) THEN + WRITE(21,14273) +14273 FORMAT(//6X,'fixed PROJECTILE ANGLE'/) + ELSE + WRITE(21,14274) +14274 FORMAT(//6X,'PROJECTILE ANGLE has 1D GAUSSIAN DISTRIBUTION '/) + ENDIF WRITE(21,1428)outnam 1428 FORMAT(1H1,//30X,'* OUTPUT DATA *',5X,A12) #if defined (OS_WIN) @@ -2754,102 +2663,72 @@ C #endif WRITE(21,1430) HLM,HLMT,SU,SUT,XC,RT,SFE,INEL,L,LJ 1430 FORMAT(//17X,'HLM',7X,'HLMT',8X,'SU',7X,'SUT',8X,'XC',8X,'RT',7X - 1 ,'SFE',6X,'INEL',9X,'L',8X,'LJ'/ - 2 10X,1F11.4,1F10.3,1F10.4,1F10.3,1F10.4,1F10.3,1F10.2,3I10) + & ,'SFE',6X,'INEL',9X,'L',8X,'LJ'/ + & 10X,1F11.4,1F10.3,1F10.4,1F10.3,1F10.4,1F10.3,1F10.2,3I10) WRITE(21,1432) NPROJ,KIB,KIT,MAXA,NALL,NPA,NSA,KIS,KIST 1432 FORMAT(//16X,'NPROJ',7X,'KIB',7X,'KIT',6X,'MAXA',6X,'NALL',7X - 1 ,'NPA',7X,'NSA',7X,'KIS',6X,'KIST'/11X,9I10) + & ,'NPA',7X,'NSA',7X,'KIS',6X,'KIST'/11X,9I10) WRITE(21,470) - 470 FORMAT(//12X,'EPS0(I)',7X,'Z2(I)',7X,'M2(I)',5X,'ARHO(I)' - 1 ,7X,'LM(I)',5X,'PDMAX(I)',5X,'ASIG(I)',7X,'SB(I)',7X,'XX(I)' - 2 ,8X,'NJ(I)') - DO 472 I=1,L - WRITE(21,473) I,EPS0(I),Z2(I),M2(I),ARHO(I),LM(I),PDMAX(I),ASIG(I) - 1 ,SB(I),XX(I),NJ(I) - 473 FORMAT(/1X,I1,6H.LAYER,1X,9E12.4,I10) - 472 CONTINUE + 470 FORMAT(//12X,'EPS0(I)',7X,'Z2(I)',7X,'M2(I)',5X,'ARHO(I)' ,7X + & ,'LM(I)',5X,'PDMAX(I)',5X,'ASIG(I)',7X,'SB(I)',7X,'XX(I)' ,8X + & ,'NJ(I)') + DO I=1,L + WRITE(21,473) I,EPS0(I),Z2(I),M2(I),ARHO(I),LM(I),PDMAX(I) + & ,ASIG(I),SB(I),XX(I),NJ(I) + 473 FORMAT(/1X,I1,6H.LAYER,1X,9E12.4,I10) + ENDDO WRITE(21,474) - 474 FORMAT(//13X, - 1 'A1(1)',3X,'A1(2)',3X,'A1(3)',3X,'A1(4)',3X,'A1(5)',3X, - 1 'A1(6)',3X,'A1(7)',3X,'A1(8)',3X,'A1(9)',2X,'A1(10)',2X, - 1 'A1(11)',2X,'A1(12)',2X,'A1(13)',2X,'A1(14)',2X,'A1(15)',2X, - 1 'A1(16)',2X,'A1(17)',2X,'A1(18)',2X,'A1(19)',2X,'A1(20)',2X, - 1 'A1(21)',2X,'A1(22)',2X,'A1(23)',2X,'A1(24)',2X,'A1(25)',2X, - 1 'A1(26)',2X,'A1(27)',2X,'A1(28)',2X,'A1(29)',2X,'A1(30)',2X, - 1 'A1(31)',2X,'A1(32)',2X,'A1(33)',2X,'A1(34)',2X,'A1(35)') -cTR 1 'A1(21)',2X,'A1(22)',2X,'A1(23)',2X,'A1(24)',2X,'A1(25)',2X, -cTR 1 'A1(26)',2X,'A1(27)',2X,'A1(28)',2X,'A1(29)',2X,'A1(30)',2X, -cTR 1 'A1(31)',2X,'A1(32)',2X,'A1(33)',2X,'A1(34)',2X,'A1(35)') + 474 FORMAT(//13X, + & 'A1(1)',3X,'A1(2)',3X,'A1(3)',3X,'A1(4)',3X,'A1(5)',3X, + & 'A1(6)',3X,'A1(7)',3X,'A1(8)',3X,'A1(9)',2X,'A1(10)',2X, + & 'A1(11)',2X,'A1(12)',2X,'A1(13)',2X,'A1(14)',2X,'A1(15)',2X, + & 'A1(16)',2X,'A1(17)',2X,'A1(18)',2X,'A1(19)',2X,'A1(20)',2X, + & 'A1(21)',2X,'A1(22)',2X,'A1(23)',2X,'A1(24)',2X,'A1(25)',2X, + & 'A1(26)',2X,'A1(27)',2X,'A1(28)',2X,'A1(29)',2X,'A1(30)',2X, + & 'A1(31)',2X,'A1(32)',2X,'A1(33)',2X,'A1(34)',2X,'A1(35)') DO 475 I=1,LJ - WRITE(21,471) A1(I) - 471 FORMAT(/1X,9X,35F8.5) - 475 CONTINUE + WRITE(21,471) A1(I) + 471 FORMAT(/1X,9X,35F8.5) + 475 CONTINUE WRITE(21,484) - 484 FORMAT(//11X, - 1 'KOR1(1)',1X,'KOR1(2)',1X,'KOR1(3)',1X,'KOR1(4)',1X,'KOR1(5)', - 2 1X,'KOR1(6)',1X,'KOR1(7)',1X,'KOR1(8)',1X,'KOR1(9)',1X,'KOR1(A)', - 3 1X,'KOR1(B)',1X,'KOR1(C)',1X,'KOR1(D)',1X,'KOR1(E)',1X,'KOR1(F)', - 4 1X,'KOR1(G)',1X,'KOR1(H)',1X,'KOR1(I)',1X,'KOR1(J)',1X,'KOR1(K)', - 5 1X,'KOR1(L)',1X,'KOR1(M)',1X,'KOR1(N)',1X,'KOR1(O)',1X,'KOR1(P)', - 6 1X,'KOR1(Q)',1X,'KOR1(R)',1X,'KOR1(S)',1X,'KOR1(T)',1X,'KOR1(U)', - 7 1X,'KOR1(V)',1X,'KOR1(W)',1X,'KOR1(X)',1X,'KOR1(Y)',1X,'KOR1(Z)') -cTR 6 1X,'KOR1(Q)',1X,'KOR1(R)',1X,'KOR1(S)',1X,'KOR1(T)',1X,'KOR1(U)', -cTR 7 1X,'KOR1(V)',1X,'KOR1(W)',1X,'KOR1(X)',1X,'KOR1(Y)',1X,'KOR1(Z)') + 484 FORMAT(//11X, + & 'KOR1(1)',1X,'KOR1(2)',1X,'KOR1(3)',1X,'KOR1(4)',1X,'KOR1(5)', + & 1X,'KOR1(6)',1X,'KOR1(7)',1X,'KOR1(8)',1X,'KOR1(9)',1X,'KOR1(A)', + & 1X,'KOR1(B)',1X,'KOR1(C)',1X,'KOR1(D)',1X,'KOR1(E)',1X,'KOR1(F)', + & 1X,'KOR1(G)',1X,'KOR1(H)',1X,'KOR1(I)',1X,'KOR1(J)',1X,'KOR1(K)', + & 1X,'KOR1(L)',1X,'KOR1(M)',1X,'KOR1(N)',1X,'KOR1(O)',1X,'KOR1(P)', + & 1X,'KOR1(Q)',1X,'KOR1(R)',1X,'KOR1(S)',1X,'KOR1(T)',1X,'KOR1(U)', + & 1X,'KOR1(V)',1X,'KOR1(W)',1X,'KOR1(X)',1X,'KOR1(Y)',1X,'KOR1(Z)') DO 486 I=1,LJ - WRITE(21,489) KOR1(I) - 489 FORMAT(/1X,9X,35F8.5) - 486 CONTINUE + WRITE(21,489) KOR1(I) + 489 FORMAT(/1X,9X,35F8.5) + 486 CONTINUE WRITE(21,476) - 476 FORMAT(//12X, - 1 'A(I,1)',2X,'A(I,2)',2X,'A(I,3)',2X,'A(I,4)',2X,'A(I,5)',2X, - 2 'A(I,6)',2X,'A(I,7)',2X,'A(I,8)',2X,'A(I,9)',1X,'A(I,10)',1X, - 3 'A(I,11)',1X,'A(I,12)',1X,'A(I,13)',1X,'A(I,14)',1X,'A(I,15)',1X, - 4 'A(I,16)',1X,'A(I,17)',1X,'A(I,18)',1X,'A(I,19)',1X,'A(I,20)',1X, - 5 'A(I,21)',1X,'A(I,22)',1X,'A(I,23)',1X,'A(I,24)',1X,'A(I,25)',1X, - 6 'A(I,26)',1X,'A(I,27)',1X,'A(I,28)',1X,'A(I,29)',1X,'A(I,30)',1X, - 7 'A(I,31)',1X,'A(I,32)',1X,'A(I,33)',1X,'A(I,34)',1X,'A(I,35)') -cTR 4 'A(I,16)',1X,'A(I,17)',1X,'A(I,18)',1X,'A(I,19)',1X,'A(I,20)',1X, -cTR 5 'A(I,21)',1X,'A(I,22)',1X,'A(I,23)',1X,'A(I,24)',1X,'A(I,25)',1X, -cTR 6 'A(I,26)',1X,'A(I,27)',1X,'A(I,28)',1X,'A(I,29)',1X,'A(I,30)',1X, -cTR 7 'A(I,31)',1X,'A(I,32)',1X,'A(I,33)',1X,'A(I,34)',1X,'A(I,35)') - DO 478 I=1,LJ - WRITE(21,477) (A(I,J),J=1,LJ) - 477 FORMAT(/1X,9X,35F8.5) - 478 CONTINUE + 476 FORMAT(//12X, + & 'A(I,1)',2X,'A(I,2)',2X,'A(I,3)',2X,'A(I,4)',2X,'A(I,5)',2X, + & 'A(I,6)',2X,'A(I,7)',2X,'A(I,8)',2X,'A(I,9)',1X,'A(I,10)',1X, + & 'A(I,11)',1X,'A(I,12)',1X,'A(I,13)',1X,'A(I,14)',1X,'A(I,15)',1X, + & 'A(I,16)',1X,'A(I,17)',1X,'A(I,18)',1X,'A(I,19)',1X,'A(I,20)',1X, + & 'A(I,21)',1X,'A(I,22)',1X,'A(I,23)',1X,'A(I,24)',1X,'A(I,25)',1X, + & 'A(I,26)',1X,'A(I,27)',1X,'A(I,28)',1X,'A(I,29)',1X,'A(I,30)',1X, + & 'A(I,31)',1X,'A(I,32)',1X,'A(I,33)',1X,'A(I,34)',1X,'A(I,35)') + DO I=1,LJ + WRITE(21,477) (A(I,J),J=1,LJ) + 477 FORMAT(/1X,9X,35F8.5) + ENDDO WRITE(21,490) 490 FORMAT(//11X, - 1 'KOR(,1)',1X,'KOR(,2)',1X,'KOR(,3)',1X,'KOR(,4)',1X,'KOR(,5)',1X, - 2 'KOR(,6)',1X,'KOR(,7)',1X,'KOR(,8)',1X,'KOR(,9)',1X,'KOR(,A)',1X, - 3 'KOR(,B)',1X,'KOR(,C)',1X,'KOR(,D)',1X,'KOR(,E)',1X,'KOR(,F)',1X, - 4 'KOR(,G)',1X,'KOR(,H)',1X,'KOR(,I)',1X,'KOR(,J)',1X,'KOR(,K)',1X, - 5 'KOR(,L)',1X,'KOR(,M)',1X,'KOR(,N)',1X,'KOR(,O)',1X,'KOR(,P)',1X, - 6 'KOR(,Q)',1X,'KOR(,R)',1X,'KOR(,S)',1X,'KOR(,T)',1X,'KOR(,U)',1X, - 7 'KOR(,V)',1X,'KOR(,W)',1X,'KOR(,X)',1X,'KOR(,Y)',1X,'KOR(,Z)') -cTR 4 'KOR(,G)',1X,'KOR(,H)',1X,'KOR(,I)',1X,'KOR(,J)',1X,'KOR(,K)',1X, -cTR 5 'KOR(,L)',1X,'KOR(,M)',1X,'KOR(,N)',1X,'KOR(,O)',1X,'KOR(,P)',1X, -cTR 6 'KOR(,Q)',1X,'KOR(,R)',1X,'KOR(,S)',1X,'KOR(,T)',1X,'KOR(,U)',1X, -cTR 7 'KOR(,V)',1X,'KOR(,W)',1X,'KOR(,X)',1X,'KOR(,Y)',1X,'KOR(,Z)') - DO 491 I=1,LJ - WRITE(21,492) (KOR(I,J),J=1,LJ) - 492 FORMAT(/1X,9X,35F8.5) - 491 CONTINUE -C WRITE(6,479) -C 479 FORMAT(//13X, -C 1 'F(I,1)',6X,'F(I,2)',6X,'F(I,3)',6X,'F(I,4)',6X,'F(I,5)',5X, -C 2 'KOR(I,1)',4X,'KOR(I,2)',4X,'KOR(I,3)',4X,'KOR(I,4)',4X, -C 3 'KOR(I,5)') -C DO 480 I=1,L -C WRITE(6,481) I,(F(I,J),J=1,5),(KOR(I,J),J=1,5) -C 481 FORMAT(/1X,I1,6H.LAYER,1X,10E12.4) -C 480 CONTINUE -C WRITE(6,483) -C 483 FORMAT(//12X, -C 1 'KL(I,1)',5X,'KL(I,2)',5X,'KL(I,3)',5X,'KL(I,4)',5X,'KL(I,5)',6X, -C 2 'K(I,1)',6X,'K(I,2)',6X,'K(I,3)',6X,'K(I,4)',6X,'K(I,5)') -C DO 482 I=1,L -C WRITE(6,485) I,(KL(I,J),J=1,5),(K(I,J),J=1,5) -C 485 FORMAT(/1X,I1,6H.LAYER,1X,10E12.4) -C 482 CONTINUE + & 'KOR(,1)',1X,'KOR(,2)',1X,'KOR(,3)',1X,'KOR(,4)',1X,'KOR(,5)',1X, + & 'KOR(,6)',1X,'KOR(,7)',1X,'KOR(,8)',1X,'KOR(,9)',1X,'KOR(,A)',1X, + & 'KOR(,B)',1X,'KOR(,C)',1X,'KOR(,D)',1X,'KOR(,E)',1X,'KOR(,F)',1X, + & 'KOR(,G)',1X,'KOR(,H)',1X,'KOR(,I)',1X,'KOR(,J)',1X,'KOR(,K)',1X, + & 'KOR(,L)',1X,'KOR(,M)',1X,'KOR(,N)',1X,'KOR(,O)',1X,'KOR(,P)',1X, + & 'KOR(,Q)',1X,'KOR(,R)',1X,'KOR(,S)',1X,'KOR(,T)',1X,'KOR(,U)',1X, + & 'KOR(,V)',1X,'KOR(,W)',1X,'KOR(,X)',1X,'KOR(,Y)',1X,'KOR(,Z)') + DO I=1,LJ + WRITE(21,492) (KOR(I,J),J=1,LJ) + 492 FORMAT(/1X,9X,35F8.5) + ENDDO C C INTEGRAL IMPLANTATION , SPUTTERING , BACKSCATTERING , TRANSMISSION C @@ -2858,40 +2737,42 @@ C HN=DBLE(NH) EMV=CVMGT(EMX/HN,E0,E0.LE.0.D0) EIM=DBLE(IIM)*EMV - DO 1550 J=1,LJ - ISPA = ISPA+IBSP(J) - 1550 ESPA = ESPA+EBSP(J) - DO 1702 J=1,LJ - ISPAT = ISPAT+ITSP(J) - 1702 ESPAT = ESPAT+ETSP(J) - WRITE(21,1500) IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,ISPAT,ESPAT, - 1 tryE,negE,epsilon,prcoeff + DO J=1,LJ + ISPA = ISPA+IBSP(J) + ESPA = ESPA+EBSP(J) + ENDDO + DO J=1,LJ + ISPAT = ISPAT+ITSP(J) + ESPAT = ESPAT+ETSP(J) + ENDDO + WRITE(21,1500) IIM,EIM,IB,EB,IT,ET,ISPA,ESPA,ISPAT,ESPAT, tryE + & ,negE,epsilon,prcoeff 1500 FORMAT(1H1,//11X,20HIMPLANTED PARTICLES=,I7,5X,7HENERGY=,E10.4, - 1 3H EV/7X,24HBACKSCATTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, - 2 3H EV/9X,22HTRANSMITTED PARTICLES=,I7,5X,7HENERGY=,E10.4, - 3 3H EV/7X,24HBACKSPUTTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, - 4 3H EV/6X,'TRANSM. SPUTT. PARTICLES=',I7,5X,7HENERGY=,E10.4, - 5 3H EV/15X,16HTRIED PARTICLES=,I7 - 6 /9X,22HPARTICLES with neg. E=,I7, - 7 /6X,25HTHOMAS FERMI RED. ENERGY=,2X,E10.4, - 8 /6X,25HSCALED PART. REFL. COEFF=,2x,E10.4) + & 3H EV/7X,24HBACKSCATTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, + & 3H EV/9X,22HTRANSMITTED PARTICLES=,I7,5X,7HENERGY=,E10.4, + & 3H EV/7X,24HBACKSPUTTERED PARTICLES=,I7,5X,7HENERGY=,E10.4, + & 3H EV/6X,'TRANSM. SPUTT. PARTICLES=',I7,5X,7HENERGY=,E10.4, + & 3H EV/15X,16HTRIED PARTICLES=,I7 + & /9X,22HPARTICLES with neg. E=,I7, + & /6X,25HTHOMAS FERMI RED. ENERGY=,2X,E10.4, + & /6X,25HSCALED PART. REFL. COEFF=,2x,E10.4) if(l.gt.1) then - do i=1,l + do i=1,l nli(i)=IDINT(xx(i)/cw+0.01) - enddo - do i=1,l + enddo + do i=1,l do j=nli(i-1)+1,nli(i) irpl(i)=irpl(i)+irp(j) enddo - enddo - WRITE(21,15001) -15001 FORMAT(/33x,'FROM BIN WIDTH',2x,'FROM LAYER THICKNESS') - IF(depth_interval_flag.EQ.0) WRITE(21,'(36x,5HWRONG)') + enddo + WRITE(21,15001) +15001 FORMAT(/33x,'FROM BIN WIDTH',2x,'FROM LAYER THICKNESS') + IF(depth_interval_flag.EQ.0) WRITE(21,'(36x,5HWRONG)') - do i=1,l + do i=1,l WRITE(21,1501) i,irpl(i),number_in_layer(i) - enddo - 1501 FORMAT(/1x,'IMPLANTED PARTICLES (LAYER ',i1,')=',i16,2x,i16) + enddo + 1501 FORMAT(/1x,'IMPLANTED PARTICLES (LAYER ',i1,')=',i16,2x,i16) endif CSUM=ICSUM CSUMS=ICSUMS @@ -2901,9 +2782,9 @@ C CST=DBLE(ICSUM-ICDI) WRITE(21,1511) AVCSUM,AVCDIS,AVCSMS 1511 FORMAT(//2X,'PROJECTILES : ', - 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, - 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F8.3/65X, - 3 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F8.3) + & 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, + & 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F8.3/65X, + & 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F8.3) IF(YH.LE.1.D0) GO TO 1502 AVNLI=ENUCLI/YH VANLI=ENL2I/YH-AVNLI*AVNLI @@ -2913,24 +2794,24 @@ C VAILI=EIL2I/YH-AVILI*AVILI SIGILI=DSQRT(VAILI) DFIILI=SIGILI/YH - CALL MOMENTS(FIX0,SEX,THX,FOX,FIX,SIX,SIGMAX,DFIX0,DSEX,DTHX, - 1 XSUM,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,YH) - CALL MOMENTS(FIR0,SER,THR,FOR,FIR,SIR,SIGMAR,DFIR0,DSER,DTHR, - 1 RSUM,R2SUM,R3SUM,R4SUM,R5SUM,R6SUM,YH) + CALL MOMENTS(FIX0,SEX,THX,FOX,FIX,SIX,SIGMAX,DFIX0,DSEX,DTHX, XSUM + & ,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,YH) + CALL MOMENTS(FIR0,SER,THR,FOR,FIR,SIR,SIGMAR,DFIR0,DSER,DTHR, RSUM + & ,R2SUM,R3SUM,R4SUM,R5SUM,R6SUM,YH) CALL MOMENTS(FIP0,SEP,THP,FOP,FIP,SIP,SIGMAP,DFIP0,DSEP,DTHP, - 1 PLSUM,PL2SUM,PL3SUM,PL4SUM,PL5SUM,PL6SUM,YH) - CALL MOMENTS(FIE0,SEE,THE,FOE,FIE,SIE,SIGMAE,DFIE0,DSEE,DTHE, - 1 EEL,EEL2,EEL3,EEL4,EEL5,EEL6,CSUM) + & PLSUM,PL2SUM,PL3SUM,PL4SUM,PL5SUM,PL6SUM,YH) + CALL MOMENTS(FIE0,SEE,THE,FOE,FIE,SIE,SIGMAE,DFIE0,DSEE,DTHE, EEL + & ,EEL2,EEL3,EEL4,EEL5,EEL6,CSUM) CALL MOMENTS(FIW0,SEW,THW,FOW,FIW,SIW,SIGMAW,DFIW0,DSEW,DTHW, - 1 EELWC,EELWC2,EELWC3,EELWC4,EELWC5,EELWC6,CSUM) - CALL MOMENTS(FII0,SEI,THI,FOI,FII,SII,SIGMAI,DFII0,DSEI,DTHI, - 1 EIL,EIL2,EIL3,EIL4,EIL5,EIL6,CSUM) - CALL MOMENTS(FIS0,SES,THS,FOS,FIS,SIS,SIGMAS,DFIS0,DSES,DTHS, - 1 EPL,EPL2,EPL3,EPL4,EPL5,EPL6,CST) + & EELWC,EELWC2,EELWC3,EELWC4,EELWC5,EELWC6,CSUM) + CALL MOMENTS(FII0,SEI,THI,FOI,FII,SII,SIGMAI,DFII0,DSEI,DTHI, EIL + & ,EIL2,EIL3,EIL4,EIL5,EIL6,CSUM) + CALL MOMENTS(FIS0,SES,THS,FOS,FIS,SIS,SIGMAS,DFIS0,DSES,DTHS, EPL + & ,EPL2,EPL3,EPL4,EPL5,EPL6,CST) WRITE(21,7117) 7117 FORMAT(/20X,' MEAN ',4X,' VARIANCE ',4X,' SKEWNESS ',4X, - 1 ' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',3X,' ERROR 2.M ', - 2 3X,' ERROR 3.M ') + & ' KURTOSIS ',5X,' SIGMA ',3X,' ERROR 1.M ',3X + & ,' ERROR 2.M ', 3X,' ERROR 3.M ') WRITE(21,7227) FIX0,SEX,THX,FOX,SIGMAX,DFIX0,DSEX,DTHX 7227 FORMAT(1X,' PENETRATION',5X,1P1E12.4,7E14.4) WRITE(21,7228) FIR0,SER,THR,FOR,SIGMAR,DFIR0,DSER,DTHR @@ -2951,10 +2832,10 @@ C WRITE(21,7234) FIS0,SES,THS,FOS,SIGMAS,DFIS0,DSES,DTHS 7234 FORMAT(1X,' SUBTHR.LOSS',5X,1P1E12.4,7E14.4) 1502 CONTINUE -c + IF(YH.LT.1.D0) GO TO 7235 - CALL MOMENT(X1SD,X2SD,X3SD,X4SD,X5SD,X6SD - 1 ,XSUM,X2SUM,X3SUM,X4SUM,X5SUM,X6SUM,YH) + CALL MOMENT(X1SD,X2SD,X3SD,X4SD,X5SD,X6SD ,XSUM,X2SUM,X3SUM,X4SUM + & ,X5SUM,X6SUM,YH) WRITE(21,7118) WRITE(21,1513) X1SD,X2SD,X3SD,X4SD,X5SD,X6SD 1513 FORMAT(1X,' PENETRATION',5X,1P1E12.4,5E14.4) @@ -2966,18 +2847,18 @@ c ACSBER=DBLE(ICSBR)/HN WRITE(21,1599) ACSUMR,ACDISR,ACSBER 1599 FORMAT(///2X,'RECOILES (PER PROJ.) : ', - 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, - 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/76X, - 3 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3) + & 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.1,3X, + & 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/76X, + & 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3) IF(NPA+NSA.EQ.0) GO TO 1453 ACSUR=CSUMR/(DBLE(NPA+NSA)) ACDIR=DBLE(ICDIR)/(NPA+NSA) ACSBR=DBLE(ICSBR)/(NPA+NSA) WRITE(21,1598) ACSUR,ACDIR,ACSBR 1598 FORMAT(/2X,'RECOILES (PER KNOCKON) : ', - 1 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.3,3X, - 2 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/,76X, - 3 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3/) + & 'MEAN NUMBER OF ELASTIC COLLISIONS = ',1F8.3,3X, + & 'MEAN NUMBER OF EL.COLL.(E > EDISPL.) = ',F10.3/,76X, + & 'MEAN NUMBER OF EL.COLL.(E > SB(1)) = ',F10.3/) IF(NJ(1).LT.2) GO TO 1453 ACDR11=DBLE(ICDIRJ(1,1))/(NPA+NSA) ACDR12=DBLE(ICDIRJ(1,2))/(NPA+NSA) @@ -2985,157 +2866,139 @@ c ACDR22=DBLE(ICDIRJ(2,2))/(NPA+NSA) WRITE(21,1451) ACDR11,ACDR12,ACDR21,ACDR22 1451 FORMAT(76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-1 = ',F10.3/ - 1 ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-2 = ',F10.3/ - 2 ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-1 = ',F10.3/ - 3 ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-2 = ',F10.3) + & ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 1-2 = ',F10.3/ + & ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-1 = ',F10.3/ + & ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-2 = ',F10.3) 1453 CONTINUE - 590 WRITE(21,600) - 600 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES,2X,10HNORM.DEPTH,1X, - 110HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X,10HELLOSS(EV), - 22X,10HDAMAGE(EV),2X,10HPHONON(EV),2X,10HCASCAD(EV),5X,3HDPA/) -C + 590 WRITE(21,600) + 600 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES,2X,10HNORM.DEPTH,1X, + & 10HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X + & ,10HELLOSS(EV), 2X,10HDAMAGE(EV),2X,10HPHONON(EV),2X + & ,10HCASCAD(EV),5X,3HDPA/) + IF(depth_interval_flag.EQ.0) THEN - WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', - & 'AGREE WITH LAYER THICKNESS' - WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', - & 'AGREE WITH LAYER THICKNESS' -c WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT AGREE WIT -c 1H LAYER THICKNESS' -c WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT AGREE WIT -c 1H LAYER THICKNESS' - ENDIF -C + WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', + & 'AGREE WITH LAYER THICKNESS' + WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', + & 'AGREE WITH LAYER THICKNESS' + ENDIF + #if defined (OS_WIN) WRITE(22,6002) - 6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES) + 6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES) #else write(22,'(a)') ' DEPTH PARTICLES' #endif IF(YH.LT.1.D0) GO TO 603 - DO 602 I=0,101 - RIRP(I) = DBLE(IRP(I))/YH - 602 CONTINUE - 603 D1=0. + DO I=0,101 + RIRP(I) = DBLE(IRP(I))/YH + ENDDO + 603 D1=0. D2=CW WRITE(21,601) D1,IRP(0),RIRP(0) - 601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4) -c DO 1441 J=1,NJ(1) - DO 1441 J=1,LJ - DO 1441 I=1,100 - ICDT(I)=ICDT(I)+ICD(I,J) - ICDTR(I)=ICDTR(I)+ICDR(I,J) - 1441 CONTINUE - DO 1442 K=1,NJ(1) - DO 1442 J=1,NJ(1) - DO 1442 I=1,100 - ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) - 1442 CONTINUE - DO 35 I=0,101 - IIRP=IIRP+IRP(I) - TRIRP=TRIRP+RIRP(I) - 35 CONTINUE - DO 36 I=1,100 - IIPL=IIPL+IPL(I) - TION=TION+ION(I) - TDENT=TDENT+DENT(I) - TDMGN=TDMGN+DMGN(I) - TELGD=TELGD+ELGD(I) - TPHON=TPHON+PHON(I) - TCASMO=TCASMO+CASMOT(I) - ICDTT=ICDTT+ICDT(I) - TIONR=TIONR+IONR(I) - TDENTR=TDENTR+DENTR(I) - TELGDR=TELGDR+ELGDR(I) - TDMGNR=TDMGNR+DMGNR(I) - TPHONR=TPHONR+PHONR(I) -C TCASMOR=TCASMOR+CASMOTR(I) - ICDTTR=ICDTTR+ICDTR(I) - 36 CONTINUE + 601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4) + DO J=1,LJ + DO I=1,100 + ICDT(I)=ICDT(I)+ICD(I,J) + ICDTR(I)=ICDTR(I)+ICDR(I,J) + ENDDO + ENDDO + DO K=1,NJ(1) + DO J=1,NJ(1) + DO I=1,100 + ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) + ENDDO + ENDDO + ENDDO + DO I=0,101 + IIRP=IIRP+IRP(I) + TRIRP=TRIRP+RIRP(I) + ENDDO + DO I=1,100 + IIPL=IIPL+IPL(I) + TION=TION+ION(I) + TDENT=TDENT+DENT(I) + TDMGN=TDMGN+DMGN(I) + TELGD=TELGD+ELGD(I) + TPHON=TPHON+PHON(I) + TCASMO=TCASMO+CASMOT(I) + ICDTT=ICDTT+ICDT(I) + TIONR=TIONR+IONR(I) + TDENTR=TDENTR+DENTR(I) + TELGDR=TELGDR+ELGDR(I) + TDMGNR=TDMGNR+DMGNR(I) + TPHONR=TPHONR+PHONR(I) + ICDTTR=ICDTTR+ICDTR(I) + ENDDO do im1=100,1,-1 -C if(ipl(im1).ne.0.or.ion(im1).ne.0.) go to 20 - if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 + if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 enddo im1=1 - 20 im1=min0(im1+2,100) + 20 im1=min0(im1+2,100) DO 11 I=1,im1 - WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I),DMGN(I) - 1 ,ELGD(I),PHON(I),CASMOT(I),ICDT(I) - Dmid=(D2-D1)/2+D1 - WRITE(22,701) Dmid,IRP(I) - 700 FORMAT(1X,F6.0,1H-,F6.0,I10,E12.4,I10,1P1E14.4,5E12.4,I8) - 701 FORMAT(1X,F6.0,2x,I10) - D1=D2 - 11 D2=D2+CW + WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I), + & DMGN(I),ELGD(I),PHON(I),CASMOT(I),ICDT(I) + Dmid=(D2-D1)/2+D1 + WRITE(22,701) Dmid,IRP(I) + 700 FORMAT(1X,F6.0,1H-,F6.0,I10,E12.4,I10,1P1E14.4,5E12.4,I8) + 701 FORMAT(1X,F6.0,2x,I10) + D1=D2 + 11 D2=D2+CW WRITE(21,604) D2-CW,IRP(101),RIRP(101) - 604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4) + 604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4) WRITE(21,710) IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO - 1 ,ICDTT - 710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8) - DO 1531 J=1,NJ(1) - DO 1531 I=1,100 - ELET(J)=ELET(J)+ELE(I,J) - ELIT(J)=ELIT(J)+ELI(I,J) - ELPT(J)=ELPT(J)+ELP(I,J) - ELDT(J)=ELDT(J)+ELD(I,J) - ICDJT(J)=ICDJT(J)+ICD(I,J) - ICDJTR(J)=ICDJTR(J)+ICDR(I,J) - ICDITR(J)=ICDITR(J)+ICDIRN(I,J) -C ELETR(J)=ELETR(J)+ELE(I,J) -C ELITR(J)=ELITR(J)+ELI(I,J) -C ELPTR(J)=ELPTR(J)+ELP(I,J) -C ELDTR(J)=ELDTR(J)+ELD(I,J) - 1531 CONTINUE -c IF(NJ(1).LT.2) GO TO 1455 + & ,ICDTT + 710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8) + DO J=1,NJ(1) + DO I=1,100 + ELET(J)=ELET(J)+ELE(I,J) + ELIT(J)=ELIT(J)+ELI(I,J) + ELPT(J)=ELPT(J)+ELP(I,J) + ELDT(J)=ELDT(J)+ELD(I,J) + ICDJT(J)=ICDJT(J)+ICD(I,J) + ICDJTR(J)=ICDJTR(J)+ICDR(I,J) + ICDITR(J)=ICDITR(J)+ICDIRN(I,J) + ENDDO + ENDDO + WRITE(21,1521) 1521 FORMAT(1H1,4X,'DEPTH(A)' - 1 ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)' - 2 ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)' - 3 ,2X,'DPA(1)',2X,'DPA(2)'/) + & ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)' + & ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)' + & ,2X,'DPA(1)',2X,'DPA(2)'/) D1=0. D2=CW do im2=100,1,-1 -C if(eli(im2,1).ne.0..or.eli(im2,2).ne.0.) go to 30 - if(.NOT.EQUAL(eli(im2,1),0.D0).or. - # (.NOT.EQUAL(eli(im2,2),0.D0))) goto 30 + if(.NOT.EQUAL(eli(im2,1),0.D0).or.(.NOT.EQUAL(eli(im2,2),0.D0)) + & ) goto 30 enddo im2=1 - 30 im2=MIN0(im2+2,100) + 30 im2=MIN0(im2+2,100) DO 1525 I=1,im2 - WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1) - 1 ,ELI(I,2),ELE(I,2),ELD(I,2),ELP(I,2),ICD(I,1),ICD(I,2) - 1523 FORMAT(1X,F6.0,1H-,F6.0,1P8E12.4,2I8) - D1=D2 - D2=D2+CW + WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1), + & ELI(I,2),ELE(I,2),ELD(I,2),ELP(I,2),ICD(I,1),ICD(I,2) + 1523 FORMAT(1X,F6.0,1H-,F6.0,1P8E12.4,2I8) + D1=D2 + D2=D2+CW 1525 CONTINUE - WRITE(21,1533) ELIT(1),ELET(1),ELDT(1),ELPT(1) - 1 ,ELIT(2),ELET(2),ELDT(2),ELPT(2),ICDJT(1),ICDJT(2) + WRITE(21,1533) ELIT(1),ELET(1),ELDT(1),ELPT(1),ELIT(2),ELET(2) + & ,ELDT(2),ELPT(2),ICDJT(1),ICDJT(2) 1533 FORMAT(/14X,1P8E12.4,2I8///) -C WRITE(21,1481) -C1481 FORMAT(1H1,2X,'D(A)' -C 1 ,2X,' DAMAGE(1)',3X,' DAMAGE(2)',1X,' DAMAGE(3)' -C 2 ,2X,' DAMAGE(4)',3X,' DAMAGE(5)'/) -C DO 1482 I=1,100 -C WRITE(6,1483) I,ELD(I,1),ELD(I,2),ELD(I,3),ELD(I,4),ELD(I,5) -C1482 CONTINUE -C1483 FORMAT(1X,I5,1P5E14.4) -C WRITE(6,1484) ELDT(1),ELDT(2),ELDT(3),ELDT(4),ELDT(5) -C1484 FORMAT(/1X,5X,1P5E14.4///) - DO 1491 I=1,L-1 - ILD(I)=IDINT(XX(I)/CW+0.01D0) - IF(ILD(I).GT.100) ILD(I)=100 - DO 1492 J=1,ILD(I) - DLI(I)=DLI(I)+DMGN(J) - 1492 CONTINUE - 1491 CONTINUE + DO I=1,L-1 + ILD(I)=IDINT(XX(I)/CW+0.01D0) + IF(ILD(I).GT.100) ILD(I)=100 + DO J=1,ILD(I) + DLI(I)=DLI(I)+DMGN(J) + ENDDO + ENDDO DLI(L)=TDMGN -C WRITE(21,*) 'L=',L,' XX=',XX,' DLI=',DLI DO 1493 I=L,2,-1 - DLI(I)=DLI(I)-DLI(I-1) + DLI(I)=DLI(I)-DLI(I-1) 1493 CONTINUE DO 1494 I=1,L - WRITE(21,1495) I,DLI(I) - 1495 FORMAT(/5X,'DAMAGE IN LAYER ',I1,' : ',1P1E12.4) + WRITE(21,1495) I,DLI(I) + 1495 FORMAT(/5X,'DAMAGE IN LAYER ',I1,' : ',1P1E12.4) 1494 CONTINUE 1455 CONTINUE if(irl.eq.0) goto 1497