dedx.f

计算能量损失,重离子经过介质后的能量沉积LET,射程

*$debug
*$declare

	program calDEDX
	include 'stpwtbl.i'

	external	range,dedx
	real		range,dedx
	real	z,a,energy,aRange,aDEDX,pRange,thick,range0,eout,eloss
	real	eps
	integer	ians,loop
	integer	n,m

	call readsptbl
	call select(z,a)

	do 1000 loop=1,10000
	  n=170
	  m=4
	  eps=0.001
	  write(*,*) '(0=Energy->Range/1=Range->Engergy/2=Change/3=Stop)'
	  read(*,'(i80)') ians
 	  if(ians.eq.2)then
		call select(z,a)
	  elseif(ians.eq.0)then
	     energy=1
1100	     if(energy.gt.0)then
	        write(*,*) 'Enter energy [MeV/n] ?? (0=exit)'
	        read(*,'(f60.0)') energy
	     if(energy.eq.0) goto 1000
	     	call ene2rng(energy,aRANGE,aDEDX)
***********	     	write(*,*) aRange,aDEDX
***********		aRange=RANGE(z,a,energy)
***********		aDEDX=DEDX(Z,A,energy)
		pRange=aRange/rho
		write(*,200) energy,aDEDX,aRange,pRange
200		format(1h ,2x,'energy=',g15.7,'MeV/n',/,
     $                  2x,'dedx  =',g15.7,'MeV/(g/cm2)',/,
     $                  2x,'range =',g15.7,'g/cm2 (',g15.7,' cm)'/)
 1110		if(.true.)then
		    write(*,*) 'Enter thickness of the material in cm (0=exit)'
		    read(*,'(f60.0)') thick
		if(thick.eq.0) goto 1100
		    thick= thick*rho
		    if(thick.ge.aRange)then
			write(*,*)'the thickness is greater than the range'
		    else
			call eneloss(energy,thick,eloss,eout)
*********			write(*,*) eloss,eout
*********			range0=aRange-thick
*********			call rng2ene(range0,eout,aDEDX)
*********			eloss=(energy-eout)*a
********			write(*,*) energy,eout,eloss
			write(*,201) energy,eout,eloss
201			format(1h ,2x,'energy=',g15.7,'MeV/n',/,
     $		               2x,'e-out=',g15.7,'MeV/n',/,
     $		               2x,'en-loss =',g15.7,'MeV',/)
		    endif
		 goto 1110
	       endif
	       goto 1100
	     endif
	  elseif(ians.eq.1)then
 1300	    if(.true.)then
		write(*,*) 'Enter Range (cm)'
		read(*,'(f100.0)') pRange
	    if(pRange.eq.0) goto 1000
		aRange=pRange*rho
		call rng2ene(aRange,energy,aDEDX)
		write(*,202) aRange,pRange,energy,aDEDX
202		format(1h ,2x,'range=',g15.7,'g/cm2 (',g15.7,' cm)',/,
     $                 2x,'energy=',g15.7,'MeV/n',/,
     $	               2x,'dedx  =',g15.7,'MeV/(g/cm2)',/)
	    goto 1300
	    endif
	  else
	  	goto 1999
	  endif
 1000	continue
 1999	continue
	END


	subroutine select(z,a)
	real	z,a

	integer	not,noi

	call askstpwin(noi,z,a)
	write(*,*) 'noi,z,a=',noi,z,a
	call setstpwin(noi,z,a)

	call askstpwtg(not)

	return
	end