using System; namespace Science.Mathematics.VectorCalculus { public class Integration { private Function.ToLastType del; private double from1, to1, result; public Integration() { } public Integration(Function.ToLastType f) { del = f; } public Integration(Function.ToLastType f, double from, double to) { del = f; from1 = from; to1 = to; } public void Compute() { Qromb qr = new Qromb(); result = qr.qromb(del,from1,to1,10e-15); } public Function.ToLastType Function { set{ del = value; } } public double From { set{from1 = value;} } public double To { set{to1 = value;} } public double Result { get{return result;} } } class Qromb { public Qromb() { } public double qromb(Function.ToLastType func, double a, double b, double eps) { int JMAX = 20, JMAXP = JMAX + 1, K = 5; double[] s = new double[JMAX]; double[] h = new double[JMAXP]; Poly_interp polint = new Poly_interp(h, s, K); h[0] = 1.0; Trapzd t = new Trapzd(func, a, b); for (int j = 1; j <= JMAX; j++) { s[j - 1] = t.next(); if (j >= K) { double ss = polint.rawinterp(j - K, 0.0); if (Math.Abs(polint.dy) <= eps * Math.Abs(ss)) return ss; } h[j] = 0.25 * h[j - 1]; } try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Too many steps in routine qromb"); } return 0.0; } } class Base_interp { public int n, mm, jsav, cor, dj; public double[] xx, yy; public double[] Y { set { yy = value; } } public Base_interp(double[] x, double[] y, int m) { n = x.Length; mm = m; jsav = 0; cor = 0; xx = x; yy = y; dj = Math.Min(1, (int)Math.Pow((double)n, 0.25)); } public virtual double interp(double x) { int jlo = (cor != 0) ? hunt(x) : locate(x); return rawinterp(jlo, x); } public int locate(double x) { int ju, jm, jl; if (n < 2 || mm < 2 || mm > n) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("locate size error"); } bool ascnd = (xx[n - 1] >= xx[0]); jl = 0; ju = n - 1; while (ju - jl > 1) { jm = (ju + jl) >> 1; if (x >= xx[jm] == ascnd) jl = jm; else ju = jm; } cor = Math.Abs(jl - jsav) > dj ? 0 : 1; jsav = jl; return Math.Max(0, Math.Min(n - mm, jl - ((mm - 2) >> 1))); } public int hunt(double x) { int jl = jsav, jm, ju, inc = 1; if (n < 2 || mm < 2 || mm > n) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("hunt size error"); } bool ascnd = (xx[n - 1] >= xx[0]); if (jl < 0 || jl > n - 1) { jl = 0; ju = n - 1; } else { if (x >= xx[jl] == ascnd) { for (; ; ) { ju = jl + inc; if (ju >= n - 1) { ju = n - 1; break; } else if (x < xx[ju] == ascnd) break; else { jl = ju; inc += inc; } } } else { ju = jl; for (; ; ) { jl = jl - inc; if (jl <= 0) { jl = 0; break; } else if (x >= xx[jl] == ascnd) break; else { ju = jl; inc += inc; } } } } while (ju - jl > 1) { jm = (ju + jl) >> 1; if (x >= xx[jm] == ascnd) jl = jm; else ju = jm; } cor = Math.Abs(jl - jsav) > dj ? 0 : 1; jsav = jl; return Math.Max(0, Math.Min(n - mm, jl - ((mm - 2) >> 1))); } public virtual double rawinterp(int jlo, double x) { return 0.0; } } class Poly_interp : Base_interp { public double dy; public Poly_interp(double[] xv, double[] yv, int m) : base(xv, yv, m) { dy = 0.0; } public override double rawinterp(int jl, double x) { int i, m, ns = 0; double y, den, dif, dift, ho, hp, w; double[] xa = xx; double[] ya = yy; double[] c = new double[mm]; double[] d = new double[mm]; dif = Math.Abs(x - xa[jl + 0]); for (i = 0; i < mm; i++) { if ((dift = Math.Abs(x - xa[jl + i])) < dif) { ns = i; dif = dift; } c[i] = ya[jl + i]; d[i] = ya[jl + i]; } y = ya[ns--]; for (m = 1; m < mm; m++) { for (i = 0; i < mm - m; i++) { ho = xa[jl + i] - x; hp = xa[jl + i + m] - x; w = c[i + 1] - d[i]; if ((den = ho - hp) == 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Poly_interp error"); } den = w / den; d[i] = hp * den; c[i] = ho * den; } y += (dy = (2 * (ns + 1) < (mm - m) ? c[ns + 1] : d[ns--])); } return y; } } class Quadrature { public Quadrature() { } public int n; public virtual double next() { return 0.0; } } class Trapzd : Quadrature { double a, b, s; Function.ToLastType func; public Trapzd(Function.ToLastType funcc, double aa, double bb) { func = funcc; a = aa; b = bb; n = 0; } public override double next() { double x, tnm, sum, del; int it, j; n++; if (n == 1) { return (s = 0.5 * (b - a) * (func(a) + func(b))); } else { for (it = 1, j = 1; j < n - 1; j++) it <<= 1; tnm = it; del = (b - a) / tnm; x = a + 0.5 * del; for (sum = 0.0, j = 0; j < it; j++, x += del) sum += func(x); s = 0.5 * (s + (b - a) * sum / tnm); return s; } } } }