Skip to content

Commit

Permalink
Refs #11515 Fix up indentation in DynamicKuboToyabe
Browse files Browse the repository at this point in the history
  • Loading branch information
Harry Jeffery committed Apr 13, 2015
1 parent 271bcfa commit cf3a069
Showing 1 changed file with 66 additions and 64 deletions.
130 changes: 66 additions & 64 deletions Code/Mantid/Framework/CurveFitting/src/DynamicKuboToyabe.cpp
Expand Up @@ -30,92 +30,94 @@ void DynamicKuboToyabe::init()

// Midpoint method
double midpnt(double func(const double, const double, const double),
const double a, const double b, const int n, const double g, const double w0) {
// quote & modified from numerical recipe 2nd edtion (page147)

const double a, const double b, const int n, const double g,
const double w0) {
// quote & modified from numerical recipe 2nd edtion (page147)

static double s;

if (n==1) {
s = (b-a)*func(0.5*(a+b),g,w0);
if (n == 1) {
s = (b - a) * func(0.5 * (a + b), g, w0);
return (s);
} else {
double x, tnm, sum, del, ddel;
int it, j;
for (it=1,j=1;j<n-1;j++) it *= 3;
tnm = it;
del = (b-a)/(3*tnm);
ddel=del+del;
x = a+0.5*del;
sum =0.0;
for (j=1;j<=it;j++) {
sum += func(x,g,w0);
x += ddel;
sum += func(x,g,w0);
x += del;
}
s=(s+(b-a)*sum/tnm)/3.0;
return s;
}
} else {
double x, tnm, sum, del, ddel;
int it, j;
for (it = 1, j = 1; j < n - 1; j++)
it *= 3;
tnm = it;
del = (b - a) / (3 * tnm);
ddel = del + del;
x = a + 0.5 * del;
sum = 0.0;
for (j = 1; j <= it; j++) {
sum += func(x, g, w0);
x += ddel;
sum += func(x, g, w0);
x += del;
}
s = (s + (b - a) * sum / tnm) / 3.0;
return s;
}
}

// Polynomial interpolation
void polint (double xa[], double ya[], int n, double x, double& y, double& dy) {
int i, m, ns = 1;
int i, m, ns = 1;
double dif;

dif = fabs(x-xa[1]);
std::vector<double> c(n+1);
std::vector<double> d(n+1);
for (i=1;i<=n;i++){
dif = fabs(x - xa[1]);
std::vector<double> c(n + 1);
std::vector<double> d(n + 1);
for (i = 1; i <= n; i++) {
double dift;
if((dift=fabs(x-xa[i]))<dif) {
ns=i;
dif=dift;
}
c[i]=ya[i];
d[i]=ya[i];
}
y=ya[ns--];
for (m=1;m<n;m++) {
for (i=1;i<=n-m;i++) {
if ((dift = fabs(x - xa[i])) < dif) {
ns = i;
dif = dift;
}
c[i] = ya[i];
d[i] = ya[i];
}
y = ya[ns--];
for (m = 1; m < n; m++) {
for (i = 1; i <= n - m; i++) {
double den, ho, hp, w;
ho=xa[i]-x;
hp=xa[i+m]-x;
w=c[i+1]-d[i];
if((den=ho-hp)==0.0){ //error message!!!
ho = xa[i] - x;
hp = xa[i + m] - x;
w = c[i + 1] - d[i];
if ((den = ho - hp) == 0.0) { // error message!!!
throw std::runtime_error("Error in routin polint");
}
den=w/den;
d[i]=hp*den;
c[i]=ho*den;
}
y += (dy=(2*(ns)<(n-m) ? c[ns+1] : d[ns--]));

}

den = w / den;
d[i] = hp * den;
c[i] = ho * den;
}
y += (dy = (2 * (ns) < (n - m) ? c[ns + 1] : d[ns--]));
}
}

// Integration
double integral (double func(const double, const double, const double),
const double a, const double b, const double g, const double w0) {
double integral(double func(const double, const double, const double),
const double a, const double b, const double g,
const double w0) {

const int JMAX = 14;
const int JMAXP = JMAX + 1;
const int K = 5;

int j;
double ss,dss;
double h[JMAXP+1], s[JMAXP];

h[1] = 1.0;
for (j=1; j<= JMAX; j++) {
s[j]=midpnt(func,a,b,j,g,w0);
if (j >= K) {
polint(&h[j-K],&s[j-K],K,0.0,ss,dss);
if (fabs(dss) <= fabs(ss)) return ss;
}
h[j+1]=h[j]/9.0;
}
double ss, dss;
double h[JMAXP + 1], s[JMAXP];

h[1] = 1.0;
for (j = 1; j <= JMAX; j++) {
s[j] = midpnt(func, a, b, j, g, w0);
if (j >= K) {
polint(&h[j - K], &s[j - K], K, 0.0, ss, dss);
if (fabs(dss) <= fabs(ss))
return ss;
}
h[j + 1] = h[j] / 9.0;
}
throw std::runtime_error("Too many steps in routine integrate");
}

Expand Down

0 comments on commit cf3a069

Please sign in to comment.