14 _lifetime ("Lifetime", 1.0, 0.0),
15 _frequency("Frequency", 0.0, 0.0),
16 _sigma ("
Sigma", 1.0, 0.0),
17 _offset ("Offset", 0.0),
23 _lifetime (right._lifetime),
24 _frequency (right._frequency),
25 _sigma (right._sigma),
26 _offset (right._offset),
67 static const double sqrtTwo = sqrt(2.0);
71 double x = argument-
offset;
75 double expG=0.0, asymm=0.0;
78 expG = exp((sigma*sigma +2*tau*(x))/(2.0*tau*tau)) *
79 erfc((sigma*sigma+tau*(x))/(sqrtTwo*sigma*tau))/(2.0*tau);
85 if (!std::isfinite(expG)) {
92 expG = exp((sigma*sigma +2*tau*(-x))/(2.0*tau*tau)) *
93 erfc((sigma*sigma+tau*(-x))/(sqrtTwo*sigma*tau))/(2.0*tau);
103 if (!std::isfinite(expG)) {
117 asymm = expG*(1/(1+tau*tau*freq*freq));
119 else if (sigma==0.0) {
121 if (x>=0) asymm= (expG*cos(freq*x));
124 if (x>=0) asymm= (expG*sin(freq*x));
128 std::complex<double> z(freq*sigma/sqrtTwo, (sigma/tau-x/sigma)/sqrtTwo);
131 asymm= 2.0*nwwerf(z).real()/tau/4.0*exp(-x*x/2.0/sigma/sigma);
134 asymm= 2.0*nwwerf(z).imag()/tau/4.0*exp(-x*x/2.0/sigma/sigma);
139 asymm= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-x*x/2.0/sigma/sigma) +
140 exp(sigma*sigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*cos(freq*x - freq/tau*sigma*sigma);
143 asymm= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-x*x/2.0/sigma/sigma) +
144 exp(sigma*sigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*sin(freq*x - freq/tau*sigma*sigma);
151 double retVal = (expG+asymm)/2.0;
154 <<
"Warning in AnalyticConvolution: negative probablity"
158 << sigma <<
' ' << tau <<
' ' << offset <<
' '
159 << freq <<
' '<< argument
162 std::cerr << retVal << std::endl;
165 else if (_type==
MIXED) {
166 double retVal = (expG-asymm)/2.0;
169 <<
"Warning in AnalyticConvolution: negative probablity"
173 << sigma <<
' ' << tau <<
' ' << offset <<
' '
174 << freq <<
' ' << argument
177 std::cerr << retVal << std::endl;
185 <<
"Unknown sign parity. State is not allowed"
193 double AnalyticConvolution::erfc(
double x)
const {
197 z = (x < 0) ? -x : x;
199 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
200 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
201 t*(-0.82215223+t*0.17087277 ))) ))) )));
202 if ( x < 0 ) ans = 2.0 - ans;
206 double AnalyticConvolution::pow(
double x,
int n)
const {
208 for(
int i=0;i<
n;i++){
214 std::complex<double> AnalyticConvolution::nwwerf(std::complex<double> z)
const {
215 std::complex<double> zh,r[38],s,t,v;
218 const double hf = z1/2;
219 const double z10 = 10;
220 const double c1 = 74/z10;
221 const double c2 = 83/z10;
222 const double c3 = z10/32;
223 const double c4 = 16/z10;
224 const double c = 1.12837916709551257;
225 const double p = pow(2.0*c4,33);
229 double xa=(x >= 0) ? x : -x;
230 double ya=(y >= 0) ? y : -y;
231 if(ya < c1 && xa < c2){
232 zh = std::complex<double>(ya+c4,xa);
233 r[37]= std::complex<double>(0,0);
235 for(
int n = 36; n>0; n--){
236 t=zh+
double(n)*std::conj(r[n+1]);
240 s=std::complex<double>(0,0);
242 for(
int k=33; k>0; k--){
249 zh=std::complex<double>(ya,xa);
250 r[1]=std::complex<double>(0,0);
252 for(
int n=9;n>0;n--){
253 t=zh+
double(n)*std::conj(r[1]);
258 if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
260 v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
261 if(x > 0) v=std::conj(v);
264 if(x < 0) v=std::conj(v);