This tutorial covers bivariate spline fitting on disc-shaped domains using FITPACK's polar coordinate types: fitpack_polar for scattered data and fitpack_grid_polar for gridded data.

Overview

Many physical problems produce data on disc-shaped or circular regions. FITPACK provides two types that handle the polar coordinate transformation internally and enforce the continuity constraints required at the origin:

Type	Input data	Boundary	Core routine
`fitpack_polar`	Scattered \( (x_i, y_i, z_i) \)	General \( r(\theta) \)	`polar`
`fitpack_grid_polar`	Grid \( z(v_j, u_i) \)	Constant radius \( r \)	`pogrid`

Always prefer fitpack_grid_polar when data lies on a polar grid – the pogrid algorithm exploits the grid structure and is significantly more efficient.

Mathematical Background

Normalized polar coordinates

Given a disc bounded by \( r(\theta) \), FITPACK maps Cartesian coordinates to normalized polar coordinates \( (u, v) \):

\[ x = u \, r(v) \cos v, \quad y = u \, r(v) \sin v, \quad 0 \leq u \leq 1, \; -\pi \leq v \leq \pi \]

The normalized radial coordinate \( u = \sqrt{x^2 + y^2} / r(\theta) \) lies in \( [0, 1] \) regardless of boundary shape, and \( v = \theta \) is the polar angle. A bicubic spline \( s(u, v) \) is fitted in the rectangular parameter domain \( [0, 1] \times [-\pi, \pi] \).

Continuity at the origin

Because all angles \( v \) map to the same physical point when \( u = 0 \), the spline must satisfy continuity constraints:

`bc_continuity_origin`	Constraint
0	\( C^0 \): \( s(0, v) = \text{const} \)
1	\( C^1 \): additionally \( \partial s / \partial u \) consistent across angles
2	\( C^2 \): second-order consistency (default for `fitpack_polar`)

See also: Dierckx, Ch. 11, §11.1 (pp. 255–263)

Scattered Data: `fitpack_polar`

fitpack_polar accepts scattered Cartesian coordinates and a user-supplied boundary function \( r(\theta) \):

use fitpack, only: fitpack_polar, fp_real, fp_flag, &
                   fitpack_success, fitpack_message, one
type(fitpack_polar) :: pol
integer(FP_FLAG) :: ierr
 
ierr = pol%new_fit(x, y, z, unit_disk, w, smoothing=real(m, fp_real))
if (.not. fitpack_success(ierr)) error stop fitpack_message(ierr)

Argument	Description
`x(:)`, `y(:)`	Cartesian coordinates of the scattered data
`z(:)`	Function values at those points
`boundary`	Procedure matching `fitpack_polar_boundary` – returns \( r(\theta) \)
`w(:)`	Optional per-point weights (default: unit weights)
`smoothing`	Optional smoothing parameter \( s \)

After the initial fit, adjust smoothing without resupplying data:

ierr = pol%fit(smoothing=10.0_fp_real)

Evaluate at Cartesian points – the \( (u, v) \) transform is internal:

z_val = pol%eval(x_val, y_val, ierr) ! single point

z_arr = pol%eval(x_arr, y_arr, ierr) ! array of points

Example: smooth function on a unit disc

integer, parameter :: m = 100
real(FP_REAL) :: x(m), y(m), z(m), w(m), r, theta
type(fitpack_polar) :: pol
integer(FP_FLAG) :: ierr
integer :: i
 
do i = 1, m
    r     = sqrt(mod(i * 0.618_fp_real, one))
    theta = 2*pi * mod(i * 0.325_fp_real, one) - pi
    x(i)  = r * cos(theta);  y(i) = r * sin(theta)
    z(i)  = (x(i)**2 + y(i)**2) / ((x(i) + y(i))**2 + half)
end do
w = one
 
ierr = pol%new_fit(x, y, z, unit_disk, w, smoothing=real(m, fp_real))
print *, pol%eval(zero, zero, ierr)
call pol%destroy()

Gridded Data: `fitpack_grid_polar`

fitpack_grid_polar accepts data on a polar grid with a constant boundary radius:

use fitpack, only: fitpack_grid_polar, fp_real, fp_flag, &
                   fitpack_success, fitpack_message, zero, one, pi
type(fitpack_grid_polar) :: gpol
integer(FP_FLAG) :: ierr
 
ierr = gpol%new_fit(u, v, r, z, z0=z0)

Argument	Description
`u(:)`	Radial grid, \( 0 < u_1 < \cdots < u_{n_u} \leq 1 \)
`v(:)`	Angular grid, \( -\pi \leq v_1 < \cdots < v_{n_v} < \pi \)
`r`	Constant boundary radius
`z(nv, nu)`	Function values (angular index first)
`z0`	Optional function value at origin \( u = 0 \)
`smoothing`	Optional smoothing parameter \( s \)

Origin boundary conditions

Use gpol\set_origin_BC to control spline behavior at the origin:

call gpol%set_origin_BC(z0=z0, exact=.true., differentiable=.true.)

ierr = gpol%fit(smoothing=real(nu * nv, fp_real))

Argument	Effect
`z0`	Prescribe the function value at the origin
`exact`	`.true.`: interpolate \( z_0 \) exactly; `.false.`: approximate
`differentiable`	`.true.`: enforce \( C^1 \) at origin; `.false.`: \( C^0 \) only

When z0 is omitted, no origin constraint is imposed.

Evaluate on \( (u, v) \) coordinates:

z_val = gpol%eval(u_val, v_val, ierr) ! single point

z_grid = gpol%eval(u_new, v_new, ierr) ! grid: shape [size(v), size(u)]

Example: data on a polar grid

integer, parameter :: nu = 10, nv = 20
real(FP_REAL) :: u(nu), v(nv), zg(nv, nu), z0
type(fitpack_grid_polar) :: gpol
integer(FP_FLAG) :: ierr
integer :: i, j
 
do i = 1, nu; u(i) = real(i, fp_real) / nu; end do
do j = 1, nv; v(j) = -pi + (j - 1) * 2*pi / nv; end do
do i = 1, nu
    do j = 1, nv
        zg(j, i) = (u(i)*cos(v(j)))**2 + (u(i)*sin(v(j)))**2
    end do
end do
z0 = zero
 
ierr = gpol%new_fit(u, v, one, zg, z0=z0)
call gpol%set_origin_BC(z0=z0, exact=.true., differentiable=.true.)
ierr = gpol%fit(smoothing=real(nu * nv, fp_real))
print *, 'Residual: ', gpol%mse()
call gpol%destroy()

Boundary Function

The boundary function for fitpack_polar must match the abstract interface fitpack_polar_boundary – a pure function receiving \( \theta \in [-\pi, \pi] \) and returning a positive radius:

pure real(FP_REAL) function my_boundary(theta) result(r)
    real(FP_REAL), intent(in) :: theta
    r = ...
end function

Unit disc (constant radius):

pure real(FP_REAL) function unit_disk(theta)
    real(FP_REAL), intent(in) :: theta
    unit_disk = one
end function

Ellipse with semi-axes \( a \) and \( b \):

pure real(FP_REAL) function ellipse(theta)
    real(FP_REAL), intent(in) :: theta
    real(FP_REAL), parameter :: a = 2.0_fp_real, b = 1.0_fp_real
    ellipse = a * b / sqrt((b * cos(theta))**2 + (a * sin(theta))**2)
end function

Complete Example

A fully worked program covering both scattered and gridded polar fitting is provided in the source tree. Build and run with:

fpm run --example example_polar

See also: examples/example_polar.f90

Error Handling

All routines return an integer(FP_FLAG) error code:

if (.not. fitpack_success(ierr)) then
    print *, 'Error: ', fitpack_message(ierr)
end if

Flag	Meaning
`FITPACK_OK` (0)	Success
`FITPACK_INTERPOLATING_OK` (-1)	Interpolating spline returned
`FITPACK_S_TOO_SMALL` (1)	Smoothing too small for the knot budget
`FITPACK_MAXIT` (2)	Maximum iterations reached
`FITPACK_INPUT_ERROR` (10)	Invalid input data