c **********************************************************************
                          PROGRAM TRIM                                  
c **********************************************************************
c **   TRIM-89 -- Ions in Solids (2-Dimensional, 4 Elements)          **
c **                                                                  **
c **       Original program by J. Biersack and J. Ziegler             **
c **                                                                  **
c **       This version rewritten & improved by P.May                 **
c **                                                                  **
c **   This version calculates substrate containing MIXED elements    **
c **   e.g. SiO2.  It no longer saves most of the data files, which   **
c **   have been commented-out, but only saves the INFO file.         **
c **   A plane potential barrier is assumed for sputtering            **
c **   whose value is given by the input file,  or assumed to be      **
c **   the same as the sublimation energy of the crystal.  This       **
c **   causes sputtered atoms to be refracted as they pass through    **
c **   the barrier, and so alters the ejection angle.                 **
c **                                                                  **
c **********************************************************************
                                                                        
c ....... Z(x), Zmass(x) = Atomic Number and Mass of elements ...........
c ...... with Z(1) and Zmass(1) for the incident ion ...................
c ...... Z2, M2, rho(g/cm3) and ATMrho are average values ..............
c ...... of atomic number, mass & density for target ...................
                                                                        
c ...... Every collision starts with ion at (x,y)=(0,0) with ...........
c ...... energy, E, and with impact parameter B ........................
c ...... Results of Magic Formula for collision are: ...................
c ...... S2 = Sin2(theta) of deflection angle theta ....................
c ...... DEN= Energy transferred during collision ......................
c ...... DEE=Electronic dE/dx ..........................................
                                                                        
c ...... Final results are in arrays: ..................................
                                                                        
c ....... RANGE(1)    =    Ion Distribution ............................
c ....... RANGE(2-5)  =    Target atom redistribution ..................
c ....... YRANGE(1)   =    Ion lateral straggling ......................
c ....... YRANGE(2-5) =    Target atom straggling ......................
c ....... xps         =    Ion x coordinate (depth into surface) .......
c ....... yps         =    Ion y coordinate (straggling) ...............
c ....... VAC(1)      =    Vacancies created by Ion ....................
c ....... VAC(2-5)    =    Vacancies of target atom #2-#5 ..............
c ....... XNOVAC      =    Replacement Collisions ......................
c ....... E2REC       =    Energy transferred to Recoils ...............
c ....... IONIZ(2)    =    Recoils energy loss by ionization ...........
c ....... PHONON(2)   =    Recoils loss by phonon creation .............
c ....... bke         =    Backscattered ion energy /eV ................
c ....... bka         =    Backscattered ion angle /degs ...............
c ....... bkx         =    Backscattered ion velocity in x dir in m/s ..
c ....... bky         =    Backscattered ion velocity in y dir in m/s ..
c ....... bkz         =    Backscattered ion velocity in z dir in m/s ..
c ....... tme         =    Transmitted ion energy /eV ..................
c ....... tma         =    Transmitted ion angle /degs .................
c ....... tmx         =    Transmitted ion velocity in x dir in m/s ....
c ....... tmy         =    Transmitted ion velocity in y dir in m/s ....
c ....... tmz         =    Transmitted ion velocity in z dir in m/s ....
c ....... spe         =    Sputtered atom energy /eV ...................
c ....... spa         =    Sputtered atom angle /degs ..................
c ....... spx         =    Sputtered atom velocity in x dir in m/s .....
c ....... spy         =    Sputtered atomvelocity in y dir in m/s ......
c ....... spz         =    Sputtered atom velocity in z dir in m/s .....
                                                                        
c .......... Other important parameters ................................
                                                                        
c ....... A=Screening length; F=LSS Epsilon/Energy; EC=Recoil factor....
c ....... rho=Density(g/cm3); ATMrho=(atoms/Ang3) ......................
c ....... natom=Total elements (ion & atoms) ; xdepth=target depth .....
c ....... Zmass2=<mass>; Z2=<Z2> .......................................
c ....... FFPATH=Free flight path; ATdist=Distance between atoms .......
c ....... P=Impact parameter; Pmax=Maximum P........................... 
c ....... SE=Electronic stopping (eV/Ang) ..............................
c **********************************************************************
                                                                        
c ............... Define variables .....................................
                                                                        
      IMPLICIT NONE                                                     
                                                                        
      LOGICAL ext                                                              
                                                                        
      INTEGER i,j,iZ,ionstop,backstop,sputstop,tmstop,counter,          
     &        ntransmit,natom,ionnum,Ion,NTGT,iZ1,iZ2,Nchan,            
     &        ncoll,ncollsum,nback,nsput1,IE,n,nsput(2:5),             
     &        nimplant,Flag1,J1,JJ1,II,Label(2000),seed           
c             Nchan1,Dchan,LLL,Iondir                                                                        

      CHARACTER Zname$(92)*2,filenm$*8,fileid$*21,dummy$*50, 
     &          layername$*10,Z$(5)*2,spError$*50,tmerror$*50,          
     &          bkError$*50                                            
                                                                        
      REAL Z(5),Zmass(5),YY(8),M2M1(5,5),F(5,5),A(5,5)                  
c     REAL BSenergy(100000),BSangle(100000)
c                                        ,BSlateral(100000)              
c     REAL BSxvel(100000),BSyvel(100000),BSzvel(100000)                    
c     REAL TMenergy(100000),TMangle(100000)
c                                        ,TMlateral(100000)              
c     REAL TMxvel(100000),TMyvel(100000),TMzvel(100000)                    
c     REAL SPenergy(2:5,100000),SPangle(2:5,100000)
c                                        ,SPlateral(100000)              
c     REAL SPxvel(100000),SPyvel(100000),SPzvel(100000)                    
      REAL ER(2000),COSXR(2000),COSYR(2000),xr(2000),yr(2000)           
c     REAL VAC(5,200),PHONON(2,200),RANGE(5,200),YRANGE(5,200)          
c     REAL XNOVAC(2000),EIONIZ(2,200),E2REC(5,200)                      
      REAL SE(1000,5),stoich(5),SE0(1000),Erec(5),EC(5,5)               
c     REAL xpos(100000),ypos(100000)                                      
                                                                        
      REAL angle,ALPHA,fuzz,Tmin,Zmass2,Z2,ATMrho,rho,ATdist, 
     &     Pmax,M1M2,VFermi,LFCTR,EE,Bragg,Esurf,EB,ED,EF,EDlow,   
     &     B,S2,C2,RN,R,RR,RZBL,E0,E0keV,FFPA,FFPF,epsDG,STbohr,KPelec, 
     &     KPvac,Epeak,E,COSX,COSY,SINX,SINY,pl,x,y,EPS,EEG,       
     &     P,Pmax1,FFPATH,PHI,SEE,DEE,STDEE,xold,IA,Eback,Esput,Etrans, 
     &     EX1,EX2,EX3,EX4,V,V1,FR,FR1,Q,ROC,SQE,CC,AA,FF,DELTA,C,CT,   
     &     ST,CU,PSI,DEN,DUM,xsum,x2sum,x3sum,x4sum,y2sum,xy2sum,
     &     x2y2sum,y4sum,plsum,pl2sum,N1,avex,avepl,avecol,varx,vary,   
     &     sigmay,sigma,xy2,x2y2,y4,vv,vv2,Gamma,Beta,sigpl,PSI2,CHI,   
     &     CPHI,SPHI,CPSI,SPSI,SINXR,SINYR,Enor,py,xp,sumstoich,DEER,   
     &     DENR,pi,xdepth,PHI2,R1,CU2,CHI2,qq,mp,Ionfactor,    
     &     dx,dy,dz,ran3
c     REAL CW,perang,Corr,VACTOT,INTERTOT,COSZ,Atomfactor                                          
                                                                        
c ............... Initialise variables .................................
                                                                        
      pi=3.1415926                                                      
      qq=1.602e-19                      ! Charge on electron /C         
      mp=1.67e-27                       ! Mass of a nucleon /kg         
      fuzz=1E-8                                                         
      counter=0                                                         
                                                                        
      DATA ionnum,nsput1,ntransmit,nback,nimplant,ncoll,ncollsum/7*0/    
      DATA xsum,x2sum,x3sum,x4sum,y2sum,xy2sum,x2y2sum,y4sum,plsum,     
     &     pl2sum,Z2,Zmass2/12*0.0/,SE/5000*0.0/,Z/5*0.0/,Zmass/5*0.0/, 
     &     YY/8*0.0/,M2M1/25*0.0/,F/25*0.0/,A/25*0.0/,                  
     &     Esput,Etrans,Eback/3*0.0/,nsput/4*0.0/                                    
                                                                        
      DATA spError$/' None'/,tmError$/' None'/,bkError$/' None'/        
                                                                        
      DATA (Zname$(i),i=1,92)/ 'H ','He','Li','Be','B ','C ','N ','O ', 
     &     'F ','Ne','Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca', 
     &     'Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge', 
     &     'As','Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru', 
     &     'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs','Ba', 
     &     'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er', 
     &     'Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg', 
     &     'Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U '/ 
                                                                        
c .................. Get input data from INPUT file ....................
                                                                        
c        PRINT*,'Please type filename - must be 8 characters long'                                  
        READ(*,'(A)')filenm$                                  
        filenm$=filenm$(1:8)                ! Ensures only 8 letters    
                                                                        
        fileid$=filenm$//'.inp'                                  
        OPEN(UNIT=13,FILE=fileid$,FORM='FORMATTED',STATUS='OLD')        
                                                                        
         READ(13,10)dummy$                  ! Dummy input for file title
10       FORMAT(A50)                                                    
                                                                        
         READ(13,11)Z$(1)                   ! Get chem symbol of ion    
11       FORMAT(A2)                                                     
                                                                        
         READ(13,*)E0keV                    ! Ion incident energy /eV
         E0kev=E0kev/1000                   ! Convert to keV  
         READ(13,*)angle                    ! Ion incident angle /degs  
         READ(13,*)Bragg                    ! Bragg correction coef     
                                                                        
         READ(13,*)ionstop                  ! No. ions before stop prog 
         READ(13,*)backstop                 ! No. backscat. before....  
         READ(13,*)sputstop                 ! No. sputtered before....  
         READ(13,*)tmstop                   ! No. transmitted before... 
                                                                        
         READ(13,10)dummy$                  ! Dummy for title           
         READ(13,13)layername$              ! Name of layer             
13       FORMAT(A10)                                                    
                                                                        
         READ(13,*)xdepth                   ! Depth of layer            
         READ(13,*)rho                      ! Density of layer g/cc     
         READ(13,*)Esurf                    ! Surface binding energy/eV 
         READ(13,*)EB                       ! Bulk binding energy/eV
         READ(13,*)ED                       ! Displacement energy/eV
         READ(13,14)Z$(2),stoich(2)         ! Atom #1 & percentage      
         READ(13,14)Z$(3),stoich(3)         ! Atom #2 & percentage      
         READ(13,14)Z$(4),stoich(4)         ! Atom #3 & percentage      
         READ(13,14)Z$(5),stoich(5)         ! Atom #4 & percentage      
14       FORMAT(A2,F5.2)                                                
         IF (stoich(3).EQ.0.0) THEN                                     
               natom=2                      ! No. atoms, ion + surface  
         ELSEIF (stoich(4).EQ.0.0) THEN                                 
               natom=3                                                  
         ELSEIF (stoich(5).EQ.0.0) THEN                                 
               natom=4                                                  
         ELSE                                                           
               natom=5                                                  
         ENDIF                                                          
                                                                        
       CLOSE(13)                                                        
                                                                        
c .... Get atomic numbers & mass of ion & atoms ........................
                                                                        
       DO 50,i=1,natom                                                  
          CALL GETATNO(Z$(i),iZ,Zname$)                                 
          Z(i)=FLOAT(iZ)                                                
          CALL SCOEF(iZ,X,Zmass(i),X,X,X,X,X,YY,X)                      
50     CONTINUE                                                         
       Ionfactor=SQRT(2.0*qq/(mp*Zmass(1)))      ! Converts E to vel.   
                                                                        
c ........ Get Binding E of surface for sputtering calculations ........
                                                                        
      IF (Esurf.le.0) THEN                       ! Haven't specified it 
         IF (natom.eq.2) THEN                    ! 1 atom-type in layer 
          i=INT(Z(2))                            ! Use Esublimation as  
          CALL SCOEF(i,X,X,X,X,X,X,X,Esurf,YY)   ! best approx to BE    
         ELSE                                                           
          Esurf=4.0                              ! Default if no better 
         ENDIF                                   ! data available       
      ENDIF                                                             
      IF (EB.lt.0) EB=0                          ! Binding E            
      IF (ED.lt.0) ED=20.0                       ! Displacment E        
      EF=Esurf                                   ! Final E              
      EDlow=Esurf                                ! Lowest E to do damage
                                                                        
c ................... Calc energy & angle of incident ion ..............
                                                                        
      E0=E0keV*1000                              ! Ion E in eV          
      ALPHA=angle*pi/180.                        ! Angle /rads          
                                                                        
c ........ Calculate channel width for file output .....................
c ........ Divide layer into 200 equal spaces (channels) for ...........
c ........ efficient filing & plotting .................................
                                                                        
c     CW=0.005*xdepth                            ! Channel width        
c     perang=1/CW                                ! Amount per channel   
                                                                        
c ....... Define minimum transferred E, this sets Pmax later ...........
                                                                        
      Tmin=5                                                            
                                                                        
c ........ Randomise seed for random number generation .................
                                                                        
      CALL tinow(seed)                                   
                                                                        
c ........... Get average atomic no. Z2 and mass, Zmass2 for surface ...
                                                                        
      DO 60,i=2,natom                                                   
        Zmass2=Zmass2+stoich(i)*Zmass(i)          ! Average mass        
        Z2=Z2+stoich(i)*Z(i)                      ! Average Z           
60    CONTINUE                                                          
                                                                        
      ATMrho=rho*0.6022/Zmass2                    ! Atomic density      
      ATdist=ATMrho**(-1/3)                       ! Atomic distance     
      Pmax=ATdist/SQRT(pi)                        ! Max impact parameter
      M1M2=Zmass(1)/Zmass2                        ! Ion:surface ratio   
                                                                        
c .......... Get 1000 values of stopping coefs for ion/atom ............
c .......... and atom/atom collisions ..................................
                                                                        
      DO 70,Ion=1,natom                                ! Ion/atom colls 
        DO 71,NTGT=2,natom                             ! Atom/atom colls
          iZ1=INT(Z(Ion))                                               
          iZ2=INT(Z(NTGT))                                              
                                                                        
          CALL  RSTOP(iZ1,Zmass(Ion),iZ2,E0keV,SE0,LFCTR,VFermi)        
                                                                        
          EE=SQRT(2*E0keV/1000/Zmass(1))                   ! Ion Vel.   
          DO 72,i=1,1000                                                
           SE(i,Ion)=SE(i,Ion)+SE0(i)*stoich(NTGT)*ATMrho               
c          IF (Bragg.ne.0) THEN                                         
c           Corr=1+(Bragg-1)/(1+EXP(1.48*(EE*SQRT(i-7))) !???           
c           SE(i,Ion)=SE(i,Ion)*Corr                       ! Compound   
c          ENDIF                                           ! stopping   
72        CONTINUE                                                      
                                                                        
71      CONTINUE                                                        
70    CONTINUE                                                          
                                                                        
c .............Calculate the mean flight path parameters .............. 
c .... STbohr=Bohr straggling, KPelec and KPvac are how much of ....... 
c .... recoil goes into electronic losses and what fraction into ...... 
c .... vacancies ...................................................... 
                                                                        
      FFPA=0.5292*0.8853/(Z(1)**0.23+Z2**0.23)      ! Univ. screening l.
      FFPF=FFPA*Zmass2/(Z(1)*Z2*14.4*(Zmass(1)+Zmass2))  ! Mean f. path 
      epsDG=Tmin*FFPF*(1+M1M2)**2/(4*M1M2)                              
      STbohr=177*Z(1)*SQRT(Z2*ATMrho)        ! Bohr straggling of dE/dx 
      KPelec=0.01*Z2**(-2.33333)                                        
      KPvac=0.1334*Z2**(0.666667)/SQRT(Zmass2)                          
                                                                        
c ....... Calculate collision parameters for atom/atom scattering ......
                                                                        
      DO 75,j=1,natom                                                   
        DO 76,i=2,natom                                                 
         M2M1(j,i)=Zmass(j)/Zmass(i)                                    
         EC(j,i)=4.0*M2M1(j,i)/(1+M2M1(j,i))**2                         
         A(j,i)=0.5292*0.8853/(Z(j)**0.23+Z(i)**0.23) ! Screening length
         F(j,i)=A(j,i)/(Z(j)*Z(i)*14.41*(1+M2M1(j,i)))                  
76      CONTINUE                                                        
75    CONTINUE                                                          
                                                                        
c ..... Peak of electronic stopping (used in straggling calculation) ...
                                                                        
      Epeak=1.e5*Z(1)**0.67*Zmass(1)                                    
                                                                        
c ......................................................................
                                                                        
c ######################################################################
c #################  Start of ion penetration loop  ####################
c ######################################################################
                                                                        
250   CONTINUE                          ! New ions start here           
      ionnum=ionnum+1                   ! Number of ions                
      E=E0                              ! Energy of ion in eV           
      COSX=COS(ALPHA)                   ! Velocity in x direction       
      SINX=SIN(ALPHA)                   ! Velocity in y direction       
      COSY=SINX                                                         
      SINY=COSX                                                         
                                                                        
      Nchan=1                           ! Channel number                
      pl=0.0                            ! Total path length of ion      
      ncoll=0                           ! No. collisions ion sees       
                                                                        
c ......... Calc initial x,y coords of ion's 1st impact ................
                                                                        
      RN=ran3(seed)                     ! Random number 0,1        
      x=ATdist*COSX*RN                                                  
      RN=ran3(seed)                                                     
      y=ATdist*COSY*RN                                                  
                                                                        
c -------------- Cycle through ion collisions --------------------------
                                                                        
280   CONTINUE                             ! Prog jumps back here after 
                                           ! each ion collision.        
c ......... Caculate free flight path of ion ...........................
                                                                        
      EPS=E*FFPF                                                        
      EEG=SQRT(EPS*epsDG)                         ! Reduced energy      
      Pmax1=FFPA/(EEG+SQRT(EEG)+0.125*EEG**0.1)   ! Max impact parameter
      RN=ran3(seed)                                                     
      P=Pmax1*SQRT(RN)                            ! Randomize imp. param
      FFPATH=1/(pi*Pmax1*Pmax1*ATMrho)            ! Free flight path    
                                                                        
c ........ Reduce FFPATH for 1st collision. We keep full ...............
c ........ Pmax since vacancies will be deposited only .................
c ........ at point of collision .......................................
                                                                        
      IF (ncoll.lt.1) THEN                     ! See if 1st collision   
        RN=ran3(seed)                                                   
        FFPATH=RN*FFPATH                       ! Random FFPATH          
      ENDIF                                                             
                                                                        
c ......... Calculate random azimuthal angle ...........................
                                                                        
      RN=ran3(seed)                                                     
      PHI=2.0*pi*RN                                                     
                                                                        
c .......Calculate Electronic energy loss for free flight path .........
                                                                        
      IE=INT((E-0.1)/E0keV)                                             
      IE=MIN0(IE,999)                             ! Ensures 1<IE<999    
      IE=MAX0(IE,1)                                                     
                                                                        
      SEE=SE(IE,1)+((E/E0keV)-IE)*(SE(IE+1,1)-SE(IE,1))                 
      IF (E.le.E0keV) SEE=SE(1,1)*SQRT(E/E0keV)                         
      DEE=FFPATH*SEE                                                    
                                                                        
c ..... Add Bohr straggling to energy loss of light ions ...............
                                                                        
      STDEE=0                                                           
      IF (Z(1).le.3) THEN                                               
       RN=ran3(seed)                                                    
       STDEE=(RN-0.5)*STbohr*SQRT(FFPATH)                               
       IF (E.le.Epeak) STDEE=STDEE*E/Epeak                              
       IF (ABS(STDEE).gt.ABS(DEE)) STDEE=DEE*STDEE/ABS(STDEE)           
      ENDIF                                                             
                                                                        
      E=E-DEE-STDEE                                                     
      IF (E.gt.E0) E=E0-fuzz                                            
      E=AMAX1(E,EF)                                                     
                                                                        
c .................. Get position of next collision ....................
                                                                        
      ncoll=ncoll+1                                 ! Update no. colls  
      pl=pl+FFPATH                                  ! Update path length
      xold=x                                                            
      x=x+FFPATH*COSX                                                   
      y=y+FFPATH*COSY                                                   
                                                                        
c ........ Calculate directions of ion velocity components .............
                                                                        
      dx=1.0                                                            
      dy=1.0                                                            
      IF (COSX.lt.0) dx=-1.0                                            
      IF (COSY.lt.0) dy=-1.0                                            
                                                                        
      dz=ran3(seed)                                                     
      IF (dz.lt.0.5) THEN                           ! Random direction  
          dz=-1.0                                   ! for z velocity    
      ELSEIF (dz.ge.0.5) THEN                                           
          dz=1.0                                                        
      ENDIF                                                             
                                                                        
c ................. Check if ion is backscattered/transmitted ..........
                                                                        
      IF (x.lt.0) THEN                              ! Backscattered     
                                                                        
       IA=(180/pi)*ABS(ATAN(SINX/(COSX+fuzz)))      ! Final angle/deg   
       nback=nback+1                                ! No. backscattered 
       n=nback                                                          
       IF (nback.gt.1000000) THEN                     ! Max array size    
         bkError$=' Array too large. Only 1000000 points saved.'          
         n=1000000                                                        
       ENDIF                                                            
       Eback=Eback+E                                ! Total E backscat. 
c      BSenergy(n)=E                                ! Energy            
c      BSangle(n)=IA                                ! Angle             
c      BSlateral(n)=y                               ! Lateral distance  
c      BSxvel(n)=dx*Ionfactor*SQRT(ABS(E*COSX))     ! x velocity        
c      BSyvel(n)=dy*Ionfactor*SQRT(ABS(E*COSY))     ! y velocity        
c      COSZ=SQRT(ABS(1-COSX**2-COSY**2))                                
c      BSzvel(n)=dz*Ionfactor*SQRT(E*COSZ)  ! z velocity                
       GOTO 660                                     ! End ion loop      
                                                                        
      ELSEIF (x.gt.xdepth) THEN                     ! Transmitted       
                                                                        
c      IA=(180/pi)*ABS(ATAN(SINX/COSX))             ! Final angle/deg   
       ntransmit=ntransmit+1                        ! No. transmitted   
       n=ntransmit                                                      
       IF (ntransmit.gt.1000000) THEN                 ! Max array size    
        tmError$=' Array too large. Only 1000000 points saved.'           
        n=1000000                                                         
       ENDIF                                                            
       Etrans=Etrans+E                              ! Total E transmit  
c      TMenergy(n)=E                                ! Energy            
c      TMangle(n)=IA                                ! Angle             
c      TMlateral(n)=y                               ! Lateral distance  
c      TMxvel(n)=dx*Ionfactor*SQRT(ABS(E*COSX))     ! x velocity        
c      TMyvel(n)=dy*Ionfactor*SQRT(ABS(E*COSY))     ! y velocity        
c      COSZ=SQRT(ABS(1-COSX**2-COSY**2))                                
c      TMzvel(n)=dz*Ionfactor*SQRT(E*COSZ)          ! z velocity        
       GOTO 660                                     ! End ion loop      
                                                                        
      ENDIF                                                             
                                                                        
c ... For next position, find target atom, J1, epsilon and B values ....
                                                                        
      Flag1=0                                                           
      sumstoich=0                                                       
      RN=ran3(seed)                                                     
      DO 360,J1=2,natom                                                 
       sumstoich=sumstoich+stoich(J1)                                   
       IF (RN.lt.sumstoich.and.Flag1.eq.0) THEN                         
         Flag1=J1                                                       
         EPS=E*F(1,J1)                           ! Reduced E of coll.   
         B=P/A(1,J1)                             ! Impact parameter     
       ENDIF                                                            
360   CONTINUE                                                          
      J1=Flag1                                                          
                                                                        
c ------------ Start scattering calculation ----------------------------
c .... We now have: epsilon (EPS) and B=Reduced impact parameter .......
c .... With these we calculate the scattering angle and the energy .....
c .... transferred to the target atom using the "Magic Formula" ........
c .... If epsilon > 10, use Rutherford scattering, else goto ...........
c .... Magic Formula. Notation:S2=Sin(theta)**2,CT=Cos(theta) ..........
c .... ST=Sin(theta) ...................................................
c ----------------------------------------------------------------------
                                                                        
      IF (EPS.gt.10) THEN                       ! Rutherford scattering 
         S2=1.0/(1.0+(1.0+B*(1.0+B))*(2.0*EPS*B)**2)                    
         C2=1.0-S2                                                      
      ELSE                                      ! Magic formula         
         R=B                                                            
         RR=-2.7*LOG(ABS(EPS*B)+fuzz)                                   
         IF (RR.ge.B) THEN                                              
           RR=-2.7*LOG(ABS(EPS*RR)+fuzz)                                
           IF (RR.ge.B) R=RR                                            
         ENDIF                                                          
                                                                        
c ........ Below is ZBL Universal Interatomic Potential ................
                                                                        
390      IF (R.gt.200) R=200                                            
         RZBL=R                                                         
         EX1=0.18175*EXP(-3.1998*RZBL)                                  
         EX2=0.50986*EXP(-0.94229*RZBL)                                 
         EX3=0.28022*EXP(-0.40290*RZBL)                                 
         EX4=0.02817*EXP(-0.20162*RZBL)                                 
         V=(EX1+EX2+EX3+EX4)/R                                          
         V1=-1.0*(V+3.1998*EX1+0.94229*EX2+0.4029*EX3+0.20162*EX4)/R    
         FR=B*B/R+V*R/EPS-R                                             
         FR1=-B*B/(R*R)+(V+V1*R)/EPS-1.0                                
         Q=FR/FR1                                                       
         R=R-Q                                                          
         IF (ABS(Q/R).gt.0.001) GOTO 390                                
                                                                        
         ROC=-2.0*(EPS-V)/V1                                            
         SQE=SQRT(EPS)                                                  
         CC=(0.011615+SQE)/(0.0071222+SQE)                              
         AA=2.0*EPS*(1+(0.99229/SQE))*B**CC                             
         FF=(SQRT(AA**2+1.0)-AA)*((9.3066+EPS)/(14.813+EPS))            
         DELTA=(R-B)*AA*FF/(FF+1)                                       
         C=(B+DELTA+ROC)/(R+ROC)                                        
                                                                        
      ENDIF                                   ! End of magic formula    
                                                                        
      C2=C**2                                 ! Scattering angle        
      C2=AMIN1(C2,1.0)                                                  
      S2=1-C2                                                           
                                                                        
c ........ Calculate ion energy and trajectory after collision .........
                                                                        
      CT=2*C2-1.0                         ! COS(theta)                  
      CT=AMIN1(CT,1.0)                                                  
      ST=SQRT(1.0-CT**2)                  ! SIN(theta)                  
      CU=CT+M2M1(1,J1)                                                  
      IF (CU.eq.0.0) CU=fuzz                                            
      PSI=ATAN(ST/CU)                                                   
      IF (PSI.lt.0) PSI=PSI+pi            ! Ion deflection angle in LAB 
      DEN=EC(1,J1)*S2*E                   ! E transferred to recoil atom
      E=E-DEN                             ! New ion energy in LAB coords
                                                                        
c ........ Store collision results .....................................
c ........ J1=Atom type hit, Nchan=Depth channel .......................
                                                                        
      IF (x.le.0.or.x.ge.xdepth) GOTO 443           ! Miss out cascades 
c     IF (xold.gt.0) DEE=DEE*x/(x-xold)             ! if outside target 
c     Nchan1=Nchan                                                      
c     Nchan1=MAX0(Nchan1,1)                                             
c     Nchan=INT((x)/CW)+1                                               
c     Nchan=MAX0(Nchan,1)                                               
                                                                        
c .......... Spread electronic energy loss over free flight path .......
                                                                        
c     Dchan=ABS(Nchan-Nchan1)+1                                         
c     IF (Dchan.le.1) THEN                                              
c       EIONIZ(1,Nchan)=EIONIZ(1,Nchan)+DEE*perang  ! Ionizations by ion
c       E2REC(1,Nchan)=E2REC(1,Nchan)+DEN*perang    ! Recoil E by ion   
c       E2REC(J1,Nchan)=E2REC(J1,Nchan)+DEN*perang  ! Recoil E by atoms 
c     ELSE                                                              
c       IF (COSX.ge.0) THEN                                             
c          Iondir=1                                                     
c       ELSE                                                            
c          Iondir=-1                                                    
c       ENDIF                                                           
c       DO 410,i=Nchan1,Nchan,iondir                                    
c         EIONIZ(1,i)=EIONIZ(1,i)+DEE*perang/Dchan                      
c         E2REC(1,i)= E2REC(1,i)+DEN*perang/Dchan                       
c         E2REC(J1,i)= E2REC(J1,i)+DEN*perang/Dchan                     
c410     CONTINUE                                                       
c      ENDIF                                                            
                                                                        
c ............ Decide if have cascade or phonons .......................
                                                                        
c     IF (DEN.ge.ED) THEN                          ! Store vacancy      
c        VACTOT=VACTOT+perang                                           
c        VAC(1,Nchan)=VAC(1,Nchan)+perang          ! Vacancies by ions  
c        VAC(J1,Nchan)=VAC(J1,Nchan)+perang        ! Vacancies by atoms 
c     ELSE                                                              
c        IF (x.gt.0.and.x.lt.xdepth) THEN    ! Dump energy into PHONONS 
c          PHONON(1,Nchan)=PHONON(1,Nchan)+DEN*perang                   
c        ENDIF                                                          
c     ENDIF                                                             
                                                                        
c ######################################################################
c ############### Calculate RECOIL CASCADES ############################
c ######################################################################
c .... From before: Label= Atom type, ER= Energy of recoil .............
c .... We follow Primary recoil until it stops. Every time it hits an ..
c .... atom we store that recoil as Atom(II) and the primary as (II+1) .
c .... When the recoil stops, the loop follows back from the last atom .
c .... hit until it get to Atom(II), the first atom hit by the recoil ..
c **********************************************************************
                                                                        
490   CONTINUE                           ! Enter here from main program 
      Label(1)=J1                        ! Atom number hit              
      ER(1)=DEN-EB                       ! Energy given to atom         
      ER(1)=AMAX1(ER(1),0.0)                                            
                                                                        
c ............. Find primary knock-on direction ........................
                                                                        
      PSI2=ATAN(ST/(1.0-CT+fuzz))                                       
      CHI=PHI-pi                              ! Recoil deflection angle 
      CPHI=COS(CHI)                                                     
      SPHI=SIN(CHI)                                                     
      CPSI=COS(PSI2)                                                    
      SPSI=SIN(PSI2)                                                    
      IF (SINX.eq.0) THEN                     ! Vertical only           
        COSXR(1)=COSX*CPSI                                              
        COSYR(1)=-SPSI*CPHI                                             
      ELSE                                                              
        DUM=-COSX*COSY/(SINX*SINY+fuzz)                                      
        DUM=AMIN1(DUM,1.0)                          ! Ensure -1<DUM<1   
        DUM=AMAX1(DUM,-1.0)                                             
        DELTA=CHI-ACOS(DUM)                                             
        COSXR(1)=CPSI*COSX+SPSI*SINX*CPHI                               
        COSYR(1)=CPSI*COSY+SPSI*SINY*COS(DELTA)                         
      ENDIF                                                             
      SINXR=SQRT(ABS(1.0-COSXR(1)**2))                                  
      SINYR=SQRT(ABS(1.0-COSYR(1)**2))                                  
                                                                        
c ............ Find starting positions, xr, yr of recoil ...............
                                                                        
      xr(1)=x                                                           
      yr(1)=y-P*(COSYR(1)-CPSI*COSY)/(SPSI+fuzz)                        
      II=1                                                              
                                                                        
c ############## Cycle loop of recoil particles ########################
                                                                        
530   J1=Label(II)                            ! J1=atom type            
                                                                        
c .. After next collision, Primary atom will be (II+1) and .............
c .. Recoil will be (II) ...............................................
                                                                        
      xr(II+1)=xr(II)+ATdist*COSXR(II)        ! x position of next atom 
      yr(II+1)=yr(II)+ATdist*COSYR(II)        ! y position of next atom 
                                                                        
c .......... If next position is outside target, .......................
c .......... calculate sputtering values ...............................
                                                                        
      IF (xr(II+1).le.0) THEN                                           
        Enor=ER(II)*COSXR(II)**2           ! Energy normal to surface   

        IF (Enor.lt.Esurf) THEN                                         
         COSXR(II)=-COSXR(II)              ! Bounces back               
         xr(II+1)=0.0001                   ! but stays close to surface 
         GOTO 530                          ! and continues onward...

        ELSE                               ! Atom sputters from surface 

c         IF (dz.lt.0.5) THEN               ! Random direction           
c          dz=-1.0                          ! for z velocity             
c         ELSEIF (dz.ge.0.5) THEN                                        
c          dz=1.0                                                        
c         ENDIF                                                          

c         dx=1.0                                                         
c         dy=1.0                                                         
c         IF (COSXR(II).lt.0) dx=-1.0       ! x velocity direction       
c         IF (COSYR(II).lt.0) dy=-1.0       ! y velocity direction       

         nsput(J1)=nsput(J1)+1              ! Add 1 to no. sputtered
         nsput1=nsput1+1     
         n=nsput(J1)                                                        
         IF (n.gt.1000000) THEN               ! Max array size             
          spError$=' Array too large. Only 1000000 points saved.'         
          n=1000000                                                       
         ENDIF                                                          

c        Atomfactor=SQRT(2.0*qq/(mp*Zmass(J1)))   ! Converts E to v      
         Esput=Esput+ER(II)-Esurf                 ! Total sputtered energy     

c        IA=(180/pi)*ABS(ATAN(SINXR/COSXR(II)))   ! Sputter angle (uncorr.)     

         IA=SQRT( (Enor-Esurf)/(ER(II)-Esurf) )
         IA=(180/pi)*ACOS(IA)               ! Convert to degs
c                                           ! Sputter angle corrected
c                                           ! for diffraction (see
c                                           ! Eckstein, pg 81)

c        SPenergy(J1,n)=ER(II)-Esurf              ! Sputter Energy      
c        SPangle(J1,n)=IA                         ! Sputter angle       
c        SPlateral(n)=yr(II)                      ! Lateral distance    
c        SPxvel(n)=dx*Atomfactor*SQRT(ABS(Esput*COSXR(II)))  ! x vel    
c        SPyvel(n)=dy*Atomfactor*SQRT(ABS(Esput*COSYR(II)))  ! y vel    
c        COSZ=SQRT(ABS(1.0-COSXR(II)**2-COSYR(II)**2))                  
c        SPzvel(n)=dz*Atomfactor*SQRT(Esput*COSZ)            ! z vel    
        ENDIF                                                           

        GOTO 650                           ! To next atom               
      ENDIF                                                             
                                                                        
c ............ If atom transmitted, stop following its trajectory ......
                                                                        
      IF (xr(II+1).ge.xdepth) GOTO 650                                  
                                                                        
c ........... Save energy of recoil ....................................
                                                                        
      Erec(J1)=ER(II)                               ! Energy in recoil  
                                                                        
c ... Don't follow recoils if they can't do damage .....................
                                                                        
      IF (ER(II).lt.EDlow) GOTO 640                 ! Recoil stopped,   
                                                    ! store results     
c ........... xp=position of primary atom (II) after collision .........
                                                                        
      SINXR=SQRT(ABS(1-COSXR(II)**2))                                   
      SINYR=SQRT(ABS(1-COSYR(II)**2))                                   
      RN=ran3(seed)                                                    
      PHI2=2.0*pi*RN                                                    
      xp=xr(II+1)-P*COS(PHI2)*SINXR                                     
      IF (xp.ge.xdepth) xp=xdepth-fuzz             ! atom stays inside  
                                                                        
c ........... Determine target species, JJ1, and EPS and B values.......
                                                                        
      Flag1=0                                                           
      sumstoich=0                                                       
      RN=ran3(seed)                                                    
      DO 590,JJ1=2,natom                                                
       sumstoich=sumstoich+stoich(JJ1)                                  
       IF (RN.lt.sumstoich.and.Flag1.eq.0) THEN                         
         Flag1=JJ1                                                      
         EPS=Erec(J1)*F(J1,JJ1)                                         
         P=SQRT(RN)*Pmax                                                
         B=P/A(J1,JJ1)                            ! Impact parameter    
       ENDIF                                                            
590   CONTINUE                                                          
      JJ1=Flag1                                                         
c ................... If Reduced E, EPS>10 use Rutherford ..............
c ................... scattering, otherwise use Magic formula ..........
                                                                        
      IF (EPS.GT.10) THEN                        ! Rutherford scattering
       S2=1.0/(1.0+(1.0+B*(1.0+B))*(2.0*EPS*B)**2)                      
       C2=1.0-S2                                                        
      ELSE                                       ! Magic for for recoils
       R=B                                                              
       RR=-2.7*LOG(ABS(EPS*B)+fuzz)                                     
       IF (RR.gt.B) THEN                                                
         RR=-2.7*LOG(ABS(EPS*RR)+fuzz)                                  
         IF(RR.GT.B) R=RR                                               
       ENDIF                                                            
                                                                        
c ........ Below is the ZBL Universal interatomic potential ............
                                                                        
600    IF (R.gt.200) R=200                                              
       R1=R                                                             
       EX1=0.18175*EXP(-3.1998*R1)                                      
       EX2=0.50986*EXP(-0.94229*R1)                                     
       EX3=0.28022*EXP(-0.40290*R1)                                     
       EX4=0.02817*EXP(-0.20162*R1)                                     
       V=(EX1+EX2+EX3+EX4)/R                                            
       V1=-(V+3.1998*EX1+0.94229*EX2+0.4029*EX3+0.20162*EX4)/R          
       FR=B*B/R+V*R/EPS-R                                               
       FR1=-B*B/(R*R)+(V+V1*R)/EPS-1.0                                  
       Q=FR/FR1                                                         
       R=R-Q                                                            
       IF (ABS(Q/R).gt.0.0001) GOTO 600                                 
       ROC=-2.0*(EPS-V)/V1                                              
       SQE=SQRT(EPS)                                                    
                                                                        
c .. The next lines are the 5-Parameter MAGIC scattering calculations ..
                                                                        
       CC=(0.011615+SQE)/(0.0071222+SQE)                                
       AA=2.0*EPS*(1.0+(0.99229/SQE))*B**CC                             
       FF=(SQRT(AA**2+1.0)-AA)*((9.3066+EPS)/(14.813+EPS))              
       DELTA=(R-B)*AA*FF/(FF+1.0)                                       
       C=(B+DELTA+ROC)/(R+ROC)                                          
       C2=C*C                                                           
       C2=AMIN1(C2,1.0)                                                 
       S2=1.0-C2                                                        
      ENDIF                                                             
                                                                        
c .................  End of recoil MAGIC .............................. 
                                                                        
      CT=2.0*C2-1.0                          ! COS(theta), scattering   
      CT=AMIN1(CT,1.0)                       ! angle in LAB frame       
      ST=SQRT(1-CT*CT)                       ! SIN(theta)               
      CU2=CT+M2M1(J1,JJ1)                                               
      IF (CU2.eq.0) CU2=fuzz                                            
      PSI2=ATAN(ST/CU2)                                                 
      IF (PSI2.lt.0) PSI2=PSI2+pi                                       
      py=-P*(COSYR(II)-CPSI*COSY)/(SPSI+fuzz)    ! y position           
                                                                        
c .............. Store atom(II) data (after collision) in (II+1) .......
                                                                        
      DUM=-COSXR(II)*COSYR(II)/(SINXR*SINYR+fuzz)                       
      DUM=AMIN1(DUM,1.0)                        ! Ensure -1<DUM<1       
      DUM=AMAX1(DUM,-1.0)                                               
      DELTA=PHI2-ACOS(DUM)                                              
      COSXR(II+1)=COSXR(II)*COS(PSI2)+SINXR*SIN(PSI2)*COS(PHI2)         
      COSYR(II+1)=COSYR(II)*COS(PSI2)+SINYR*SIN(PSI2)*COS(DELTA)        
      Label(II+1)=J1                                                    
                                                                        
c .......Calculate electronic energy loss of atom(J1), DEER ............
                                                                        
      IE=INT((Erec(J1)-fuzz)/E0keV)                                     
      IF (IE.gt.999) THEN                                               
       IE=999                                                           
      ELSEIF (IE.ne.0) THEN                                             
       DEER=ATdist*(SE(IE,J1)+((Erec(J1)/E0keV)-IE)*(SE(IE+1,J1)        
     &      -SE(IE,J1)))                                                
      ENDIF                                                             
      ER(II)=AMAX1(ER(II),0.0)                                          
      IF (IE.eq.0) DEER=ATdist*SE(1,J1)*SQRT(ER(II)/E0keV)              
                                                                        
c ........ Calculate nuclear energy loss of recoil, DENR ...............
                                                                        
      DENR=Erec(J1)*S2*EC(J1,JJ1)                                       
      ER(II+1)=Erec(J1)-DENR-DEER         ! Energy primary recoil keeps 
                                                                        
c ........ Store electronic energy-loss ................................
                                                                        
c     IF (xp.gt.0.and.xp.lt.xdepth) THEN                                
c       LLL=INT(xr(II+1)/CW)+1                                          
c       LLL=MIN0(LLL,200)                                               
c       EIONIZ(2,LLL)=EIONIZ(2,LLL)+DEER*perang                         
c     ENDIF                                                             
                                                                        
c ............... Store data for NEW recoil ............................
c ..... Secondary Knock-on=atom(II) since Primary is now atom(II+1) ....
                                                                        
      Label(II)=JJ1                       ! Secondary recoil atom type   
      ER(II)=DENR-EB                      ! Energy secondary gets       
                                                                        
c ......... Origin of new recoil, xp, was previously determined ........
                                                                        
      xr(II)=xp                           ! x position of secondary     
      yr(II)=yr(II+1)+py                  ! y position of secondary     
                                                                        
c ........ Find new direction of secondary atom (II) recoil ............
                                                                        
      PSI2=ATAN(ST/(1-CT+fuzz))                                         
      CHI2=PHI2-pi                                                      
      CPHI=COS(CHI2)                                                    
      SPHI=SIN(CHI2)                                                    
      CPSI=COS(PSI2)                                                    
      SPSI=SIN(PSI2)                                                    
      DUM=-COSXR(II)*COSYR(II)/(SINXR*SINYR+fuzz)                       
      DUM=AMIN1(DUM,1.0)                           ! Ensures -1<DUM<1   
      DUM=AMAX1(DUM,-1.0)                                               
      DELTA=CHI-ACOS(DUM)                                               
      COSXR(II)=COSXR(II)*CPSI+SINXR*SPSI*CPHI                          
      COSYR(II)=COSYR(II)*CPSI+SINYR*SPSI*COS(DELTA)                    
                                                                        
c .........  Continue following primary recoil which is now (II+1) .....
                                                                        
620   II=II+1                                                           
                                                                        
c ........ Check if recoil exceeds storage array .......................
                                                                        
      IF  (II.gt.1999) THEN                                             
       II=1999                                                          
       spError$='Recoil Cascade too BIG. Sputter data wrong.'           
      ENDIF                                                             
                                                                        
c ............. Cycle recoil if ER>EDlow ...............................
                                                                        
      IF (ER(II).gt.EDlow) GOTO 530                  ! To next collision
                                                                        
c .............. Current recoil stopped. Store results .................
                                                                        
640   CONTINUE                                                          
c     IF (xr(II).le.0.or.xr(II).ge.xdepth) GOTO 650  ! Atom left surface
c     J1=Label(II)                                                      
c     LLL=INT(xr(II)/CW)+1                                              
c     LLL=MIN0(LLL,200)                                                 
c     PHONON(2,LLL)=PHONON(2,LLL)+ER(II)*perang                         
                                                                        
c ............... Last Atom in Cascade? ................................
                                                                        
c     IF (II.le.1) THEN                    ! Atom stops, store position 
c       RANGE(J1,LLL)=RANGE(J1,LLL)+perang                              
c       INTERTOT=INTERTOT+perang                                        
c     ELSEIF (ER(II-1).gt.ED.and.J1.eq.Label(II-1)) THEN                
c       XNOVAC(LLL)=XNOVAC(LLL)+perang      ! Atom ends as replacement  
c       VAC(JJ1,LLL)=VAC(JJ1,LLL)-perang    ! collision.                
c       VACTOT=VACTOT-perang                                            
c     ELSE                                                              
c       RANGE(J1,LLL)=RANGE(J1,LLL)+perang                              
c       INTERTOT=INTERTOT+perang                                        
c     ENDIF                                                             
                                                                        
c .............. Choose next recoil from array .........................
                                                                        
650   II=II-1                                 ! End of recoils when I=0 
      IF (II.gt.0) THEN                                                 
                                                                        
c ........... See if new recoil has enough energy to move ..............
                                                                        
        IF (ER(II).lt.EDlow) GOTO 650       ! Not enough E, so next atom
        J1=Label(II)                                                    
        IF (xr(II).gt.0.and.xr(II).lt.xdepth) THEN                      
c         LLL=INT(xr(II)/CW)+1                                          
c         LLL=MIN0(LLL,200)                                             
c         VAC(J1,LLL)=VAC(J1,LLL)+perang    ! Recoil moved, make vacancy
c         VACTOT=VACTOT+perang                                          
          GOTO 530                          ! To cascade calculation    
        ENDIF                                                           
        GOTO 650                            ! Next recoil               
      ENDIF                                 ! If I=0, ends recoils here 
                                                                        
c ######################################################################
c ################# End of Recoil section ##############################
c ######################################################################
                                                                        
c ................ Get new ion direction ...............................
c ...... PHI=Random azimuthal angle; PSI= Nuclear scattering angle .....
                                                                        
443   DUM=-COSX*COSY/(SINX*SINY+fuzz)                                   
      DUM=AMIN1(DUM,1.0)                          ! Ensures -1<DUM<1    
      DUM=AMAX1(DUM,-1.0)                                               
      DELTA=PHI-ACOS(DUM)                                               
      COSX=COSX*COS(PSI)+SINX*SIN(PSI)*COS(PHI)                         
      COSY=COSY*COS(PSI)+SINY*SIN(PSI)*COS(DELTA)                       
      SINX=SQRT(ABS(1-COSX*COSX))                                       
      SINY=SQRT(ABS(1-COSY*COSY))                                       
                                                                        
c ................. See if energy > Final energy, EF .................. 
                                                                        
      IF (E.gt.EF) GOTO 280               ! Go to next collison of ion  
                                                                        
c ######################################################################
c ###################### End Cycling Ions ##############################
c ######################################################################
                                                                        
c ............Ion has stopped, save data arrays ....................... 
                                                                        
      IF (x.gt.0.and.x.lt.xdepth) THEN           ! Ion inside surface   
c        Nchan=MIN0(100,INT(x/CW)+1)                                    
c        RANGE(1,Nchan)=RANGE(1,Nchan)+perang    ! Save range in array  
c        YRANGE(1,Nchan)=YRANGE(1,Nchan)+1       ! Straggling, by ion   
c        YRANGE(2,Nchan)=YRANGE(2,Nchan)+Y*Y     ! & by atoms           
         nimplant=nimplant+1                     ! Total no. implanted  
c        n=MIN0(nimplant,100000)                                         
c        xpos(n)=x                               ! x when ion stops     
c        ypos(n)=y                               ! y when ion stops     
      ENDIF                                                             
                                                                        
c ... Store sums for values of current range and straggling ............
                                                                        
      xsum=xsum+x                                ! Range                
      x2sum=x2sum+x**2                           ! Straggling           
      x3sum=x3sum+x**3                           ! Skewness             
      x4sum=x4sum+x**4                           ! Kurtosis             
      y2sum=y2sum+y**2                                                  
      xy2sum=xy2sum+x*y**2                                              
      x2y2sum=x2y2sum+X**2*y**2                                         
      y4sum=y4sum+y**4                                                  
      plsum=plsum+pl                             ! Total path length    
      pl2sum=pl2sum+pl**2                                               
      ncollsum=ncollsum+ncoll                    ! Total no. collisions 
                                                                        
c ................ See if program terminator flags exceeded ............
                                                                        
660   counter=counter+1                                                 
c      IF (counter.gt.INT(ionstop/10)) THEN    ! Message on screen to let
c        PRINT*,'Working...'                  ! a real-time user know   
c         counter=0                            ! prog hasn't crashed.    
c      ENDIF                                                             
                                                                        
      IF ((ionnum.lt.ionstop).and.(ntransmit.lt.tmstop).and.            
     &    (nback.lt.backstop).and.(nsput1.lt.sputstop)) GOTO 250         
                                                ! Loop back to next ion 
                                                                        
c ######################################################################
c ################# End of TRIM calculation ############################
c ######################################################################
                                                                        
c ............. Calculate statistical parameters .......................
                                                                        
      Flag1=0                                                           
      nimplant=ionnum-ntransmit-nback           ! No. ions implanted    
      IF (nimplant.le.0) THEN                                           
        Flag1=1                                                         
        GOTO 797                                ! Miss out stats calcs  
      ENDIF                                                             
      N1=FLOAT(nimplant)                        ! Made a REAL number    
      avex=xsum/N1                              ! Average depth of ion  
      avepl=plsum/N1                            ! Average path length   
      avecol=ncollsum/FLOAT(ionnum)             ! Average no. collisions
      varx=x2sum/N1-avex**2                     ! Variance of depth dist
      vary=y2sum/N1                             ! Variance of y dist    
      sigmay=SQRT(vary)                         ! Standard dev of y dist
      xy2=xy2sum/(N1*sigmay**3)                                         
      x2y2=x2y2sum/(N1*vary**2)                                         
      y4=y4sum/(N1*vary**2)                                             
      sigma=SQRT(varx)                          ! Standard dev of x dist
      vv=avex/sigma                                                     
      vv2=vv**2                                                         
      Gamma=x3sum/(N1*sigma**3)-vv*(3.0+vv2)                            
      Beta=x4sum/N1/varx**2-4.0*vv*Gamma-vv2*(6.0+vv2)                  
      sigpl=SQRT(pl2sum/N1-avepl**2)            ! Standard dev of pl dis
                                                                        
c .......... Save all data files .......................................
                                                                        
797   IF (nimplant.ne.0) THEN                                           
c        n=MIN0(nimplant,100000)                                          
c        fileid$=filenm$//'.xps'                                   
c        CALL SAVEDATA(fileid$,xpos,n)                                   
c        fileid$=filenm$//'.yps'                                   
c        CALL SAVEDATA(fileid$,ypos,n)                                   
      ENDIF                                                             
                                                                        
c      IF (nback.ne.0) THEN                                              
c        n=MIN0(nback,100000)                                             
c        fileid$=filenm$//'.bse'                               
c        CALL SAVEDATA(fileid$,BSenergy,n)                               
c        fileid$=filenm$//'.bsa'                                
c        CALL SAVEDATA(fileid$,BSangle,n)                                
c        fileid$=filenm$//'.bsl'                                  
c        CALL SAVEDATA(fileid$,BSlateral,n)                              
c        fileid$=filenm$//'.bsx'                                 
c        CALL SAVEDATA(fileid$,BSxvel,n)                                 
c        fileid$=filenm$//'.bsy'                                 
c        CALL SAVEDATA(fileid$,BSyvel,n)                                 
c        fileid$=filenm$//'.bsz'                                 
c        CALL SAVEDATA(fileid$,BSzvel,n)                                 
c      ENDIF                                                             
                                                                        
c     IF (ntransmit.ne.0) THEN                                          
c       n=MIN0(ntransmit,100000)                                         
c       fileid$=filenm$//'.tme'                               
c       CALL SAVEDATA(fileid$,TMenergy,n)                               
c       fileid$=filenm$//'.tma'                                
c       CALL SAVEDATA(fileid$,TMangle,n)                                
c       fileid$=filenm$//'.tms'                                 
c       CALL SAVEDATA(fileid$,TMlateral,n)                              
c       fileid$=filenm$//'tmx'                                 
c       CALL SAVEDATA(fileid$,TMxvel,n)                                 
c       fileid$=filenm$//'.tmy'                                 
c       CALL SAVEDATA(fileid$,TMyvel,n)                                 
c       fileid$=filenm$//'.tmz'                                 
c       CALL SAVEDATA(fileid$,TMzvel,n)                                 
c     ENDIF                                                             
                                                                        
c      IF (nsput1.ne.0) THEN                                              
c        n=MIN0(nsput1,100000)                                             
c        fileid$=filenm$//'.spe'                               
c        CALL SAVEDATA(fileid$,SPenergy,n)                               
c        OPEN(UNIT=13,FILE=fileid$,FORM='UNFORMATTED',STATUS='NEW')
c         WRITE(13)nsput1,natom
c         WRITE(13)((SPenergy(J1,I),J1=2,natom),I=1,nsput1)
c        CLOSE(13)

c        fileid$=filenm$//'.spa'                                
c       CALL SAVEDATA(fileid$,SPangle,n)                                
c        OPEN(UNIT=13,FILE=fileid$,FORM='UNFORMATTED',STATUS='NEW')              
c         WRITE(13)nsput1                                                        
c         WRITE(13)((SPangle(J1,i),J1=2,natom),I=1,nsput1)
c        CLOSE(13)

c       fileid$=filenm$//'spl'                                 
c       CALL SAVEDATA(fileid$,SPlateral,n)                              
c        fileid$=filenm$//'.spx'                                 
c        CALL SAVEDATA(fileid$,SPxvel,n)                                 
c        fileid$=filenm$//'.spy'                                 
c        CALL SAVEDATA(fileid$,SPyvel,n)                                 
c        fileid$=filenm$//'.spz'                                 
c        CALL SAVEDATA(fileid$,SPzvel,n)                                 
c      ENDIF                                                             
                                                                        
c ...... Note: Other data arrays are not saved in this version of ......
c ...... TRIM to save time. These include all the recoil data, .........
c ...... phonons, straggling, vacancies and ionizations ................
c ...... Arrays:- RANGE(i,j), YRANGE(i,j), PHONONS(i,j), VAC(i,j) ......
c ...... ER(i), E2rec(i,j), Eioniz(i,j) etc ............................
                                                                        
c ............. Save info file .........................................
                                                                        
      fileid$=filenm$//'.inf'                                     
      INQUIRE(file=fileid$,EXIST=ext)                                   
      IF (ext) THEN                                                     
         PRINT*,fileid$,' already exists. Cannot save new file'         
         GOTO 850                                                       
      ENDIF                                                             
                                                                        
      OPEN(UNIT=13,FILE=fileid$,FORM='FORMATTED',STATUS='NEW')          
                                                                        
       WRITE(13,833)                                                    
c       WRITE(13,801)fileid$
       WRITE(13,831)                                             
       WRITE(13,833)                                                    
       WRITE(13,802)Z$(1),Zmass(1)                                      
       WRITE(13,803)E0keV*1000                                          
       WRITE(13,804)angle                                               
       WRITE(13,805)LFCTR                                               
       WRITE(13,800)                                                    
       WRITE(13,806)layername$,rho                                      
       WRITE(13,807)Esurf,EB,ED                                               
       WRITE(13,808)Bragg                                               
       WRITE(13,800)                                                    
       WRITE(13,809)                                                    
       DO 799,i=2,natom                                                 
         WRITE(13,810)Z$(i),NINT(Z(i)),Zmass(i),stoich(i),nsput(i),
     &        nsput(i)/FLOAT(ionnum)               
799    CONTINUE
       WRITE(13,832)                                                         
       WRITE(13,800)                                                    
       WRITE(13,811)Z2,Zmass2                                           
       WRITE(13,812)ionnum                                              
       WRITE(13,813)nsput1,nsput1/FLOAT(ionnum)                           
       WRITE(13,814)ntransmit,ntransmit/FLOAT(ionnum)                   
       WRITE(13,815)nback,nback/FLOAT(ionnum)                           
       WRITE(13,816)nimplant,nimplant/FLOAT(ionnum)                     
       WRITE(13,800)                                                    
       IF (nsput1.ne.0) THEN                                             
         WRITE(13,817)Esput/FLOAT(nsput1)                                
       ENDIF                                                            
       IF (ntransmit.ne.0) THEN                                         
         WRITE(13,818)Etrans/FLOAT(ntransmit)                           
       ENDIF                                                            
       IF (nback.ne.0) THEN                                             
         WRITE(13,819)Eback/FLOAT(nback)                                
       ENDIF                                                            
       WRITE(13,800)                                                    
       IF (Flag1.eq.1) GOTO 798                 ! Miss out stats if all 
       WRITE(13,820)avex                        ! ions transmitted      
       WRITE(13,821)avepl                                               
       WRITE(13,822)avecol                                              
       WRITE(13,823)sigma                                               
       WRITE(13,824)sigpl                                               
       WRITE(13,825)Gamma                                               
       WRITE(13,826)Beta                                                
798    WRITE(13,800)                                                    
       WRITE(13,827)                                                    
       WRITE(13,828)bkError$                                            
       WRITE(13,829)tmError$                                            
       WRITE(13,830)spError$                                            
       WRITE(13,800)                                                    
                                                                        
      CLOSE(13)                                                         
                                                                        
800    FORMAT(1H ,'<HR>')                                               
c 801    FORMAT(1H ,T20,A21)                                              
802    FORMAT(1H ,'Ion beam = ',A2,' of mass ',F6.2,' amu <BR>')             
803    FORMAT(1H ,'Initial energy = ',F11.3,' eV <BR>')                      
804    FORMAT(1H ,'Initial angle = ',F4.1,' degs <BR>')                      
805    FORMAT(1H ,'Screening length of ion = ',F7.2,'<BR>')                    
806    FORMAT(1H ,'Target is <B>',A10,'</B> of density ',F7.3,
     &' g/cc<BR>')         
807    FORMAT(1H ,'Energies, Esurf = ',F5.2,' Eb = ',F5.2,' Ed = '
     &,F5.2,' eV<BR>')          
808    FORMAT(1H ,'Bragg correction factor = ',F4.2,'<BR>')
809    FORMAT(1H ,'<TABLE BORDER><TR ALIGN=CENTER><TH>Atom<TH>At.No.
     &<TH>Mass<TH>%<TH>No. Sputtered<TH>Sputter Yield')     
810    FORMAT(1H ,'<TR ALIGN=CENTER><TD>',A2,'<TD>',I7,'<TD>',F12.3,
     &'<TD>',F10.3,'<TD>',I7,'<TD>',F10.4)                         
811    FORMAT(1H ,'Mean values for layer:  Z = ',F5.2,' and M = ',F6.2,
     & '<BR>') 
812    FORMAT(1H ,'Total no. ions used =   ',I10,'<BR>')                       
813    FORMAT(1H ,'Total no. atoms sputtered =   ',I7,' : Sputter yield      
     & = ',F6.4,'<BR>')                                                        
814    FORMAT(1H ,'No. ions transmitted =  ',I10,' : Transmit yield     
     & = ',F6.4,'<BR>')                                                        
815    FORMAT(1H ,'No. ions backscattered =',I10,' : Backscatter yield  
     & = ',F6.4,'<BR>')                                                        
816    FORMAT(1H ,'No. ions implanted =    ',I10,' : Implant yield      
     & = ',F6.4,'<BR>')                                                        
817    FORMAT(1H ,'Average energy of sputtered atoms      =',F12.3,' eV
     &<BR>')                                                                 
818    FORMAT(1H ,'Average energy of transmitted ions     =',F12.3,' eV
     &<BR>')                                                                 
819    FORMAT(1H ,'Average energy of backscattered ions   =',F12.3,' eV
     &<BR>')                                                                 
820    FORMAT(1H ,'Average depth of ions                  =',F10.2,' A
     &<BR>')
821    FORMAT(1H ,'Average path length of ions            =',F10.2,' A
     &<BR>')
822    FORMAT(1H ,'Average no. collisions per ion         =',F10.2,
     &'<BR>')
823    FORMAT(1H ,'Standard deviation of ion distribution =',F10.2,' A
     &<BR>')
824    FORMAT(1H ,'Standard deviation of path length dist =',F10.2,' A
     &<BR>')
825    FORMAT(1H ,'Skewness (Gamma) of ion distribution   =',F10.2,
     &'<BR>')     
826    FORMAT(1H ,'Kurtosis (Beta) of ion distribution    =',F10.2,
     &'<BR>')     
827    FORMAT(1H ,'Errors:<BR>')                                            
828    FORMAT(1H ,'Backscattered ions : ',A50,'<BR>')                          
829    FORMAT(1H ,'Transmitted ions   : ',A50,'<BR>')                          
830    FORMAT(1H ,'Sputtered atoms    : ',A50,'<BR>')                          
831    FORMAT(1H ,'<H1 Align=CENTER>TRIM Output Data</H1>')
832    FORMAT(1H ,'</TABLE><P>')
833    FORMAT(1H ,'<CENTER>
     &<IMG SRC="http://www.tlchm.bris.ac.uk/clrline.gif"></CENTER><P>')
                                                                        
c ................ End program .........................................
                                                                        
850   STOP                                                              
      END                                                               
                                                                        
c **********************************************************************
c *********** End of main Program **************************************
c **********************************************************************
                                                                        
                                                                        
c **********************************************************************
                                                                        
      SUBROUTINE GETATNO(symbol$,atno,Zname$)                           
                                                                        
       INTEGER atno                                                     
       CHARACTER Zname$(92)*2,symbol$*2                                 
                                                                        
       DO 100,i=1,92                                                    
         IF (Zname$(i).eq.symbol$) atno=i                               
100    CONTINUE                                                         
                                                                        
      RETURN                                                            
      END                                                               
                                                                        
c **********************************************************************
                                                                        
       SUBROUTINE RSTOP(Z1,M1,Z2,EE,SE,LFCTR,VFERMI)                    
                                                                        
c ......................................................................
                                                                        
c     "RSTOP" : ION STOPPING CROSS-SECTION SUBROUTINE USED BY T 
c                  COMPILED BY E. HOWARD  1990                          
                                                                        
c     This program calculates 1000 values of Electronic Stopping        
c     in (eV-Ang**2). Nuclear Stopping is NOT calculated.               
                                                                        
c ......................................................................
                                                                        
c            Z1 = Ion atomic number                                     
c            MM1 = Ion atomic mass                                      
c            M1 = Ion atomic weight (amu)                               
c            Z2 = Target atomic number                                  
c            M2 = Target atomic weight (amu)                            
c            EE = Ion energy  (keV)                                     
c            RHO = Target density (g/cm3)                               
c            ATRHO = Target density (atoms/cm3)                         
c            VFERMI = (Fermi velocity of solid) / V0                    
c            LFCTR = Lambda screening factor for ions                   
c            PCOEF = Stopping coefficients for protons                  
c            SE = Calculated electronic stopping (eV-A2)                
                                                                        
c ......................................................................
                                                                        
c       Data file 'SCOEF DATA' contains the following coefficients:     
                                                                        
c            SCOEF(*,1) = atomic number of element                      
c            SCOEF(*,2) = atomic mass of most abundant isotope          
c            SCOEF(*,3) = atomic weight of m.a.i.                       
c            SCOEF(*,4) = target mass (amu)                             
c            SCOEF(*,5) = g/cm3 density (normal abundance)              
c            SCOEF(*,6) = atoms/cm3*1e22 (normal abundance)             
c            SCOEF(*,7) = (fermi velocity of solid) / V0                
c            SCOEF(*,8) = lambda screening factor for ions              
                                                                        
c       Data file 'SCOEF PSTOP' contains the atomic no. of the element  
c       followed by the 8 proton stopping coefficients.                 
                                                                        
c ......................................................................
                                                                        
      IMPLICIT NONE                                                     
                                                                        
      REAL Y1,MM1,Y2,Y3,Y4,Y5,LFCTR,ATRHO,VFERMI,M1,E0,E,A,B,HE0,       
     &     POWER,EEE,VMIN,L0,L,Q1,L1,Q,YR,YRMIN,V,SP,VR,HE,HEH,         
     &     PCOEF(8),X(8),SE(1000),VRMIN,ZETA,EE,Esubl                   
                                                                        
      INTEGER Z1,Z2,I                                                   
                                                                        
c ........... Obtain coefficients from subroutines .....................
                                                                        
c ..First call is for data on the ion,the second is about the target ...
c ..LFACTOR is the screening length of the ion, VFERMI is of the target.
c .... MM1 is the mass (amu) of the most abundant isotope of the ion ...
c .... Y1-5 and X() are dummy variables.................................
                                                                        
      CALL SCOEF(Z1,Y1,MM1,Y2,Y3,Y4,Y5,LFCTR,Y1,X)                      
      CALL SCOEF(Z2,Y1,Y2,Y3,Y4,ATRHO,VFERMI,Y5,Esubl,PCOEF)            
                                                                        
c ................. Calculate energy increment (ev/amu) ................
                                                                        
      IF (M1.EQ.0.) M1=MM1                                              
      E0=0.001*EE/M1                          ! energy increment        
                                                                        
c ..Loop to get 1000 values of stopping up to the initial ion Energy ...
                                                                        
      DO 110,I=1,1000                                                   
       E=E0*I                                                           
                                                                        
c .... Branch to the part of the program for ions:Z1=1,Z1=2,OR Z1>2 ....
                                                                        
       IF (Z1.EQ.1) THEN                                                
                                                                        
c ........ Get electronic stopping powers for Hydrogen .................
c ........ in eV/(1e15 atoms/cm2) ......................................
                                                                        
100     CALL RPSTOP(Z2,E,PCOEF,SE(I))                                   
                                                                        
       ELSEIF (Z1.EQ.2) THEN                                            
                                                                        
c..................Helium electronic stopping powers....................
c......... Velocity-proportional stopping below HE0 keV/amu ............
                                                                        
140     HE0=1.0             ! velocity proportional stopping threshold. 
        HE=AMAX1(HE0,E)                                                 
        B=ALOG(HE)                                                      
                                                                        
c ... Below is formula for the ratio of He to proton stopping powers ...
                                                                        
        A=0.2865+0.1266*B-0.001429*B*B+0.02402*B*B*B-0.01135*B**4       
     &   +0.001475*B**5                                                 
        HEH=1.0-EXP(-AMIN1(30.0,A))         ! He:H stopping power ratio 
                                                                        
c ........Add Z1**3 effect to He/H stopping power ratio, HEH ...........
                                                                        
        A=(7.6-AMAX1(0.0,ALOG(HE)))**2                                  
        A=(1.0+(0.007+0.00005*Z2)*EXP(-A))                              
                                                                        
        CALL RPSTOP(Z2,HE,PCOEF,SP)                                     
                                                                        
        SE(I)=SP*HEH*4.0                                                
                                                                        
c.....Calc. He velocity proportional stopping ..........................
                                                                        
        IF (E.LE.HE0) THEN                                              
         SE(I)=SE(I)*SQRT(E/HE0)                                        
        ENDIF                                                           
                                                                        
       ELSE                                                             
                                                                        
c ................ Heavy ion electronic stopping powers ................
c ... Use velocity stopping for (YRMIN=VR/Z1**.67) <= 0.13 .............
c ... or for VR <= 1.0 .................................................
                                                                        
150     YRMIN=0.13                                 ! Effective ion vel  
        VRMIN=1.0                                  ! Relative ion vel   
        V=SQRT(E/25)/VFERMI                                             
                                                                        
        IF (V.GE.1.0) THEN                                              
         VR=V*VFERMI*(1+1/(5*V*V))                                      
        ELSE                                                            
         VR=(3*VFERMI/4)*(1+(2*V*V/3)-V**4/15)                          
        ENDIF                                                           
                                                                        
c ... Set YR = maximum of (VR/Z1**.67),(VRMIN/Z1**.67) or YRMIN ........
                                                                        
180     YR=AMAX1(YRMIN,VR/Z1**0.6667)                                   
        YR=AMAX1(YR,VRMIN/Z1**0.6667)                                   
                                                                        
c ......... Calculate relation of ion's fractional charge state,Q ......
c ......... to its velocity, YR ........................................
                                                                        
        A=-0.803*YR**0.3+1.3167*YR**0.6+0.38157*YR+0.008983*YR*YR       
        Q=AMIN1(A,50.0)                                                 
        Q=AMIN1(1.0,AMAX1(0.0,1.0-EXP(-Q)))                             
                                                                        
c ... Calculate screening length, L and also convert ionization ........
c ... level to the effective charge ....................................
                                                                        
        B=(AMIN1(0.43,AMAX1(0.32,0.12+0.025*Z1)))/Z1**0.3333            
        L0=(0.8-Q*(AMIN1(1.2,0.6+Z1/30.0)))/Z1**0.3333                  
                                                                        
        IF (Q.LT.0.2) THEN                                              
          L1=0.0                                                        
        ELSEIF (Q.LT.(AMAX1(0.0,0.9-0.025*Z1))) THEN                    
          Q1=0.2                                                        
          L1=B*(Q-0.2)/ABS(AMAX1(0.0,0.9-0.025*Z1)-0.2000001)           
        ELSEIF (Q.LT.(AMAX1(0.0,1.0-0.025*AMIN1(16.0,1.0*Z1)))) THEN    
          L1=B                                                          
        ELSE                                                            
          L1=B*(1.0-Q)/(0.025*AMIN1(16.0,1.0*Z1))                       
        ENDIF                                                           
                                                                        
        L=AMAX1(L1,L0*LFCTR)                   ! Screening length       
                                                                        
c ...... Calculate the effective charge of the heavy ion, ZETA .........
c ...... This is called gamma in chapter 3 .............................
                                                                        
        ZETA=(4*L*VFERMI/1.919)**2                                      
        ZETA=Q+(1.0/(2.0*VFERMI**2))*(1.0-Q)*ALOG(1+ZETA)               
                                                                        
c ... Add Z1**3 effect as shown in ref. 779 ............................
                                                                        
        A=(7.6-AMAX1(0.0,ALOG(E)))**2                                   
        ZETA=ZETA*(1.0+(1.0/Z1**2)*(0.18+0.0015*Z2)*EXP(-A))            
                                                                        
c ... Calculate velocity stopping for YR greater than YRMIN ............
                                                                        
        IF (YR.GT.AMAX1(YRMIN,VRMIN/Z1**0.6667)) THEN                   
                                                                        
240      CALL RPSTOP (Z2,E,PCOEF,SP)                                    
         SE(I)=SP*(ZETA*Z1)**2                                          
                                                                        
        ELSE                                                            
                                                                        
c ... Calculate velocity stopping for YR less than YRMIN ...............
                                                                        
250      VRMIN=AMAX1(VRMIN,YRMIN*Z1**0.666667)                          
         VMIN=0.5*(VRMIN+SQRT(AMAX1(0.0,VRMIN**2-0.8*VFERMI**2)))       
         EEE=25.0*VMIN**2                                               
                                                                        
         CALL RPSTOP(Z2,EEE,PCOEF,SP)                                   
                                                                        
         POWER=0.5                                                      
         IF((Z2.EQ.6).OR.(((Z2.EQ.14).OR.(Z2.EQ.32))                    
     &   .AND.(Z1.LE.19))) POWER=0.375                                  
         SE(I)=(SP*(ZETA*Z1)**2)*(E/EEE)**POWER                         
                                                                        
       ENDIF                                                            
                                                                        
      ENDIF                                                             
                                                                        
c .......... Convert to eV-Angstom**2 ..................................
                                                                        
       SE(I)=SE(I)*10.0                                                 
                                                                        
110   CONTINUE                                                          
                                                                        
c ........... Converts stopping to eV/Angstrom .........................
                                                                        
      A=SE(1000)*ATRHO*1E-24                      ! Stopping in eV/A    
      B=SE(1000)/10.0                             ! Stopping in keV-cm2 
                                                                        
c ......... Print out calculated values on screen ......................
                                                                        
c     WRITE(6,120) Z1,LFCTR,Z2,VFERMI                                   
c     PRINT*,'Ion energy = ',EE                                         
c     PRINT*,'Stopping = ',A,' eV/A'                                    
c     PRINT*,'Stopping = ',B,'keV-cm2'                                  
                                                                        
120   FORMAT(' Z1=',I4,' LAMBDA =',F5.2,' ; Z2=',I4,' VFERMI=',F5.2)    
                                                                        
260   RETURN                                                            
      END                                                               
                                                                        
c ********* Proton electronic stopping powers **************************
                                                                        
      SUBROUTINE RPSTOP(Z2,E,PCOEF,SP)                                  
                                                                        
      IMPLICIT NONE                                                     
      REAL E,SP,PE0,PE,SL,SH,VELPWR,PCOEF(8)                            
      INTEGER Z2                                                        
                                                                        
c .... Velocity proportional stopping below velocity=PE ................
                                                                        
      PE0=25.0                                 ! Bohr velocity 25keV/amu
      PE=AMAX1(PE0,E)                                                   
      SL=(PCOEF(1)*PE**PCOEF(2))+PCOEF(3)*PE**PCOEF(4)                  
      SH=PCOEF(5)/PE**PCOEF(6)*ALOG((PCOEF(7)/PE)+PCOEF(8)*PE)          
      SP=SL*SH/(SL+SH)                         ! Proton stopping factor 
                                                                        
c .... VELPWR is the power of velocity stopping below PE0 ..............
                                                                        
      IF (E.LE.PE0) THEN                                                
        VELPWR=0.45                                                     
        IF (Z2.LE.6) VELPWR=0.25                                        
        SP=SP*(E/PE0)**VELPWR                                           
      ENDIF                                                             
                                                                        
100   RETURN                                                            
      END                                                               
                                                                        
c **********************************************************************
                                                                        
      SUBROUTINE SCOEF(Z1,MM1,M1,M2,RHO,ATRHO,VFERMI,LFCTR,Esubl,PCOEF) 
                                                                        
c ............. This routine reads the SCOEF data ......................
                                                                        
      IMPLICIT NONE                                                     
                                                                        
      INTEGER i,z,j,Z1                                                  
                                                                        
      REAL A,MM1,M1,M2,RHO,ATRHO,VFERMI,LFCTR,Esubl,PCOEF(8),M0         
                                                                        
      CHARACTER G$*70                                                   
                                                                        
c .............. Define atomic number of ion ...........................
                                                                        
      z=Z1                                                              
                                                                        
c .............. Read in ion data from file ............................
                                                                        
      OPEN(UNIT=13,FILE='scoef.dat',FORM='FORMATTED',STATUS='OLD') 
                                                                        
       READ(13,103)G$                            ! Dummy string for     
       READ(13,103)G$                            ! file title etc       
       DO 201,i=1,z                                                     
          READ(13,*) A,M0,M1,M2,RHO,ATRHO,VFERMI,LFCTR,Esubl            
201    CONTINUE                                                         
                                                                        
      CLOSE(13)                                                         
                                                                        
      MM1=M0                                                            
                                                                        
103   FORMAT(A70)                                                       
                                                                        
c ................. Read in Proton stopping data .......................
                                                                        
      OPEN(UNIT=13,FILE='scoef.psp',FORM='FORMATTED',STATUS='OLD') 
                                                                        
       READ(13,103)G$                                                   
       DO 202,i=1,z                                                     
          READ(13,*)A,(PCOEF(j),j=1,4)                                  
          READ(13,*)(PCOEF(j),j=5,8)                                    
202    CONTINUE                                                         
                                                                        
      CLOSE(13)                                                         
                                                                        
c ......... Calculate target density atoms/cm3..........................
                                                                        
      ATRHO=ATRHO*1e22                                                  
                                                                        
c ................. Return to main program .............................
                                                                        
      RETURN                                                            
      END                                                               
                                                                        
c **********************************************************************
                                                                        
c      SUBROUTINE SAVEDATA(fileid$,array,size)                           
c                                                                        
c       REAL array(100000)                                                
c       INTEGER size,i                                                   
c       CHARACTER fileid$*21                                             
c       LOGICAL ext                                                      
c                                                                        
c       INQUIRE (FILE=fileid$,EXIST=ext)                                 
c       IF (ext) THEN                                                    
c          PRINT*,fileid$,' already exists, cannot save data file.'      
c       ELSE                                                             
c          OPEN(UNIT=13,FILE=fileid$,FORM='UNFORMATTED',STATUS='NEW')    
c           WRITE(13)size                                                
c           WRITE(13)(array(i),i=1,size)                                 
c          CLOSE(13)                                                     
c          PRINT*,fileid$,' succesfully saved.'                          
c       ENDIF                                                            
c                                                                        
c      RETURN                                                            
c      END                                                               
c                                                                        
c.............Randomize seed for random number generation...............

	SUBROUTINE tinow(seed)
	 IMPLICIT REAL (a-h,o-z)
	 CHARACTER*8 t
	 INTEGER seed
	 seed=0
	 CALL TIME(t)
	 seed = (ichar(t(7:7))-48) * 10
	 seed = seed + ichar(t(8:8)) - 48
	RETURN
	END

c ........................................................................
c ......................Random Number Generator (Steve Wilson's)..........
c ........................................................................

      REAL FUNCTION RAN3(IDUM)
	  save
c         IMPLICIT REAL*4(M)
c         PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5E-7)
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.E-9)
      DIMENSION MA(55)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        MJ=MSEED-IABS(IDUM)
        MJ=MOD(MJ,MBIG)
        MA(55)=MJ
        MK=1
        DO 11 I=1,54
          II=MOD(21*I,55)
          MA(II)=MK
          MK=MJ-MK
          IF(MK.LT.MZ)MK=MK+MBIG
          MJ=MA(II)
11      CONTINUE
        DO 13 K=1,4
          DO 12 I=1,55
            MA(I)=MA(I)-MA(1+MOD(I+30,55))
            IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG
12        CONTINUE
13      CONTINUE
        INEXT=0
        INEXTP=31
        IDUM=1
      ENDIF
      INEXT=INEXT+1
      IF(INEXT.EQ.56)INEXT=1
      INEXTP=INEXTP+1
      IF(INEXTP.EQ.56)INEXTP=1
      MJ=MA(INEXT)-MA(INEXTP)
      IF(MJ.LT.MZ)MJ=MJ+MBIG
      MA(INEXT)=MJ
      RAN3=MJ*FAC
      RETURN
      END

c ........................................................................
