=head1  NAME

V4 -- a 4-vector handling package

=head1 DESCRIPTION

Implements an object called a "4-vector" that keeps track of a
location in spacetime.  4-vectors can be converted between a variety
of known coordinate systems.  The system is designed to be easily
extensible -- you can add a new coordinate system to the soup by
writing a single pair of interconversion routines to/from any known
system and the new one.  The package automatically finds the shortest
conversion route from that system to any of the desired ones.

In keeping with the PDL philosophy, an individual V4 object can hold a 
whole list of points -- just add extra dimensions for threading.

=head1 SYNOPSIS

    use V4;

    $a = v4('Solar coords',['time','arcsec','arcsec','solar-radii'],$coordslist)
    $b = $a->convert('Carrington');
    $c = $a->convert('Carrington',['time','radian','radian','furlong'])

    Conversions are handled behind-the-scenes.  You can get the default
    units or specify your own (provided that they're interconvertible).

=head1 DATA STRUCTURE

V4's are implemented as PDLs, but usage of crossways slicing is deprecated,
because that may lead to difficulties later.  In particular, I reserve
the right to store times as strings and not numbers, which will require
that the single PDL inside the V4 be broken out into separate variables.
That will break any code that uses PDL::slice() instead of V4::x() to 
access the raw structures.

Ancillary information about the coordinate system is stored in the V4 --
currently in the PDL header.  The values of interest are:

  coords - the name of the coordinate system used for the V4.
  
  units  - an array ref with the names of the units for each coord (0-3).

  sys_params - An optional hash ref containing ancillary information
        that's required for some of the coordinate systems, in particular
        derived systems.  You can pass this in through the v4() 
        constructor method and out via the parms() method.

        Sys_params is used by the offset (derived) coordinate systems
        to store the required derivation information:

        
        REF_SYSTEM  The name of the reference system (unused at present,
	            as heliographic coords are the only supported 
		    offset-reference system just now)

        ORIGIN      The origin of the offset system, in HC coordinates
                    (this is a string to be eval'ed to get the V4 back --
                    a workaround to the PDL-in-header problem).

        Z_ORIG      If set, then the reference angle for the ROT
                    keyword points directly away from the origin of the 
                    original system (to match Heliographic Cartesian -- Z
                    toward the viewer!).  If it exists, the ROT matrix is 
                    applied to the origin-facing system.

        ROT         A 3-list containing the three rotation angles about 
		    X, Y, Z (in that order) to get from orig. to offset coords.


=head1 Shortcomings & caveats

This code was written long before the PDL::Transform module.  It
reallyl needs to be overhauled to take advantage of the newer
functionality.

For now, time handling is severely limited for the scientific world.
Everything is converted to UNIX time using Time::ParseDate.  More
effort should be given to time handling -- in particular, a routine
is needed that understands the difference between UT and TAI.
It would be terrific to augment Math::Units::convert to handle 
TAI, UT, Carrington second, Carrington hour, Julian second, Julian day,
and other numerical time formats.  Needless to say, parsedate and
strftime are a little underpowered for this stuff.

All of the existing conversions are non-relativistic.

More support is needed for relatively specified coordinate systems:
that is to say, for a general spherical system at some arbitrary location
in one of the other named systems. 

Currently, only 1-D lists of vectors are supported.  This makes some
of the slicing and converting easier, at the cost of some hassle in 
the interface if you want 2-D lists of points.  

Currently, the only derived systems are referred to the Heliographic
system.  

This is a port of the SSW/IDL code V4.pro, which offers much the same
functionality.

=head1 Author, license, no warranty

Copyright 2000, Craig DeForest.  This code may be modified and distributed
under the same terms as PDL itself.  The code comes with NO WARRANTY.

=head1 FUNCTIONS

=cut

  $PDL::V4::VERSION = 0.6;

######################################################################
package PDL::V4;

BEGIN {
package PDL::V4;
$VERSION = 0.3;

use Exporter;

@ISA = ("PDL");
@EXPORT = qw ( v4 xform canonicalize copy append replicate);
@EXPORT_OK = qw ( v4 xform canonicalize resolve unit copy append );
%EXPORT_TAGS = (Func=>[ qw ( v4 xform canonicalize copy appnd ) ]);

sub PDL::V4::append (;@);
sub PDL::V4::index (;@): lvalue;
sub PDL::V4::xform;
sub PDL::V4::x  : lvalue;
sub PDL::V4::x0 : lvalue;
sub PDL::V4::x1 : lvalue;
sub PDL::V4::x2 : lvalue;
sub PDL::V4::x3 : lvalue;
}

######################################################################
## Initialize constants!
#

BEGIN { 
  *debug = \0;

  use PDL;
  print __PACKAGE__,"loading...\n" if($debug);
  use PDL::V4::Xforms;

  use Math::Units;
  use Time::CTime;
  use Time::ParseDate;

  *convert = \(&{Math::Units:convert});
  *PI = \3.141592653589793238462643383279502;
  *NaN = \PDL::asin(pdl(2)); # E Pluribus unum -- non-equality doesn't
                             # guarantee non-nanness!
			     # (use (($a * 0) != 0) to test for NaNness)
}
                      

######################################################################
=pod

=head2 v4 

=for ref

4-vector constructor -  takes a PDL and some other info, returns a V4.

=for usage
    
    $a = v4(<coord-name>, [<hash-with-hdr>], [<units array ref>], <data>)

=for example
  
  $a = v4( "Coord-system-name", {hdr}, [<units array ref>], <PDL> )
  $a = v4( "coord-system-name", [<units>], {hdr}, ($x0,$x1,$x2,[$x3]) )
  $a = v4( "coord-system-name", [<units>], {hdr}, [1],[p2],[p3],...)
  $a = v4( <Hash-with-name-and-units>, [p1],[p2],[p3],...)
  $a = v4( <Hash-with-name-and-units>, {hdr}, <PDL> )

You feed in the name of a coordinate system, optional ancillary header
information, an optional array of unit names, and the data to V4ize. You
get back a V4 of undef.

Note that when supplying data points as array refs, a units array ref is
REQUIRED to avoid ambiguity!  If the hash syntax is used, then the hash
ref should contain keys "coords" (with the coord system name) and 
"units" (an array ref with the units to be used for each coordinate).  If 
the "units" key is missing, default units are extracted from the coordinate
system.

If you pass in a PDL, it must be at most two-dimensional.  The first coordinate
loops across dimensions and the second coordinate loops across points.

=cut

sub v4 {
  local $_;

  if($_[0] eq 'PDL::V4') {
    Carp::confess("V4 Funky constructor flag tripped...  someone's been using V4s as PDLs...\n");
  }

### Parse out unit names and coordinate system from the leading args.
  my($units);
  my($coords) = shift;
  if(ref $coords eq "HASH") {
      ($coords,$units) = ($coords->{'coords'},$coords->{'units'});
  } else {
      $units = shift if(ref $_[0] eq 'ARRAY');
  }

### Check that the coordinate system really exists and that the
##   units are compatible with the defined ones in the coordinate system.
  
  my($type) = $systems{$coords};
  Carp::confess "Couldn't find coordinate system named '$coords'\n" unless ($type);
  
  if(! defined($units)) {
      $units = $type->{'units'};
  } elsif( ((ref $units) ne 'ARRAY') || (@{$units} != 4)) {
      barf "Mangled units sent to V4...";
  } 
  
  my($sys_params) = (ref$_[0] eq 'HASH') ? shift : {};

## Kludge until we handle timestamps right - skip units containing "time".
## Check that units agree with the defaults!
  for(my ($i)=0;$i<4;$i++) {
      if(!($systems{$coords}->{'units'}->[$i]) =~ m/time/) {
	  barf "Unit $units[$i] don't agree with canonical $units[$i] ".
	      "(coord $i)\n" 
		  unless (defined 
			  convert(1,
				  $units[$i],
				  $systems{$coords}->{'units'}->[$i]
				  ));
      }
  }


## Parse the remaining arguments and stuff 'em into a PDL.
  my $a;
  while(scalar(@_)) {
      $a = shift;

      if('' eq ref $a) {
	  $a = [$a,@_];
	  @_ = ();
      }
      
      if('ARRAY' eq ref $a) {
	  ## Kludge for timestamp handling...  zeroth coord is always time
	  ## (so far).  parsedate will ignore all timezone specifiers after
	  ## the first, so we append a "UT" before converting, to make sure
	  ## we default to the right zone.
	  
	  ## This is complex as we might have been passed a list of 
	  ## times 4-vectors and not just a single 4-vector.

	  if('ARRAY' eq ref($a->[0])) { # Got a list of 4-vectors...
	      for(@{$a}) { 
		  $_->[0] = Time::ParseDate::parsedate(($_->[0]).' UT')
		      if(($_->[0]).0 ne ($_->[0]));
	      }
	  } elsif((($a->[0]+1)-1) ne $a->[0]) { # If $a->[0] is a string...
	      $a->[0] = Time::ParseDate::parsedate($a->[0].' UT'); 
	  }
	  $a = pdl($a);
      }elsif('PDL' eq ref $a) {
	  barf "v4: v4's can only have up to 2 dimensions" if($a->dims > 2); 
	  if($a->dims <= 1) {
	    ## Assume one PDL per coordinate
	    my($out) = zeroes(4,$a->nelem);
	    my($foo) =  $out->slice('(0)');
	    $foo .= $a;
	    
	    $foo = $out->slice('(1)');
	    $foo .= shift;
	    
	    $foo = $out->slice('(2)');
	    $foo .= shift;
	    
	    $foo = $out->slice('(3)');
	    $foo .= shift;

	    $a = $out;
	  }
	}
      
      
      if((($a->dims)[0]) < 4) {
	$a = PDL::append($a,zeroes(4-(($a->dims)[0])));
      }
    }
  
  my $self = {'coords'=> $coords
		,'units'=> ((defined $units) ? 
			    $units : 
			    $systems{$coords}->{'units'}
			    )
		  ,'sys_params' => $sys_params
		    ,'PDL' => $a
		    };
  
  bless $self,__PACKAGE__;
}

#*initialize = \&v4;
*main::v4 = \&v4;

######################################################################
=pod

=head2 append

=for ref

Like PDL::append, but for V4's -- appends in the 2nd dimension,
after making sure that all the coords are the same.

=for usage

    $a = $v->append($v2);

=cut
sub append(;@) {
  my($a,@victims) = @_;
  local(@_);
  $a = $a->copy;
  foreach $_(@victims) {
    my($foo) = $_->xform($a->{coords},$a->{sys_params});
    my($z);
    if($a->{PDL}->dims == 1) {
      $a->{PDL} = $a->{PDL}->dummy(1,1);
    }
    if($_->{PDL}->dims == 1) {
      $a->{PDL} = (PDL::append ($a->{PDL}->xchg(0,1),$_->{PDL}->dummy(1,1)->xchg(0,1)))->xchg(0,1);
    } else {
      $a->{PDL} = (PDL::append ($a->{PDL}->xchg(0,1),$_->{PDL}->xchg(0,1)))->xchg(0,1);
    }
    
  }
  return $a;
}

######################################################################
=pod

=head2 replicate

=for ref

Given a V4, makes <n> copies of it in a nice array

=cut
sub replicate(;@) {
  my($a,$n) = @_;
  
  if($a->{PDL}->is_inplace) {$b=$a;} else {$b=$a->copy;}
  $b->{PDL} = $b->{PDL}->dummy(1,1) unless ($b->{PDL}->dims>1);
  $b->{PDL} = $b->{PDL}->dummy(2,$n)->copy->xchg(2,0)->clump(2)->xchg(1,0);
  $b->set_inplace(0);
  $b;
}
    

######################################################################
=pod

=head2 index

=for ref

Like PDL::index, but for V4's -- index out individual points.

=for usage

    $a = $v->index(pdl(3,4,5));

=cut

sub index (;@) : lvalue {
    my($a,$index) = @_;
    $index = pdl($index) unless (ref $index eq 'PDL');
    my($out) = {PDL=>1};
    my $foo;
    foreach $foo(keys %{$a}){
      $out->{$foo} = main::deep_copy($a->{$foo}) unless($foo eq 'PDL');
    }

    print "V4 index: pdl dims: ",$a->{PDL}->dims if($debug);
    print "\nindex dim is ",pdl($index)->nelem,"; max is ",max($index),"; min is ",min($index),"\n" if($debug);
    $out->{PDL} = $a->{PDL}->dice((xvals zeroes 4),$index);
    print "v4::index returning" if($debug);
    bless($out,'V4');
    return $out;
}

	
######################################################################
#
=pod 

=head2 x, x0, x1, x2, x3

=for ref 

Access to individual coordinate slices of a V4

=for usage
    
    $v = v4('Heliographic',0,1,2,3)
    $r = $v->x3;

Provide a way to get out individual coordinates of a system.  These
are implemented as direct calls to slice() but are preferred over them
for extensibility reasons.  These return slice-like connected references --
you can .= them like you could PDLs (because, after all, that's what they are!

x0...x3 are syntactical sugar for x(0), x(1), etc.

Just now the x0...x3 work by calling x -- that's because of the extra
hassle involved in the ->dummy(0,1) workaround for the glitches with
0-dimensional piddles in PDL 2.005.

=cut

sub x0 : lvalue { return (my $a = x($_[0],0)) }
sub x1 : lvalue { return (my $a = x($_[0],1)) }
sub x2 : lvalue { return (my $a = x($_[0],2)) }
sub x3 : lvalue { return (my $a = x($_[0],3)) }
sub xlist {
    my $self = shift;
    local $_;
    my @out;
    if(!@_) { @_ = (0,1,2,3); }
    while(defined ($_=shift)) {
	my($a) = $self->{PDL}->slice("($_)");
	bless($a,"PDL");
	if($a->dims) {	
	    push @out,$a; 
	} else {
	    push @out,$a->dummy(0,1);
	}
    }
    return @out if(@out > 1);
    return $out[0];
}

sub x : lvalue {
    my $self = shift;
    local $_;
    my @out;
    if(!@_) { @_ = (0,1,2,3); }
    while(defined ($_=shift)) {
	my($a) = $self->{PDL}->slice("($_)");
	bless($a,"PDL");
	if($a->dims) {	
	    push @out,$a; 
	} else {
	    push @out,$a->dummy(0,1);
	}
    }
    my $out = (@out > 1 ? cat(@out) : $out[0]);
    return $out;
}

######################################################################
#
=pod

=head2 system

=for ref 

Returns the canonical name of the coordinate system associated with a v4.
Add additional parameters to set the name.

=for usage

    $a = system($name)->{'coords'};

=cut
sub system {
    if(@_>1) {
	$_[0]->{'coords'} = $_[1];
    }

    return $systems{$_[0]->{'coords'}}->{'coords'}->[0];
}

######################################################################
#
=pod

=head2 coords

=for ref

Return the standard name of the coordinate system in a V4

=for usage

    <out> = $v4->coords;

=for example

    $a = v4('HS',undef,0,0,1);
    print $a->coords;

Pass in a list of coords to retrieve the names of.  If you pass in
just one, you get back a scalar, otherwise a list.

=cut

sub coords {
    my $self = shift;

    
    if(!@_) {
	push(@_,0,1,2,3);
    }

    my @out;
    local $_;
    while($_=shift) {
	push(@out,systems->{$self->system}->{udesc}->[$_]);
    }
    return @out if(@out > 1);
    return $out[0];
}

######################################################################
=pod

=head2 units

=for ref 

Access to the units strings in a V4

=for usage 

    $a = $v->units;

=for example

    @a = units($v4)                               - Retrieve all units
    @a = units($v4,1,2)                           - Retrieve units for x1 & x2
    units($v4,[undef,undef,'arcsec','minutes']);  - Set units 2 & 3

Pass in an array ref to set units; a list of numbers to retrieve
specific ones.

=cut
    
sub units {
    my $self = shift;
    @_=(0,1,2,3) if(!@_);
    my @out;
    local $_;
	
    foreach $_(@_) {
	if(ref $_ eq 'ARRAY') {
	    print "UNIT SETTING NOT READY\n";
	}
	push @out,$self->{'units'}->[$_];
    }
    return @out if(@out > 1);
    return $out[0];
}


######################################################################
#
=pod

=head2 PDL::V4::string

=for ref

Stringify a V4 

This routine takes a V4 and returns a string.  It overlays the PDL
stringify-overlay method, in order to return the V4 header and unit
information.

=cut

sub string {
    my @out;

    my($self,$format)=@_;
    if($PDL::_STRINGIZING) {
	return "ALREADY_STRINGIZING_NO_LOOPS_V4";
    }
    local $PDL::_STRINGIZING = 1;
    my $ndims = $self->getndims;
    if($self->nelem > 500) {
	return "Too many vectors to print";
    }
    if($ndims == 0) {
	my @x = $self->at();
	return ($format ? sprintf($format,$x[0]) : "$x[0]");
    }
    return "Null" if $self->isnull;
    return "Empty" if $self->isempty;
    local $PDL::Core::sep  = $PDL::use_commas ? "," : " ";
    local $PDL::Core::sep2 = $PDL::use_commas ? "," : "";
    my $coord =  $systems{$self->{'coords'}} || {};
		
    push(@out,sprintf("4-vector:  %s\n",
		      ( (defined $coord) ? ($coord->{'coords'}->[0]) : 
			("invalid coordinate system '$self->{'coords'}'")
			)));
    
    local($i);
    for($i=0;$i<scalar(@{$coord->{'units'}});$i++){
	push(@out,sprintf("   x$i: %15s (%s)\n"
			  ,$coord->{'udesc'}->[$i]
			  ,$self->{'units'}->[$i]
	     ));
    }
    
    if($ndims == 1){
	return join("",@out,PDL::Core::str1D($self,$format));
    } else {
	return join("",@out,PDL::Core::strND($self,$format));
    }
}

use overload (
	      ( "\"\"" => \&string )
	      );


######################################################################
#
=pod

=head2 V4::resolve

=for ref

Find the best known conversion path between two coordinate systems.

=for usage
    @chain = &resolve($from,$to);

    $from and $to are names of coordinate systems.
    @chain gets a list of coordinate system names separated by 
    "atomic" transitions.  It's returned as an array so that
    callers have a convenient way of getting the number of steps
    (with $chain = &resolve($from,$to)).
    
    Meanwhile, the list of subroutines gets stored in the V4 global
    data structures. 

=for description

This is an internal method used by xform to find the proper sequence
of atomic conversions to call for a particular conversion.

The code works by traversing the existing @conversions array, looking for
the conversion in question.  If it's not there then a treewalker explores
possible lines of conversion.  The first path found is the shortest
possible, and that sequence is returned.  If there's no such path, then 
the empty list is returned.

resolve starts by testing the trivial chain (an atomic conversion).  Then
it starts building a list of all non-repeating conversion chains 
that start from A.  The central loop goes from L -> L+1 by scanning
the entire conversions database for conversions from the last item in
each existing chain.  The possible chain extensions are filtered -- no circular
chain is allowed.  The process repeats until there are NO chains of the
current length that start with A; or a chain is found that reaches B. 

The speed of the algorithm is something like O(N^2.5).  The fastest
method is faster, but more of a hassle to implement.  TIMTOWTDI.

History: Ported from V4_resolve_xform.pro, CED 1-Apr-2000

=cut

sub resolve {
  my ($from,$to) = @_;
  
  # Regularize coordinate systems and barf if nonexistent
  if(ref $from =~ m/V4/) {
      $from = $from->{'coords'};
  }
  $from = ($systems{$from}->{'coords'}->[0] || barf "V4::resolve: Unknown system `$from'");
  $to   = ($systems{$to}  ->{'coords'}->[0] || barf "V4::resolve: Unknown system `$to'");
  
  ##############################
  # Check for the trivial cases
  #
  return $from if($from eq $to);

  my $conv = $systems{$from}->{conversions};

  if(defined($conv->{$to})) {
    my $steps = (scalar(@{$conv->{$to}}));
    barf "V4::resolve: Impossible to convert from `$from' to `$to'!" unless ($steps);

    return ($to) if($steps == 1);
    
    my @route = @{$systems{$from}->{conv_route}->{$to}};
    barf "V4::resolve: Inconsistency in the databases! Help! (route has ".@route." steps; conv has $steps)" 
      if(scalar(@route) != $steps);

    return @route;
  }
  
  ######################################################################
  # Build chains of conversions that start from $from and fan out until
  # either we find no new ones or we reach maximal length.  We keep track
  # of where we've been with %beenthere to avoid loops.  The chains 
  # consist of linked lists back from the endpoint to the source.  
  # They're hashes with "val" and "back" tags.  
  #
  # We deliberately avoid taking advantage of previously discovered
  # routes -- this ensures we'll always find the shortest route through
  # the graph.

  my (%beenthere,$root,$current,$prev,@paths,$i);
  $beenthere{$from} = 1;
  $root = {val=>$from};
  $current = [$root];
  
  for($i=0; !$beenthere{$to} && scalar(@{$current}); $i++) {
    $prev=$current;
    $current = [];
    my $newfrom_ref, $newto;
    foreach $newfrom_ref(@{$prev}) {

      foreach $newto(keys %{$systems{$newfrom_ref->{'val'}}->{conversions}}) {
	next if($beenthere{$newto});

	my $newto_ref = $systems{$newfrom_ref->{'val'}}->{conversions}->{$newto};
	if(scalar(@{$newto_ref}) == 1) {  # Simple conversion found

	  push (@{$current},{ 'val' => $newto 
				,'back'=> $newfrom_ref } );
	  $beenthere{$newto} = 1;
	}
      }
    }
    push(@paths,@{$current});  # Cache discovered paths -- may come in handy.
  }


  ## Now we've discovered all paths FROM the originating system up to 
  ## the length of the first path that contains the desired system.
  ## Cache 'em all in the data structures! (saves another trip here.)

  my $to_ref,@route,@conv;
  my $begin_to_ref;

  foreach $begin_to_ref(@paths) {

    $to_ref = $begin_to_ref;
    for(@route=@conv=();defined $to_ref->{'back'};$to_ref=$to_ref->{'back'}) {
      unshift(@route,$to_ref->{'val'});
      unshift(@conv,
	      @{ $systems{ $to_ref->{'back'}->{'val'} } 
		 ->{conversions}->{ $to_ref->{'val'} } 
	       }
	      );
    }

    barf "V4::resolve: HELLoo -- conv has ".@conv." elements; route has ".@route 
      unless( @route == @conv );  # scalar comparison

    $systems{$from}->{conversions}->{$begin_to_ref->{'val'}} = [@conv];
    $systems{$from}->{conv_route} ->{$begin_to_ref->{'val'}} = [@route];
  }

  my $out = $systems{$from}->{conv_route}->{$to};
  if(!defined $out) {
    $systems{$from}->{conversions}->{$to}=[];
    $systems{$from}->{conv_route}->{$to}=[];
    barf "V4::resolve: No route found from `$from' to `$to'!\n\t(Marked it impossible.)\n";
  }
  
  return @{$out};
} 

    
######################################################################
=pod

=head2 copy

=for ref

Copies a 4-vector.  Overlay for the PDL ->copy.  It uses deep_copy 
to copy the structure.

=cut

sub copy {
    print "V4 copy...\n" if($debug);
    my(%out);
    my($key);
    my($in) = shift;
    foreach $key(keys %{$in}) {
	if(ref($in->{$key})) {
	    $out{$key} = main::deep_copy($in->{$key});
	} else{
	    $out{$key} = $in->{$key};
	}
    }
    return bless( \ %out, 'V4');
}
    
######################################################################
=pod

=head2 xform

=for ref

Convert a 4-vector to another coordinate system.

=for usage
  
    $a = xform( $b, 'SRS', [{Coord-system-parms}] );

=for description

xform just walks the conversion tables given to it by resolve(), 
calling the conversion routines listed in the main %V4::systems data
structure.

Note that because only one set of system_args is currently
supported -- and that is passed to each of the conversion
routines, whether it wants 'em or not! -- only one transition
(out of the whole chain) that supports arbitrary parameters is
currently allowed -- or, I suppose, multiple ones may be allowed,
provided that their parameters are appropriately selected to 
avoid collision!

=cut
sub xform {
    my($vect,$to,$system_args) = @_;
    
    # Make sure it's not the identity transform!
    if($systems{$vect->{'coords'}} == 
       $systems{$to}) {
	return $vect;
    }

    # Make sure path exists
    if(!scalar(resolve($vect->{coords},$to))) {
	barf "You can't get here!  (resolve should've barfed)";
    }

    $out = $vect->copy unless($vect->{PDL}->is_inplace);

    local $conv_sub;
    print "Converting: " if($debug);
    foreach $conv_sub(@{$systems{$vect->{coords}}->{conversions}->
			{$systems{$to}->{coords}->[0]}}) {
      $out = &{$conv_sub}((inplace $out),$system_args);
    }
    $out->{coords} = $to;
    $out->{sys_params} = $system_args if (defined($system_args));
    bless $out,'V4';
    return $out;
}


######################################################################
=pod

=head2 canonicalize

=for ref

Put a 4-vector into "standardized" units for that coordinate system

=for usage

    $a = $a->canonicalize();

=cut

sub canonicalize {
    my $a = @_[0];

    my ($type) = $systems{$a->{coords}};
    barf "canonicalize:  unknown system!" unless($type);

    $a->{'coords'} = $type->{coords}->[0];

    for(my $i=0;$i<4;$i++) {
	my $foo;
	$foo = $a->{PDL}->slice("($i)");
	$foo .= PDL::V4::Xforms::unit($foo,$a->{'units'}->[$i]
			      , $type->{'units'}->[$i]
		     );

	$a->{'units'}->[$i] = $type->{'units'}->[$i];
    }
    bless $a;
    return $a;
}
  

######################################################################
# Data structures follow.
#

=pod

=head2 @V4::systems, %V4::systems.

Array containing system definition data structures.  The array is 
defined in the source code; but the hash is filled by the package
constructor, so that all the official synonyms can be used to access
a package's information.


Each element is a hash ref.  

 Tags are:

       coords       A list of names for each system.  The 0th element
                    is a verbose name; the 1st is the canonical abbreviation.
                    Alternate names follow.

       units        A list of four unit-names containing the canonical
                    scientific units for the system (e.g. "km").

       udesc        A list of strings describing (briefly!) the meaning
                    of each coordinate.

       desc         A string containing a (one-line!) description of the 
	            coordinate system.
			
       metric       A string containing the name of the subroutine used
	            to measure distance in the coordinate system.  The metric
		    returns values in the canonical units for the coordinate
		    system.

       conversions  A hash of other coordinate systems to which we know
                    the conversions.  The value of each of these tags is
		     an array ref containing a sequence of subroutines 
		     to be called in sequence to get to the target.
		     An empty array ref signifies that it's impossible
		     to get there from here!

       conv_route   A hash of routes to other coordinate systems to which
		    we know the conversions.  Each of these tags points
		    to an array ref containing the names of the 
		    coordinate systems in the conversion.  Only nontrivial
		    routes are listed.  (The conv_route tag is written
		    to by the resolve routine and is not present in the
                    static data structure declared here.)

       cartesian    A string containing the name of the associated Cartesian
                    coordinate system, if any.  This is used for conversions
                    (such as derived or offset systems) that are generic but
                    only work between cartesian systems.  If the string
                    doesn't exist, there is no known associated Cartesian
                    system; and if it exists but is the empty string, then 
                    this system *is* Cartesian.

=head1 COORDINATE SYSTEMS

=cut

BEGIN {

=pod

=head2 Earth Range Solar ("ERS","Solar","S")

(Spherical) Earth-centered "observing coordinates".  These are used 
for historical reasons, though observing coordinates are probably best
treated as an EQA projection rather than true spherical.  This system
is centered on Earth, with argument along the Earth-Sun line and pole in
the plane defined by the E-S line and the Earth's north pole.  It is a
LEFT HANDED system -- longitude gets POSITIVE to the west.  Again, this
is for historical reasons -- the longitude is like the traditional 
"solar-X" and latitude is like "solar-Y".  This is the coordinate system used
at the focal plane of near-Earth solar observations.

=head2 Solar Range Solar ("SRS")

(Modified spherical) These are Earth Range Solar coordinates, with the
radius replaced by radius from Sun center.  This is not a complete coordinate
system, as *two* points satisfy any given 4-tuple (near side and far side
of the Sun).  These coordinates are used primarily as an intermediate 
step for heliographic conversion, and the conversion routines take the 
closer branch of the ambiguity (because you can't see the far side of the
surface of the Sun).

=cut

push(@systems,(
{coords=> ["Earth Range Solar","ERS","Solar","S"]
,units=> ["unixtime",	"arcsec",	"arcsec",	"solar-radii"]
,udesc=> ["time",	"W angle",	"Solar-N angle","distance"]
,desc=>  "Spherical coords; 2-D origin at Disk Center; P=0; arg. thru Sun; LEFT-HANDED!"
,metric=>"V4MSPH", cartesian=>"Solar Cartesian"
,conversions => {
    "Solar Cartesian" => [ \&PDL::V4::Xforms::ERS2SCart ]
   ,"Heliographic Cartesian" => [ sub { &PDL::V4::Xforms::ERS2SCart( $_[0], $_[1], 1); } ]
}}

,{coords=> ["Solar Range Solar","SRS"]
,units=> ["unixtime",	"arcsec",	"arcsec",	"solar-radii"]
,udesc=> ["time",	"W angle",	"Solar-N angle","altitude"]
,desc=> "Observing coords: Earth-cent. sph.; radius rel. to Sun center; P=0"
,metric=>"V4MSRS"
,conversions => {
}}
));

=pod

=head2 Solar Cartesian ("SC", "Earth Solar Cartesian", "ESC")

This is the cartesian coordinate system that corresponds best to the 
Earth-Range Solar and Solar-Range Solar coordinates.  It is a cartesian
system centered on the Sun's center and oriented to match ERS:  "Z"
is directly away from the Earth on the E-S line; "X" is positive-west,
and "Y" lies in the plane defined by the E-S line and the (Earth's) north
pole.

=cut

push(@systems,(
{coords=> ["Solar Cartesian","SC","Earth Solar Cartesian","ESC"]
,units=> ["unixtime",	"solar-radii", 	"solar-radii", 	"solar-radii"]
,udesc=> ["time", 	"W fr. Sun ctr","up fr. Sun ctr","toward Earth"]
,desc=>  "Sun-centered Earth-N-oriented cartesian (approx. SRS)"
,metric=>"V4MCART", cartesian=>""
,conversions => {
    "Earth Range Solar" => [ \&SCart2ERS ]
   ,"Heliographic Cartesian" => [ sub { SCart2HC($_[0],$_[1],1); } ]
}}
));

=pod

=head2 Heliographic Cartesian ("HC")

This is the cartesian system matching Heliographic coordinates.  It
is centered on the Solar center, +Y is true Solar north, +Z is approximately
away from the Earth, and +X is west. 

=cut

push(@systems,(
{coords=> ["Heliographic Cartesian","HC"]
,units=> ["unixtime",   "solar-radii",  "solar-radii", "solar-radii"]
,udesc=> ["time",  "W fr. Sun ctr.", "Solar-N fr. Sun ctr.", "toward Earth"]
,desc=> "Sun-centered solar-N-oriented cartesian"
,metric=>"V4MCART", cartesian=>""
,conversions => {
   "Heliographic"       => [sub { Cart2Sph($_[0],{ix=>[0,3,1,2],%{$_[1]}}); }]
   ,"Offset Cartesian"  => [  \&PDL::V4::Xforms::offset_cartesian  ]
   ,"Solar Cartesian"   => [  sub {SCart2HC($_[0],$_[1],-1); } ]
   ,"Earth Range Solar" => [  sub { PDL::V4::Xforms::SCart2ERS($_[0],$_[1],1); } ]
}}
));

=pod

=head2 Heliographic ("H")

These are Solar spherical coordinates -- *Solar* latitude and longitude,
with the prime meridian intersecting the Earth-Sun line.

=cut

push(@systems,(
{coords=> ["Heliographic","H"]
,units=> ["unixtime",	"deg", 		"deg", 		"solar-radii"]
,udesc=> ["time", 	"H-Longitude", 	"Latitude", 	"above Sun ctr"]
,desc=>  "Spherical, with pole at Solar North and prime meridian at disk ctr"
,metric=>"V4MSPH", cartesian=>"HC"
,conversions => {
    "Heliographic Cartesian"=> [ sub { Sph2Cart($_[0],{ox=>[0,3,1,2],coords=>'Heliographic Cartesian'}); } ]
   ,"Carrington"            => [ sub { barf "H -> C not here yet"; } ]
}}
));

=pod

=head2 Carrington (C)
 
The traditional Carrington solar spherical coordinates, in the frame
of the approximate overall solar rotation.  North is true Solar north;
longitude is the traditional Carrington longitude.

=head2 Carrington Cartesian (CC)

Cartesian system matching the Carrington system.  +Z rolls around with 
the Carrington prime meridian.

=cut

push(@systems,(
{coords=> ["Carrington", "C"]
,units=> ["unixtime", "deg", "deg", "solar-radii"]
,udesc=> ["time", "C-Longitude", "Latiude", "above Sun ctr"]
,desc=> "Pole at Solar North; prime meridian spinning with average Sun"
,metric=> "V4MSPH", cartesian=>"CC"
,conversions => {
   "Heliographic"            => [ sub {barf "C -> H not here yet"; } ]
  ,"Carrington Cartesian"   => [ sub {Sph2Cart($_[0],{ox=>[0,3,1,2],coords=>'Carrington Cartesian'}); } ]
}}

,{coords=> ["Carrington Cartesian", "CC"]
,units=> ["unixtime",	"km", 		"km", 		"km"]
,udesc=> ["time",	"Solar-E",	"Solar-N",	"Out pr. meridian"]
,desc=> "Cartesian coords matching the Carrington system"
,metric=> "V4MCART", cartesian=>""
,conversions => {
    "Carrington"         => [ sub {Cart2Sph($_[0],{ix=>[0,3,1,2],coords=>"Carrington"}); } ]
}}
));

=pod

=head2 Earth Solar Polar ("ESP", "zunwrap")

Spherical coordinates based at Earth; but the pole goes through the 
Earth-Sun line, and the argument is North.  The "longitude" is azimuth
(CW) around the Sun, and the colatitude measures impact parameter of each
line of sight.  (x2 is in fact colatitude rather than latitude.)  Useful
for interpreting radially-symmetric objects in solar observations.  Called
"zunwrap" after an earlier coordinate transform routine.

=cut

push(@systems,(
{coords=> ["Earth Solar Polar", "ESP", "zunwrap"]
,units=>  ["unixtime",	"deg",		"solar-radii",	"solar-radii"]
,udesc=>  ["time",	"Azimuth",	"degrees",      "Dist along E-S line"]
,desc=> "Zunwrapped solar coords -- azimuth and impact parameter of L.O.S."
,metric=> "V4MCYL", 
,conversions => {
}}

,{coords=>["Heliocentric Ecliptic Cartesian", "HEC", "HE", "Heliocentric Ecliptic"]
,units => ["unixtime",  "km",		"km",		"km"]
,udesc => ["time",	"toward Aries",	"perpl",	"Ecliptic North"]
,desc  => "Cartesian coords centered on Sun, aligned with Ecliptic"
,metric => "V4MCART", cartesian=>""
,conversions => {
}}
));

=pod

=head2 Offset Cartesian ("OC", "Perspective Cartesian", "PC")

Cartesian coordinates with parametrized origin and Euler angle rotation
from another coordinate system.

=head2 Offset Spherical ("Perspective Spherical", "PS", "OS")

Spherical coordinates with parametrized origin and Euler angle rotation
from another coordinate system (currently, the only parent system
that is supported is the Heliographic one).  [This system is derived
from OC]

=head2 Perspective ("Offset Equi-Angular Azimuthal", "P", "OEQA")

Equi-Angular Azimuthal Perspective coordinates with parametrized origin
and Euler angle rotation from another coordinate system.  x1 and x2 define
points in an image plane such that the distance to the origin in the image
plane is proportional to the angle between that line of sight and the 
"Z" axis, and azimuth is preserved.  x3 is radius from the origin.
[This system is derived from OC]

=cut

push(@systems,(
{coords=>["Offset Spherical", "Perspective Spherical", "PS", "OS"]
,units=> ["unixtime","deg","deg","solar-radii"]
,udesc=>["time","horiz. angle","vert. angle","range"]
,desc=>"Perspective:  Sph. coords w/ (lat,lon) origin pointed at major origin"
,metric => "V4MSPH", cartesian=>"Offset Cartesian"
,conversions => {
    "Offset Cartesian"          => [ sub {Sph2Cart($_[0],{ox=>[0,3,1,2],coords=>"Offset Cartesian",%{$_[1]}});} ]
    }
}

,{ coords=>[ "Offset Cartesian", "OC", "Perspective Cartesian", "PC" ]
,units=>["unixtime","solar-radii","solar-radii","solar-radii"]
,udesc=>["time","arb. X","arb. Y","arb. Z"]
,desc=>"Cartesian system offset at arb. angle from Heliographic Cartesian"
,metric=> "V4MCART", cartesian=>""
,conversions => { 
    "Heliographic Cartesian" => [ \&PDL::V4::Xforms::unoffset_cartesian ]
   ,"Offset Spherical"       => [ sub{Cart2Sph($_[0],{ix=>[0,3,1,2],coords=>"Offset Spherical"});} ]
   ,"Perspective"   => [ \&PDL::V4::Xforms::Cart2EQA ]
   }
}

,{ coords=>[ "Perspective", "Offset Equi-Angular Azimuthal", "P","OEQA" ]
,units=>["unixtime","deg","deg","solar-radii"]
,desc =>["Spherical projection preserving angular dist. fr. origin line"]
,udesc=>["time","proj. X", "proj. Y", "range"]
,metric=>"to OC", cartesian=>"Offset Cartesian"
,conversions=> {
	"Offset Cartesian"   => [ \&PDL::V4::Xforms::EQA2Cart ]
}}
));


}

BEGIN {
  local ($a);

  print "V4 initializing systems list...\n" if ($debug);
  foreach $a(@systems) {
    foreach $b(@{$a -> {'coords'}}) {
      $systems{$b}=$a;
    } 
  }
}

1;