Skip to content

Commit

Permalink
Checkpointing work. Refs #7653.
Browse files Browse the repository at this point in the history
  • Loading branch information
wdzhou committed Aug 14, 2013
1 parent 180ecf2 commit 1342c46
Show file tree
Hide file tree
Showing 2 changed files with 2 additions and 73 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -124,9 +124,6 @@ namespace CurveFitting

};

/// Integral for Gamma
std::complex<double> E1(std::complex<double> z);


} // namespace CurveFitting
} // namespace Mantid
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -533,8 +533,8 @@ namespace CurveFitting
const double SQRT_H_5 = sqrt(H)*.5;
std::complex<double> p(alpha*x, alpha*SQRT_H_5);
std::complex<double> q(-beta*x, beta*SQRT_H_5);
double omega2a = imag(exp(p)*E1(p));
double omega2b = imag(exp(q)*E1(q));
double omega2a = imag(exp(p)*Mantid::API::E1(p));
double omega2b = imag(exp(q)*Mantid::API::E1(q));
omega2 = -1.0*N*eta*(omega2a + omega2b)*TWO_OVER_PI;
}
const double omega = omega1+omega2;
Expand All @@ -557,73 +557,5 @@ namespace CurveFitting
return omega;
}

//------------------------- External Functions ---------------------------------------------------
// FIXME : This is the same function used for ThermalNeutron... ...
/** Implementation of complex integral E_1
*/
std::complex<double> E1(std::complex<double> z)
{
const double el = 0.5772156649015328;

std::complex<double> exp_e1;

double rz = real(z);
double az = abs(z);

if (fabs(az) < 1.0E-8)
{
// If z = 0, then the result is infinity... diverge!
complex<double> r(1.0E300, 0.0);
exp_e1 = r;
}
else if (az <= 10.0 || (rz < 0.0 && az < 20.0))
{
// Some interesting region, equal to integrate to infinity, converged
// cout << "[DB] Type 1" << endl;

complex<double> r(1.0, 0.0);
exp_e1 = r;
complex<double> cr = r;

for (size_t k = 1; k <= 150; ++k)
{
double dk = double(k);
cr = -cr * dk * z / ( (dk+1.0)*(dk+1.0) );
exp_e1 += cr;
if (abs(cr) < abs(exp_e1)*1.0E-15)
{
// cr is converged to zero
break;
}
} // ENDFOR k

// cout << "[DB] el = " << el << ", exp_e1 = " << exp_e1 << endl;

exp_e1 = -el - log(z) + (z*exp_e1);
}
else
{
// Rest of the region
complex<double> ct0(0.0, 0.0);
for (int k = 120; k > 0; --k)
{
complex<double> dk(double(k), 0.0);
ct0 = dk / (10.0 + dk / (z + ct0));
} // ENDFOR k

exp_e1 = 1.0 / (z + ct0);
exp_e1 = exp_e1 * exp(-z);
if (rz < 0.0 && fabs(imag(z)) < 1.0E-10 )
{
complex<double> u(0.0, 1.0);
exp_e1 = exp_e1 - (PI * u);
}
}

// cout << "[DB] Final exp_e1 = " << exp_e1 << "\n";

return exp_e1;
}

} // namespace CurveFitting
} // namespace Mantid

0 comments on commit 1342c46

Please sign in to comment.