=head1 NAME

PDL::Interpolate::Slatec - simple interface to SLATEC interpolation routines

=head1 SYNOPSIS

 use PDL::Interpolate::Slatec;
 use PDL::Math;

 # somewhat pointless way to estimate cos and sin,
 # but is shows that you can thread if you want to
 #
 my $x   = pdl( 0 .. 45 ) * 4 * 3.14159 / 180;
 my $y   = cat( sin($x), cos($x) );
 #
 my $obj = new PDL::Interpolate::Slatec( x => $x, y = $y );
 #
 my $xi  = pdl( 0.5, 1.5, 2.5 );
 my $yi  = $obj->interpolate( $xi );
 #
 print "cos( $xi ) equals ", $yi->slice(':,(0)'), "\n";
 cos( [0.5 1.5 2.5] ) equals  [0.87759844 0.070737667 -0.80115622]
 #
 print "sin( $xi ) equals ", $yi->slice(':,(1)'), "\n";
 sin( [0.5 1.5 2.5] ) equals  [ 0.4794191 0.99768655 0.59846449]
 #
 print cos($xi), "\n", sin($xi), "\n";
 [0.87758256 0.070737202 -0.80114362]
 [0.47942554 0.99749499 0.59847214]

=head1 DESCRIPTION

Use the interface defined by L<PDL::Interpolate|PDL::Interpolate>
to provide a simple way to use the SLATEC interpolation
routines (e.g. see L<PDL::Slatec|PDL::Slatec>).
Hence the name for this library - as returned by the C<library>
method - is C<"Slatec">.

Currently, only the 
L<piecewise cubic Hermite polynomial routines|PDL::Slatec/Piecewise_cubic_Hermite_interpol> 
are available (C<type == "pch">).

=head2 Attributes

The following changes are made to the attributes 
of L<PDL::Interpolate|PDL::Interpolate>:

 Attribute  Flag  Description
 bc         sgr   boundary conditions
 g          g     estimated gradient at x positions

 Attribute  Default    Allowed values
 bc         "simple"   see Boundary conditions section
 type       "pch"

Given the initial set of points C<(x,y)>, the C<"pch"> 
library estimates the gradient at these points using the
given boundary conditions (as specified by the
C<bc> attribute).
The estimated gradient can be obtained using

 $gradient = $obj->get( 'g' );

As described in the L<interpolate|/interpolate> method,
the C<"pch"> routines can also estimate the gradient,
as well as the function value, for a set of C<$xi>.

=head2 Boundary conditions for the pch routines

If your data is monotonic, and you are not too bothered about
edge effects, then the default value of C<bc> of C<"simple"> is for you.
Otherwise, take a look at the description of
L<PDL::Slatec::chic|PDL::Slatec/chic> and use a hash reference
for the C<bc> attribute, with the following keys:

=over 3

=item monotonic

0 if the interpolant is to be monotonic in each interval (so
the gradient will be 0 at each switch point), 
otherwise the gradient is calculated using a 3-point difference
formula at switch points. 
If E<gt> 0 then the interpolant is forced to lie close to the 
data, if E<lt> 0 no such control is imposed.
Default = B<0>.

=item start

A perl list of one or two elements. The first element defines how the
boundary condition for the start of the array is to be calculated;
it has a range of C<-5 .. 5>, as given for the C<ic> parameter
of L<chic|PDL::Slatec/chic>. 
The second element, only used if options 2, 1, -1, or 2
are chosen, contains the value of the C<vc> parameter.
Default = B<[ 0 ]>.

=item end

As for C<start>, but for the end of the data.

=back

An example would be

 $obj->set( bc => { start => [ 1, 0 ], end => [ 1, -1 ] }

which sets the first derivative at the first point to 0, 
and at the last point to -1.

=head2 Errors

The C<status> method provides a simple mechanism to check if
the previous method was successful. The C<err> attribute
contains the C<$ierr> piddle returned by the Slatec
routine if a more precise diagnostic is required.
To find out which routine was called, use the
C<routine> method.

=cut

package PDL::Interpolate::Slatec;

use PDL::Interpolate;
use PDL::Slatec;
use Carp;

use strict;
use vars qw( @ISA );

@ISA     = qw ( PDL::Interpolate );

####################################################################
#

####################################################################
# 
## Public routines:

sub new {
    my $this  = shift;
    my $class = ref($this) || $this;
    my $self  = $class->SUPER::new();  # note: do not send in values

    # change from PDL::Interpolate to PDL::Interpolate::Slatec
    bless ($self, $class);

    # change class attributes
    $self->_change_attr( 
			 bc => { required => 1, settable => 1 }, # already gettable
			 );
    $self->_set_value( bc => "simple", type => "pch" );
    $self->_add_attr( 
		      g => { gettable => 1 },
		     );

    $self->{flags}->{library} = "Slatec";
    $self->{flags}->{routine} = "none";

    # set variables
    $self->set( @_ );

    return $self;         
                                               
} # sub: new

####################################################################

# set up the interpolation
#
sub _initialise {
    my $self = shift;

    # set up error flags
    $self->{flags}->{status} = 0;
    $self->{flags}->{routine} = "none";

    # get values in one go
    my ( $x, $y, $g, $bc ) = $self->_get_value( qw( x y g bc ) );

    # check 1st dimention of x and y are the same
    #  ie allow the possibility of threading
    my $xdim = $x->getdim( 0 );
    my $ydim = $y->getdim( 0 );
    croak "ERROR: x and y piddles must have the same first dimension.\n"
	unless $xdim == $ydim;

    # if a gradient has been specified, then we don't need to do anything
    # - other than check the dimensions
    if ( defined $g ) {
	croak "ERROR: gradient piddle must have the same first dimension as x and y.\n"
	    unless $g->getdim( 0 ) == $xdim;
	$self->{flags}->{status} = 1;
	return;
    }

    my $ierr;
    if ( ref($bc) eq "HASH" ) {
	my $monotonic = $bc->{monotonic} || 0;
	my $start     = $bc->{start}     || [ 0 ];
	my $end       = $bc->{end}       || [ 0 ];

	my $ic = $x->short( $start->[0], $end->[0] );
	my $vc = $x->float( 0, 0 ); # it will get promoted if required

	if ( $#$start == 1 ) { $vc->set( 0, $start->[1] ); }
	if ( $#$end   == 1 ) { $vc->set( 1, $end->[1] ); }

	my $wk = $x->zeroes( $x->float, 2*$xdim );
	( $g, $ierr ) = chic( $ic, $vc, $monotonic, $x, $y, $wk );

	$self->{flags}->{routine} = "chic";

    } elsif ( $bc eq "simple" ) {
	# chim
	( $g, $ierr ) = chim( $x, $y );

	$self->{flags}->{routine} = "chim";

    } else {
	# Unknown boundary condition
	croak "ERROR: unknown boundary condition <$bc>.\n";
	# return;
    }

    $self->_set_value( g => $g, err => $ierr );

    if ( all $ierr == 0 ) {
	# everything okay
	$self->{flags}->{status} = 1;
    } elsif ( any $ierr < 0 ) {
	# a problem
	$self->{flags}->{status} = 0;
    } else {
	# there were switches in monotonicity
	$self->{flags}->{status} = -1;
    }

} # sub: _initialise

####################################################################

=head2 interpolate

=for usage

 my $yi          = $obj->interpolate( $xi );
 my ( $yi, $gi ) = $obj->interpolate( $xi );

=for ref

Returns the interpolated function and derivative
at a given set of points.

If evaluated in scalar mode, it returns only 
the interpolated function values.

=cut

sub interpolate {
    my $self = shift;
    my $xi   = shift;

    croak 'Usage: $obj->interpolate( $xi )' . "\n"
	unless defined $xi;

    # check everything is fine
    $self->_check_attr();

    # get values in one go
    my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );

    my ( $yi, $gi, $ierr );

    if ( wantarray ) {
	( $yi, $gi, $ierr ) = chfd( $x, $y, $g, 0, $xi );
	$self->{flags}->{routine} = "chfd";
    } else {
	( $yi, $ierr ) = chfe( $x, $y, $g, 0, $xi );
	$self->{flags}->{routine} = "chfe";
    }

    # set err/status info
    $self->_set_value( err => $ierr );

    if ( all $ierr == 0 ) {
	# everything okay
	$self->{flags}->{status} = 1;
    } elsif ( all $ierr > 0 ) {
	# extrapolation was required
	$self->{flags}->{status} = -1;
    } else {
	# a problem
	$self->{flags}->{status} = 0;
    }
	
    return wantarray ? ( $yi, $gi ) : $yi;

} # sub: interpolate

=head2 integrate

=for usage

 my $ans = $obj->integrate( index => pdl( 2, 5 ) );
 my $ans = $obj->integrate( x => pdl( 2.3, 4.5 ) );

=for ref

Integrate the function stored in the PDL::Interpolate::Slatec
object. 

The integration can either be between points of
the original C<x> array (C<index>), or arbitrary x values
(C<x>). For both cases, a two element piddle
should be given,
to specify the start and end points of the integration.

=over 7

=item index

The values given refer to the indices of the points
in the C<x> array.

=item x

The array contains the actual values to integrate between.

=back

If the C<status> method returns a value of -1, then
one or both of the integration limits did not
lie inside the C<x> array. I<Caveat emptor> with the
result in such a case.

The reason for using piddles, rather than arrays, is that
it allows for threading.

=cut

sub integrate {
    my $self = shift;

    croak 'Usage: $obj->integrate( $type => $limits )' . "\n"
	unless $#_ == 1;

    # check everything is fine
    $self->_check_attr();

    $self->{flags}->{status} = 0;
    $self->{flags}->{routine} = "none";

    my ( $type, $indices ) = ( @_ );

    croak "Unknown type ($type) sent to integrate method.\n"
	unless $type eq "x" or $type eq "index";

    my $fdim = $indices->getdim(0);
    croak "Indices must have a first dimension of 2, not $fdim.\n"
	unless $fdim == 2;

    my $lo = $indices->slice('(0)');
    my $hi = $indices->slice('(1)');

    my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );
    my ( $ans, $ierr );

    if ( $type eq "x" ) {
	( $ans, $ierr ) = chia( $x, $y, $g, 0, $lo, $hi );
	$self->{flags}->{routine} = "chia";

	if ( all $ierr == 0 ) {
	    # everything okay
	    $self->{flags}->{status} = 1;
	} elsif ( any $ierr < 0 ) {
	    # a problem
	    $self->{flags}->{status} = 0;
	} else {
	    # out of range
	    $self->{flags}->{status} = -1;
	}

    } else {
	( $ans, $ierr ) = chid( $x, $y, $g, 0, $lo, $hi );
	$self->{flags}->{routine} = "chid";

	if ( all $ierr == 0 ) {
	    # everything okay
	    $self->{flags}->{status} = 1;
	} elsif ( all $ierr != -4 ) {
	    # a problem
	    $self->{flags}->{status} = 0;
	} else {
	    # out of range (ierr == -4)
	    $self->{flags}->{status} = -1;
	}

    }

    $self->_set_value( err => $ierr );
    return $ans;
}

=head1 AUTHOR

Copyright (C) 2000 Doug Burke (burke@ifa.hawaii.edu).
All rights reserved. There is no warranty. 
You are allowed to redistribute this software / documentation as 
described in the file COPYING in the PDL distribution.                                    

=head1 SEE ALSO

L<PDL::Interpolate>, L<PDL>, perltoot(1).

=cut

####################################################################
# End with a true
1;