OPTIONS /EXTEND_SOURCE

	SUBROUTINE chargeStateYields(E,masse,Yield_plus,Yield_zero)
c       ===========================================================

c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
c Die Funktion sowie die Parameter sind uebernommen aus:
c
c	M.Gonin, R.Kallenbach, P.Bochsler: 'Charge exchange of hydrogen atoms
c	in carbon foils at 0.4 - 120 keV', Rev.Sci.Instrum. 65 (3), March 1994
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

	IMPLICIT NONE

	real E		! kinetische Energie in keV
	real masse	! in keV / c**2

	real a_zero,a_minus
	real k_Fermi,k_zero,k_minus
	real zwo_k_Fermi
	real k_Fermi_Quad,k_zero_Quad,k_minus_Quad
	real vc_minus,vc_plus,v_Bohr,v_rel

	parameter ( a_zero   = 0.953,  a_minus = 0.029 )
	parameter ( k_Fermi  = 1.178                   )	! [v_Bohr]
	parameter ( k_Fermi_Quad = k_Fermi * k_Fermi   )
	parameter ( zwo_k_fermi = 2. * k_Fermi         )
	parameter ( k_zero   = 0.991*k_Fermi           )	! [v_Bohr]
	parameter ( k_zero_Quad = k_zero * k_zero      )
	parameter ( k_minus  = 0.989*k_Fermi           )	! [v_Bohr]
	parameter ( k_minus_Quad = k_minus * k_minus   )
	parameter ( vc_minus = 0.284,  vc_plus = 0.193 )	! [v_Bohr]
	parameter ( v_Bohr   = 7.2974E-3               )	! [c]

	real Q_zero,Q_minus,D
	real Yield_minus,Yield_zero,Yield_plus

	real help1,help2,help3


c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

	if (E.LT.0) then
	    write(*,*)
	    write(*,*) 'error in subroutine ''chargeStateYields'':'
	    write(*,*) 'E = ',E,' < 0!'
	    write(*,*) '-> STOP'
	    write(*,*)
	    STOP
	endif


c Energie in Geschwindigkeit umrechnen (in Einheiten von v_Bohr):

c - klassisch:

	v_rel = SQRT(2.*E/masse) / v_Bohr

c - relativistisch:

c	help1 = 1. + E/masse
c	v_rel = SQRT(1. - 1./(help1*help1)) / v_Bohr


c Die geladenen Anteile berechnen (vgl. obige Referenz):

	help1 = v_rel*v_rel
	help2 = zwo_k_Fermi*v_rel
	Q_zero  = 1. + (k_zero_Quad  - k_Fermi_Quad - help1) / help2
	Q_minus = 1. + (k_minus_Quad - k_Fermi_Quad - help1) / help2


	help1 = a_zero * Q_zero
	help2 = a_minus * Q_minus
	help3 = (1.-Q_zero)*(1.-Q_minus)
	D = help1*(help2 + (1.-Q_minus)) + help3

	Yield_minus = help1*help2 / D
	Yield_plus  = help3 / D

	Yield_minus = Yield_minus * exp(-vc_minus/v_rel)
	Yield_plus  = Yield_plus  * exp(-vc_plus /v_rel)

	Yield_zero  = 1. - (Yield_minus + Yield_plus)

c	write(6,*) 'E	vrel	Neutral	    Plus	 Minus'
c	write(6,*) E, v_rel, yield_zero, yield_plus, yield_minus

	END