      program ELO
      
      real ion(92), iongs(92), con(92), aveb(40), stp(40,101), edens(92)
      real fr(40), pres(40), th(40), aab(40), zab(40), uv(40)
      integer ne(10), mn(40,40), na(40,40)
      integer mnp(40), mzp(40)
      integer*2 an(40,40)
      character*1 st, st1, st2, st3, st4
      character*2 gsym, ssym, psym(40)
      real iave,cfact,spwhr,eres(101,40),ene(101), ewe(101),erem(101,40)
      real*8 r0, rnge(5001), ere(101,40), ers(101), ert(101)
      
      common /ro/r0,e0,emn,der
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH
      COMMON /STOPPING/IAVE,CFACT

      
      DATA EDENS/ 8.99E-5, 17.85E-5, 0.534,  1.85,    2.34,    2.25, 
     2 1.25E-3, 1.43E-3, 1.69E-3, 9.00E-4, 0.97,    1.74, 
     3 2.702,   2.329,   1.82,    2.07,    3.214E-3,1.78E-3,
     4 0.86,    1.55,    2.5,     4.5,     5.96,    7.20, 
     5 7.20,    7.86,    8.9,     8.90,    8.92,    7.14, 
     6 5.904,   5.35,    5.727,   4.82,    2.928,   3.71E-3,
     7 1.532,   2.6,     5.51,    6.4,     8.55,    10.2, 
     8 -1.,     12.06,   12.4,    11.40,   10.5,    8.642,
     9 7.30,    7.28,    6.684,   6.25,    4.93,    5.85E-3,
     1 1.873,   3.5,     6.15,    6.7,     6.5,     6.9,
     2 -1.,     7.7,     5.22,    7.95,    8.33,    8.56, 
     3 8.76,    9.16,    9.35,    7.01,    9.74,    13.3, 
     4 16.6,    19.3,    20.53,   22.48,   22.42,   21.45,
     5 19.3,    13.55,   11.85,   11.34,   9.80,    9.24, 
     6 -1.,     9.73E-3, -1.,     5.,      -1.,     11.2, 
     7 15.4,    18.7/
      
      DATA ION/17.1,45.2,47.,63.,75.,79.,84.4,104.8,126.4,150.9,
     1 141.,149.,162.,159.,168.9,179.2,187.2,200.,189.4,195.,
     2 215.,228.,237.,257.,275.,284.,304.,314.,330.,323.,
     3 335.4,323.,354.7,343.4,360.5,368.2,349.7,353.3,365.,382.,
     4 391.3,393.,416.2,428.6,436.4,456.,470.,466.,479.,511.8,
     5 491.9,491.3,491.8,489.5,484.8,485.5,493.8,512.7,520.2,540.,
     6 537.,545.9,547.5,567.,577.2,578.,612.2,583.3,629.2,637.,
     7 655.1,662.9,682.,695.,713.6,726.6,743.7,760.,742.,768.4,
     8 764.8,761.,762.9,765.1,761.7,764.2,762.3,760.1,767.9,776.4,
     9 807.,808./

      DATA IONGS/19.,42.,0.,0.,0.,66.2,86.,99.,118.8,135.,

     1 0.,0.,0.,0.,0.,0.,170.3,180.,0.,0.,
     2 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     3 0.,0.,0.,0.,339.3,347.,0.,0.,0.,0.,
     4 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     5 0.,0.,452.4,459.,0.,0.,0.,0.,0.,0.,
     6 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     7 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     8 0.,0.,0.,0.,0.,733.1,0.,0.,0.,0.,
     9 0.,0./

      DATA CON/1.674,6.647,11.52,14.97,17.95,19.95,23.26,26.57,31.55
     1  ,33.52,
     1 38.15,40.37,44.8,46.64,51.41,53.24,58.87,66.34,64.89,66.56,
     2 74.62,79.54,84.53,86.34,91.96,92.74,97.85,97.49,105.5,108.6,
     3 115.7,120.5,124.3,131.1,132.6,139.2,141.9,145.5,147.6,151.5,
     4 154.2,159.3,164.3,167.8,170.9,176.6,179.1,186.6,191.3,197.0,
     5 202.2,211.8,210.7,218.0,220.6,228.1,230.6,232.6,234.8,239.4,
     6 244.0,250.0,252.3,261.0,263.8,269.8,273.8,277.6,280.3,287.3,
     7 290.4,296.2,300.5,305.3,309.0,315.7,319.0,324.,327.1,332.9,
     8 339.2,344.1,347.,348.6,348.6,368.5,370.1,375.1,376.8,385.1,
     9 383.4,395.3/

      const = 1./(22.4*760.)            
      der = 50.
      emn = 10.

c      write(*,*) ' Enter EMN'
c      write(*,*) ' '
c      read(*,*) emn

C***********************************************************************

C     INPUT FILE

c      open(unit = 1, file = 'elo.inp', status = 'old')
      
C     OUTPUT FILE      

      open(unit = 2, file = 'elo.out')
      
C***********************************************************************      

c      write(*,11) emn
c      write(2,11) emn
c  11  format(1x,'EMN =',f10.2)
      
C     READ INPUT DATA

C     ABSORBERS SPECIFICATIONS

C     NUMBER OF  ABSORBERS

      write(*,*) '  Notes to ELO users:'
      write(*,*) ' '
      write(*,*) '  1) Maximum of energy steps = 100'
      write(*,*) '  2) Maximum of absorbers number = 40'
      write(*,*) '  3) Maximum of projectiles number = 40'
      write(*,*) '  4) For energies less than 10 keV/u uncertainty of ra
     &nge (and energy loss)'
      write(*,*) '     may be greater then 10%, especially for lighter p
     &rojectiles or absorbers'
      write(*,*) ' '

      write(*,*) '                     *** Input data ***'
      write(*,*) ' '
      write(2,*) '                     *** Input data ***'
      write(2,*) ' '

      lala = 1

66    continue

      write(*,*) '              *** Absorbers specification ***'
      write(2,*) '              *** Absorbers specification ***'
      write(*,*) ' '
      write(2,*) ' '

      if(lala.gt.1) go to 8070

      write(*,*)'Enter number of absorbers '
      
 1011 continue
       
      write(*,*) ' '
      read(*,*,err=1011) Nab

 8070 continue
      
      i = 0
      
   77 continue
   
      if(i.ge.Nab) go to 1350

C     ABSORBER NUMBER AND ITS STATE

        write(*,101)
101     format(1x,'Choose option:'/'    G - Gas, S - Solid, C - CsI(Tl),
     & N - NaI(Tl), K - Silicon, M - Mylar')

 272  continue

      write(*,*) ' '
      read(*,'(a1)') st

      if(st.eq.'C'.or.st.eq.'c') go to 6077
      if(st.eq.'S'.or.st.eq.'s') go to 1150
      if(st.eq.'K'.or.st.eq.'k') go to 6577
      if(st.eq.'M'.or.st.eq.'m') go to 6677
      if(st.eq.'N'.or.st.eq.'n') go to 7000
      if(st.eq.'G'.or.st.eq.'g') go to 95      
      
      write(*,*) 'Invalid absorber state specification, try again '
      go to 272

 6077 continue

      i = i + 1
      
      write(*,*) 'Enter thickness of CsI detector ( < 0 - mm, > 0 - mg/
     &cm**2)'
      write(*,*) ' '
      read(*,*) thck 
      zab(i) = 54.
      aab(i) = 130.
      aveb(i) = 488.
      lz = int(zab(i))
      la = int(aab(i))
      th(i) = thck
      thck = - thck
      if(thck.gt.0.) th(i) = thck*451.
      if(thck.lt.0.) thck = - thck/451.

C     FOR SOLID ABSORBER PRESSURE IS EQUAL 0

      pres(i) = 0.
      uv(i) = 451.
      write(2,404) i, thck, th(i)
      write(*,404) i, thck, th(i)
  404 format(/,' Absorber #',i3,': CsI ',f12.4,' mm (',f12.4,' mg/cm**2)
     &')
      write(*,405) lz, la, aveb(i) 
      write(2,405) lz, la, aveb(i) 
  405 format(' Equivalent Values: Z =',i3,', A =',i4,', Av.Ion.En. =',
     &f10.3,' eV')      

C     NEXT ABSORBER

      go to 77            

 6577 continue

      i = i + 1

      write(*,*) 'Enter thickness of silicon detector ( < 0 - mm, > 0 - 
     &mg/cm**2)'
      write(*,*) ' '
      read(*,*) thck 
      zab(i) = 14.
      aab(i) = 28.
      lz = int(zab(i))
      la = int(aab(i))
      aveb(i) = 159.
      th(i) = thck
      thck = - thck
      if(thck.gt.0.) th(i) = thck*232.9
      if(thck.lt.0.) thck = - thck/232.9

C     FOR SOLID ABSORBER PRESSURE IS EQUAL 0

      pres(i) = 0.
      uv(i) = 232.9
      write(2,406) i, thck, th(i)
      write(*,406) i, thck, th(i)
  406 format(/,' Absorber #',i3,': Si ',f12.4,' mm (',f12.4,' mg/cm**2)'
     &)
C     NEXT ABSORBER

      go to 77            

 6677 continue

      i = i + 1

      write(*,*) 'Enter thickness of Mylar detector ( < 0 - mm, > 0 - mg
     &/cm**2)'
      write(*,*) ' '
      read(*,*) thck 
      zab(i) = 6.
      aab(i) = 12.
      lz = int(zab(i))
      la = int(aab(i))
      aveb(i) = 79.
      th(i) = thck
      thck = - thck
      if(thck.gt.0.) th(i) = thck*139.0
      if(thck.lt.0.) thck = - thck/139.0

C     FOR SOLID ABSORBER PRESSURE IS EQUAL 0

      pres(i) = 0.
      uv(i) = 139.
      write(2,476) i, thck, th(i)
      write(*,476) i, thck, th(i)
  476 format(/,' Absorber #',i3,': Mylar ',f12.4,' mm (',f12.4,' mg/cm**
     &2)')

C     NEXT ABSORBER

      go to 77            

 7000 continue

      i = i + 1

      write(*,*) 'Enter thickness of NaI(Tl) detector ( < 0 - mm, > 0 - 
     &mg/cm**2)'
      write(*,*) ' '
      read(*,*) thck 
      ani = 0.32
      bni = 0.68
      zab(i) = 11.*ani + 53.*bni
      aab(i) = (con(11)*ani + con(53)*bni)*.6023
      aveb(i) = 141.*ani + 491.8*bni
      th(i) = thck
      thck = - thck
      if(thck.gt.0.) th(i) = thck*367.
      if(thck.lt.0.) thck = - thck/367.

C     FOR SOLID ABSORBER PRESSURE IS EQUAL 0

      pres(i) = 0.
      uv(i) = 367.
      write(2,408) i, thck, th(i)
      write(*,408) i, thck, th(i)
  408 format(/,' Absorber #',i3,': NaI(Tl) ',f12.4,' mm (',f12.4,' mg/cm
     &**2)')
      write(*,435) zab(i), aab(i), aveb(i) 
      write(2,435) zab(i), aab(i), aveb(i) 
  435 format(' Equivalent Values: Z =',f8.2,', A =',f8.2,', Av.Ion.En. =
     &',f10.3)      
 
C     NEXT ABSORBER

      go to 77            


   95 continue

      i = i + 1

C     GAS ABSORBERS SPECIFICATIONS

C     NUMBER OF COMPONENTS IN GAS ABSORBER # I
      
      write(*,*)' Enter number of gas components'
      write(*,*) ' '
      read(*,*) nc
      write(2,532) i, nc
  532 format(/,' Absorber #',i3,':',' Gas of',i3,' components')

C     THICKNESS AND PRESSURE OF GAS ABSORBER # I
      
      write(*,104)
104   format(' Enter thickness ( < 0 - mm, > 0 - mg/cm**2 )',/,' and pre
     &ssure ( in Torr. ) of gas absorber ')
      write(*,*) ' '
      read(*,*) thck, pres(i)
      
C     GAS COMPONENTS SPECIFICATIONS

      do 1000 k = 1, nc
      
C     FRACTION AND NUMBER OF ELEMENTS IN COMPONENT # K

        write(*,102)k
102     format(' Enter fraction and number of elements in component # ',
     &i2)
      write(*,*) ' '
      read(*,*) fr(k), ne(k)
      write(2,534) k, fr(k), ne(k)
  534 format(1x,'Component #',i3,', Fraction =',g8.3,', Number of elemen
     &ts =',i3)      
      
C     SPECIFICATIONS OF ELEMENTS IN GAS MOLECULE # K

      do 1000 j = 1, ne(k)
      
C     MASS AND ATOMIC NUMBERS AND NUMBER OF ATOMS OF ELEMENT # J

      write(*,103)j
103   format(' Enter mass, atomic numbers and number of atoms',
     1' of element # ',i2)
      write(*,*) ' '
      read(*,*) mn(k,j), an(k,j), na(k,j)

      call CONV(gsym,an(k,j),ngr)
      
      if(ngr.eq.0) go to 168
      
      write(*,164) j, mn(k,j), gsym, na(k,j)
      write(2,164) j, mn(k,j), gsym, na(k,j)
  164 format(' Element #',i3,':',i4,a2,2x,i2,' atoms')
  
      go to 1000
      
  168 continue      
      
      write(*,163) j, mn(k,j), an(k,j), na(k,j)
      write(2,163) j, mn(k,j), an(k,j), na(k,j)
  163 format(' Element #',i3,': A =',i4,', Z =',i3,2x,i2,' atoms')
      
 1000 continue                              
      
      zab(i) = 0.
      aab(i) = 0.
      avg = 0.            
      avf = 0.
   
      do 1600 k = 1, nc

      zmo = 0.
      amo = 0.
      avemo = 0.
      
      do 1700 j = 1, ne(k)
      
      nzg = an(j,k)
      
      vion =iongs(nzg)      
      if(vion.eq.0.) vion = ion(nzg)
      if(mn(k,j).eq.0) mn(k,j) = con(nzg)*0.6023
      
C     EFFECTIVE Z, A AND AVERAGE IONISATION ENERGY FOR GAS MOLECULE # K     

      zmo = zmo + nzg*na(k,j)
      amo = amo + mn(k,j)*na(k,j)
      avemo  =  avemo + vion*na(k,j)*fr(k)            
      avf  =  avf + fr(k)*na(k,j)
      
 1700 continue      
                                  
      zgab = zgab + zmo*fr(k)
      agab = agab + amo*fr(k)
      avg = avg + avemo
      
 1600 continue

C     EFFECTIVE Z, A, AVERAGE IONISATION ENERGY FOR GAS ABSORBER # I 

      zab(i) = zgab
      aab(i) = agab
      aveb(i) = avg/avf
      th(i) = thck 
      thck = - thck
      if(thck.gt.0.) th(i) = aab(i)*thck*pres(i)*const
      if(thck.lt.0.) thck = - thck/(aab(i)*pres(i)*const)

      uv(i) = aab(i)*pres(i)*const
      
      write(2,300)  zab(i), aab(i), aveb(i)
      write(*,300)  zab(i), aab(i), aveb(i)
  300 format(' Equivalent Values : Z =',f8.2,', A =',f8.2,', Av.Ion.En. 
     &=',f10.3,' eV')
      write(2,310) thck, th(i), pres(i)      
  310 format(' Thickness =',f12.4,' mm (',f12.4,' mg/cm**2 )','Pressure 
     &=',f8.2,' Torr.')
  
C     NEXT ABSORBER

      go to 77      
      
 1150 continue
  
C     SOLID ABSORBERS SPECIFICATIONS 

      i = i + 1
      
C     NUMBER OF COMPONENTS IN SOLID ABSORBER # I

      write(*,105)i
 105  format(' Enter number of components in solid absorber # ',i2)
      write(*,*) ' '
      read(*,*) nc
      write(*,165) i, nc
      write(2,165) i, nc
 165  format(1x,/,' Absorber #',i3,': Solid of',i3,' components')
      


C     SPECIFICATIONS OF COMPONENTS OF SOLID ABSORBER
      
      do 1200 k = 1, nc
      
C     FIXED PARAMETER J = 1 FOR SOLID ABSORBER
      
C     MASS AND ATOMIC NUMBERS AND NUMBER OF ATOMS IN COMPONENT # K

      write(*,106)k
  106 format(' Enter mass and atomic numbers and number of atoms in comp
     &onent # ',i3)
      write(*,*) ' '
      read(*,*) mn(k,1), an(k,1), na(k,1)            
        
      call CONV(ssym,an(k,1),nsr)
        
      if(nsr.eq.0) go to 171

      write(*,174) k, mn(k,1), ssym, na(k,1)
      write(2,174) k, mn(k,1), ssym, na(k,1)
  174 format(' Element #',i2,':',i4,a2,2x,i2,' atoms')
  
      go to 1200
      
  171 continue      
      
      write(*,173) k, mn(k,1), an(k,1), na(k,1)
      write(2,173) k, mn(k,1), an(k,1), na(k,1)
  173 format(' Element #',i3,': A =',i4,', Z =',i3,2x,i2,' atoms')

 1200 continue
 
C     THICKNESS OF SOLID ABSORBER # I

      write(*,107)
107   format(' Enter thickness of solid absorber ( < 0 - mm, > 0 - mg/cm
     &**2)')
      write(*,*) ' '
      read(*,*) thck
      zsab = 0.
      asab = 0.
      avs = 0.
      fs = 0.
      
      do 1900 k = 1, nc
      
      nzs = an(k,1)
      
      if(mn(k,1).eq.0.) mn(k,1) = con(nzs)*0.6023
      
C     EFFECTIVE Z, A, AVERAGE IONISATION ENERGY FOR SOLID COMPONENT # K

      zsab = zsab + nzs*na(k,1)      
      asab = asab + mn(k,1)*na(k,1)
      avs = avs + ion(nzs)*na(k,1)
      fs = fs + na(k,1)           

 1900 continue
      
C     EFFECTIVE Z, A, AVERAGE IONISATION ENERGY FOR SOLID ABSORBER # I       

      zab(i) = zsab/fs
      aab(i) = asab/fs
      aveb(i) = avs/fs

C     FOR SOLID ABSORBER PRESSURE IS EQUAL 0

      pres(i) = 0.
      uv(i) = edens(nzs)*100.
      
C     THICKNESS FOR SOLID ABSORBER # I (IN MG/CM**2)

      th(i) = thck
      thck = - thck
      if(thck.gt.0.) th(i) = thck*edens(nzs)*100. 
      if(thck.lt.0.) thck = - thck/edens(nzs)*0.01
    
 2000 continue
 
      write(2,400)  zab(i), aab(i), aveb(i)
      write(*,400)  zab(i), aab(i), aveb(i)
  400 format(' Equivalent Values : Z =',f8.2,', A =',f8.2,', Av.Ion.En. 
     &=',f10.3,' eV')
      write(2,410) thck, th(i)
      write(*,410) thck, th(i)
  410 format(' Thickness =',f12.4,' mm (',f12.4,' mg/cm**2 )') 
 
C     NEXT ABSORBER

      go to 77            

 1350 continue

      if(kp.eq.1.or.lala.eq.1) go to 1352
 
      write(*,*) ' New projectile ? (y or n)'
      
 1351 continue      
 
      write(*,*) ' '
      read(*,'(a1)') st3
      if(st3.eq.'y') go to 1352
      if(st3.eq.'n') go to 1444
      write(*,*) ' Wrong letter, try again !'
 
      go to 1351
      
 1352 continue
 
C     PROJECTILES SPECIFICATION
      write(*,*) ' '
      write(*,*) '                  *** Projectile specification ***'
      write(2,*) ' '
      write(2,*) '                  *** Projectile specification ***'

      if(lala.gt.1) go to 8080
      
      write(*,*) 'Enter number of projectiles'

C     NUMBER OF PROJECTILES
c      read(1,*) Np
      write(*,*) ' '
      write(2,*) ' '
      read(*,*) Np

 8080 continue

      do 1400 i = 1, Np

C     MASS AND ATOMIC NUMBERS OF PROJECTILE # I
      write(*,241) i
 241  format(1x,' Mass and atomic numbers of projectile #',i2)
      write(*,*) ' '
      read(*,*) mnp(i), mzp(i)
      
      call CONV(psym(i),mzp(i),npr)
      
      if(npr.eq.0) go to 911
      
      write(*,912) i, mnp(i), psym(i)
      write(2,912) i, mnp(i), psym(i)
  912 format(' Projectile #',i3,':',i4,a2)
  
      go to 1400
            
  911 continue
        
      write(*,247) i, mnp(i), mzp(i)
      write(2,247) i, mnp(i), mzp(i)
 247  format(1x,'Projectile #',i3,': A =',i4,', Z =',i3)
  
 1400 continue
 
C     PROJECTILE ENERGY

 1444 continue

      if(ke.eq.1.or.lala.eq.1) go to 1447

      write(*,*) ' New energy ? (y or n)'
      
 1446 continue
       
      write(*,*) ' '
      read(*,'(a1)') st4
      if(st4.eq.'y') go to 1447
      if(st4.eq.'n') go to 99
      write(*,*) ' Wrong letter, try again !'
      
      go to 1446
      
 1447 continue
 
      write(*,*) ' '
      write(2,*) ' '
      write(*,*) '                     *** Energy specification ***'
      write(2,*) '                     *** Energy specification ***'
      write(2,*) ' '

      if(lala.gt.1) go to 8090
     
      write(*,*) 'Enter indeks, if indeks = 0 - total energy, if indeks 
     &= 1 - energy per nucleon' 

C     READ INDEKS
C     INDEKS = 0 - TOTAL ENERGY
C     INDEKS = 1 - ENERGY PER NUCLEON
c      read(1,*) ind

  924 continue
  
      write(*,*) ' '
      read(*,*) ind

 8090 continue
      
      if(ind.eq.0) go to 921
      if(ind.eq.1) go to 922
      
      write(*,923) ind
  923 format(' Wrong indeks value, ind =',i2,', don`t worry, try again !
     &')      
      
      go to 924
      
  921 continue
  
      write(2,*) ' Total energy [ MeV ]'
      
      go to 925
            
  922 write(2,*) ' Energy per nucleon [ MeV/u ]'
  
  925 continue
  
C     MINIMUM, MAXIMUM AND STEP OF ENERGY

   82 continue

      write(*,*) 'Enter minimum, maximum and step of energy'
      write(*,*) ' '
      read(*,*) Emin, Emax, Estep
      write(2,*) ' '
      write(2,930) Emin, Emax, Estep
  930 format(' Emin =',f12.3,'  Emax =',f12.3,'  Estep =',f12.3)      
      if(Estep.eq.0.) go to 88
      Nen = int((Emax-Emin)/Estep)+1      
      enk = Emin + (Estep - 1)*real(Nen)
      if(enk.lt.Emax) Nen = Nen + 1
      if(Nen.le.101) go to 84
      write(*,*) ' Number of energy steps greater than 101 !'
      write(*,*) ' Define energy once again, please'
      go to 82
 84   continue
      go to 99
 88   continue      
      Nen = 1
 99   continue      

      if(lala.gt.1) go to 8232
      
C     OUTPUT SPECIFICATION
      write(*,*) 'Choose output specification - enter iout'
      write(*,*) 'iout = 1 - only energy losses '
      write(*,*) 'iout = 2 - energy losses and ranges '
      write(*,*) 'iout = 3 - energy losses, ranges and stopping powers '
       
C     READ IOUT
C     IOUT = 1 - ONLY ENERGY LOSSES
C     IOUT = 2 - ENERGY LOSSES AND RANGES
C     IOUT = 3 - ENERGY LOSSES, RANGES AND STOPPING POWERS
c      read(1,*) iout
      write(*,*) ' '
      read(*,*) iout
      
C     END OF INPUT DATA READING
      write(*,*) ' End of input data reading'

      write(*,*) 'Output data will be written on file ELO.OUT'
      write(*,*) '              '
      write(*,*) '                   *** Results ***'
      write(2,*) '              '
      write(2,*) '                   *** Results ***'
      write(2,*) ' '

 8232 continue
                   
C***********************************************************************

      if (iout.ne.3) go to 4000
      
C     CALCULATE STOPPING POWERS

C     FOR PROJECTILE # I

      write(*,*) '                *** Stopping powers ***'
      write(2,*) '                *** Stopping powers ***'

      do 3900 i = 1, Np
      
      write(2,3500) mnp(i),psym(i)       
      write(*,3500) mnp(i),psym(i)       
 3500 format(/,15x, ' Stopping powers for projectile: ',i3,a2)
      
      zp = float(mzp(i))
      ap = float(mnp(i))
      izp = mzp(i)
      
C     FOR ENERGY # L

      if(ind.eq.0) go to 3501
      
      write(2,3300)
      write(*,3300)
      
      go to 3502

 3501 continue
 
      write(2,3301)
      write(*,3301)
     
 3502 continue
        
 3300 format(/, ' Energy [MeV/u]       Stopping power [MeV/(mg/cm**2)] i
     &n absorber:')      

 3301 format(/, ' Energy [MeV]         Stopping power [MeV/(mg/cm**2)] i
     &n absorber:')      

      do 3000 l = 1, Nen

      ese = Emin + float(l-1)*Estep
      
C     CONVERT ENERGY TO keV/u

      if (ind.eq.0) ese = ese/ap
      
      en = ese*1000.                  

C     AND FOR ABSORBER # K      

       
      do 3600 k = 1, Nab
            
      zt = zab(k)
      at = aab(k)
      izt = int(zt+.5)
      iave = aveb(k)*1.e-06
      
      stp(k,l) = stopp(en)      
      
 3600 continue      

 3000 continue
 
      do 3900 k = 1, Nab

      do 3900 l = 1, Nen

      ese = Emin + float(l-1)*Estep
      
      if (ind.eq.0) ese = ese/ap

      uk = real(k)
      wk = uk/5.
      iwk = int(wk)
      zk = real(iwk)

      n5 = Nab/5
      nd = Nab - n5*5
      nn = Nab - nd + 1

      if(nd.eq.0) nn = k - 4

      if(l.eq.1.and.k.eq.Nab) go to 3615
      
      if(l.eq.1.and.zk.eq.wk) go to 3614

      go to 3613
      
 3615 continue

      write(2,3704) (ki , ki = nn, k)      
      write(*,3704) (ki , ki = nn, k)            
 3704 format(1x,/,12x,5(4x,i3,5x:))           
      write(*,*) ' '      
      write(2,*) ' '      
      
      go to 3613
      
 3614 continue
       
      write(2,3700) (ki , ki = k-4, k)      
      write(*,3700) (ki , ki = k-4, k)            
 3700 format(1x,/,12x,5(4x,i3,5x))           
      write(2,*) ' ' 
      write(*,*) ' ' 
      
 3613 continue
      
      if(k.ne.Nab) go to 3632

      if(ind.eq.0) go to 3631
      
      write(2,3804) ese, (stp(jk,l), jk = nn, k)   
      write(*,3804) ese, (stp(jk,l), jk = nn, k)   
      
      go to 3602
      
 3631 continue
    
      ese = ese*ap      

      write(2,3804) ese, (stp(jk,l), jk = nn, k)   
      write(*,3804) ese, (stp(jk,l), jk = nn, k)   
 3804 format(1x,g11.4,2x,5(g11.4,1x:))           
      
 3632 continue      
      
      if(zk.ne.wk) go to 3900

      if(ind.eq.0) go to 3601
      
      write(2,3800) ese, (stp(jk,l), jk = k-4, k)   
      write(*,3800) ese, (stp(jk,l), jk = k-4, k)   
      
      go to 3602
      
 3601 continue
    
      ese = ese*ap      

      write(2,3800) ese, (stp(jk,l), jk = k-4, k)   
      write(*,3800) ese, (stp(jk,l), jk = k-4, k)   
      
 3602 continue      
      
 3800 format(1x,g11.4,2x,5(g11.4,1x:))
 
 3900 continue
       
C***********************************************************************

 4000 continue
 
      if (iout.lt.2) go to 5000
      
C     CALCULATE RANGES

C     IN ABSORBER # I

      do 4700 i = 1, Nab
      
      write(2,4100) i
      write(*,4100) i
 4100 format(/,20x,'Ranges for absorber #',i2)
 
      zt = zab(i)
      at = aab(i)
      izt = int(zt+.5)
      iave = aveb(i)*1.e-06

      if(ind.eq.0) go to 4224
             
      write(2,4200)
      write(*,4200)
 4200 format(/, '   Energy [MeV/u]    Range [mg/cm**2] [-mm] for project
     &ile:')      
     
      go to 4225
     
 4224 continue
 
      write(2,4207)
      write(*,4207)
 4207 format(/, '   Energy [MeV]      Range [mg/cm**2]  [-mm] for project
     &ile:')      
      
 4225 continue
       

C     FOR PROJECTILE # K

      do 4600 k = 1, Np
      
      zp = mzp(k)
      ap = mnp(k)
      izp = mzp(k)

      r0 = 0.d0
      cap = ap/1000.
      el = 0.0

C     AND FOR ENERGY # L      
       
      do 4600 l = 1, Nen
      
      ene(l) = Emin + float(l-1)*Estep
      
C     CONVERT ENERGY TO keV/u

      ewe(l) = ene(l)

      if (ind.eq.0) ene(l) = ene(l)/ap
      
      en = ene(l)*1000.                   
      de = (en - el)/4.
      
      eres(l,k) = range(en,de)*cap
      erem(l,k) = - eres(l,k)/uv(i)
      
      el = en

 4600 continue      

       do 4700 k = 1, Np
 
       uk = real(k)
       wk = uk/3.
       iwk = int(wk)
       zk = real(iwk)
       
       do 4700 l = 1, Nen
       
       if(l.eq.1.and.zk.eq.wk) go to 4721
       
       if(k.eq.Np) go to 4730
       
       go to 4722
       
 4721  continue

       write(2,4300) (mnp(kx), psym(kx), kx = k-2, k)      
       write(*,4300) (mnp(kx), psym(kx), kx = k-2, k)      
       write(2,*) ' ' 
 4300  format(1x,/,12x,3(8x,i3,a2,9x))           
       
 4722  continue       
 
       if(zk.ne.wk) go to 4700
        
       write(2,4800) ewe(l), (eres(l,ks), erem(l,ks), ks = k-2, k)   
       write(*,4800) ewe(l), (eres(l,ks), erem(l,ks), ks = k-2, k)   
 4800 format(g11.4,2x,3(g10.3,1x,g10.3,1x))
 
       go to 4700
      
 4730  continue
 
       n3 = Np/3
       nd = Np - n3*3
       nn = Np - nd + 1

       if(nd.eq.0) nn = k - 2

       if(l.ne.1) go to 4740

       write(2,4350) (mnp(kx), psym(kx), kx = nn, k)      
       write(*,4350) (mnp(kx), psym(kx), kx = nn, k)      
       write(2,*) ' ' 
 4350  format(1x,/,12x,3(8x,i3,a2,9x:))           

 4740  continue
       
       write(2,4803) ewe(l), (eres(l,ks), erem(l,ks), ks = nn, k)   
       write(*,4803) ewe(l), (eres(l,ks), erem(l,ks), ks = nn, k)   
 4803 format(g11.4,2x,3(g10.3,1x,g10.3,1x:))
       
 4700 continue

C***********************************************************************

 5000 continue
 
C     CALCULATE ENERGY LOSSES

C     FOR PROJECTILE # I      
      
      do 6900 i = 1, Np
      
      ap = float(mnp(i))
      ca = ap/1000.
      zp = float(mzp(i))
      izp = mzp(i)
      
      do 5100 in = 1, Nen
      
      ene(in) = Emin + float(in-1)*Estep
      est = Estep
      
C     CONVERT ENERGY TO keV/u

      if (ind.eq.0) ene(in) = ene(in)/ap
      if (ind.eq.0) est = Estep/ap
      
      ene(in) = ene(in)*1000.                   
      
      est = est*1000.
      
      eres(in,1) = ene(in)
      
 5100 continue
 
      eh = ene(Nen)
      de = amax1(50.,est/10.) 
      de = amin1(de,ene(1)/10.)
      nes = ifix(eh/de)
      nes = min0(nes,5000)
      de = eh/nes
      el = de
      
C     IN ABSORBER # J
      
      do 5900 j = 1, Nab
      
      j1 = j + 1
      zt = zab(j)
      at = aab(j)
      izt = int(zt+.5)
      iave = aveb(j)*1.e-06
      thk =  th(j)
      r0 = 0.d0
      r = range(el/2.,de)
      en = el
      
      do 5400 k = 1, nes+1
      
      rnge(k) = range(en,de)*ca
      
      en = en + de
      
 5400 continue
 
      do 5900 k = 1, Nen
      
      en = eres(k,j)
      
      dkk = (en-el)/de+1.
      kk = ifix(dkk)
      
      if(kk.gt.0) go to 5600
      
      kk = 1
      rng = 0.
      
      go to 5700
      
 5600 continue
      
      rng = rnge(kk)
      
 5700 continue
 
      dif = dkk - kk
      rng = (rnge(kk+1) - rng)*dif + rng
      rngl = rng - thk
       
      if(rngl.gt.0.) go to 5800
      
      eres(k,j1) = 0.
       
      go to 5900
      
 5800 continue
 
      kkl = kk
      
 6000 continue
 
      if(rngl.gt.rnge(kkl)) go to 6100
      
      kkl = kkl - 1
      
      if(kkl.gt.1) go to 6000
      
      rngo = 0.
      
      go to 6200
      
 6100 continue
 
      rngo = rnge(kkl)
      
 6200 continue
   
      difr = rnge(kkl+1) - rngo
      
      eres(k,j1) = (kkl-1.+(rngl - rngo)/difr)*de + el
      
 5900 continue
 
      if(ind.eq.1) go to 5732
 
      write(2,5750) mnp(i), psym(i)
      write(*,5750) mnp(i), psym(i)
 5750 format(/,18x,'Energy losses [MeV] for projectile:',i3,a2,/)      
      write(2,5950)
      write(*,5950)
 5950 format(1x,' E0 [MeV]    ',18x,'   DE in absorbers #:')      
      write(2,*) ' '
      write(*,*) ' '
      
      go to 5634
      
 5732 continue      
 
      write(2,5754) mnp(i), mzp(i)
      write(*,5754) mnp(i), mzp(i)
 5754 format(/,18x,'Energy losses [MeV] for projectile:',i3,',',i3,/)      
      write(2,5954)
      write(*,5954)
 5954 format(1x,' E0 [MeV/u] ',18x,'   DE in absorbers #:')      
      write(2,*) ' '
      write(*,*) ' '
      
 5634 continue      

      do 6400 k = 1, Nen
      
      ers(k) = eres(k,1)*ca

      ert(k) = eres(k,Nab+1)*ca
      
      if(ind.eq.1) ers(k) = ers(k)/ap

      do 6400 l = 1, Nab
      
      ere(k,l) = (eres(k,l) - eres(k,l+1))*ca
      
 6400 continue
 
      do 6900 k = 1, Nab

      uk = real(k)
      wk = uk/5.
      iwk = int(wk)
      zk = real(iwk)

      n5 = Nab/5
      nd = Nab - n5*5
      nn = Nab - nd + 1

      if(nd.eq.0) nn = k - 4

      do 6900 l = 1, Nen

      if(k.eq.Nab) go to 6910
      
      if(l.eq.1.and.zk.eq.wk) go to 6920
      
      go to 6930
      
 6920 continue
 
      write(2,5854) (kk, kk = k-4, k)
      write(*,5854) (kk, kk = k-4, k)
 5854 format(12x,5(4x,i3,4x))      
      write(*,*) ' '
      write(2,*) ' '
 
 6930 continue
  
      if(zk.ne.wk) go to 6900
        
      write(2,6500) ers(l), (ere(l,kl), kl = k-4, k)
      write(*,6500) ers(l), (ere(l,kl), kl = k-4, k)
 6500 format(1x,g11.4,2x,5(g10.4,1x))
 
      go to 6900
      
 6910 continue      

      if(l.ne.1) go to 6970
      
      write(2,5859) (kk, kk = nn, k+1)
      write(*,5859) (kk, kk = nn, k+1)
 5859 format(12x,6(4x,i3,4x:))      
      write(*,*) ' '
      write(2,*) ' '

 6970 continue
 
      write(2,6508) ers(l), (ere(l,kl), kl = nn, k), ert(l)
      write(*,6508) ers(l), (ere(l,kl), kl = nn, k), ert(l)
 6508 format(1x,g11.4,2x,6(g10.4,1x:))
 
 6900 continue
 
      write(*,*) ' Enter option:'
      write(*,*) '              a - new absorber'
      write(*,*) '              p - new projectile'
      write(*,*) '              e - new energy'
      write(*,*) '              q - quit program'

      lala = lala + 1
      ke = 0
      kp = 0
      
 8307 continue
       
      write(*,*) ' '
      read(*,'(a1)') st2
      if(st2.eq.'a') go to 66
      if(st2.eq.'p') go to 7350
      if(st2.eq.'e') go to 7444
      if(st2.eq.'q') go to 9000
      write(*,*) ' Wrong letter, try again !'

      go to 8307

 7444 ke = 1
      go to 1444
 
 7350 kp = 1
      go to 1350

 9000 continue      

      stop
      end

      FUNCTION STOPP(EA)
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH
      REAL*4 NUSPWR
      SPWRH=ELSPWR(EA)
      STOPP=SPWRH
      IF (IZP.GT.1) STOPP=SPWRHI(EA)
      IF (EA.GT.2000.) RETURN
      STOPP=STOPP+NUSPWR(EA)
      RETURN
      END

      FUNCTION ELSPWR(EA)
      
      REAL*4 IAVE,C2S(92)
      REAL*4 A0S(10),A1S(10),A2S(10),A3S(10),A4S(10)
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH
      COMMON /STOPPING/IAVE,CFACT
      
      DATA C2S/1.44,1.397,1.6,2.59,2.815,2.989,3.35,3.,2.352,2.199,
     1 2.869,4.293,4.739,4.7,3.647,3.891,5.714,6.5,5.833,6.252,
     2 5.884,5.496,5.055,4.489,3.907,3.963,3.535,4.004,4.175,4.75,
     3 5.697,6.3,6.012,6.656,6.335,7.25,6.429,7.159,7.234,7.603,
     4 7.791,7.248,7.671,6.887,6.677,5.9,6.354,6.554,7.024,7.227,
     5 8.48,7.81,8.716,9.289,8.218,8.911,9.071,8.444,8.219,8.,
     6 7.786,7.58,7.38,7.592,6.996,6.21,5.874,5.706,5.542,5.386,
     7 5.505,5.657,5.329,5.144,5.851,5.704,5.563,5.034,5.46,4.843,
     8 5.311,5.982,6.7,6.928,6.979,6.954,7.82,8.448,8.609,8.679,
     9 8.336,8.204/

      DATA A0S/5.05,4.41,7.65,11.57,14.59,19.07,17.11,15.0,17.8,20.4/
      DATA A1S/2.05,1.88,2.99,4.39,5.46,7.02,6.10,5.15,6.18,6.97/
      DATA A2S/.304,.281,.419,.598,.733,.931,.778,.626,.766,.874/
      DATA A3S/.020,.018,.025,.035,.042,.053,.043,.032,.041,.048/
      DATA A4S/4.7,4.2,5.5,7.5,9.0,11.2,8.6,6.0,7.9,9.5/
      
      DATA TMC2,PI,E4/1.02195,3.14159,2.074E-05/     
      
C       ROUTINE TO CALCULATE THE ELECTRONIC STOPPING POWER SE(ZH*)
C      OF HYDROGEN IN ANY MATERIAL USING EQ'S OF ZIEGLER.  STOPPING
C      POWER FOR OTHER IONS ARE THEN GIVEN BY SE(HI)=SE(H)*((ZHI*)/
C      (ZH*))**2 WHERE ZH* AND ZHI* ARE THE EFFECTIVE CHARGES OF
C      HYDROGEN AND THE H.I. RESPECTIVELY.
C       IZT=Z OF TARGET
C       EA=KEV/AMU OF PARTICLE
C       TMC2=2 MC**2 TWICE ELECTRON REST MASS
C       BETA2=(V/C)**2 SQUARE OF LAB. VELOCITY.
C       E4=E**4  E=ELECTRON CHARGE.
C       AI,I=0,1,2,3,4 = APPROXIMATE BETHE SHELL CORRECTION PARAMETERS.
C       CZT=APPROXIMATE BETHE SHELL CORRECTION.
C       IAVE=AVERAGE IONIZATION ENERGY.

C       INITIALIZATION

      C2=C2S(IZT)
      INDX=(IZT+15)/10
      A0=-A0S(INDX)
      A1= A1S(INDX)
      A2=-A2S(INDX)
      A3= A3S(INDX)
      A4=-A4S(INDX)*1.E-04
C       CONVERSION FACTOR TO MEV/(MG/CM**2) (=.6023/A, A=TGT MASS #)
      CFACT=.6023/AT
  100 CONTINUE

C       LOW ENERGY STOPPING USING VARELAS-BIERSACK FORMULA
C      (PAGES 10,16 VOL.3 ANDERSON & ZIEGLER).  FOR HYDROGEN
C      PROJECTILE OF 10-1000 KEV.  ALL COEFFICIENS BUT C2 ARE
C      APPROXIMATED.

      IF(EA-1000.) 150,150,200
  150 SLOW=C2*EA**.45
      ear = 1.+500./EA+2.195E-06*EA/IAVE
      if(ear.gt.0.0) go to 155
      write(*,*)'Alog argument less than 0'
      write(*,*)'ear =', ear
      write(*,*)'ea =', ea
      pause
  155 continue      
      SHIGH=(243.-.375*ZT)*ZT/EA*ALOG(ear)
      ELSPWR=SLOW*SHIGH/(SLOW+SHIGH)*CFACT
      RETURN
  200 CONTINUE

C       CALCULATION OF STOPPING POWER USING BETHE STOPPING (SEE
C      VOL. 5 OF ZIEGLER).  STOPPING POWER FOR IDEAL UNIT CHARGE.

C       APPROXIMATE CALCULATION OF BETHE SHELL CORRECTION. MAKES
C      <10% CORRECTION OVER 1-30MEV/A.

      CZT=0.
      IF (EA-40000.) 240,240,250
  240 continue
      if(ea.gt.0.) go to 35
      write(*,*)' Alog argument less than 0'
      write(*,*)'ea =', ea
      pause
   35 continue     
      ALE=ALOG(EA) 
      CZT=A0+(A1+(A2+(A3+A4*ALE)*ALE)*ALE)*ALE
  250 CONTINUE
      BETA2=1.-1./(1.+EA/931189)**2
      betka = BETA2*TMC2/((1.-BETA2)*IAVE)
      if(betka.gt.0.) go to 255
      write(*,*)'Alog argument less than 0'
      write(*,*)'betka =', betka
      write(*,*)'ea =', ea
      pause
  255 continue      
      SPWR=8.*PI*E4*ZT/(BETA2*TMC2)*(ALOG(betka)-BETA2-CZT)

C       EFFECTIVE CHARGE FOR HYDROGEN.
      if(ea.ge.0.) go to 285
      write(*,*)'An argument of sqrt less than 0'
      write(*,*)'ea =', ea
      pause
  285 continue      
      EX=0.2*SQRT(EA)+.0012*EA+1.443E-05*EA*EA
      ZH=1.
      IF (EX.LT.20.) ZH=1.-EXP(-EX)

C       CORRECTION FOR TARGET DEPENDENCE OF EFFECTIVE CHARGE.
C      MAKES <5% CORRECTION OVER 1-2 MEV/A
      ZFACT=1.
      IF (EA-1999.) 290,290,300
  290 B=1.
      IF (IZT.LT.35) B=(ZT-1.)/34.
      ZFACT=1.+B*(0.1097-5.561E-05*EA)
  300 CONTINUE

C       STOPPING POWER OF HYDROGEN
      ELSPWR=SPWR*ZH*ZH*CFACT*ZFACT
      RETURN
      END
      
      FUNCTION SPWRHI(EA)
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH

C       ROUTINE TO CALCULATE THE ELECTRONIC STOPPING POWER OF A HEAVY
C      ION USING THE HI'S EFFECTIVE CHARGE AND THE STOPPING POWER FOR
C      HYDROGEN.

      SPWRHI=SPWRH
      IF (IZP.LT.2) RETURN
      IF (IZP-3) 30,30,100
   30 continue
      if(ea.gt.0.) go to 35
      write(*,*)' Alog argument less than 0'
      write(*,*)'ea =', ea
      pause
   35 continue      
      X=ALOG(EA)
      EX=(7.6-X)**2
      GAMMA=1.
      IF (EX.LT.20.) GAMMA=1.+(.007+.00005*ZT)*EXP(-EX)
      IF (IZP-3) 40,50,100
   40 EX=.7446+.1429*X+.01562*X*X-.00267*X**3+1.325E-06*X**8
      ZHEZH=2.*GAMMA
      IF (EX.LT.20.) ZHEZH=ZHEZH*(1.-EXP(-EX))

C       STOPPING POWER FOR HE ION.
      SPWRHI=SPWRH*ZHEZH*ZHEZH
      RETURN

   50 ZLIZH=3.*GAMMA
      EX=.7138+.002797*EA+1.348E-06*EA*EA
      IF (EX.LT.20.) ZLIZH=ZLIZH*(1.-EXP(-EX))

C       STOPPING POWER FOR A LITHIUM ION.
      SPWRHI=SPWRH*ZLIZH*ZLIZH
      RETURN

  100 B=.886/5.*SQRT(EA)*ZP**(-.666666)
      A=B+.0378*SIN(1.5708*B)
      ZHIZH=ZP
      IF (A.LT.20.) ZHIZH=ZP*(1.-EXP(-A)*(1.034-.1777*EXP(-.08114*ZP)))

C       STOPPING POWER OF A HEAVY ION.
      SPWRHI=SPWRH*ZHIZH*ZHIZH
      RETURN
      END
      
      FUNCTION NUSPWR(EA)
      REAL*4 NUSPWR,IAVE
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH
      COMMON /STOPPING/IAVE,CFACT

C       ROUTINE TO CALCULATE THE NUCLEAR (COULOMB) STOPPING
C      POWER SN(E).  MAKES <5% CORRECTION AT 1-2MEV/A.
C       REDE=REDUCED ENERGY
C       SN=REDUCED NUCLEAR STOPPING POWER
C       NUSPWR=NUCLEAR STOPPING POWER IN MEV/(MG/CM**2)

      axe = zt**.666666 + zp**.666666
      
      if(axe.gt.0.) go to 60

      write(*,*) 'zt =', zt      
      write(*,*) 'zp =', zp      
      write(*,*) 'axe =', axe
      write(*,*) 'An attempt to calculate sqrt for negative number !'
      pause
      
   60 continue
         
      AZ=(AT+AP)*SQRT(axe)

      REDE=32.53*AT*AP*EA/(ZT*ZP*AZ)
      
      if(rede.gt.-1.) go to 70

      write(*,*) 'at =', at      
      write(*,*) 'ap =', ap      
      write(*,*) 'ea =', ea      
      write(*,*) 'az =', az      
      write(*,*) 'rede =',rede
      write(*,*) 'An attempt to calculate log for negative number !'
      pause
      
   70 continue
         
      SN=0.5*ALOG(1.+REDE)/(REDE+0.10718*REDE**.37544)
   
      NUSPWR=SN*(8.462*ZT*ZP*AP)*CFACT/AZ
   
      RETURN
      END
      
      FUNCTION RANGE(EA,DE)
      REAL*8 K1,K2,K4
      REAL*8 R0,R

      common /ro/r0,e0,emn,der      
      COMMON /PARTICLE/IZT,IZP,ZT,ZP,AT,AP,SPWRH
            

C       ROUTINE TO CALCULATE THE REDUCED RANGE R BY SOLVING THE
C      DIFF. EQ. DR/DEA=1/S(E) USING THE 4TH ORDER RUNGE-KUTTA
C      METHOD. THE QUANTITIES ARE DEFINED AS:
C       R=RANGE/AP*1000 (1000 SINCE W IS IN KEV BUT S(E) IS IN MEV)
C       AP=PROJECTILE MASS
C       EA=E/AP ENERGY PER NUCLEON IN KEV
C       S(E)=STOPPING POWER FOR PROJECTILE IN MEV/MG/CM**2

C       CALL FIRST TIME WITH R0=0. TO OBTAIN RANGE AT E0.  CALL
C      SUBSEQUENTLY WITH R0=R FROM PREVIOUS CALL.

      R=R0

      IF (R0) 40,40,100

   40 if(ea.lt.100..and.ea.ge.10.) em = 0.4
      if(ea.ge.100..and.ea.lt.500.) em = 2.5
      if(ea.ge.500.) em = 25.
      if(ea.lt.10..and.ea.ge.1.) em = .04
      if(ea.lt.1.) em = ea/4.
c	w linii ponizej byl komentarz 'clp'. Bez komentarza czasem pada, z
c	komentarzem zle liczy dla duzych energii.
      EM=AMIN1(EA,EMN)
      E0=EM/2.
      R=(E0/2.)/STOPP(E0)

  100 DDE=AMIN1(DER,DE)
      N=(EA-E0)/DDE
      IF(N.LT.4) N=4
      W=(EA-E0)/N
      K1=1./STOPP(E0)
      DO 200 I=1,N
      K2=1./STOPP(E0+W/2.)
      K4=1./STOPP(E0+W)
      R=R+(K1+4.*K2+K4)*W/6.
      K1=K4
  200 E0=E0+W
      R0=R
      RANGE=R
      RETURN
      END

      subroutine CONV(symb,znum,ner)
      character*2 symb, nucn(103)
      integer*2 znum
      
      DATA NUCN/
     1'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na','Mg'
     2,'Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Tl','V ','Cr'
     3,'Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr'
     4,'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd'
     5,'In','Sn','Sb','Te','I ','Xe','Cs','Ba','La','Ce','Pr','Nd'
     6,'Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf'
     7,'Ta','W ','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po'
     8,'At','Rn','Fr','Ra','Ac','Th','Pa','U ','Np','Pu','Am','Cm'
     9,'Bk','Cf','Es','Fm','Md','No','Lw'/
     
      iz = 0
      
      do 10 i = 1, 103
 
      if(znum.eq.i) iz = i
            
   10 continue

      if(iz.ne.0) go to 15
        
      write(*,*) '                 NO SUCH SYMBOL IN NUCLEI ARRAY'
 
      ner = 0
      symb = '  '
 
      go to 30
      
   15 continue
   
      symb = nucn(iz)
      ner = 1
      
   30 continue      
      
      return
      end
      
