diff --git a/trimsp/src/trimsp7l.F b/trimsp/src/trimsp7l.F index 889748e..1ef405a 100644 --- a/trimsp/src/trimsp7l.F +++ b/trimsp/src/trimsp7l.F @@ -34,7 +34,7 @@ C C MAJOR CHANGES: C C Sep 1998: to create an executable for PC running under WIN95 -C using the DIGITAL VISUAL FORTRAN compiler (ed. Aug. 97) +C using the DIGITAL VISUAL FORTRAN compiler (ed. Aug. 97) C Insertion of !DEC$REAL:8 (all REAL are REAL*8) C Conversion of function to double precession C Conversion of function to integer*4 @@ -46,7 +46,7 @@ C UNIT 21 (output.data) changed to file name AUSGABE1.out C Insertion of UNIT 22 (only range data), file name AUSGABE1.rge C Dec 1998: Introduction of gaussian distributed projectile energies C new variable RI2 = random number initializer for gaussian energy distribution -C new variable ISEED2 = random number for gaussian energy distribution +C new variable ISEED2 = random number for gaussian energy distribution C new variable Esig = sigma of the energy distribution C if Esig=0. then fixed projectile energy C new variable Epar = projectile energy @@ -69,7 +69,7 @@ C Unit 99 for file ausgabe1.err inserted, if and IF... (see below C then a message is included into this file C IF's inserted - can be found by using find WRITE(99, C DABS inserted for DLOG and DLOG10 arguments -C Mai 1999: Version h +C Mai 1999: Version h C for TRIMSP simulations running in batch mode (i.e. to calculate C a serie of different range profiles using the same target an output C file FOR33 is created. In this file the following parameters are @@ -87,8 +87,8 @@ C Range file (output to unit 22) slightly changed. Now the midth o C is calculated and only the number of particles is written to the file C C -C Jun 1999: Version j1 -C changed calculation of Firsoc screening length according to the IRCU 49 report. +C Jun 1999: Version j1 +C changed calculation of Firsoc screening length according to the IRCU 49 report. C the Andersen-Ziegler tables (file stopping.dat) are revised,i.e. the coefficients C from the ICRU 49 report are included for all elements (file stopicru.dat) C @@ -103,8 +103,8 @@ C Nov 1999: Version j1a C UNIT33: now Etrans, sigmaEtrans, Eback, sigmaEback included C Bug fixed in line 2084, location of 'write to unit33' changed. C Dec 1999: Version k -C inclusion of variables -C tryE :how often a random energy is calculated by the random generator +C inclusion of variables +C tryE :how often a random energy is calculated by the random generator C for large E0 this number should be the same as the number of projectiles nh C negE :how often a negative energy is calculated by the random generator C tryE - negE should be nh @@ -118,8 +118,8 @@ C the options alpha = -1 or alpha = -2 (see file TRVMC95-v4k.txt) C inclusion of RI3: random seed for the calculation of the gaussian distribution C of alpha. C Header line for file for33 included. -C Jan 2000: no new version but -C included energy scaling of particle reflection coefficients after Thomas et al. +C Jan 2000: no new version but +C included energy scaling of particle reflection coefficients after Thomas et al. C NIM B69 (1992) 427 C the calculated particle reflection coefficients prcoeff are written to the fort.33 file C prc is calculated if 1st layer consists only one element ! @@ -149,47 +149,40 @@ C in the 100 depth intervals agrees wit C number of particles stopped in the different layers C IF 0 not, message will be written in the range file C UNIT 21 and UNIT22 -c November 2000. Tanya Riseman -c Starts porting program to Open VMS and DEC UNIX. -c Compile options: f90/list/warn/ext +c Nov 2000: Tanya Riseman +c Starts porting program to Open VMS and DEC UNIX. +c Compile options: f90/list/warn/ext c Minor problems with continuation characters in wrong column and c with a few variables undeclared, mostly functions. -c Make code fit on 72 columns, because -extend_source +c Make code fit on 72 columns, because -extend_source c does not appear to work on DEC Unix. -c Start changing strings in format and write statments so that +c Start changing strings in format and write statments so that c strings don't straddle continuation character in c column 6. Straddling can be non-portable! -c Comment out !DEC$REAL:8 because all +c Comment out !DEC$REAL:8 because all c reals are already declared REAL*8. Add implicit none. c -c June 2002. Thomas Prokscha PSI -c Stopped porting to f90, use f77 only. -c replaced all non-standard f77 function by standard functions. -c use ranlux random number generator from the CERN library (libmathlib.a must -c be installed) to get rid of machine specific random number generators. -c Add pre-compiler instructions for making different output for the .rge files -c on Windows and Unix. +c Jun 2002: Thomas Prokscha PSI +c Stopped porting to f90, use f77 only. +c replaced all non-standard f77 function by standard functions. +c use ranlux random number generator from the CERN library (libmathlib.a must +c be installed) to get rid of machine specific random number generators. +c Add pre-compiler instructions for making different output for the .rge files +c on Windows and Unix. c -c October 2002 Thomas Prokscha PSI -c corrected error in the calculation of the Thomas-Fermi reduced energy: -c it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3)) -c it must be: -c 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3))) +c Oct 2002: Thomas Prokscha PSI +c corrected error in the calculation of the Thomas-Fermi reduced energy: +c it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3)) +c it must be: +c 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3))) c -c August 2009 Zaher Salman PSI -c changed input file reading to non-formatted. This still works -c with old input files, but makes it much easier to create new -c input files without the hassle of formatting. -c Included the source of ranlux for random numbers generation -c into trimsp7l source. No need for cern libraries to be installed. +c Aug 2009: Zaher Salman PSI +c Changed input file reading to non-formatted. This still works +c with old input files, but makes it much easier to create new +c input files without the hassle of formatting. +c Included the source of ranlux for random numbers generation +c into trimsp7l source. No need for cern libraries to be installed. c -CDIR$ NOLIST -C -cTR !DEC$REAL:8 -C -C IMPLICIT INTEGER (i-j) -C IMPLICIT REAL*8 (a-h,k-z) -C c------------------------------------------- c check OS c @@ -237,9 +230,7 @@ c # ,MEAGS(102,12,22,30) # ,MEASL(75,21,30) INTEGER*4 MEAST(102,22,30),MAGST(62,22,30) -C # ,MEAGST(102,36,22,10) von Eckstein herauskommentiert # ,MEASTL(75,21,30) -C REAL*8 MEAGSL(75,36,21),EAGSL(75) von Eckstein herauskommentiert INTEGER*4 NJ(7),JT(7),ILD(7) INTEGER*4 LLL(64),JJJ(64),IK(64) INTEGER*4 me(5000),nli(0:7),irpl(7) @@ -257,7 +248,7 @@ C REAL*8 MEAGSL(75,36,21),EAGSL(75) von Eckstein herauskommentiert INTEGER*4 IPB1,KIB,IPT,IE,IERLOG,IAGB,KIT,IMA,IIM INTEGER*4 im1,im2,im3,IG2,ies,ias INTEGER*4 JE,JA,JG,JTJ,JTK,JTL -C REAL Variablen +C REAL Variables REAL*8 CVMGT REAL*8 X(64),Y(64),Z(64),E(64),PL(64) # ,COSX(64),COSY(64),COSZ(64),SINE(64) @@ -320,13 +311,11 @@ C REAL Variablen # ,KL(35,35),KOR(35,35),KLM(7,35) REAL*8 MU1(35),EC1(35),A1(35),F1(35),KL1(35),KOR1(35) # ,DI(35),EP(35),ZZ(35),TM(35) -C REAL*8 SL(64,5),SM(64,5),SH(64,5),CH(92,12),CHE(92,9) REAL*8 CH1(7,5),CH2(7,5),CH3(7,5),CH4(7,5),CH5(7,5) REAL*8 CHM1(7) REAL*8 SM(64),SH(64,5) REAL*8 FIESB(30),SEESB(30),THESB(30),FOESB(30) # ,SGMESB(30),DFIESB(30),DSEESB(30),DTHESB(30) -C REAL*8 ESVDL(2000) REAL*8 pi,c,E0,de,alfa,z1,rtheta,cu,enot,esb,est,esp REAL*8 E2,AB,FP,AN REAL*8 Esig,Epar @@ -406,7 +395,7 @@ c fuer part. reflec. coeff. Berechnung REAL*8 PL1S,PL2S,PL3S,PL4S,PL5S,PL6S REAL*8 YH,HN,CST,BI,BIL,YSP,YSPL,EEE REAL*4 random,ran2(2) -C CHARACTER Variablen +C CHARACTER Variables CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2 CHARACTER filein*8,inext*4,fileout*8,outext*4,innam*12,outnam*12 CHARACTER rgenam*12,rgeext*4,errnam*12,errext*4 @@ -416,10 +405,10 @@ C CHARACTER Variablen CHARACTER month_start*4,month_stop*4,day_start*2,day_stop*2 CHARACTER year_start*4,year_stop*4,hour_start*2,hour_stop*2 CHARACTER min_start*2,min_stop*2,sec_start*2,sec_stop*2 -C + COMMON /A/ M1,VELC,ZARG COMMON /B/ TI,SHEATH,CALFA -C + DATA PI/3.14159265358979D0/, ICW/100/, E2/14.399651D0/ DATA AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/ DATA inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/ @@ -477,7 +466,6 @@ C DATA ICDR/3500*0/,ICDIRN/3500*0/,IONR/100*0.D0/ DATA DENTR/100*0.D0/,DMGNR/100*0.D0/ DATA IPL/100*0/,IPLB/100*0/,IPLT/100*0/ -c DATA IRP/102*0/ gibt witzigweise einen Fehler, aber warum???? DATA IRPL/7*0/ DATA ICDJT/35*0/,ICDJTR/35*0/,ICDITR/35*0/ DATA ICD/3500*0/,ELP/3500*0.D0/,ELD/3500*0.D0/ @@ -498,21 +486,18 @@ c DATA IRP/102*0/ gibt witzigweise einen Fehler, aber warum???? DATA ELETR/35*0.D0/,ELITR/35*0.D0/,ELPTR/35*0.D0/ DATA ELDTR/35*0.D0/ DATA Epar/0.D0/ -c part. refl. coeff. from Thomas et al. +C part. refl. coeff. from Thomas et al. DATA PRC/0.825D0,21.41D0,8.606D0,0.6425D0,1.907D0,1.927D0/ DATA number_in_layer /7*0/ -C -C EXTERNAL CVMGT,ILLZ,SCOPY,ISRCHEQ,ISRCHFGE,ISRCHFGT -C + innam=filein//inext outnam=fileout//outext rgenam=fileout//rgeext errnam=fileout//errext -C + OPEN(UNIT=99,file=errnam,STATUS='new') OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=13591) -C OPEN(UNIT=31,NAME='energ.dat',STATUS='new') -C + READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 # ,IPOT,IPOTR,IRL @@ -546,24 +531,15 @@ C WRITE(*,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) WRITE(*,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) WRITE(*,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) -c WRITE(*,103) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) -c WRITE(*,103) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) -c WRITE(*,103) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) -c WRITE(*,107) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) -c WRITE(*,107) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) -c WRITE(*,107) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) -c WRITE(*,107) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) -c WRITE(*,107) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) 1359 CONTINUE -C + 100 FORMAT(2F7.2,1F12.2,7F9.2) 101 FORMAT(I9,3F8.0,1F7.2,1F7.0,2F7.2,6I4,I3) 102 FORMAT(7F9.2,14F7.2) 103 FORMAT(5F9.4) 107 FORMAT(5F12.6) -C + C open statement for output files, removed from line 2449 ff to here -C OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000) GOTO 6001 6000 WRITE(*,*)' File schon vorhanden, Gib neue Ausgabedatei an (A8)' @@ -574,13 +550,11 @@ C 6001 OPEN(UNIT=22,FILE=rgenam,STATUS='new') WRITE(21,1000) 1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - V TrimSP7L 17.Oct.02 TP *') -c TP 1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - Vers. TrimSP7L 22.Nov.00 *') -c TR 1000 FORMAT(1H1/,6X,'** PROGRAM TRVMC95 - V TrimSP7L 26.06.00 **') -C + C 1st CALL DATE_AND_TIME CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3), # Date_Time) -C + IF(Date_Time(2).EQ.1) THEN month_start='Jan.' days_start_total=Date_Time(3) @@ -621,13 +595,13 @@ C C in seconds from beginning of year seconds_start_total=Date_Time(7)+(Date_Time(6)*60)+ # (Date_Time(5)*3600)+(days_start_total-1)*86400 -C + READ(Real_Clock(1)(1:4),'(A4)')year_start READ(Real_Clock(1)(7:8),'(A2)')day_start READ(Real_Clock(2)(1:2),'(A2)')hour_start READ(Real_Clock(2)(3:4),'(A2)')min_start READ(Real_Clock(2)(5:6),'(A2)')sec_start -C + WRITE(21,*) WRITE(21,10050)day_start,month_start,year_start, # hour_start,min_start,sec_start @@ -635,7 +609,6 @@ C # 1x,A2,':',A2,':',A2) C SET INTERVAL CONSTANTS FOR OUTPUT -C DE = 1.D0 DA = 3.D0 DG = 3.D0 @@ -656,18 +629,17 @@ C DAW = BW/DA DGW = BW/DG DGIK = BW/DGI -C + C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS -C PI2=2.D0*PI ABC=AB*FP LMAX=7 JMAX=5 L=ISRCHEQ(LMAX,DX(1),1,0.D0)-1 -C + C Checks wether depth interval is an integer denominator of layer thickness or not C If not, calculated implantation profile is not correct. -C + depth_interval_flag = 1 LOOP_Check_layer_thick : DO K=1,L-1 IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN @@ -679,7 +651,6 @@ C DO 165 I=1,L DO 155 J=1,JMAX IF(EQUAL(CO(I,J),0.D0)) GOTO 156 -C IF(CO(I,J).D0EQ.D00.D000) GO TO 156 155 CONTINUE J=JMAX+1 156 NJ(I)=J-1 @@ -819,9 +790,7 @@ C ZBL POTENTIAL (IPOTR=3) C KLM(LL,I) = KLM(LL,I)+CK(LL)*CO(LL,J)*KL(I,J+JT(LL)) KLM(LL,I) = KLM(LL,I)+CO(LL,J)*KL(I,J+JT(LL)) 193 CONTINUE -C -C ALPHA = CVMGT( .001, ALPHA, ALPHA.EQ.0. ) -C ALPHA = CVMGT( 179.999, ALPHA, ALPHA.EQ.180.) + ALPHA = CVMGT( .001D0, ALPHA, EQUAL(ALPHA,0.D0)) ALPHA = CVMGT( 179.999D0, ALPHA, EQUAL(ALPHA,180.D0)) @@ -831,9 +800,9 @@ C ALPHA = CVMGT( 179.999, ALPHA, ALPHA.EQ.180.) 8883 FORMAT(1X,'ERROR : IF ALPHA.GE.90. THEN IT MUST BE X0.LE.0.') GO TO 8000 8882 CONTINUE -C + C SET CONSTANT DISTANCES -C + TT = XX(L) INEL = 0 HLM = CVMGT( 0.D0, -.5D0*LM(1), INEL.EQ.0) @@ -848,7 +817,7 @@ C SUT = DMAX1(SUTR,DMAX1(SUT1,SUT2)) XC = CVMGT( X0, -SU, X0.GT.0.D0) RT = TT-RD -C + IF(E0.GE.0.D0) GO TO 51 C C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION @@ -868,13 +837,10 @@ C IY = IDINT(RI) IY2 = IDINT(RI2) IY3 = IDINT(RI3) -CC ANFANG = RANSET(IY) -CC ANFANG = SRAND48(IY) ISEED = IY ISEED2 = IY2 ISEED3 = IY3 -C WRITE(*,*)ISEED -C + IF ( E0.GT.0.D0 ) GO TO 47 IF ( ALPHA.GE.0.D0 ) THEN C @@ -884,10 +850,10 @@ C DO 49 IV=1,NUM EMX = EMX+E(IV) 49 CONTINUE - DO iv=1,num - ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) - me(ne) = me(ne)+1 - ENDDO + DO iv=1,num + ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) + me(ne) = me(ne)+1 + ENDDO c GO TO 56 C @@ -907,30 +873,26 @@ C IF(EQUAL(Esig,0.D0)) THEN C FIXED PROJECTILE ENERGY DO IV=1,NUM - E(IV) = E0 -C WRITE(*,*)' Da Esig=0 ist E=E0' - ENDDO - ELSE + E(IV) = E0 + ENDDO + ELSE C GAUSSIAN ENERGY DISTRIBUTION - DO IV=1,NUM -5200 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) - tryE = tryE+1 - IF(Epar.LE.0.0D0) THEN - negE = negE+1 - GO TO 5200 - ENDIF - E(IV)=Epar -C WRITE(*,*)E(IV),Epar,E0 - ENDDO + DO IV=1,NUM + 5200 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) + tryE = tryE+1 + IF(Epar.LE.0.0D0) THEN + negE = negE+1 + GO TO 5200 + ENDIF + E(IV)=Epar + ENDDO ENDIF C C die nachfolgende Zeile wurden von Linie 633 hier hin verschoben C SFE = DMIN1(SB(1),SB(L)) -C IF ( ALPHA.GE.0.D0 ) THEN -C - IF(EQUAL(ALPHASIG,0.D0))THEN + IF(EQUAL(ALPHASIG,0.D0))THEN C C FIXED PROJECTILE ANGLE C @@ -967,22 +929,16 @@ C C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL) C CDIR$ IVDEP - DO 53 IV=1,NUM -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) = DSQRT(RTHETA) - SINE(IV) = DSQRT(1.D0-RTHETA) - COSY(IV) = SINE(IV)*DCOS(RPHI) - COSZ(IV) = SINE(IV)*DSIN(RPHI) - 53 CONTINUE -C + DO IV=1,NUM + call ranlux(ran2,2) + RPHI = PI2*DBLE(ran2(1)) + RTHETA = DBLE(ran2(2)) + COSX(IV) = DSQRT(RTHETA) + SINE(IV) = DSQRT(1.D0-RTHETA) + COSY(IV) = SINE(IV)*DCOS(RPHI) + COSZ(IV) = SINE(IV)*DSIN(RPHI) + ENDDO + ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN C ELSEIF ( ALPHA.EQ.-1. AND. X0.GT.0. ) THEN C @@ -993,16 +949,16 @@ CDIR$ IVDEP CC RPHI = PI2*RANF() CC RPHI = PI2*DRAND48() CC RPHI = PI2*DBLE(RAN(ISEED)) - call ranlux(ran2,2) - RPHI = PI2*DBLE(ran2(1)) + 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) + 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) 50 CONTINUE C C ELSEIF ( ALPHA.EQ.-1. AND. X0.LE.0. ) THEN @@ -1013,48 +969,48 @@ CDIR$ IVDEP CC RPHI = PI2*RANF() CC RPHI = PI2*DRAND48() CC RPHI = PI2*DBLE(RAN(ISEED)) - call ranlux(ran2,2) - RPHI = PI2*DBLE(ran2(1)) + 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) + 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) 55 CONTINUE C - ENDIF + ENDIF C - 56 IF ( X0.GT.0.D0 ) GO TO 59 + 56 IF ( X0.GT.0.D0 ) GO TO 59 C C EXTERNAL START C DO 57 IV=1,NUM - 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 57 CONTINUE C C LOCUS OF FIRST COLLISION C - 59 JL = ISRCHFGT(L,XX(1),1,X0) + 59 JL = ISRCHFGT(L,XX(1),1,X0) C WRITE(*,*)X0 DO 58 IV=1,NUM CC RA = CVMGT(RANF(),1.0,X0.LE.0.0) CC RA = CVMGT(DRAND48(),1.0,X0.LE.0.0) CC RA = CVMGT(DBLE(RAN(ISEED)),1.D0,X0.LE.0.0D0) - call ranlux(random, 1) - RA = CVMGT(DBLE(random),1.D0,X0.LE.0.0D0) - X(IV) = XC + LM(JL) *RA *COSX(IV) - Y(IV) = LM(JL) *RA *COSY(IV) - Z(IV) = LM(JL) *RA *COSZ(IV) - PL(IV) = CVMGT(0.D0,LM(JL)*RA,X0.LE.0.0) + call ranlux(random, 1) + RA = CVMGT(DBLE(random),1.D0,X0.LE.0.0D0) + X(IV) = XC + LM(JL) *RA *COSX(IV) + Y(IV) = LM(JL) *RA *COSY(IV) + Z(IV) = LM(JL) *RA *COSZ(IV) + PL(IV) = CVMGT(0.D0,LM(JL)*RA,X0.LE.0.0) 58 CONTINUE C DO 199 IV=1,NUM @@ -1063,18 +1019,18 @@ C C C PROJECTILE LOOP C - 1 CONTINUE + 1 CONTINUE C NPROJ=NPROJ+1 C DO 63 IV=1,IH1 - CX(IV)=COSX(IV) - CY(IV)=COSY(IV) - CZ(IV)=COSZ(IV) - SX(IV)=SINE(IV) - DEES(IV)=0.D0 - DENS(IV)=0.D0 - DEN(IV)=0.D0 + CX(IV)=COSX(IV) + CY(IV)=COSY(IV) + CZ(IV)=COSZ(IV) + SX(IV)=SINE(IV) + DEES(IV)=0.D0 + DENS(IV)=0.D0 + DEN(IV)=0.D0 63 CONTINUE KK1=KK0 C @@ -1085,40 +1041,40 @@ C C CHOICE OF COLLISION PARTNERS C DO 298 IV=1,IH1 - call ranlux(random, 1) - JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1 + call ranlux(random, 1) + JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1 CC # ,RANF())+JT(LLL(IV)) CC # ,DRAND48())+JT(LLL(IV)) CC # ,DBLE(RAN(ISEED)))+JT(LLL(IV)) - # ,DBLE(random))+JT(LLL(IV)) + # ,DBLE(random))+JT(LLL(IV)) 298 CONTINUE DO 67 IV=1,IH1 - EPS(IV)=E(IV)*F1(JJJ(IV)) + EPS(IV)=E(IV)*F1(JJJ(IV)) 67 CONTINUE C CDIR$ IVDEP DO 64 IV=1,IH1 -C +C C RANDOM AZIMUTHAL ANGLE AND IMPACT PARAMETER C CC PHIP=PI2*RANF() CC PHIP=PI2*DRAND48() CC PHIP=PI2*DBLE(RAN(ISEED)) - call ranlux(ran2, 2) - PHIP=PI2*DBLE(ran2(1)) - CPHI(IV)=DCOS(PHIP) - SPHI(IV)=DSIN(PHIP) + call ranlux(ran2, 2) + PHIP=PI2*DBLE(ran2(1)) + CPHI(IV)=DCOS(PHIP) + SPHI(IV)=DSIN(PHIP) CC P(IV)=PDMAX(LLL(IV))*DSQRT(RANF()+KK) CC P(IV)=PDMAX(LLL(IV))*DSQRT(DRAND48()+KK) CC P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(RAN(ISEED))+KK) - P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(ran2(2))+KK) + P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(ran2(2))+KK) C C POSITION OF TARGET ATOM C - X1(IV)=X(IV)-P(IV)*CPHI(IV)*SX(IV) - P(IV)=CVMGT(1.D10,P(IV),X1(IV).LT.0.D0.OR.X1(IV).GT.TT) + X1(IV)=X(IV)-P(IV)*CPHI(IV)*SX(IV) + P(IV)=CVMGT(1.D10,P(IV),X1(IV).LT.0.D0.OR.X1(IV).GT.TT) C IF(A1(JJJ(IV)).EQ.0.) WRITE(99,'(A50)')' A1 vor Label 64 ' - B(IV)=P(IV)/A1(JJJ(IV)) + B(IV)=P(IV)/A1(JJJ(IV)) 64 CONTINUE CALL SCOPY(IH1,B,1,R,1) C WRITE(99,*)IH1,B(IV),R(IV) @@ -1133,10 +1089,10 @@ C KRYPTON-CARBON POTENTIAL C CALL MAGICKRC(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 104 DO 105 IV=IVMIN,IVMAX - IF(R(IV).LT.1.D-20)THEN + IF(R(IV).LT.1.D-20)THEN WRITE(99,'(A70)')'in DO 104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' R(IV)=0.00001D0 - ENDIF + ENDIF EX1(IV)=DEXP(-.278544D0*R(IV)) EX2(IV)=DEXP(-.637174D0*R(IV)) EX3(IV)=DEXP(-1.919249D0*R(IV)) @@ -1159,14 +1115,14 @@ C GET MAX AND MIN INDEX OF TEST FAILURES IF(IVMIN.GT.IVMAX) GO TO 106 GO TO 104 106 DO 108 IV=1,IH1 - ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) - SQE=DSQRT(DABS(EPS(IV))) - CC=(.235809D0+SQE)/(.126000D0+SQE) - AA=2.D0*EPS(IV)*(1.D0+(1.0144D0/SQE))*B(IV)**CC - FF=(DSQRT(AA*AA+1.)-AA)*((69350.D0+EPS(IV))/(83550.D0+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) + ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) + SQE=DSQRT(DABS(EPS(IV))) + CC=(.235809D0+SQE)/(.126000D0+SQE) + AA=2.D0*EPS(IV)*(1.D0+(1.0144D0/SQE))*B(IV)**CC + FF=(DSQRT(AA*AA+1.)-AA)*((69350.D0+EPS(IV))/(83550.D0+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) 108 S2(IV)=1.D0-(1.D0*C2(IV)) GO TO 4103 C @@ -1174,10 +1130,10 @@ C C MOLIERE POTENTIAL C CALL MAGICMOL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 4104 DO 4105 IV=IVMIN,IVMAX - IF(R(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 4104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - R(IV)=0.00001D0 - ENDIF + IF(R(IV).LT.1.D-20)THEN + WRITE(99,'(A70)')'in DO 4104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + R(IV)=0.00001D0 + ENDIF EX1(IV)=DEXP(-.3D0*R(IV)) EX2(IV)=DEXP(-1.2D0*R(IV)) EX3(IV)=DEXP(-6.0D0*R(IV)) @@ -1199,14 +1155,14 @@ C GET MAX AND MIN INDEX OF TEST FAILURES GO TO 4104 4106 DO 4108 IV=1,IH1 C IF((EPS(IV)-V(IV)).EQ.0.D0)WRITE(99,'(A50)')' nach Label 4106 ' - ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) - SQE=DSQRT(EPS(IV)) - CC=(.009611D0+SQE)/(.005175D0+SQE) - AA=2.D0*EPS(IV)*(1.D0+(0.6743D0/SQE))*B(IV)**CC - FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPS(IV))/(10.D0+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) + ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) + SQE=DSQRT(EPS(IV)) + CC=(.009611D0+SQE)/(.005175D0+SQE) + AA=2.D0*EPS(IV)*(1.D0+(0.6743D0/SQE))*B(IV)**CC + FF=(DSQRT(AA*AA+1.D0)-AA)*((6.314D0+EPS(IV))/(10.D0+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) 4108 S2(IV)=1.D0-(1.D0*C2(IV)) GO TO 4103 C @@ -1214,7 +1170,7 @@ C C ZBL POTENTIAL C CALL MAGICZBL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 5104 DO 5105 IV=IVMIN,IVMAX - IF(R(IV).LT.1.D-20)THEN + 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 @@ -1257,68 +1213,68 @@ C C END OF MAGIC C DO 65 IV=1,IH1 - DEN(IV)=EC1(JJJ(IV))*E(IV)*S2(IV) + 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 + 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)) + 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))) - CT(IV)=C2(IV)+C2(IV)-1.D0 - ST(IV)=DSQRT(DABS(1.D0-CT(IV)*CT(IV))) - CU=CT(IV)+MU1(JJJ(IV)) - CU=CVMGT(CU,1.0D-8,DABS(CU).GE.1.0D-8) - TA=ST(IV)/CU - TA2=1.D0/DSQRT(DABS(1.D0+TA*TA)) - 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, + 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))) + CT(IV)=C2(IV)+C2(IV)-1.D0 + ST(IV)=DSQRT(DABS(1.D0-CT(IV)*CT(IV))) + CU=CT(IV)+MU1(JJJ(IV)) + CU=CVMGT(CU,1.0D-8,DABS(CU).GE.1.0D-8) + TA=ST(IV)/CU + TA2=1.D0/DSQRT(DABS(1.D0+TA*TA)) + 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) - DENS(IV)=DENS(IV)+DEN(IV) - DEES(IV)=DEES(IV)+DEEOR + DENS(IV)=DENS(IV)+DEN(IV) + DEES(IV)=DEES(IV)+DEEOR 65 CONTINUE 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) - 245 CONTINUE + 245 CONTINUE C C END OF COLLISION LOOP C C INELASTIC ENERGY LOSS( 5 POSSIBILITIES) C DO 14 IV=1,IH1 - ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV)) - TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV)) + ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV)) + TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV)) 14 CONTINUE 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)), + DEE(IV)=CVMGT(0.D0,KLM1(LLL(IV))*ASIGT(IV)*DSQRT(E(IV)), # X(IV).LT.HLM.OR.X(IV).GT.HLMT) 26 CONTINUE GO TO 40 16 DO 21 IV=1,IH1 - DEE(IV)=DEES(IV) + DEE(IV)=DEES(IV) 21 CONTINUE GO TO 40 17 DO 22 IV=1,IH1 - DEE(IV)=CVMGT(DEES(IV),0.5D0*(KLM1(LLL(IV))*ASIGT(IV)* + 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 GO TO 40 18 DO 23 IV=1,IH1 - SM(IV)=0.D0 - EM(IV)=E(IV)*0.001D0/M1 + 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)) + 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))) @@ -1329,112 +1285,113 @@ C 66 CONTINUE DO 73 IV=1,IH1 DO 73 J=1,NJ(LLL(IV)) - SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) + SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) 73 CONTINUE - DO 78 IV=1,IH1 - DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10.D0) - 78 CONTINUE - DO 69 IV=1,IH1 - DEE(IV)=10.D0*ASIGT(IV)* - # CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT) - 69 CONTINUE + DO IV=1,IH1 + DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10 + & .D0) + ENDDO + DO IV=1,IH1 + DEE(IV)=10.D0*ASIGT(IV)*CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR + & .X(IV).GT.HLMT) + ENDDO GO TO 40 - 19 FHE=CVMGT(1.3333D0,1.D0,M1.LT.4.00D0) - DO 25 IV=1,IH1 - SM(IV)=0.D0 - EM(IV)=E(IV)*0.001D0*FHE - 25 CONTINUE - DO 74 IV=1,IH1 - DO 74 J=1,NJ(LLL(IV)) - SH(IV,J)=CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)* - # (CH3(LLL(IV),J)/(EM(IV)*0.001D0)) - # *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/(EM(IV)*0.001D0))+ - # CH5(LLL(IV),J)*EM(IV)*0.001D0)) - # /(CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)+ - # (CH3(LLL(IV),J)/(EM(IV)*0.001D0)) - # *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/(EM(IV)*0.001D0))+ - # CH5(LLL(IV),J)*EM(IV)*0.001D0))) - 74 CONTINUE - DO 92 IV=1,IH1 - DO 92 J=1,NJ(LLL(IV)) - SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) - 92 CONTINUE - DO 79 IV=1,IH1 - DEE(IV)=10.D0*ASIGT(IV)* - # CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT) - 79 CONTINUE - 40 CONTINUE + 19 FHE=CVMGT(1.3333D0,1.D0,M1.LT.4.00D0) + DO IV=1,IH1 + SM(IV)=0.D0 + EM(IV)=E(IV)*0.001D0*FHE + ENDDO + DO IV=1,IH1 + DO J=1,NJ(LLL(IV)) + SH(IV,J)=CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)* (CH3(LLL(IV) + & ,J)/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J) + & /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0)) + & /(CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)+ (CH3(LLL(IV),J + & )/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J) + & /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0))) + ENDDO + ENDDO + DO IV=1,IH1 + DO J=1,NJ(LLL(IV)) + SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J) + ENDDO + ENDDO + DO IV=1,IH1 + DEE(IV)=10.D0*ASIGT(IV)* CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV + & ).GT.HLMT) + ENDDO + 40 CONTINUE C - DO 44 IV=1,IH1 - DEL=DMAX1(1.0D-20,DENS(IV)+DEE(IV)) - DENS(IV)=CVMGT(E(IV)*DENS(IV)/DEL,DENS(IV),DEL.GT.E(IV)) - DEE(IV)=CVMGT(E(IV)*DEE(IV)/DEL,DEE(IV),DEL.GT.E(IV)) - 44 CONTINUE + DO IV=1,IH1 + DEL=DMAX1(1.0D-20,DENS(IV)+DEE(IV)) + DENS(IV)=CVMGT(E(IV)*DENS(IV)/DEL,DENS(IV),DEL.GT.E(IV)) + DEE(IV)=CVMGT(E(IV)*DEE(IV)/DEL,DEE(IV),DEL.GT.E(IV)) + ENDDO C C INCREMENT OF DAMAGE, CASCADE AND PHONON ENERGY C DO 70 IV=1,IH1 C IF(X(IV).LT.0.OR.X(IV).GT.TT) GO TO 70 - I=MAX0(MIN0(IDINT(X1(IV)/CW+1.D0),100),1) - DENT(I)=DENT(I)+DENS(IV) - DMGN(I)=DMGN(I)+DEN(IV) - ION(I)=ION(I)+DEE(IV) - ELE(I,JJJ(IV))=ELE(I,JJJ(IV))+DEN(IV) - ELI(I,JJJ(IV))=ELI(I,JJJ(IV))+DEE(IV) - IF(DEN(IV).LE.DI(JJJ(IV))) GO TO 28 - EPS(IV)=F1(JJJ(IV))*DEN(IV) - G=EPS(IV)+.40244D0*EPS(IV)**.75D0+3.4008D0*EPS(IV)**.16667D0 - MOT=DEN(IV)/(1.D0+K2(LLL(IV))*G) - CASMOT(I)=CASMOT(I)+MOT - ELGD(I)=ELGD(I)+DEN(IV) - ELD(I,JJJ(IV))=ELD(I,JJJ(IV))+DEN(IV) - ICD(I,JJJ(IV))=ICD(I,JJJ(IV))+1 - GO TO 70 - 28 PHON(I)=PHON(I)+DEN(IV) - ELP(I,JJJ(IV))=ELP(I,JJJ(IV))+DEN(IV) + I=MAX0(MIN0(IDINT(X1(IV)/CW+1.D0),100),1) + DENT(I)=DENT(I)+DENS(IV) + DMGN(I)=DMGN(I)+DEN(IV) + ION(I)=ION(I)+DEE(IV) + ELE(I,JJJ(IV))=ELE(I,JJJ(IV))+DEN(IV) + ELI(I,JJJ(IV))=ELI(I,JJJ(IV))+DEE(IV) + IF(DEN(IV).LE.DI(JJJ(IV))) GO TO 28 + EPS(IV)=F1(JJJ(IV))*DEN(IV) + G=EPS(IV)+.40244D0*EPS(IV)**.75D0+3.4008D0*EPS(IV)**.16667D0 + MOT=DEN(IV)/(1.D0+K2(LLL(IV))*G) + CASMOT(I)=CASMOT(I)+MOT + ELGD(I)=ELGD(I)+DEN(IV) + ELD(I,JJJ(IV))=ELD(I,JJJ(IV))+DEN(IV) + ICD(I,JJJ(IV))=ICD(I,JJJ(IV))+1 + GO TO 70 + 28 PHON(I)=PHON(I)+DEN(IV) + ELP(I,JJJ(IV))=ELP(I,JJJ(IV))+DEN(IV) 70 CONTINUE DO 80 IV=1,IH1 - ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV)))) - ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1))) - ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0)) + ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV)))) + ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1))) + ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0)) 80 CONTINUE DO 72 IV=1,IH1 - DEN2=DEN(IV)*DEN(IV) - DEN3=DEN2*DEN(IV) - EEL=EEL+DEN(IV) - EEL2=EEL2+DEN2 - EEL3=EEL3+DEN3 - EEL4=EEL4+DEN2*DEN2 - EEL5=EEL5+DEN3*DEN2 - EEL6=EEL6+DEN3*DEN3 - DEE2=DEE(IV)*DEE(IV) - DEE3=DEE2*DEE(IV) - EIL=EIL+DEE(IV) - EIL2=EIL2+DEE2 - EIL3=EIL3+DEE3 - EIL4=EIL4+DEE2*DEE2 - EIL5=EIL5+DEE3*DEE2 - EIL6=EIL6+DEE3*DEE3 - EPL=EPL+CVMGT(DEN(IV),0.D0,DEN(IV).LT.DI(JJJ(IV))) - EPL2=EPL2+CVMGT(DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) - EPL3=EPL3+CVMGT(DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV))) - EPL4=EPL4+CVMGT(DEN2*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) - EPL5=EPL5+CVMGT(DEN3*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) - EPL6=EPL6+CVMGT(DEN3*DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV))) - ENUCL(IV)=ENUCL(IV)+DENS(IV) - EINEL(IV)=EINEL(IV)+DEE(IV) + DEN2=DEN(IV)*DEN(IV) + DEN3=DEN2*DEN(IV) + EEL=EEL+DEN(IV) + EEL2=EEL2+DEN2 + EEL3=EEL3+DEN3 + EEL4=EEL4+DEN2*DEN2 + EEL5=EEL5+DEN3*DEN2 + EEL6=EEL6+DEN3*DEN3 + DEE2=DEE(IV)*DEE(IV) + DEE3=DEE2*DEE(IV) + EIL=EIL+DEE(IV) + EIL2=EIL2+DEE2 + EIL3=EIL3+DEE3 + EIL4=EIL4+DEE2*DEE2 + EIL5=EIL5+DEE3*DEE2 + EIL6=EIL6+DEE3*DEE3 + EPL=EPL+CVMGT(DEN(IV),0.D0,DEN(IV).LT.DI(JJJ(IV))) + EPL2=EPL2+CVMGT(DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) + EPL3=EPL3+CVMGT(DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV))) + EPL4=EPL4+CVMGT(DEN2*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) + EPL5=EPL5+CVMGT(DEN3*DEN2,0.D0,DEN(IV).LT.DI(JJJ(IV))) + EPL6=EPL6+CVMGT(DEN3*DEN3,0.D0,DEN(IV).LT.DI(JJJ(IV))) + ENUCL(IV)=ENUCL(IV)+DENS(IV) + EINEL(IV)=EINEL(IV)+DEE(IV) 72 CONTINUE IF(KK0.EQ.0) GO TO 89 DO 71 IV=1,IH1 - DEWC=DENS(IV)-DEN(IV) - DEWC2=DEWC*DEWC - DEWC3=DEWC2*DEWC - EELWC=EELWC+DEWC - EELWC2=EELWC2+DEWC2 - EELWC3=EELWC3+DEWC3 - EELWC4=EELWC4+DEWC2*DEWC2 - EELWC5=EELWC5+DEWC3*DEWC2 - EELWC6=EELWC6+DEWC3*DEWC3 + DEWC=DENS(IV)-DEN(IV) + DEWC2=DEWC*DEWC + DEWC3=DEWC2*DEWC + EELWC=EELWC+DEWC + EELWC2=EELWC2+DEWC2 + EELWC3=EELWC3+DEWC3 + EELWC4=EELWC4+DEWC2*DEWC2 + EELWC5=EELWC5+DEWC3*DEWC2 + EELWC6=EELWC6+DEWC3*DEWC3 71 CONTINUE 89 CONTINUE C