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 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.3,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.3,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.3,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.3,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