function [A,b,x]=geomag(n,example,a,b,h)

A=[];x=[];
[nargout,nargin] = argn(0)
//GEOMAG Test problem: 1-D geomagnetic prospecting model problem
// 
// [A,b,x] = geomag(n,example,a,b,h)
// 
// Discretization of a 1-D model problem in geomagnetic prospecting,
// in which a deposit with magnetid field f(t) is located at depth h,
// while the magnetic field g(s) is measured at the surface.
// 
// The resulting problem is a first-kind Fredholm integral equation
// with kernel
//    K(s,t) = h*(h^2 + (s-t)^2)^(-3/2) .
// The following three examples are implemented (example = 1 is default):
//    1: f(t) = sin(pi*t) + 0.5*sin(2*pi*t),
//    2: f(t) = piecewise linear function,
//    3: f(t) = piecewise constant function.
// The problem is discretized by means of the midpoint quadrature
// rule with n points, leading to the matrix A and the vector x.
// Then the right-hand side is computed as b = A*x.
// 
// The t integration interval is fixed to [0,1], while the s
// integration interval [a,b] can be specified by the user.
// The default interval is [0,1], leading to a symmetric positive
// definite Toeplitz matrix.
// 
// The parameter h is the depth at which the magnetic deposit
// is located, and the default value is h = 0.25. The larger the h,
// the faster the decay of the singular values.
 
// Reference: G. M. Wing and J. D. Zahrt, """"A Primer on Integral
// Equations of the First Kind"""", SIAM, Philadelphia, 1991; p. 17.
 
// Per Christian Hansen, IMM, 07/27/98.
 
// Initialization.
if nargin<2 then
  example = 1;
end
if nargin<4 then
  a = 0;
  b = 1;
end
if nargin<5 then
  h = 0.25;
end
 
// Set up abscissas and matrix.
dt = 1/n;
ds = (b-a)/n;
t = dt*(((1:n)')-0.5);
s = a+ds*(((1:n)')-0.5);
// inline translation of meshgrid(t,s)
T=matrix(t,1,n)
S=matrix(s,n,1)
T=T(ones(T),:)
S=S(:,ones(S)')
A = dt*h*ones(n,n) ./ ((h^2+(S-T).^2).^(3/2));
 
// Set up solution vector and right-hand side.
nt = round(n/3);
nn = round(n*7/8);
x = ones(n,1);
select example
case 1 then
  x = sin(%pi*t)+0.5*sin(2*%pi*t);
case 2 then
  %v1 = 2/nt*((1:nt)')
  x(1:nt,1) = %v1(:);
  %v1 = (2*nn-nt-((nt+1:nn)'))/(nn-nt)
  x(nt+1:nn,1) = %v1(:);
  %v1 = (n-((nn+1:n)'))/(n-nn)
  x(nn+1:n,1) = %v1(:);
case 3 then
  x(1:nt) = 2*ones(nt,1)
else
  error('Illegal value of example');
end
b = A*x;