SUB RSTOP(Z1,M1,Z2,EE,SE(1),ENUMB,EFLAG) '#### If EFLAG=0, then return 1000 stopping values up to energy EE. '#### If EFLAG=1, then use ENUMB energies hidden in SE. LOCAL EKEV(),PCOEF() DIM DYNAMIC EKEV(1000), PCOEF(8) '********************************************************************* '**** "RSTOP" : ION STOPPING CROSS-SECTION SUBROUTINE USED BY TRIM * '**** VERSION : 1987 (C) 1985, J.F.ZIEGLER, YORKTOWN,NY,10598 * '**** Converted from FORTRAN to BASIC G.Cuomo J.Ziegler 5-1987 * '********************************************************************* '******** (OUTPUT IS 1000 VALUES OF SE (EV-ANG.2)) ******** '******** (NUCLEAR STOPPING IS NOT CALCULATED.) ******** '******************************************************************* '******** Z1 = ION ATOMIC NUMBER ******** '******** MM1= ION ATOMIC MASS ******** '******** M1 = ION ATOMIC WEIGHT (AMU) ******** '******** Z2 = TARGET ATOMIC NUMBER ******** '******** M2 = TARGET ATOMIC WEIGHT (AMU) ******** '******** EE = ION ENERGY (KEV) ******** '******** RHO = TARGET DENSITY (G/CM3) ******** '******** ATRHO=TARGET DENSITY (ATOMS/CM3) ******** '******** VFERMI = (FERMI VELOCITY OF SOLID) / V0 ******** '******** LFCTR = LAMBDA SCREENING FACTOR FOR IONS ******** '******** PCOEF = STOPPING COEFFICIENTS FOR PROTONS ******** '******** SE = CALCULATED ELECTRONIC STOPPING (EV-A2) ******** '******************************************************************* '******************************************************************* '******** NEED DATA FILE : SCOEF DATA (184,80) ******** '******** SCOEF(*,1) = ATOMIC NUMBER OF ELEMENT ******** '******** SCOEF(*,2) = ATOMIC MASS OF MOST ABUNDANT ISOTOPE ***** '******** SCOEF(*,3) = ATOMIC WEIGHT OF M.A.I. ******** '******** SCOEF(*,4) = TARGET MASS (AMU) ******** '******** SCOEF(*,5) = G/CM3 DENSITY (NORMAL ABUNDANCE) ******** '******** SCOEF(*,6) = ATOMS/CM3*1E22 (NORMAL ABUNDANCE) ******** '******** SCOEF(*,7) = (FERMI VELOCITY OF SOLID) / V0 ******** '******** SCOEF(*,8) = LAMBDA SCREENING FACTOR FOR IONS ******** '******** SCOEF(92,*) TO SCOEF(184,*) = COEF.OF PROTON STOPPING * '******************************************************************* '******************************************************************* '----------------------------------------------------------------------------- ' CALL DATA FILE: SCOEF DATA TO GET VALUES FOR STOPPING CALCULATION ' FIRST CALL IS FOR DATA ON THE ION, THE SECOND IS ABOUT THE TARGET. ' LFACTOR IS THE SCREENING LENGTH OF THE ION, VFERMI IS OF THE TARGET. CALL SCOEF(Z1,Y1,MM1,Y2,Y3,Y4 ,Y5 ,LFCTR,PCOEF()) CALL SCOEF(Z2,Y1,Y2 ,Y3,Y4,ATRHO,VFERMI,Y5 ,PCOEF()) IF M1=0 THEN M1=MM1 ' MM1 = Mass (amu) of most abundant isotope. '----------------------------------------------------------------------------- IF EFLAG=0 GOTO RS00 'EFLAG=0 (get 1000 Se), EFLAG=1 (ENUMB values only) EMAX=0 FOR I=1 TO ENUMB '---- Use Preset Energies (keV) EKEV(I)=SE(I)/M1 ' Energies passed in array SE. IF EMAX < EKEV(I) THEN EMAX = EKEV(I) 'Maximum energy of array. NEXT I GOTO RS0 RS00: ENUMB=1000 '--- 1000 Equally Spaced Energies EMAX=EE/M1 ' Maximum Energy of array E=EE/1000/M1 FOR I=1 TO 1000 EKEV(I)=E*I 'Array of Energies NEXT I '------------------------------------------------------------------------ '--------------- High Energy Stopping Powers ( 10 MeV - 2 GeV ) ------------- RS0: IF EMAX < 10000 GOTO RS1 OPEN "HICOEF.DAT" FOR INPUT AS #4 'High Energy Stopping Powers FOR I=1 TO Z2 INPUT#4,A$ NEXT I CLOSE #4 HI1=VAL(MID$(A$,1 ,10)):HI2=VAL(MID$(A$,11,10)) HI3=VAL(MID$(A$,21,11)):HI4=VAL(MID$(A$,32,11)) '----------------------------------------------------------------------------- RS1: '------- LOOP to get all Stopping values. FOR I=1 TO ENUMB IF (I MOD 50) <> 0 GOTO RS50 LOCATE 18,50: PRINT I 'Print out every 50th Calcualtion. RS50: E=EKEV(I) '----- Branch for Stopping of : Z1=1, Z1=2 OR Z1>2 . IF (Z1=1) THEN GOSUB RS100: GOTO RS10 IF (Z1=2) THEN GOSUB RS200: GOTO RS10 IF (Z1>2) THEN GOSUB RS300: GOTO RS10 RS10: '########################## SE(I)=SE(I)*10 'Stopping in eV/(1E15-cm2). This converts to eV-A2 NEXT I GOTO RS999 'End Subroutine '********** Proton Stopping Powers ****************************** RS100: CALL RPSTOP(Z2,E,PCOEF(),SE(I),HI1,HI2,HI3,HI4) RETURN '******************************************************************* '******* HELIUM Stopping Powers ****************************** '---- VELOCITY PROPORTIONAL STOPPING BELOW KEV/AMU >> HE0 <<. RS200: HE0=1 IF HE0>E THEN HE=HE0 ELSE HE=E B=LOG(HE) A=.2865+.1266*B-.001429*B*B+.02402*B*B*B-.01135*B^4+.001475*B^5 IF A>30 THEN A=30 HEH=1.-EXP(-A) '**** ADD Z1^3 EFFECT TO HE/H STOPPING POWER RATIO** HEH **. IF HE<1 THEN HE=1 A=(1.+(.007+.00005*Z2)*EXP(-(7.6-LOG(HE))^2)) HEH=HEH*A*A CALL RPSTOP(Z2,HE,PCOEF(),SP,HI1,HI2,HI3,HI4) SE(I)=SP*HEH*4. IF E>HE0 GOTO RS299 '*** CALC. HE VELOCITY PROPORTIONAL STOPPING SE(I)=SE(I)*SQR(E/HE0) RS299: RETURN '******************************************************************* '********** HEAVY ION ELECTRONIC STOPPING POWERS ******************* '**** USE VELOCITY STOPPING FOR (YRMIN=VR/Z1**.67) <= 0.13 '**** OR FOR VR <= 1.0 RS300: YRMIN=0.13 VRMIN=1.0 V=SQR(E/25)/VFERMI 301 IF V>=1 GOTO 302 VR=(3*VFERMI/4)*(1+(2*V*V/3)-V^4/15) GOTO 303 302 VR=V*VFERMI*(1+1/(5*V*V)) '**** SET YR = MAXIMUM OF (VR/Z1^.67),(VRMIN/Z1^.67) OR YRMIN. ''303 YR=FNAMAX1(YRMIN,VR/Z1^.6667) 303 YR=VR/Z1^.6667: IF YR50 THEN A=50 'Prevents Underflow Q=1-EXP(-A): IF Q<0 THEN Q=0 ELSE IF Q>1 THEN Q=1 '**** Q = IONIZATION LEVEL OF THE ION AT VELOCITY * YR *. '**** NOW WE CONVERT IONIZATION LEVEL TO EFFECTIVE CHARGE. ''''''LL=(FNAMIN1(0.43,FNAMAX1(.32,.12+.025*Z1)))/Z1^.3333 LL=.12+.025*Z1: IF LL<.32 THEN LL=.32 ELSE IF LL>.43 THEN LL=.43 LL=LL/Z1^.3333 '---- The following calculates the screening length of the ion, using '---- Brandt & Kitagawa eq.(15) with constant a = .343. '---- It includes s-shell correction shown in ZBL, figure (3-24). ''''''L0=(.8-Q*(FNAMIN1(1.2,0.6+Z1/30.)))/Z1^.3333 A=0.6+Z1/30.: IF A>1.2 THEN A=1.2 L0=(.8-Q*A)/Z1^.3333 IF Q<0.2 GOTO 307 ''''''IF (Q<(FNAMAX1(0.,.9-.025*Z1))) GOTO 306 A=.9-.025*Z1: IF A<0 THEN A=0 IF Q16 THEN A=16 B=1-.025*A: IF B<0 THEN B=0 IF Q < B GOTO 305 ''304 L1=LL*(1-Q)/(.025*FNAMIN1(16.,1*Z1)) 304 A=Z1: IF A>16 THEN A=16 L1=LL*(1-Q)/(.025*A) GOTO 308 305 L1=LL GOTO 308 306 Q1=0.2 ''''' L1=LL*(Q-.2)/ABS(FNAMAX1(0,.9-.025*Z1)-.2000001) A=.9-.025*Z1: IF A<0 THEN A=0 L1=LL*(Q-.2)/ABS(A-.2000001) GOTO 308 307 L1=0. ''308 L=FNAMAX1(L1,L0*LFCTR) 308 L=L0*LFCTR: IF L1E4) GOTO RP20 '---- High Energy Stopping '**** VELOCITY PROPORTIONAL STOPPING BELOW VELOCITY ** PE0 **. PE0=25. ''''''PE=FNAMAX1(PE0,E) PE=PE0: IF PE < E THEN PE=E SL=(PCOEF(1)*PE^PCOEF(2))+PCOEF(3)*PE^PCOEF(4) SH=PCOEF(5)/PE^PCOEF(6)*LOG((PCOEF(7)/PE)+PCOEF(8)*PE) SP=SL*SH/(SL+SH) IF (E>PE0) GOTO RP10 '**** VELPWR IS THE POWER OF VELOCITY STOPPING BELOW PE0. VELPWR=0.45 IF (Z2<=6) THEN VELPWR=0.25 SP=SP*(E/PE0)^VELPWR GOTO RP10 RP20: X=LOG(E)/E '********* High Energy Stopping (6/87) SP=HI1+(HI2*X)+(HI3*X*X)+(HI4/X) RP10: END SUB '***************************************************************** '***************************************************************** '********** GET STOPPING COEFFICIENTS ****************************** SUB SCOEF(Z1,MM1,M1,M2,RHO,ATRHO,VFERMI,LFCTR,PCOEF(1)) LOCAL C$(),I DIM C$(11) 'When make up SCOEF, need to add 2 extra chars. for E2 data files !! OPEN "SCOEF.88" AS #1 LEN=88 FIELD #1,4 AS C$(1),5 AS C$(2),9 AS C$(3),9 AS C$(4),10 AS C$(5),8 AS C$(6),9 AS C$(7),7 AS C$(8),5 AS C$(9) GET#1,Z1+2 ' Get Record from data file, Add 2 for text in data file Z1=VAL(C$(1)) :MM1=VAL(C$(2)) :M1=VAL(C$(3)) :M2=VAL(C$(4)) RHO=VAL(C$(5)) :ATRHO=VAL(C$(6)) :VFERMI=VAL(C$(7)) :LFCTR=VAL(C$(8)) SUBLIM=VAL(C$(9)) 'For proton stopping coefficients, add 2+1 for text in data file. FIELD #1,2 AS C$(1),11 AS C$(2),12 AS C$(3),10 AS C$(4),10 AS C$(5),10 AS C$(6),10 AS C$(7),10 AS C$(8),11 AS C$(9) GET#1,Z1+92+2+1 FOR I=2 TO 9: PCOEF(I-1)=VAL(C$(I)): NEXT I ATRHO= ATRHO*1E+22 CLOSE #1 END SUB