(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *)
(* Version 1, Update 2 *)

MODULE CalcConvolution;   (** AUTHOR "adf"; PURPOSE "Accurate computation of a convolution integral"; *)

(* Based on Algorithm 2 in: Diethelm, K. and Freed, A.D., "An Efficient Algorithm for the Evaluation of
	Convolution Integrals."  In review.  The level of effort required by this algorithm is O(N ln N).  *)

IMPORT NbrInt, NbrRe, DataErrors, MathInt, MathRe, CalcFn;

	PROCEDURE LogP( p, x: NbrRe.Real ): NbrRe.Real;
	BEGIN
		RETURN MathRe.Ln( x ) / MathRe.Ln( p )
	END LogP;

	(** Returns the smallest dimension (number of rows) of the solution array returned in  Solve  that is allowed. *)
	PROCEDURE GetArrayDimension*( x, t: NbrRe.Real;  n: NbrInt.Integer ): NbrInt.Integer;
	BEGIN
		RETURN NbrRe.Ceiling( n * x / t ) + 1
	END GetArrayDimension;

	(** Solves the convolution integral  y =  x0x k(x-u) f(u,x) du.
	The input parameters are:
		k	is a univariate kernel function, which must be positive and decay monotonically,
			and  k(0) = %  is allowed,
		f	is a bivariate forcing function,
		x	is the upper limit of integration, the lower limit is set to zero,  x > 0,
		t	is the characteristic time affiliated with kernel  k,  t > 0,
		n	is the number of integration steps applied to each interval of characteristic length  t, n > 3,
		p	is a parameter that balances speed (smaller values) against accuracy (larger values)
			with typical values ranging between 4 and 6,  p > 1.
	The return parameter is:
		soln	is an array of dimension  hn*x/ti + 1 4 2  or greater containing pairs  {xi , yi}  where:
			xi = i*x/J,  yi =  x0xi k(xi- u) f(u,xi) du,  i = 0,1, 2, ..., J,  J = hn*x/ti.
	*)
	PROCEDURE Solve*( k: CalcFn.ReArg;  f: CalcFn.Re2Arg;  x, t: NbrRe.Real;  n, p: NbrInt.Integer;
						   VAR soln: ARRAY OF ARRAY OF NbrRe.Real );
	VAR i, j, l, s, z, gamma, iBar, jj, mu, power, sigma: NbrInt.Integer;  h, prod, sum, xi: NbrRe.Real;
		kernel: POINTER TO ARRAY OF ARRAY OF NbrRe.Real;
		laplaceArray: ARRAY 5 OF NbrRe.Real;
		laplaceMatrix: ARRAY 8, 8 OF NbrRe.Real;
	BEGIN
		IF k # NIL THEN
			IF f # NIL THEN
				IF x > 0 THEN
					IF t > 0 THEN
						IF LEN( soln, 0 ) > NbrRe.Ceiling( n * x / t ) THEN
							IF LEN( soln, 1 ) > 1 THEN
								IF n < 4 THEN n := 4 END;
								IF p < 2 THEN p := 2 END;
								(* Initialize variables. *)
								h := t / n;  mu := NbrInt.Max( 1, NbrRe.Ceiling( LogP( p, x / t ) ) );
								(* Create the table of required kernel values. *)
								NEW( kernel, mu + 1, p * n + 1 );
								FOR i := 1 TO p * n DO kernel[1, i] := k( h * (i - 0.5) ) END;
								FOR j := 2 TO mu DO
									FOR i := 1 TO (p - 1) * n DO kernel[j, i] := k( (t + h * (i - 0.5)) * MathInt.Power( p, j - 1 ) ) END
								END;
								(* Determine the weights of Laplace quadrature when there are 4, 5, 6 or 7 nodes. *)
								FOR j := 4 TO 7 DO
									FOR i := 1 TO j DO laplaceMatrix[i, j] := 1 END;
									laplaceMatrix[1, j] := laplaceMatrix[1, j] + 703 / 5760;
									laplaceMatrix[2, j] := laplaceMatrix[2, j] - 463 / 1920;
									laplaceMatrix[3, j] := laplaceMatrix[3, j] + 101 / 640;
									laplaceMatrix[4, j] := laplaceMatrix[4, j] - 223 / 5760;
									laplaceMatrix[j, j] := laplaceMatrix[j, j] + 703 / 5760;
									laplaceMatrix[j - 1, j] := laplaceMatrix[j - 1, j] - 463 / 1920;
									laplaceMatrix[j - 2, j] := laplaceMatrix[j - 2, j] + 101 / 640;
									laplaceMatrix[j - 3, j] := laplaceMatrix[j - 3, j] - 223 / 5760;
								END;
								(* Determine the weight of Laplace quadrature when there are 8 or more nodes. *)
								laplaceArray[1] := 1 + 703 / 5760;  laplaceArray[2] := 1 - 463 / 1920;
								laplaceArray[3] := 1 + 101 / 640;  laplaceArray[4] := 1 - 223 / 5760;
								(* The main part of the algorithm. *)
								jj := NbrRe.Ceiling( x / h );  soln[0, 0] := 0;  soln[0, 1] := 0;
								(* Solve at the grid points in sub-interval  [0,  p*t]. *)
								IF mu = 1 THEN iBar := jj ELSE iBar := p * n END;
								FOR i := 1 TO NbrInt.Min( 3, iBar ) DO
								(* MacLaurin approximation for  x0xi k(xi- u) f(u,xi) du. *)
									xi := i * h;  soln[i, 0] := xi;
									soln[i, 1] :=
										i * h *
										(13 * (k( 7 * xi / 8 ) * f( xi / 8, xi ) + k( xi / 8 ) * f( 7 * xi / 8, xi )) +
										  11 * (k( 5 * xi / 8 ) * f( 3 * xi / 8, xi ) + k( 3 * xi / 8 ) * f( 5 * xi / 8, xi ))) / 48
								END;
								FOR i := 4 TO NbrInt.Min( 7, iBar ) DO
								(* Laplace approximation for  x0xi k(xi- u) f(u,xi) du. *)
									xi := i * h;  soln[i, 0] := xi;  soln[i, 1] := 0;
									FOR j := 1 TO i DO
										prod := laplaceMatrix[j, i] * kernel[1, j] * f( (i - j + 0.5) * h, xi );  soln[i, 1] := soln[i, 1] + prod
									END;
									soln[i, 1] := h * soln[i, 1]
								END;
								FOR i := 8 TO iBar DO
								(* Continue with Laplace approximation for  x0xi k(xi- u) f(u,xi) du. *)
									xi := i * h;  soln[i, 0] := xi;  soln[i, 1] := 0;
									FOR j := 1 TO 4 DO
										prod := laplaceArray[j] * kernel[1, j] * f( (i - j + 0.5) * h, xi );  soln[i, 1] := soln[i, 1] + prod
									END;
									FOR j := 5 TO i - 4 DO prod := kernel[1, j] * f( (i - j + 0.5) * h, xi );  soln[i, 1] := soln[i, 1] + prod END;
									FOR j := i - 3 TO i DO
										prod := laplaceArray[i + 1 - j] * kernel[1, j] * f( (i - j + 0.5) * h, xi );
										soln[i, 1] := soln[i, 1] + prod
									END;
									soln[i, 1] := h * soln[i, 1]
								END;
								(* Solve at the remaining grid points. *)
								FOR i := iBar + 1 TO jj DO
								(* Determine the contribution from  x0pst k(xi- u) f(u,xi) du. *)
									xi := i * h;  soln[i, 0] := xi;  soln[i, 1] := 0;
									FOR j := 1 TO 4 DO
										prod := laplaceArray[j] * kernel[1, j] * f( (i - j + 0.5) * h, xi );  soln[i, 1] := soln[i, 1] + prod
									END;
									FOR j := 5 TO p * n - 4 DO
										prod := kernel[1, j] * f( (i - j + 0.5) * h, xi );  soln[i, 1] := soln[i, 1] + prod
									END;
									FOR j := p * n - 3 TO p * n DO
										prod := laplaceArray[p * n + 1 - j] * kernel[1, j] * f( (i - j + 0.5) * h, xi );
										soln[i, 1] := soln[i, 1] + prod
									END;
									sigma := NbrRe.Ceiling( LogP( p, i / n ) ) - 1;
									FOR s := 1 TO sigma - 1 DO
										sum := 0;  power := MathInt.Power( p, s );
										FOR j := 1 TO 4 DO
											sum := sum + laplaceArray[j] * kernel[s + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
										END;
										FOR j := 5 TO n * (p - 1) - 4 DO
											sum := sum + kernel[s + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
										END;
										FOR j := n * (p - 1) - 3 TO n * (p - 1) DO
											sum :=
												sum +
												laplaceArray[n * (p - 1) + 1 - j] * kernel[s + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
										END;
										soln[i, 1] := soln[i, 1] + power * sum
									END;
									soln[i, 1] := h * soln[i, 1];
									(* Add to this the contribution from  xpstXi k(xi- u) f(u,xi) du. *)
									power := MathInt.Power( p, sigma );  l := i - n * power;
									IF l < 4 * power THEN
										(* Use a MacLaurin approximation. *)
										soln[i, 1] :=
											soln[i, 1] +
											l * h *
						  (13 * (k( (i - l / 8) * h ) * f( l * h / 8, xi ) + k( (i - 7 * l / 8) * h ) * f( 7 * l * h / 8, xi )) +
						    11 * (k( (i - 3 * l / 8) * h ) * f( 3 * l * h / 8, xi ) + k( (i - 5 * l / 8) * h ) * f( 5 * l * h / 8, xi ))) / 48
									ELSE
										(* Use a Laplace approximation for  xpstXps(n+g) k(xi- u) f(u,xi) du. *)
										sum := 0;  gamma := NbrRe.Floor( l / power );
										IF gamma < 8 THEN
											FOR j := 1 TO gamma DO
												sum :=
													sum +
													laplaceMatrix[j, gamma] * kernel[sigma + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
											END;
											soln[i, 1] := soln[i, 1] + h * power * sum
										ELSE
											FOR j := 1 TO 4 DO
												sum :=
													sum + laplaceArray[j] * kernel[sigma + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
											END;
											FOR j := 5 TO gamma - 4 DO
												sum := sum + kernel[sigma + 1, j] * f( xi - (t + h * (j - 0.5)) * power, xi )
											END;
											FOR j := gamma - 3 TO gamma DO
												sum :=
													sum +
													laplaceArray[gamma + 1 - j] * kernel[sigma + 1, j] *
										  f( xi - (t + h * (j - 0.5)) * power, xi )
											END;
											soln[i, 1] := soln[i, 1] + h * power * sum
										END;
										(* Use a MacLaurin approximation for  xXps(n+g)Xi k(xi- u) f(u,xi) du. *)
										z := i - power * (n + gamma);
										soln[i, 1] :=
											soln[i, 1] +
											h * z *
						  (13 * (k( (i - z / 8) * h ) * f( z * h / 8, xi ) + k( (i - 7 * z / 8) * h ) * f( 7 * z * h / 8, xi )) +
						    11 * (k( (i - 3 * z / 8) * h ) * f( 3 * z * h / 8, xi ) + k( (i - 5 * z / 8) * h ) * f( 5 * z * h / 8, xi ))) / 48
									END
								END
							ELSE DataErrors.Error( "The second dimension of the solution array must be at least 2." )
							END
						ELSE DataErrors.IntError( LEN( soln, 0 ), "The solution array is too short." )
						END
					ELSE DataErrors.ReError( t, "The characteristic time is required to be positive." )
					END
				ELSE DataErrors.ReError( x, "The upper limit of integration is required to be positive." )
				END
			ELSE DataErrors.Error( "The supplied forcing function was NIL." )
			END
		ELSE DataErrors.Error( "The supplied kernel function was NIL." )
		END
	END Solve;

END CalcConvolution.