=head1 NAME

V4::Xforms - Transformation formulae for 4-vectors

=head1 DESCRIPTION

Methods and functions for converting 4-vectors between coordinate
systems.  These routines all take a standard set of inputs and make
specific, direct conversions between coordinate systems.  They are
intended to be used primarily within V4, and hence aren't fully
protected against malicious input.  Furthermore, they don't keep proper
track of all the coordinate-system labeling that makes v4 useful -- only
coordinate units are updated.  Look in V4 for actual user-level 
methods.

=head1 NOTES

Because many of these routines are on the inside of tight loops,
extra ->copy operations are deprecated.  In particular, conversion
routines will clobber the input vector unless you explicitly copy it.
They should all recognize the ->inplace method.

=head1 SYNOPSIS

  use V4;  # Normal routines
use Xform; # magic happens here

=head1 METHODS

=cut


package PDL::V4::Xforms;
use PDL::V4;
*constants = \&main::constants;
constants();
print __PACKAGE__." loading...\n";

use PDL;
use Data::Dumper;

*rot_matrix = \ &main::rot_matrix;

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

=head2 unit

Local version of Math::Units::convert just puts the Math::Units version
in an eval to hide the results of Croak (returns undef on
error).

=cut

use Math::Units;
sub unit {
    my $out;
    eval q#$out = eval '&Math::Units::convert(@_)'#;
}

@EXPORT = qw( Sph2Cart Cart2Sph offset_cartesian, unoffset_cartesian,
		ERS2SCart, SCart2ERS, SCart2HC, PER2Cart,Cart2PER);
@EXPORT_OK = qw( indices unit );
%EXPORT_TAGS = (
     Func=>[qw(Sph2Cart Cart2Sph offset_cartesian unoffset_cartesian,
		ERS2SCart, SCart2ERS, SCart2HC )],
     Internal=>[qw/indices/],
);
use Exporter;
use DynaLoader;


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

=head2 V4::Xform::indices

Grab coordinate indirection indices or defaults from a hash ref.  Used
internally by the various conversions.

=cut

sub indices {
  $opt = $_[0];
  local($ix,$ox) = ($opt->{'ix'},$opt->{'ox'});
  $ix = [0,1,2,3] unless defined($ix);
  $ox = [0,1,2,3] unless defined($ox);
  print ":: ix=",join(",",@{$ix}),"; ox=",join(",",@{$ox})," ::" if($V4::debug);
  return ($ix,$ox);
}

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

=head2 Sph2Cart

Convert generic spherical to Cartesian coordinates
  
    $sph = Sph2Cart($coords,{options});


   Coord:       Interpretation (input)        Interpretation (output)
   0:          time  (left alone)             time 
   1:          longitude (angle)              X    [ ALL 3 
   2:          latitude  (angle)              Y      are the
   3:          radius    (arb.)               Z      same... ]

The pole points along the Z axis; the central meridian (0 longitude)
is through the X axis; and the Y axis forms a right-handed coordinate
system.  

History: Ported from v4xsph2cart.pro, CED 30-Mar-2000
    V. 0.5 -- Modified to use v4 hashes, 23-Apr-2000
   
=cut

sub Sph2Cart {
    print "Sph2Cart..." if($V4::debug);
    
    my ($in,$opt) = @_;
    my ($ix,$ox) = &indices($opt);
    
    ### We always have to copy the PDL, because we scrozzle it during
    ### the conversion.  But we can avoid the eval in v4 copying.
    ### (It would be conceptually cleaner but codily messier to 
    ### modify the original PDL using the copy as a buffer...)
    ### The way this is done, now the original PDL is left by the wayside 
    ### I'm really trusting the perl GC here.
    my($outpdl,$out);
    if($in->is_inplace) {
      $in->{PDL}->set_inplace(0);
      $outpdl = $in->{PDL}->copy;
      $in->{PDL}->set_inplace(1);
      $out = $in;
    } else {
      $out = $in->copy;
    }
# Cut input into longitude, latitude, and radius; and output into x,y,z
    my ($lon,$lat,$r) = ($in->xlist(@{$ix}[1..3]));

### Have to move the PDL AFTER taking the 'x', if we're doing inplace.
    if($in->is_inplace) {
	$out->{PDL} = $outpdl;
    }

    my ($x,$y,$z) = $out->xlist(@{$ox}[1..3]);

# Detach and convert to radians if necessary
    my($a,$b) =
	(unit(1,$in->units($ix->[1]),'rad')
	,unit(1,$in->units($ix->[2]),'rad'));

    barf("Lat and lon (theta and phi) must have angular units") unless($a && $b);
    $lon = $lon * $a if($a != 1);  # Sic -- detach
    $lat = $lat * $b if($b != 1);  # Sic -- detach
# Do the conversion.
    $x .= $r * cos($lat) * cos($lon);
    $y .= $r * cos($lat) * sin($lon);
    $z .= $r * sin($lat);
    
    if($ix->[0] != $ox->[0]) {
	my $foo = $out->slice("($ox->[0])");
	$foo .= $in->slice("($ix->[1])");
    }
	 
# Set the output units.
    #time

    local(@a);
    $a[$ox->[0]] = $in->units($ix->[0]);
    $a[$ox->[1]]=( $a[$ox->[2]]=( $a[$ox->[3]]= ($in->units($ix->[3]))));
    $out->{'units'} = \ @a;
    $out->{'coord'} = $opt->{coords};
# Return.
    return $out;
}
######################################################################

=pod

=head2 Cart2Sph

Convert generic Cartesian to spherical coordinates

    $sph = Cart2Sph($coords,{options});

    Coord:	Interpretation (input)		Interpretation (output)
    0:		time (left alone)		(left alone)
    1:		X [ All 3			longitude (angle)
    2:		Y   must be			latitude (angle)
    3:		Z   interconvertible ]          radius (units of X)

The pole points along the Z axis; the central meridian (0 longitude)
is through the X axis; and the Y axis forms a right-handed coordinate
system.  

See Sph2Cart for the inverse transform.

History: Ported from V4xcart2sph.pro, CED 31-Mar-2000
    Modified to use V4 hashes, 23-Apr-2000

=cut
    
sub Cart2Sph {
    print "Cart2Sph..." if($V4::debug);
    my($in,$opt) = splice(@_,0,2);
    my ($ix,$ox) = &indices($opt);
  
    my($out,$outpdl);
    $in->set_inplace(0);
    $out = $in->copy;

# Cut input into x,y,z; and output into lon, lat, rad.
    my ($x,$y,$z) = ($in->xlist(@{$ix}[1..3]));
    my($lon,$lat,$r) = ($out->xlist(@{$ox}[1..3]));
    
# If necessary, detach and convert y and z into X's units
    my $xu = $in->{'units'}->[$ix->[1]];
    my $yu = $in->{'units'}->[$ix->[2]];
    my $zu = $in->{'units'}->[$ix->[3]];

    if($yu ne $xu) {
	my($yy) = unit(1,$yu,$xu);  
	$y = $y * $yy if($yy != 1);
	barf("V4::Xforms::Cart2Sph: X and Y coords don't agree! ('$xu','$yu'") if(!$yy);
	
    }

    if($zu ne $xu) {
	my($zz) = unit(1,$zu,$xu);
	$z = $z * $zz if($zz != 1);
	barf("V4::Xforms::Cart2Sph: X and Z coords don't agree!") if(!$zz);
    }

# Do the conversion.
    my($r12) = ($x*$x+$y*$y);
    my ($r1) = sqrt($r12);
    my ($radeg) = unit(1,'rad','deg');
    
    $lon .= atan2($y,$x) * $radeg;
    $lat .= atan2($z,$r1)* $radeg;;
    print "post-lat z=$z\n" if($V4::debug);
    $r   .= sqrt($r12 + $z*$z);

    if($ix->[0] != $ox->[0]) {
	my $foo = $out->slice("($ox->[0])");
	$foo .= $in->slice("($ix->[0])");
    }
    
    
# Set output units
    my(@a);
    $a[$ox->[0]] = $in->units($ix->[0]);
    $a[$ox->[1]] = $a[$ox->[2]] = "degrees";
    $a[$ox->[3]] = $xu;
    $out->{'units'} = \ @a;

    return $out;
}

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

=head2 offset_cartesian


Convert between Heliographic Cartesian and Offset Cartesian
coordinates. 


    $foo = offset_cartesian($coords,{parameters});

    Coord:	Interpretation (input)		Interpretation (output)
    0:		time (left alone)		(left alone)
    1:		X [ All 3			longitude (angle)
    2:		Y   must be			latitude (angle)
    3:		Z   interconvertible ]          radius (units of X)

System parameters are indexed by hash:

    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.

=cut

$cartesian_identity = zeroes(4,4);
my($foo) = $cartesian_identity->diagonal(0,1);
$foo++;
    
sub offset_cartesian {
    print "offset_cartesian" if($V4::debug);
    my($in,$parameters) = @_;
    my ($matrix) = $cartesian_identity->copy;

    if( !defined($parameters)) {
      use Carp;
      Carp::confess("Hey!  offset_cartesian needs a parameter set!");
    };
    my ($t0,$x0,$y0,$z0) = $parameters->{ORIGIN}->xform('HC')->list;
    
    if(defined $parameters->{'Z_ORIG'}) {
      ## Note: the Z_ORIG matrix includes an implicit about-face
      ## to match observed coordinates.  I couldn't bring myself to 
      ## make it a left-handed system, so "X" ends up reversed. --CED
      my $xzr0 = sqrt($x0*$x0 + $z0*$z0);
      my $r0   = sqrt($x0*$x0 + $y0*$y0 + $z0*$z0);

      # Rotate:
      $matrix x= rot_matrix(4,pdl(1,3), -$z0/$xzr0,-$x0/$xzr0)if($xzr0); # Y 
      $matrix x= rot_matrix(4,pdl(2,3), $xzr0/$r0, -$y0/$r0)if($r0);     # X
      $matrix x= rot_matrix(4,pdl(1,2),$parameters->{'Z_ORIG'});         # Z
    }
    
    if (!defined $parameters->{'ROT'}) {$parameters->{'ROT'}=[0,0,0]}
    else {
      $matrix x= rot_matrix(4,pdl(2,3),$parameters->{'ROT'}->[0]);
      $matrix x= rot_matrix(4,pdl(1,3),$parameters->{'ROT'}->[1]);
      $matrix x= rot_matrix(4,pdl(1,2),$parameters->{'ROT'}->[2]);
    }

    if($in->is_inplace) {
      $out = $in;
    } else {
      $out = $in->copy;
    }

    if(defined($parameters->{'ORIGIN'})) {
	$out->{PDL} -= ($parameters->{'ORIGIN'}->xform('HC'));
      }


###################################################
# Because the time specification is often "nan" and doesn't
# interact with the rest of the coordinates, I treat NaNs
# in a way that's arguably Wrong -- they're treated as 0's 
# and then renanified after the matrix multiply.
    my $foo,$w;
    $w = which (0 * ($out->{PDL}->clump(-1)));
    $foo = ($out->{PDL}->clump(-1)->index($w));
    $foo .= 0;
    $out->{PDL} x= $matrix;
    $foo .= $V4::NaN;

    $out->{'sys_params'} = $parameters;
    return $out;
}


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

=head2 unoffset_cartesian

Converts from a derived cartesian coordinate system to the cartesian
system from which it is derived.  First rotates using the Euler angles
in the sys_params entry in the V4's header, then offsets back to the
origin of the original coordinate system.  The information in the 
sys_params field of the input V4's header is used for the Euler
angles and offset. 

If the sys_params field includes "Z_ORIG" then an initial rotation is 
applied that points the Z axis directly AWAY from the origin of the 
original system, reversing the X axis from its usual sense as well.
This is the Right Thing to do for perspective transformations (where
you want +Z to be toward the object and away from the observer in the 
observer's coordinates, but where +Z in the object's coordinates generally
points toward the observer).

At the moment, unoffset_cartesian ONLY takes you back to the 
Heliographic Cartesian system; but the mathematics are general
and it should be expanded to handle other systems as well.  (The
output ought to be in the coordinate system used by the ORIGIN field).

   $cart = unoffset_cartesian($coords,{options});
  

   Coord:      Interpretation (input)      Interpretation (output)
   0:          time (left alone)           time
   1-3:        input "X","Y","Z"           Output "X","Y","Z"

   Input header information:

    $coords->{sys_params}->{Z_ORIG} :
           If this field is *defined*, then the Euler angles are
           referred to an initial orientation where the -Z axis of the
           derived system points directly toward the origin of the 
           main system, X is nearly reversed, and Y is nearly the same 
           as its original configuration.

    $coords->{sys_params}->{ROT} :
           Contains the Euler angles used to generate the rotated
           derived system.  They're in degrees, and (because this
           is the UNoffset routine) applied in reverse order and 
           direction (i.e. [30,20,10] yields a 10 degree CCW rotation 
           about the Z axis, followed by 20 degrees CCW about the 
           *new* Y axis, and 20 degrees CCW about the *newer* X axis.
              0: Rotation about the X axis
              1: Rotation about the Y axis
              2: Rotation about the Z axis
           
    $coords->{sys_params}->{ORIGIN} :
           A V4 containing the origin of the derived system.  Because
           header-nested PDLs are currently broken, this is currently
           a string that evaluates to the correct V4; it needs fixing
           once the reference problem is fixed.

=cut

sub unoffset_cartesian {
    print "unoffset_cartesian..." if($V4::debug);
    my($in) = @_;
    my($origin_xformed);

    $parameters = $in->{'sys_params'};
    my $matrix;

    $matrix = $cartesian_identity->copy;
    
    $origin_xformed = $parameters->{ORIGIN}->xform('HC');
    my ($t0,$x0,$y0,$z0) = $origin_xformed->list;

    my $xzr0= sqrt($x0*$x0 + $z0*$z0);
    my $r0  = sqrt($x0*$x0 + $y0*$y0 + $z0*$z0);     
    
    if(defined $parameters->{'ROT'}) {
      $matrix x= rot_matrix(4,pdl(1,2),-$parameters->{'ROT'}->[2]);
      $matrix x= rot_matrix(4,pdl(1,3),-$parameters->{'ROT'}->[1]);
      $matrix x= rot_matrix(4,pdl(2,3),-$parameters->{'ROT'}->[0]);
      }	    
    
    if( (!defined($_uc_cache_origin))
	|| (sum ($_uc_cache_origin != $parameters->{ORIGIN}->{PDL}))
	|| ($_uc_cache_parameters->{ROT}->[0] != $parameters->{ROT}->[0])
	|| ($_uc_cache_parameters->{ROT}->[1] != $parameters->{ROT}->[1])
	|| ($_uc_cache_parameters->{ROT}->[2] != $parameters->{ROT}->[2])
	|| ($_uc_cache_parameters->{Z_ORIG}   != $parameters->{Z_ORIG}) ) {
	
      if(($r0 != 0) && (defined $parameters->{'Z_ORIG'})) {
	$matrix x= rot_matrix(4,pdl(1,2),-$parameters->{'Z_ORIG'});       #  Z
	$matrix x= rot_matrix(4,pdl(2,3),$xzr0/$r0, $y0/$r0)   if($r0);   #  X
	$matrix x= rot_matrix(4,pdl(1,3),-$z0/$xzr0,$x0/$xzr0) if($xzr0); #  Y
      }
      
      
      $_uc_cache_matrix = $matrix;
      $_uc_cache_parameters = main::deep_copy($parameters);
      $_uc_cache_origin = $parameters->{ORIGIN}->{PDL}->copy;
    }
    else {
      $matrix = $_uc_cache_matrix;
    }

    if($in->is_inplace) {
      $in->{PDL} x= $matrix;
      $in->{PDL} += $origin_xformed;
      $in->{coords}="Heliographic Cartesian";
      return $in;
    }
    
    $in->set_inplace(0);
    my($out) = $in->copy;
    $out->{PDL} = (($in->{PDL} x $matrix) + $origin_xformed);
    $out->{coords}="Heliographic Cartesian";
    
    return $out;
  }

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

=pod

=head2 Cart2EQA

Converts Cartesian coordinates to the equi-angular azimuthal spherical 
projection. (in EQA projection, the pole points along the "Z" axis, 
the [x0,x1]-azimuth is the same as the [X,Y]-azimuth, and radial
distance in the [x0,x1]-plane is the same as the subtended angle between
the point and the "Z" axis.  This is the system used by fisheye lenses.)
The first two coordinates are angular, in degrees, and the last coordinate
is the radius.

     $eqa = cart2eqa($coords,[{options}]);
	  
   Coord:       Interpretation (input)        Interpretation (output)
   0:           time  (left alone)            time 
   1:  	        "X"                           projection-x
   2:           "Y"                           projection-y
   3:           "Z"                           radius from origin

History:  Written 21-Apr-2000, CED, v. 0.1 (not yet tested!)
    0.5 Modified to use v4 hashes, 23-Apr-2000

=cut

constants('V4::Xforms');

sub Cart2EQA {
    print "Cart2EQA..." if($V4::debug);
    my($in,$opt) = @_;
    my($ix,$ox) = &indices($opt);

    my($out,$outpdl);
    if($in->is_inplace) {
	$out = $in;
	$in->set_inplace(0);
	$outpdl = $in->{PDL}->copy;
	$in->set_inplace(1);
    } else {
	$out = $in->copy;
    }

    my($x,$y,$z) =   $in->xlist(@{$ix}[1..3]);

    $out->{PDL} = $outpdl   if($in->is_inplace);

    my($xp,$yp,$r) = $out->xlist(@{$ox}[1..3]);
    
    # if necessary, detach and convert X and Y to same units as Z.
    my $xu = $in->units($ix->[1]);
    my $yu = $in->units($ix->[2]);
    my $zu = $in->units($ix->[3]);
    $y = units($y,$yu,$xu) if($yu ne $xu);
    $z = units($z,$zu,$xu) if($zu ne $xu);

    $r .= sqrt($x*$x + $y*$y + $z*$z);
    my($angle)   = PDL::acos($z/$r) * $RADEG;
    my($xyr) = pdl(sqrt($x*$x + $y*$y));  # = sin($angle)*$r
    $w = which($xyr == 0);

    $xp .= $angle * $x / $xyr;
    $yp .= $angle * $y / $xyr;

    ($foo = $xp->index($w)) .= 0;
    ($foo = $yp->index($w)) .= 0;

    local($_);
    $a[$ox->[0]] = $in->units($ix->[0]);
    $a[$ox->[1]] = "degrees";
    $a[$ox->[2]] = "degrees";
    $a[$ox->[3]] = $in->units($ix->[1]);
    $out->{'units'} = \@a;

    return $out;
}

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

=pod

=head2 EQA2Cart

Converts equi-angular azimuthal coordinates to Cartesian coordinates.
(See Cart2EQA for a description of the EQA projection).

     $eqa = cart2eqa($coords,[{options}]);
	  
   Coord:       Interpretation (input)        Interpretation (output)
   0:           time  (left alone)            time 
   1:  	        projection-x                  X
   2:           projection-y                  Y
   3:           radius                        Z

 History:  Written 21-Apr-2000, CED, v. 0.1
    Modified to use v4 hashes, 23-Apr-2000, v. 0.5

=cut

sub EQA2Cart {
    print "EQA2Cart..." if($V4::debug);
    my($in,$opt) = @_;
    my($ix,$ox) = &indices($opt);

    my($out,$outpdl);
    if($in->is_inplace) {
	$out = $in; 
	$in->set_inplace(0);
	$outpdl = $in->{PDL}->copy;
	$in->set_inplace(1);
    } else {
	$out = $in->copy;
    }
    
    my($xp,$yp,$r) = $in->xlist(@{$ix}[1..3]);

    $out->{PDL} = $outpdl if($in->is_inplace);

    my($x,$y,$z) = $out->xlist(@{$ox}[1..3]);

    my($xxp,$yyp) = 
	(unit($xp,$in->units($ix->[1]),'rad')
	,unit($yp,$in->units($ix->[2]),'rad')
	 );

    $angle = sqrt($xxp*$xxp + $yyp*$yyp);

    my($w) = which($angle == 0);
    if( $w->nelem ) {
      $angle->index($w) .= 2*3.1415926535897;
    }

    my($xypr) = $r * sin($angle);
    $x .= $xypr * ($xxp/$angle);
    $y .= $xypr * ($yyp/$angle) ;
    $z .= $r * cos($angle);

    local(@a);
    local($_);
    $a[$ox->[0]] = $in->units($ix->[0]);
    $a[$ox->[1]] = $a[$ox->[2]] = $a[$ox->[3]] = $in->units($ix->[3]);
    $out->{units} = \ @a;
    return $out;

}
			     
######################################################################
=pod

=head2 ERS2SCart

Convert from Earth Range Solar coordinates to Solar Cartesian
coordinates.  (ERS is spherical based at Earth, with Z axis along the
Earth-Sun line, and is a LEFT handed spherical coordinate system, by
observing convention (west is positive-X; north is positive-Y).


   Coord:       Interpretation (input)        Interpretation (output)
   0:          time  (left alone)             time 
   1:          W angle (angle)                W fr. Sun ctr 
   2:          N angle (angle)                N fr. Sun ctr 
   3:          radius  (arb.)                 Toward Earth  

 History:  Written 23-Apr-2000, CED
=cut

sub ERS2SCart {
    print "ERS2SCart..." if($V4::debug);
    my($in,$opt,$to_HC) = @_;

    my($i2,$i2pdl);
    if($in->is_inplace) {
	$i2 = $in;
    } else {
	$i2 = $in->copy;
    }
    
# Calculate the Earth-Sun distance...
    my($pb0rr) = main::pb0r($in->x0);
    
# Place the "plane of the sky" at the Sun if radius is undefined.
    my($foo);

    my($w) = which(0*($i2->x3));
    ($foo = $i2->x3->index($w)) .= $pb0rr->slice('(3)')->index($w);
    
# Convert from Earth L.H. spherical to Earth-centered L.H. cartesian...
    $i2 = Sph2Cart($i2,{ox=>[0,3,2,1]});
    
# Translate the origin from the Earth to the Sun & switch to R.H. cart.
    (($foo = $i2->x3) *= -1) += $pb0rr->slice('(3)');
    
    if($to_HC) {
	my($x2) = $i2->x2->copy;
	my($x3) = $i2->x3->copy;
	my($b) = $pb0rr->slice('(1)');
	
	$i2 x= rot_matrix(4,pdl(2,3),$b);

	$i2->{coords} = "Heliographic Cartesian";
    } else {

	$i2->{coords} = "Solar Cartesian";

    }

    return $i2;
    
}

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

=head2 SCart2ERS

Convert from Solar Cartesian coordinates to Earth Range Solar 
coordinates.  (See ERS2Cart for details)  Because the "B" angle 
is handy, we also allow conversion from HC->ERS with a special,
magic extra parameter.

=cut

sub SCart2ERS {
    print "SCart2ERS..." if($V4::debug);
    my($in,$opt,$from_HC) = @_;
    my($foo);

    my($i2);
    if($in->is_inplace) {
	$i2 = $in;
    } else {
	$i2 = $in->copy;
    }
    
# Calculate the Earth-Sun distance...
    $pb0rr = main::pb0r($in->x0);

# If we're going from HC rather than ERS, rotate by the B angle:
    if($from_HC) {
	my($b) = $pb0rr->slice('(1)');
	$i2->{PDL} x= rot_matrix(4,pdl(2,3),$b);
    }

# Translate the origin from the Sun to the Earth & switch to L.H. cart...

    (($foo = $i2->x3) -= $pb0rr->slice('(3)')) *= -1;

# Convert from Earth-centered L.H. cart to L.H. spherical
    $i2 = Cart2Sph($i2,{ix=>[0,3,2,1]});
    
    $i2->{coords} = "Earth Range Solar";

    return $i2;
}

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

=head2 SCart2HC

Convert from Solar Cartesian to Heliographic Cartesian coordinates.
This is a simple rotation by the "B" angle.  To save namespace, the
routine works both directions -- you feed it a multiplier on the 
angle.  To go from HC to SCart, you feed in a -1.

Because the conversion is time-dependent, if you feed in a 'nan' time,
you get all 'nan's back out.

=cut

sub SCart2HC {
    print "SCart2HC..."    if($V4::debug);

    my ($in,$opt,$mult) = @_;
    
    $mult = 1 unless($mult);
## Find the B angle...
    my($b) = $mult * 
      main::pb0r($in->x0)->slice('(1)');
## rotate
    if($in->is_inplace) {
	$out = $in;
    } else {
	$out = $in->copy;
    }
    
    $out->{PDL} x= rot_matrix(4,pdl(2,3),$b);

    if($mult == 1) {
	$out->{coords} = "HC";
    } else {
	$out->{coords} = "SC";
    }
    return bless($out,"V4");
}

1;