2005-11-11 12:35:21 +00:00

103 lines
2.7 KiB
C++

#include "yields.h"
#include <iomanip>
#include <fstream>
#include <iostream>
#include <stdlib.h>
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Die Funktion sowie die Parameter sind uebernommen aus:
M.Gonin, R.Kallenbach, P.Bochsler: 'Charge exchange of hydrogen atoms
in carbon foils at 0.4 - 120 keV', Rev.Sci.Instrum. 65 (3), March 1994
fIRST IMPLEMENTATION BY ANLSEM,H. IN FORTRAN
C++ CONVERSION T.K.PARAISO 04-2005
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
Yields::Yields(){;}
Yields::~Yields(){;}
void Yields::GetYields(
double E, // kinetische Energie in keV
double masse,// in keV / c**2
double yvector[3])// pointer to yields table
{
// YIELDS FUNCTIONS PRAMETERS
double a_zero,a_minus;
double k_Fermi,k_zero,k_minus;
double zwo_k_Fermi;
double k_Fermi_Quad,k_zero_Quad,k_minus_Quad;
double vc_minus,vc_plus,v_Bohr,v_rel;
// YIELDS FUNCTIONS PARAMETERs VALUES
a_zero = 0.953 ;
a_minus = 0.029 ;
k_Fermi = 1.178 ;// ! [v_Bohr]
k_Fermi_Quad = k_Fermi * k_Fermi ;
zwo_k_Fermi = 2. * k_Fermi ;
k_zero = 0.991*k_Fermi ; //! [v_Bohr]
k_zero_Quad = k_zero * k_zero ;
k_minus = 0.989*k_Fermi ; // [v_Bohr]
k_minus_Quad = k_minus * k_minus ;
vc_minus = 0.284 ;
vc_plus = 0.193 ; // [v_Bohr]
v_Bohr = 7.2974E-3 ; // [c]
// ABORT IF ENERGY NEGaTIVE-------------------------------------
if (E < 0)
{
std::cout<< "Error in method ''Yields'':"<<std::endl;
std::cout<< "E = "<< E <<" < 0!"<<std::endl;
std::cout<< "-> ABORT!"<<std::endl;
return;
}
//---------------------------------------------------------------
// VARIABLES DEFINITION
//Energie in Geschwindigkeit umrechnen (in Einheiten von v_Bohr):
// - klassisch:
v_rel = sqrt(2.*E/masse)/ v_Bohr;
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);
if(Yield_minus > exp(-vc_minus/v_rel)) Yield_minus=exp(-vc_minus/v_rel);
if(Yield_plus > exp(-vc_plus/v_rel)) Yield_plus=exp(-vc_plus/v_rel);
Yield_zero = 1. - (Yield_minus + Yield_plus);
yvector[0]=Yield_plus;
yvector[1]=Yield_zero;
yvector[2]=Yield_minus;
}