# -*-perl-*-
#
# The contents depend on whether we have bad-value support in PDL.
# (there's only a small amount of code that uses bad values)
#
# - this is called by Makefile.PL to ensure the file exists
#   when the makefile is written
#

use strict;

use Config;
use File::Basename qw(&basename &dirname);

# check for bad value support
use PDL::Config;

my $bvalflag = $PDL::Config{WITH_BADVAL};

# This forces PL files to create target in same directory as PL file.
# This is so that make depend always knows where to find PL derivatives.
chdir(dirname($0));
my $file;
($file = basename($0)) =~ s/\.PL$//;
$file =~ s/\.pl$//
        if ($Config{'osname'} eq 'VMS' or
            $Config{'osname'} eq 'OS2');  # "case-forgiving"

if ( $bvalflag ) {
    print "Extracting $file (WITH bad value support)\n";
} else {
    print "Extracting $file (NO bad value support)\n";
}

open OUT,">$file" or die "Can't create $file: $!";
chmod 0644, $file;

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

print OUT <<"!WITH!SUBS!";

# This file is automatically created by NDF.pm.PL
# - bad value support = $bvalflag
!WITH!SUBS!

print OUT <<'!NO!SUBS!';

package PDL::IO::NDF;

=head1 NAME

PDL::IO::NDF - PDL Module for reading and writing Starlink
N-dimensional data structures as PDLs.

=head1 SYNOPSIS

 use PDL::IO::NDF;

 $a = PDL->rndf($file);

 $a = rndf('test_image');
 $a = rndf('test_image', 1);

 $a->wndf($file);
 wndf($a, 'out_image');

 propndfx($a, 'template', 'out_image');

=head1 DESCRIPTION

This module adds the ability to read and write Starlink N-dimensional data
files as N-dimensional PDLs.

You must have the Starlink NDF library installed to use it.  The library
is distributed under the GPL and is available from "http://www.starlink.ac.uk".

=head1 FUNCTIONS

=cut

@EXPORT_OK = qw/rndf wndf/;
%EXPORT_TAGS = (Func=>[@EXPORT_OK]);

@ISA    = qw( PDL::Exporter );

use PDL::Core;
use PDL::Types;
use PDL::Options;
use Carp;
use strict;

use Astro::FITS::Header::Item;

# Starlink data type conversion
use vars qw/%pdltypes %startypes $ndf_loaded $VERSION $EXTNAME/;

$VERSION = '1.06';

# Set PDL -> Starlink data types
%pdltypes = ("$PDL_B"  => "_UBYTE",     # was "_BYTE"
	     "$PDL_S"  => "_WORD",
	     "$PDL_US" => "_UWORD",
	     "$PDL_L"  => "_INTEGER",
	     "$PDL_F"  => "_REAL",
	     "$PDL_D"  => "_DOUBLE"
	    );

# Set Starlink -> PDL data types
%startypes = ("_BYTE"    =>$PDL_B,
	      "_UBYTE"   =>$PDL_B,
	      "_WORD"    =>$PDL_S,
	      "_UWORD"   =>$PDL_US,
	      "_INTEGER" =>$PDL_L,
	      "_REAL"    =>$PDL_F,
	      "_DOUBLE"  =>$PDL_D
	     );

# This is the name I use in the PDL header hash that contains the
# NDF extensions so that they can be recreated

$EXTNAME = 'NDF_EXT';

# load NDF library
sub _init_NDF {
    return if $ndf_loaded++; # a bit pointless increasing each time?
    # Sometimes $@ can give the wrong answer for test. Be explicit.
    my $retval = eval 'use NDF; 1;';
    croak "Cannot use NDF library: $@" if $retval != 1;
}

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';

use PDL::Bad;

# internal functions for checking bad values

my $starbad = 0;  # flag

use vars qw( %starbad );

sub _init_starbad {
    return if $starbad;
    _init_NDF();  # just to be safe

    %starbad = (
        $PDL_S => &NDF::VAL__BADW(),
        $PDL_US => &NDF::VAL__BADUW(),
        $PDL_L => &NDF::VAL__BADI(),
        $PDL_F => &NDF::VAL__BADR(),
        $PDL_D => &NDF::VAL__BADD(),
        );

    $starbad = 1;
}

sub starlink_bad_value ($) {
    _init_starbad();
    return $starbad{$_[0]} if (exists $starbad{$_[0]});
    croak "Unable to generate a bad Starlink value for PDL data type $_[0]\n";
}

# given a piddle, check if the Starlink bad value matches
# the bad value used for this type
sub compare_bad_values ($) {
    _init_starbad();
    return ($starbad{$_[0]} == badvalue($_[0]));
}

!NO!SUBS!

} # if: $bvalflag

print OUT <<'!NO!SUBS!';

=head2 rndf()

=for ref

Reads a piddle from a NDF format data file.

=for example

 $pdl = rndf('file.sdf');
 $pdl = rndf('file.sdf',1);

The '.sdf' suffix is optional. The optional second argument turns off
automatic quality masking and returns a quality array as well.

A hash options version of rndf() also exists

 $pdl = rndf( 'file.sdf', nomask => 1, quickread => 1 );

where the supported options are (defaults in square brackets):

  nomask => (boolean) Disable automatic quality masking. [false]
  quickread => (boolean) If true, only read the primary data array and
            FITS header [false]

Header information and NDF Extensions are stored in the piddle as a hash
which can be retreived with the C<$pdl-E<gt>gethdr> command.
Array extensions are stored in the header as follows:

 $a - the base DATA_ARRAY

If C<$hdr = $a-E<gt>gethdr>;

then:

 %{$hdr}        contains all the FITS headers plus:
 $$hdr{Error}   contains the Error/Variance PDL
 $$hdr{Quality} The quality byte array (if reqeusted)
 @{$$hdr{Axis}} Is an array of piddles containing the information
                for axis 0, 1, etc.
 $$hdr{NDF_EXT} Contains all the NDF extensions
 $$hdr{Hist}    Contains the history information
 $$hdr{NDF_EXT}{_TYPES}
                Data types for non-PDL NDF extensions so that
                wndf can reconstruct a NDF.

All extension information is stored in the header hash array.
Extension structures are preserved in hashes, so that the PROJ_PARS
component of the IRAS.ASTROMETRY extension is stored in
C<$$hdr{NDF_EXT}{IRAS}{ASTROMETRY}{'PROJ_PARS'}>. All array structures are
stored as arrays in the Hdr: numeric arrays are stored as PDLs,
logical and character arrays are stored as plain Perl arrays. FITS
arrays are a special case and are expanded as scalars into the header.

PDL does not have a signed byte datatype, so any '_BYTE' data
is read into a C<byte> (unsigned) piddle and a warning is printed
to C<STDOUT>.

=cut

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';
=for bad

If the starlink bad flag is set, then the bad flag on the output
piddle is also set.  Starlink bad values are converted to the
current bad value used by the piddle type (if they are different).

=cut

!NO!SUBS!
} # if: $bvalflag

print OUT <<'!NO!SUBS!';
# This is one form of the new command

sub rndf {PDL->rndf(@_)}

# And this is the real form
# Allows the command to be called in OO form or as a function
sub PDL::rndf {  # Read a piddle from a NDF file

  my $class = shift;
  barf 'Usage: $a = rndf($file,[%args|$nomask]]); $a = PDL->rndf(...) '
    if (@_ == 0);

  # ensure we have the NDF module
  _init_NDF();

  # Setup the Header Hash
  my $header = {};

  # Read in the filename
  # And setup the new PDL
  my $file = shift;
  my $pdl = $class->new;

  # only one argument suggests we are using the old API, two
  # or more indicates the hash API
  my %opts = (nomask => 0, quickread => 0);
  if (@_ == 1) {
    $opts{nomask} = shift;
  } else {
    my %args = @_;
    %opts = iparse( \%opts, \%args );
  }

  my ($infile, $status, $indf, $entry, @info, $value);

  # Strip trailing .sdf if one is present
  # remove the trailing extension names
  # (anything after a `.' UNLESS it's part of a directory
  # name, hence the check for ^/) and any trailing (...)
  # specifier (NDF sections, a la PDL slices)
  #
  $file =~ s/\.sdf$//;

  $infile = $file;
  $infile =~ s!\.[^/]*$!!;
  $infile =~ s/\(.*\)$//;

  # If file is not there
  croak "Cannot find $infile.sdf" unless -e "$infile.sdf";

  # Set status
  $status = &NDF::SAI__OK;

  # Begin NDF and ERR session
  err_begin($status);
  ndf_begin();

  # Open the file
  ndf_find(&NDF::DAT__ROOT(), $file, $indf, $status);

  # unset automatic quality masking if nomask is true
  ndf_sqmf(!$opts{nomask}, $indf, $status);

  # Read the data array (need to pass the header as well)
  my $badflag;
  ( $status, $badflag ) = rdata($indf, $pdl, \%opts, $header, $class, $status);

  # Read the axes
  $status = raxes($indf, $pdl, $header, $class, $status) unless $opts{quickread};

  # Read the header
  $status = rhdr($indf, $pdl, $header, \%opts, $class, $status);

  # Read history information
  $status = rhist($indf, $pdl, $header, $status) unless $opts{quickread};

  # Read labels
  @info = ('Label', 'Title', 'Units');
  for $entry (@info) {
    $value = 'NULL';
    ndf_cget($indf, $entry, $value, $status);
    $$header{"$entry"} = $value if $value ne 'NULL';
    print "$entry is $value\n" if ($PDL::verbose && $value ne 'NULL');
    undef $value;
  }

  # Tidy up
  ndf_annul($indf, $status);
  ndf_end($status);

  if ($status != &NDF::SAI__OK) {
    my $starstat = $status;
    my $starerr = err_flush_to_string($status);
    err_end($status);
    croak("rndf: Starlink NDF error reading $file ($starstat):\n$starerr\n");
  }

  err_end($status);

  # Put the header into the main pdl
  $pdl->sethdr($header);

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';
  # set the bad flag if it was set in the input file
  print "Bad flag status is: $badflag\n" if $PDL::verbose;
  $pdl->badflag(1) if $badflag;

!NO!SUBS!

} # if; $bvalflag

  print OUT <<'!NO!SUBS!';
  return $pdl;
}

=head2 wndf()

=for ref

Writes a piddle to a NDF format file:

=for example

   $pdl->wndf($file);
   wndf($pdl,$file);

wndf can be used for writing PDLs to NDF files.
The '.sdf' suffix is optional. All the extensions
created by rndf are supported by wndf.  This means that error, axis
and quality arrays will be written if they exist. Extensions are also
reconstructed by using their name (ie FIGARO.TEST would be expanded as
a FIGARO extension and a TEST component). Hdr keywords Label, Title
and Units are treated as special cases and are written to the label,
title and units fields of the NDF.

Header information is written to corresponding NDF extensions.
NDF extensions can also be created in the {NDF} hash by using a key
containing '.', ie {NDF}{'IRAS.DATA'} would write the information to
an IRAS.DATA extension in the NDF. rndf stores this as
C<$$hdr{NDF}{IRAS}{DATA}> and the two systems are interchangeable.

rndf stores type information in {NDF}{'_TYPES'} and below so that
wndf can reconstruct the data type of non-PDL extensions. If no
entry exists in _TYPES, wndf chooses between characters, integer and
double on a best guess basis.  Any perl arrays are written as CHAR
array extensions (on the assumption that numeric arrays will exist as
PDLs).

=cut

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';
=for bad

The bad flag is written to the file, if set.

=cut

!NO!SUBS!

} # if: $bvalflag

print OUT<<'!NO!SUBS!';
# This is one form of the new command
# OO version.

*wndf = \&PDL::wndf;

# And this is the real form
# Allows the command to be called in OO form or as a function
sub PDL::wndf {  # Write a PDL to an NDF format file

  barf 'Usage: wndf($pdl,$file)' if $#_!=1;

  # ensure we have the NDF module
  _init_NDF();

  my ($indf, $place, $status, $outndf);
  my (@lbnd, @ubnd);

  my ($pdl, $outfile) = @_;

  # Check that we are dealing with a PDL
  croak 'Argument is not a PDL variable' unless $pdl->isa('PDL');

  # Set status
  $status = &NDF::SAI__OK;

  # Strip trailing .sdf if one is present
  $outfile =~ s/\.sdf$//;

  # Begin context
  ndf_begin();
  err_begin($status);

  # Open the new file
  ndf_place(&NDF::DAT__ROOT(), $outfile, $place, $status);

  # Need to create data_array component
  # Change the bounds to match the PDL

  @ubnd = $pdl->dims;
  @lbnd = ();
  @lbnd = map { 1 } (0..$#ubnd);

  ndf_new($pdltypes{$pdl->get_datatype}, $#ubnd+1, \@lbnd, \@ubnd, $place,
          $outndf, $status);

  # Set the application name
  ndf_happn('PDL::wndf', $status);

  # Write the data to this file
  $status = wdata($outndf, $pdl, $status);

  # Write the header and title
  $status = whdr($outndf, $pdl, $status);

  # Create a history extension
  ndf_hcre($outndf, $status);

  # Annul the identifier
  ndf_annul($outndf, $status);

  # End context
  ndf_end($status);

  if ($status != &NDF::SAI__OK) {
    my $starstat = $status;
    my $starerr = err_flush_to_string($status);
    err_end($status);
    croak "wndf: Starlink error writing NDF $outfile ($starstat)\n$starerr\n";
  }

  err_end($status);

  return 1;
}

=head2 propndfx()

Routine to write a PDL to an NDF by copying the extension information
from an existing NDF and writing DATA,VARIANCE, QUALITY and AXIS info
from a PDL (if they exist).

Extensions, labels and history are propogated from the old NDF.
No new extension information is written.

This command has been superseded by L<wndf()|/wndf()>.

=cut

*propndfx = \&PDL::propndfx;

sub PDL::propndfx {  # Write a PDL to a NDF format file

  barf 'Usage: propndfx($pdl, $infile, $outfile)' if $#_!=2;

  # ensure we have the NDF module
  _init_NDF();

  my ($indf, $in_place, $status, $outplace, $outndf);

  my ($pdl, $infile, $outfile) = @_;

  # Set status
  $status = &NDF::SAI__OK;

  # Strip trailing .sdf if one is present
  $outfile =~ s/\.sdf$//;
  $infile =~ s/\.sdf$//;
  # File is the first thing before a .
  my $file = (split(/\./,$infile))[0];

  croak "$file does not exist\n" unless -e "$file.sdf";

  # Begin context
  ndf_begin();
  err_begin($status);

  # Open the old file
  ndf_find(&NDF::DAT__ROOT(), $infile, $indf, $status);

  # Set the application name
  ndf_happn('PDL::propndfx', $status);

  # Open the new file
  ndf_place(&NDF::DAT__ROOT(), $outfile, $outplace, $status);

  # Copy the extension information to outfile
  ndf_scopy($indf, ' ', $outplace, $outndf, $status);

  # Annul the input file identifier
  ndf_annul($indf, $status);

  # Write the data to this file
  $status = wdata($outndf, $pdl, $status);

  # Annul the identifier
  ndf_annul($outndf, $status);

  # End context
  ndf_end($status);

  if ($status != &NDF::SAI__OK) {
    my $starstat = $status;
    my $starerr = err_flush_to_string($status);
    err_end($status);
    croak "propndfx: Starlink error propagating from $infile to $outfile ($starstat)\n$starerr\n";
  }

  err_end($status);

  return 1;
}

=head1 NOTES

The perl NDF module must be available. This is available from
Starlink (http://www.starlink.ac.uk).

If an NDF is read which contains AST World Coordinate information
(a .WCS component) this information is currently ignored. Currently
WCS information can only be written and stored using standard FITS headers.
See L<http://www.starlink.ac.uk/ast> for more information on AST.

=head1 AUTHOR

This module was written by Tim Jenness E<lt>tjenness@cpan.orgE<gt>.

PDL bad value support was added by Douglas Burke E<lt>dburke@cfa.harvard.eduE<gt>.

Copyright (C) 1998 - 2012 Tim Jenness. All Rights Reserved.  This program is free
software; you can redistribute it and/or modify it under the same
terms as Perl itself.

=head1 SEE ALSO

L<PDL::FAQ> for general information on the Perl Data language,
L<NDF> for information on the NDF module.

=cut


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

# These are the generic routines used by the rndf command

#  RDATA
#    Read Data, Quality and Variance

# returns ( $status, $badflag ), where $badflag is:
#
#   1 - the piddle may contain bad data
#   0 - the piddle contains no bad data
#
# This reports starlink's bad-data flag (ie it's valid
# even in bad value support is not available in PDL)
#
sub rdata {
  my ($indf, $pdl, $optref, $header, $class, $status) = @_;

  return ( $status, undef ) if $status != &NDF::SAI__OK;

  my ($maxdims, $ndim, @dim, $dcomp, $tcomp, $exist);
  my ($type, $data_pntr, $el, $temppdl, $nbytes, $badbit, $dref);

  ####################################################################
  #                      D  I  M  S                                  #
  ####################################################################

  # Find the dimensions of the data array
  $maxdims = 100;		# Maximum number of dimensions allowed

  &ndf_dim($indf, $maxdims, \@dim, $ndim, $status);
  @dim = @dim[0..$ndim-1];

  print "Dims = @dim\n" if $PDL::verbose;

  $$header{Extensions} = [];

  ####################################################################
  #    D A T A  +  V  A  R  I  A  N  C  E  +  Q U A L I T Y          #
  ####################################################################

  my $badflag = 0;

  my @comps = ( "Data" );
  if (!$optref->{quickread}) {
    # we can read additional components
    push(@comps, "Error");
    push(@comps, 'Quality') if $optref->{nomask};
  }

  for $dcomp (@comps) {

    $tcomp = $dcomp;
    $tcomp = 'Variance' if $tcomp eq 'Error';

    ndf_state($indf, $tcomp, $exist, $status);

    if ($exist && ($status == &NDF::SAI__OK)) {

      # Find the type of the data array
      ndf_type($indf, $tcomp, $type, $status);

      # Map the data array
      ndf_map($indf, $dcomp, $type, 'READ', $data_pntr, $el, $status);

      # Check the bad-data status of this component
      my $this_badflag = 0;
      my $ecomp = ( $dcomp =~ /Err/ ? "Variance" : $dcomp);
      ndf_bad($indf, $ecomp, 0, $this_badflag, $status );

      if ($status == &NDF::SAI__OK) {

	print "Reading $dcomp array...\n" if $PDL::verbose;

	# combine the bad data flag with the "global" one for
	# the file
	$badflag |= ($this_badflag != 0);

	# Create a temporary PDL to deal with separate components
	if ($dcomp eq 'Data') {
	  $temppdl = $pdl;
	} else {
	  $dcomp = 'Errors' if $dcomp eq 'Error';
	  $$header{"$dcomp"} = $class->new;
	  $temppdl = $$header{"$dcomp"};
	}

	# Set up the PDL
	$temppdl->set_datatype($startypes{$type});
	$temppdl->setdims([@dim]);
	$dref = $temppdl->get_dataref;

	# warning about data loss
	if ( $type eq "_BYTE" ) {
	    warn "Warning: NDF type (signed byte) converted to PDL's unsigned byte.\n";
	}

	# How many bytes in this data type?
	$type =~ s/^_|^_U//;
	$nbytes = &NDF::byte_size($type);

	mem2string($data_pntr, $nbytes*$el, $$dref);
        $temppdl->upd_data();

	push(@{$$header{Extensions}}, $dcomp) if $dcomp ne 'Data';

	if ($dcomp eq 'Quality') {
	  # Quality bad bit mask
	  ndf_bb($indf, $badbit, $status);

	  my $qhdr = {};
	  $$qhdr{BADBIT} = $badbit;
	  $temppdl->sethdr($qhdr);
	}

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';
	# if badflag is true, convert any STARLINK bad values
	# to PDL values
	if ( $badflag and
	     !compare_bad_values($temppdl->get_datatype) ) {
	    my $star_bad = starlink_bad_value($temppdl->get_datatype);
	    print "Converting from STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
	    $temppdl->inplace->setvaltobad( $star_bad );
	}

!NO!SUBS!

} # if: $bvalflag

print OUT <<'!NO!SUBS!';
	# Free temporary PDL
	undef $temppdl;
      }

      ndf_unmap($indf, $tcomp, $status);
    }
  }

  return ( $status, $badflag );
}



  ####################################################################
  #                      A  X  I  S                                  #
  ####################################################################

# Axes are stored as follows:
#  Array of PDL axes as an array in header of main pdl (called 'Axis')
#  Header of this PDL contains label and errors and units etc.

sub raxes {
  my ($indf, $pdl, $header, $class, $status) = @_;

  return $status if $status != &NDF::SAI__OK;

  my ($there, $naxis, @dims, $axcomp, $exist, $axtype, $axpntr, $el);
  my ($nbytes, $temppdl, $tcomp, $daxref, $dref);
  my ($axhdr, $ndims);

  # Read in axis information
  ndf_state($indf, 'Axis', $there, $status);

  $ndims = $pdl->getndims; # Find number of dimensions

  if ($there) {

    # Label from 0..$#dims to match array index handling
    # Want to put the axes into an array of axes
    # And update the extensions field
    $$header{Axis} = [];
    push(@{$$header{Extensions}}, "Axis");

    for $naxis (0..$ndims-1) {
      # Does a CENTRE component exist
      $axcomp = 'CENTRE';
      ndf_astat($indf, $axcomp, $naxis+1, $exist, $status);

      if ($exist && ($status == &NDF::SAI__OK)) {

	print "Reading Axis $naxis...\n" if $PDL::verbose;
	ndf_atype($indf, $axcomp, $naxis+1, $axtype, $status);

	ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
		 $el, $status);

	# Set up new PDL for axis info if map was okay
	if ($status == &NDF::SAI__OK) {
	  print "Number of elements: $el\n" if $PDL::verbose;
	  $$header{Axis}[$naxis] = $class->new;
	  $temppdl = $$header{Axis}[$naxis];
	  $temppdl->set_datatype($startypes{$axtype});
	  $temppdl->setdims([$el]);
	  $dref = $temppdl->get_dataref;


	  # How many bytes?
	  $axtype =~ s/^_|^_U//;
	  $nbytes = &NDF::byte_size($axtype);

	  mem2string($axpntr, $nbytes*$el, $$dref);
	  $temppdl->upd_data();

	  # Header info for axis
	  $axhdr = {};   # somewhere to put the pdl

	  # Read VARIANCE array if there
	  $axcomp = 'Error';
	  $tcomp = 'Variance';
	  ndf_astat($indf, $tcomp, $naxis+1, $exist, $status);

	  if ($exist && ($status == &NDF::SAI__OK)) {
	    print "Reading Axis $naxis errors...\n" if $PDL::verbose;
	    ndf_atype($indf, $tcomp, $naxis+1, $axtype, $status);

	    ndf_amap($indf, $axcomp, $naxis+1, $axtype, 'READ', $axpntr,
		     $el, $status);

	    # Set up new PDL for axis info if map was okay
	    if ($status == &NDF::SAI__OK) {

	      $$axhdr{Errors} = $class->new;
	      $$axhdr{Extensions} = ["Errors"];

	      $$axhdr{Errors}->set_datatype($startypes{$axtype});
	      $$axhdr{Errors}->setdims([$el]);
	      $daxref = $$axhdr{Errors}->get_dataref;

	      # How many bytes?
	      $axtype =~ s/^_|^_U//;
	      $nbytes = &NDF::byte_size($axtype);

	      mem2string($axpntr, $nbytes*$el, $$daxref);
	      $$axhdr{Errors}->upd_data();

	    }
	  }

	  # Get label and units
	  for my $entry ('Units', 'Label') {
	    ndf_astat($indf, $entry, $naxis+1, $exist, $status);
	    if ($exist) {
	      my $value = '';
	      ndf_acget($indf,$entry, $naxis+1, $value, $status);
	      $axhdr->{"$entry"} = $value;
	    }
	  }
	}
	ndf_aunmp($indf,'*',$naxis+1,$status);

	# Now store this header into temppdl
	$temppdl->sethdr($axhdr);

      }

    }

  }

 return $status;
}


  ####################################################################
  #                      E  X  T  E  N  S  I  O  N  S                #
  ####################################################################

sub rhdr {
  my ($indf, $pdl, $header, $optref, $class, $status) = @_;

  return $status if $status != &NDF::SAI__OK;

  my ($nextn, $xname, $xloc, $extn, $size, $entry, $value);
  my ($comment, $item, $el);

  # Read in any header information from the extensions

  ndf_xnumb($indf, $nextn, $status);
  print "Reading in header information...\n"
    if (($nextn > 0) && $PDL::verbose);

  # Read in extensions one at a time
  for $extn (1..$nextn) {

    $xname = 'NULL';
    ndf_xname($indf, $extn, $xname, $status);

    undef $xloc;
    ndf_xloc($indf, $xname, 'READ', $xloc, $status);

    # We only read a FITS header if quickread is enabled
    if ($optref->{quickread} && $xname ne 'FITS') {
      dat_annul( $xloc, $status );
      next;
    }

    # FITS is a special case and must be expanded
    if ($xname eq 'FITS') {
      # Read in the data
      my @fits;
      dat_size($xloc, $size, $status);
      &dat_getvc($xloc, $size, \@fits, $el, $status);

      if ($status == &NDF::SAI__OK) {
	print "Reading FITS header...\n" if $PDL::verbose;
        # unfortunately we can not store a tied Astro::FITS::Header object
        # because we are storing other things in the header
        for my $card (@fits) {
          my $item = Astro::FITS::Header::Item->new( Card => $card );
          $header->{$item->keyword} = $item->value;
          $header->{'_COMMENTS'}{$item->keyword} = $item->comment;
        }
      }
    } else {

      # Read in extensions to $EXTNAME
      $status = &find_prim($xloc, \%{$$header{$EXTNAME}}, $class, $status);
    }

    dat_annul($xloc, $status);
  }
  return $status;
}

  ####################################################################
  #                      H  I  S  T  O  R  Y                         #
  ####################################################################

sub rhist {
  my ($indf, $pdl, $header, $status) = @_;

  return $status if $status != &NDF::SAI__OK;

  my $has_hist;
  ndf_state( $indf, "HISTORY", $has_hist, $status );
  return $status unless $has_hist;

  my ($nrec, $app, $date, $i);

  # Read HISTORY information into a Hist header

  ndf_hnrec($indf,$nrec,$status);

  if ($status == &NDF::SAI__OK && ($nrec > 0)) {
    print "Reading history information into 'Hist'...\n" if $PDL::verbose;

    $$header{Hist}{'Nrecords'} = $nrec;

    # Loop through the history components and find last "APPLICATION"
    for $i (1..$nrec) {

      ndf_hinfo($indf, 'APPLICATION', $i, $app, $status);
      ndf_hinfo($indf, 'DATE', $i, $date, $status);

      $$header{Hist}{"Application$i"} = $app;
      $$header{Hist}{"Date$i"} = $date;

      print "  $app on $date\n" if $PDL::verbose;
    }
  }
  return $status;
}



################ Generic routines ###########################
#################### INTERNAL ###############################

# Find primitive components below a given HDS locator
# FITS information is stored in Hdr
# NDF extensions stored in NDF

sub find_prim {

  my ($loc, $href, $class, $status) = @_;
  my ($prim, $type, $size, @dim, $ndims, $struct, $name, $nloc, $el);
  my ($value, @values, $item, $entry, $maxdims);
  my ($packtype, $ncomp, $packed, $comp);
  my ($pntr, $nbytes, $comment, $dref);

  # Return if bad status
  return  $status if $status != &NDF::SAI__OK;

  # Most dimensions
  $maxdims = 7;

  # Is it a primitive component
  dat_prim($loc, $prim, $status);

  # Type, size and shape
  dat_type($loc, $type, $status);
  dat_size($loc, $size, $status);
  dat_shape($loc, $maxdims, \@dim, $ndims, $status);
  dat_name($loc, $name, $status);

  if ($prim) {
    print "\tReading component: $name ($type) " if $PDL::verbose;

    # Store type of information
    $$href{"_TYPES"}{"$name"} = $type;

    # Read it into the header
    if (($size == 1) && ($ndims == 0)) {
      dat_state( $loc, my $isdefined, $status );
      if (!$isdefined) {
          $value = undef;
      } elsif ($type  eq '_DOUBLE') {
	# dat_get0c doesn't return full precision for doubles
	dat_get0d($loc, $value, $status);
        err_rep("", "Error reading $type $name scalar as double",
                $status ) if $status != &NDF::SAI__OK;
      } else {
	dat_get0c($loc, $value, $status);
        err_rep("", "Error reading scalar $name with type $type as string",
                $status ) if $status != &NDF::SAI__OK;
      }

      $$href{"$name"} = $value if ($status == &NDF::SAI__OK);
      undef $value;

    } else {
      @dim = @dim[0..$ndims-1];

      # Inform the user of dimensions
      print "[@dim]" if $PDL::verbose;

      # Map arrays (don't need to unpack them)
      if ($type =~ /CHAR|LOG/) {

	# Read in the data
	dat_getvc($loc, $size, \@values, $el, $status);

	if ($status != &NDF::SAI__OK) {
	  err_rep("", "Error reading $name of type $type as string.",
                  $status );
	} else {
	  $$href{"$name"} = [@values];
	  undef @values;
	}

      } else {

	# map the data
	dat_mapv($loc, $type, 'READ', $pntr, $el, $status);

	if ($status != &NDF::SAI__OK) {
          err_rep("", "Error mapping $name with type $type", $status);
	}


	if ($status == &NDF::SAI__OK) {
	  # Create the pdl
	  $$href{"$name"} = $class->new;
	  $$href{"$name"}->set_datatype($startypes{$type});
	  $$href{"$name"}->setdims([@dim]);
	  $dref = $$href{"$name"}->get_dataref();

	  # How many bytes in this data type?
	  $type =~ s/^_|^_U//;
	  $nbytes = &NDF::byte_size($type);

	  mem2string($pntr, $nbytes*$el, $$dref);
	  $$href{"$name"}->upd_data();
	}

	dat_unmap($loc,$status);
      }
    }
    print "\n" if $PDL::verbose;

  } else {

    my ($ndimx, @dim, $ndim);
    dat_shape($loc, 7, \@dim, $ndim, $status);

    # If we have a single dimension ($ndim=0) then
    # proceed. Cant yet cope with multi-dimensional
    # extensions
    if ($ndim == 0) {
      dat_ncomp($loc, $ncomp, $status);
      print "$name ($type) has $ncomp components\n" if $PDL::verbose;

      # Store type of extension
      $$href{$name}{"_TYPES"}{STRUCT}{"$name"} = $type;

      for $comp (1..$ncomp) {
	dat_index($loc, $comp, $nloc, $status);
	$status = &find_prim($nloc, \%{$$href{$name}}, $class, $status);
	dat_annul($nloc, $status);
      }
    } else {
      print "Extension $name is an array structure -- skipping\n"
	if $PDL::verbose;
    }

  }
  return $status;
}


###### Routines for WRITING data ######################################

sub wdata {

  my ($outndf, $pdl, $status) = @_;
  my (@bounds, $ndims, @lower, @comps, $dcomp, $temppdl, $type);
  my ($pntr, $el, $nbytes, $axis, $axcomp, $axpntr, %hdr, %axhdr);
  my ($entry, $value, $i, $axndims, $axtype, $tcomp, @acomps);


  # Return if bad status
  return $status if $status != &NDF::SAI__OK;

  ####################################################################
  #                      D  I  M  S                                  #
  ####################################################################
  # Change the bounds to match the PDL
  @bounds = $pdl->dims;
  $ndims = $#bounds;
  @lower = map { 1 } (0..$ndims);

  &ndf_sbnd($ndims+1, \@lower, \@bounds, $outndf, $status);


  ####################################################################
  #    D A T A  +  V  A  R  I  A  N  C  E  +  Q U A L I T Y          #
  ####################################################################
  # Map data, variance, quality...
  @comps = ('Data','Errors','Quality');

  # Retrieve header and check that the header is a hash ref
  my $hdr = $pdl->gethdr;
  if (ref($hdr) eq 'HASH') {
    %hdr = %$hdr;
  } else {
    %hdr = ();
  }

  for $dcomp (@comps) {

    if ($dcomp eq 'Data') {
      $temppdl = $pdl;
    } else {
      if (exists $hdr{"$dcomp"}) {
        $temppdl = $hdr{"$dcomp"};
        # Check that we have a PDL
        next unless UNIVERSAL::isa($temppdl, 'PDL');
      } else {
        next;			# Skip this component
      }
    }

    # Check that we have some data
    if ($temppdl->nelem > 0) {

      if (($dcomp eq 'Quality') && ($temppdl->get_datatype != $PDL_B)) {
        # Quality must be bytes so convert to BYTE if this is not true
        $temppdl = byte ($$pdl{"$dcomp"});
      }

      $type = $pdltypes{$temppdl->get_datatype};
      $type = '_UBYTE' if $dcomp eq 'Quality'; # No choice for QUALITY

      # Set the output type of the NDF
      $tcomp = $dcomp;
      $tcomp = 'Variance' if $tcomp eq 'Errors';
      ndf_stype($type, $outndf, $tcomp, $status) unless $dcomp eq 'Quality';

      print "Mapping $dcomp , $type\n" if $PDL::verbose;
      # Map NDF
      $dcomp = 'Error' if $dcomp eq 'Errors';
      ndf_map($outndf, $dcomp, $type, 'WRITE', $pntr, $el, $status);

      if ($status == &NDF::SAI__OK) {
        # Number of bytes per entry
        $nbytes = PDL::Core::howbig($temppdl->get_datatype);

        # Convert to 1D data stream
        my $p1d = $temppdl->clump(-1);
        # Total number of bytes
        $nbytes *= $p1d->getdim(0);

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';
	# if badflag is true, convert any PDL bad values to STARLINK ones
	# note: we have to create a copy of the data to ensure that the
	# bad data value of the original piddle does not get changed.
	# per-piddle bad values would obviate the need for this.
	if ( $temppdl->badflag() ) {
	    print "Setting bad flag for component $dcomp ...\n" if $PDL::verbose;
	    ndf_sbad( 1, $outndf, $dcomp, $status );
	    if ( !compare_bad_values($temppdl->get_datatype) ) {
		my $star_bad = starlink_bad_value($temppdl->get_datatype);
		print "Converting from PDL to STARLINK bad value ($star_bad)...\n" if $PDL::verbose;
		$temppdl = $temppdl->setbadtoval( $star_bad ); # do NOT do inplace
	    }
	}

!NO!SUBS!

} # if: $bvalflag

print OUT <<'!NO!SUBS!';
        # Copy to disk
        string2mem(${$temppdl->get_dataref}, $nbytes, $pntr);

        # Write badbit mask
        if ($dcomp eq 'Quality') {
          ndf_sbb($$temppdl{Hdr}{BADBIT}, $outndf, $status)
            if defined $$temppdl{Hdr}{BADBIT};
        }
      }
      # Unmap Data
      ndf_unmap($outndf, $tcomp, $status);

    }

    # free the temporary pdl
    undef $temppdl;

  }


  ####################################################################
  #                      A  X  I  S                                  #
  ####################################################################
  # A X I S information is expected as an array in the header
  # called 'Axis'. These PDLs contain the CENTRE data and any further
  # info is stored in their header.
  # @bounds is accessed to find expected shape


  # Simply look in the header for a Axis
  if (exists $hdr{Axis}) {
     # Check that we have an array
     if (ref($hdr{Axis}) eq 'ARRAY') {
       # Now loop over axes
       for my $i (0..$#{$hdr{Axis}} ) {
	 # Loop unless status is bad
	 last unless $status == &NDF::SAI__OK;

	 # Extract the ith axis PDL from the array
         my $axis = $hdr{Axis}->[$i];

         # If we have a PDL
         if (UNIVERSAL::isa($axis, 'PDL')) {
            # We now want to copy the data and if necessary the
            # Error array. Since there are only two I will do it the
            # long way by explcitly storing data and then error

            # Set data type
            $axtype = $pdltypes{$axis->get_datatype};
            ndf_astyp($axtype, $outndf, 'Centre', $i+1, $status);

            # Okay we can now map this pdl
            ndf_amap($outndf, 'Centre', $i+1, $axtype, 'WRITE', $axpntr,
                   $el, $status);

            # Check that we have the correct number of entries
	    my $nelem = $axis->nelem;
            if ($el == $nelem) {
	      print "Mapping axis " , $i+1 , "\n"  if $PDL::verbose;

              # Number of bytes per entry
              $nbytes = PDL::Core::howbig($axis->get_datatype) * $el;

              # Copy to disk
              string2mem( $ { $axis->get_dataref }, $nbytes, $axpntr)
                 if ($status == &NDF::SAI__OK);

            } else {
              carp "Axis ",$i+1 .
                   " is the wrong size ($el values required but got ".
                   $nelem . ")- ignoring";
            }
            # Unmap
            ndf_aunmp($outndf, '*', $i+1, $status);

            # Errors
            my $axhdr = $axis->gethdr; # Retrieve and check header
            if (ref($axhdr) eq 'HASH') {
              %axhdr = %$axhdr;
            } else {
              %axhdr = ();
            }

	    # Look for an Errors component in the header hash
            if (exists $axhdr{Errors}) {
               my $axerr = $axhdr{Errors};
               if (UNIVERSAL::isa($axerr, 'PDL')) {

                 # Set data type
                 $axtype = $pdltypes{$axerr->get_datatype};
                 ndf_astyp($axtype, $outndf, 'Variance', $i+1, $status);

                 # Okay we can now map this pdl
                 ndf_amap($outndf, 'Errors', $i+1, $axtype, 'WRITE',
                   $axpntr, $el, $status);

                 # Check that we have the correct number of entries
		 my $nelem = $axerr->nelem;
                 if ($nelem == $el) {
                   print "Mapping errors for axis " . $i+1 . "\n"
                     if $PDL::verbose;

                   # Number of bytes per entry
                   $nbytes = PDL::Core::howbig($axerr->get_datatype) * $el;

                   # Copy to disk
                   string2mem($ {$axerr->get_dataref}, $nbytes, $axpntr)
                     if ($status == &NDF::SAI__OK);

                 } else {
                    carp "Error PDL for Axis ",$i+1,
                       " is the wrong size ($el values required but got ".
                        $axerr->nelem . ")- ignoring";
                 }
                 # Unmap
                 ndf_aunmp($outndf, '*', $i+1, $status);

               }

            }

            # Now character components
            foreach my $char (keys %axhdr) {
              if ($char =~ /LABEL|UNITS/i) {
                $value = $axhdr{"$char"};
                ndf_acput($value, $outndf, $char, $i+1, $status)
                  if length($value) > 0;
              }

            }

         }

       }

     }

  }

  return $status;
}


# Write header information to NDF extension
# This routine does not reconstruct an NDF identically to one that
# was read in. FITS extensions are made. All non-integer numeric types
# are written as doubles (apart from PDLs which carry there own type
# information) unless _TYPES information exists in {NDF}.


sub whdr {

  my ($outndf, $pdl, $status) = @_;

  my (%header, @fitsdim, $fitsloc, $value);
  my (%unused, @fits, $hashref, $hdr);

  # Return if bad status
  return $status if $status != &NDF::SAI__OK;

  print "Writing header information...\n" if $PDL::verbose;

  # Write FITS header from {Hdr}

  # Retrieve and check header from PDL
  $hdr = $pdl->gethdr;
  if (ref($hdr) eq 'HASH') {
    %header = %$hdr;
  } else {
    %header = ();
  }

  %unused = ();
  @fits = ();

  foreach my $key (sort keys %header) {

    next if $key eq '_COMMENTS';

    if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
      # This is not extension info
      ndf_cput($header{$key}, $outndf, $key, $status)
        if length($header{$key}) > 0;
      next;
    }

    # Only write scalars
    if (not ref($header{"$key"})) {

       my $item = Astro::FITS::Header::Item->new( Keyword => $key,
                                                  Value => $header{$key},
                                                  Comment => $header{'_COMMENTS'}{$key});
       if (defined $item) {
         push(@fits, $item->card() );
       } else {
         $unused{$key} = $header{$key};
       }
    } else {
      # Deal with later
      $unused{$key} = $header{$key};
    }

  }

  # Write the FITS extension
  if ($#fits > -1) {
    push(@fits, "END");
    # Write it out
    @fitsdim = ($#fits+1);
    ndf_xnew($outndf, 'FITS', '_CHAR*80', 1, @fitsdim, $fitsloc, $status);
    dat_put1c($fitsloc, $#fits+1, \@fits, $status);
    dat_annul($fitsloc, $status);
  }

  # Loop through all NDF header info and array Hdr information
  # Note that %unused still contains our NDF hash so
  # we must remove that before proceeding (else we end up with
  # an 'NDF_EXT' extension!

  delete $unused{$EXTNAME};

  foreach $hashref (\%unused, $header{$EXTNAME}) {
     $status = whash($hashref, $outndf, '', '', $status);
  }
  return $status;
}


#-----------------------------------------------------
# This routine writes a hash array to a NDF extension
# Keys that have leading underscores are skipped
# $hashname contains the extension name to use by default.
# This is extended  by any '.' in the key.

sub whash {

  my ($hash, $outndf, $hashname, $stypes, $status) = @_;

  my ($comp, $extension, $xtype, @dim, $ndims, $loc);
  my (@locs, $there, $type, $packed, $length);
  my (@bounds, $nbytes, $pntr, $maxsize, $which, $el);
  my ( $structs);

  my %header = (defined $hash ? %$hash : () );
  my $root = 1 if length($hashname) == 0; # Mark a ROOT structure

  # Return if bad status
  return $status if $status != &NDF::SAI__OK;

  foreach my $key (sort keys %header) {

    next if $key =~ /^_/;

    # If the key matches the Extensions (Axis or errors)
    # then skip. Should probably use the Extensions array
    # itself to work out which components to skip

    next if $key eq 'Axis';
    next if $key eq 'Errors';
    next if $key eq 'Extensions';
    next if $key eq 'Hist';

    # Deal with HASH arrays early
    ref($header{"$key"}) eq 'HASH' && do {
      my $nhashname = $hashname . ".$key";
      $nhashname =~ s/^\.//; # remove leading dot
      my $nstypes = $stypes . ".".$header{"$key"}{'_TYPES'}{STRUCT}{"$key"};
      $nstypes =~ s/^\.//; # remove leading dot
      $status = whash(\%{$header{"$key"}}, $outndf, $nhashname, $nstypes, $status);
      next;
    };

    my $path = $hashname . ".$key";
    $path =~ s/^\.//; # remove leading dot

    if ($key =~ /^TITLE$|^UNITS$|^LABEL$/i) {
      # This is not extension info
      ndf_cput($header{$key}, $outndf, $key, $status);
    } else {

      # Find nested structures
      my @structures = split(/\./, uc($path));

      # Last entry in structure is the important name
      $comp = pop(@structures);

      # Find list of structures
      $structs = join(".", @structures);
      $structs = 'PERLDL' unless (defined $structs && $structs =~ /./);
      $stypes = 'PERLDL_HDR' unless (defined $stypes && $stypes =~ /./);

      ($status, $loc) = mkstruct($outndf, $structs, $stypes, $status);

      # Now write the component
      # $loc is a locator to the last structure component
      $type = $header{"_TYPES"}{"$key"};

      # What is the data type?
      unless ($type =~ /./) {
        # number or character or array or pdl
        # All arrays are chars - should be PDLs if nums
        ref($header{"$key"}) eq 'ARRAY' && ($type = "_CHAR");

        ref($header{"$key"}) eq 'PDL' &&
          ($type = $pdltypes{$header{"$key"}{"Datatype"}});

        ref($header{"$key"}) || do {
          if ($header{"$key"} =~ /^(-?)(\d*)(\.?)(\d*)([Ee][-\+]?\d+)?$/) {
            $header{"$key"} =~ /\.|[eE]/
              && ($type = "_DOUBLE")
                ||   ($type = "_INTEGER");
          } else {
            ($type = '_CHAR')   # Character
          }
        };

      }

      # Create and write component
      # PDLs

      if (ref($header{"$key"}) eq 'PDL') {
        my $pdl = $header{"$key"};
        @bounds = $pdl->dims;
	my $n = $#bounds + 1;
        dat_new($loc, $comp, $type, $n, \@bounds, $status);
        $nbytes = PDL::Core::howbig($pdl->get_datatype) *
            $pdl->nelem;
        cmp_mapv($loc, $comp, $type, 'WRITE', $pntr, $el, $status);
        string2mem(${$pdl->get_dataref}, $nbytes, $pntr)
          if ($status == &NDF::SAI__OK);
        cmp_unmap($loc, $comp, $status);
      }

      # SCALARS
      ref($header{"$key"}) || do {

        if ($type =~ /^_CHAR/) {
          $length = length($header{"$key"});
          dat_new0c($loc, $comp, $length, $status);
          cmp_put0c($loc, $comp, $header{"$key"}, $status)
            if $length > 0;
	} elsif ($type =~ /^_LOG/) {
	  # In this case, add a check for FALSE or TRUE as strings
	  dat_new0l($loc, $comp, $status);
	  my $value = $header{$key};
          if (defined $value) {
            if ($value eq 'FALSE') {
                $value = 0;
            } elsif ($value eq 'TRUE') {
                $value = 1;
            }
            cmp_put0l($loc, $comp, $value, $status);
          }
        } else {
          # HDS does not have b,ub,w or uw types for new0 and put0
          # so we have to do create the item explicitly and put a bigger
          # type.
          my @dims = (0);
          dat_new( $loc, $comp, $type, 0, \@dims, $status );
          if ( $type =~ /_(BYTE|UBYTE|UWORD|WORD|INTEGER)$/) {
            cmp_put0i($loc, $comp, $header{$key}, $status);
          } elsif ($type =~ /_(REAL|DOUBLE)$/) {
            cmp_put0d($loc, $comp, $header{$key}, $status);
          } else {
            if ($status == &NDF::SAI__OK) {
              $status = &NDF::SAI__ERROR;
              err_rep("", "Unrecognized HDS type $type when writing component $comp",
                      $status);
            }
          }
        }
      };

      # CHARACTER ARRAYS
      ref($header{"$key"}) eq 'ARRAY' && do {
        # First go through the array and remove any
        # entries that are references. I draw the line at storing
        # nested arrays
        my @norefs = grep { not ref($_) } @{$header{"$key"}};

        # Determinine maximum string length.
        $maxsize = 0;
        foreach (@norefs) {
          my $size = length($_);
          $maxsize = $size if $size > $maxsize;
        }

        $ndims = 1;
        @dim = ($#norefs+1);
        dat_newc($loc, $comp, $maxsize, $ndims, @dim, $status);
        cmp_putvc($loc, $comp, $dim[0], \@norefs, $status);
      };

      # Annul all the locators
      dat_annul($loc, $status);

    }
  }

  return $status;
}

# Subroutine to make extension structures of arbritrary depth
# Just pass dot separated string to this function
# A locator to the final structure is returned

sub mkstruct {

  my ($outndf, $structs, $types, $status) = @_;
  my ($extension, @locs, @dim, $ndims, $xtype, $loc, $retloc);
  my ($set, $prmry, $there);

  my (@structures) = split('\.', $structs);
  my (@types) = split('\.', $types);

  return ($status, &NDF::DAT__NOLOC()) if $status != &NDF::SAI__OK();

  if ($#structures < 0) {
    $status = &NDF::SAI__ERROR;
    err_rep("", "Too few structures provided for mkstruct function (possible programming error)",
            $status);
    return ($status, &NDF::DAT__NOLOC());
  }

  # Does extension exist

  $extension = shift(@structures);
  $xtype = shift(@types);

  ndf_xstat($outndf, $extension, $there, $status);

  # Make it if necessary
  unless ($there) {
    print "Writing $extension extension ($xtype)...\n" if $PDL::verbose;
    @dim = (0);         # No pdl extensions!
    $ndims = 0;
    ndf_xnew($outndf, $extension, $xtype, $ndims, @dim, $loc, $status);
  } else {                      # Get the locator to the extension
    ndf_xloc($outndf, $extension, 'UPDATE', $loc, $status);
  }


  @locs = ();           # store the locators so that they can be annulled later
  push(@locs, $loc);

  # Make the structures, end up with a locator to the correct component
  for (0..$#structures) {
    dat_there($loc, $structures[$_], $there, $status);
    unless ($there) {   # Make it if it isnt there
      my $type = $types[$_];
      @dim = (0);
      dat_new($loc, $structures[$_], $type, 0, @dim, $status);
    }
    dat_find($loc, $structures[$_], $loc, $status);
    push(@locs, $loc);  # Store the new locator
  }

  # Clone the locator
  dat_clone($loc, $retloc, $status);
  $set = 1;
  $prmry = 1;
  dat_prmry($set, $retloc, $prmry, $status);

  # Annul all except the last locator
  foreach (@locs) {
    dat_annul($_, $status);
  }

  return ($status, ($status == &NDF::SAI__OK ? $retloc : &NDF::DAT__NOLOC() ));
}

# end of module
1;

!NO!SUBS!