view release on metacpan or  search on metacpan

Astro/Coord.pm  view on Meta::CPAN

                     );
  @EXPORT_OK   = qw ( fk4fk5r fk5fk4r fk4galr galfk4r
                      ephem_vars nutate precsn $bepoch );
  @EXPORT_FAIL = qw ( );

  use Carp;
  use POSIX qw( asin acos fmod tan );

  use Astro::Time qw( $PI rad2turn turn2rad mjd2lst );
}

$bepoch = 1950.0;

use constant JULIAN_DAY_J2000 => 2451545.0;
use constant JULIAN_DAYS_IN_CENTURY => 36525.0;

# The E-terms vector for FK4 <--> other coordinate system transforms
# (used in fk4fk5 fk5fk4 fk4gal galfk4)
my @eterm = (-1.62557E-06, -0.31919E-06, -0.13843E-06);

## The precession matrix for FK4 <--> FK5 conversions (used in
## fk4fk5 and fk5fk4)
#my @btoj = ([+0.999925678186902,-0.011182059642247,-0.004857946558960],
#	    [+0.011182059571766,+0.999937478448132,-0.000027176441185],
#	    [+0.004857946721186,-0.000027147426498,+0.999988199738770]);

# The precession matrix for FK4 <--> Galactic conversions (used in
# fk4gal and galfk4)
my @etog = ([-0.066988739415,-0.872755765852,-0.483538914632],
	    [+0.492728466075,-0.450346958020,+0.744584633283],
	    [-0.867600811151,-0.188374601723,+0.460199784784]);

# Values used in SLALIB routines

use constant D2PI => 6.283185307179586476925287;

#  Radians per year to arcsec per century
use constant PMF => 100*60*60*360/D2PI;

#  Small number to avoid arithmetic problems
use constant TINY => 1e-30;

#  Km per sec to AU per tropical century
#  = 86400 * 36524.2198782 / 149597870
use constant  VF => 21.095;

#  Vectors A and Adot, and matrix M
my @a = ( -1.62557e-6,  -0.31919e-6, -0.13843e-6,
	  +1.245e-3, -1.580e-3, -0.659e-3);

my @ad =(+1.245e-3,    -1.580e-3,   -0.659e-3);

my @em = ([+0.9999256782, -0.0111820611, -0.0048579477],
	  [+0.0111820610, +0.9999374784, -0.0000271765],
	  [+0.0048579479, -0.0000271474, +0.9999881997],
	  [-0.000551,	    -0.238565,     +0.435739],
	  [+0.238514,     -0.002667,     -0.008541],
	  [-0.435623,     +0.012254,     +0.002117]);

my @emi = ([+0.9999256795, +0.0111814828, +0.0048590039,
	    -0.00000242389840, -0.00000002710544, -0.00000001177742],
	   [-0.0111814828, +0.9999374849, -0.0000271771,
	    +0.00000002710544, -0.00000242392702, +0.00000000006585],
	   [-0.0048590040, -0.0000271557, +0.9999881946,
	    +0.00000001177742, +0.00000000006585, -0.00000242404995],
	   [-0.000551,     +0.238509,     -0.435614,
	    +0.99990432,       +0.01118145,       +0.00485852],
	   [-0.238560,     -0.002667,     +0.012254,
	    -0.01118145,       +0.99991613,       -0.00002717],
	   [+0.435730,     -0.008541,     +0.002117,
	    -0.00485852,       -0.00002716,       +0.99996684]);

=item B<pol2r>

  ($x, $y, $z) = pol2r($polar1, $polar2);

 Converts a position in polar coordinates into rectangular coordinates
   $polar1, $polar2   The polar coordinates to convert (turns)
   $x, $y, $z         The rectangular coordinates

=cut

sub pol2r ($$) {
  my ($p1, $p2) = @_;

  # Converts polar coordinates into rectangluar
  my @rect;
  $rect[0] = cos(turn2rad($p1))*cos(turn2rad($p2));
  $rect[1] = sin(turn2rad($p1))*cos(turn2rad($p2));
  $rect[2] = sin(turn2rad($p2));
  return(@rect);
}

=item B<r2pol>

  ($polar1, $polar2) = r2pol($x, $y, $z);

 Converts a position in rectangular coordinates into polar coordinates
   $x, $y, $y         The rectangular coordinates to convert
   $polar1, $polar2   The polar coordinates (turns);
 Returns undef if too few or too many arguments are passed.

=cut

sub r2pol (@) {
  # First check that we have 3 arguments
  if (scalar @_ != 3) {
    carp 'Inconsistent arguments';
    return undef ;
  }
  my ($x, $y, $z) = @_;

  # Converts rectangular coordinates to polar
  my ($tmp, $left, $right);
  $tmp = atan2($y, $x)/(2.0*$PI);

  if (ref($tmp) =~ /PDL/ ) {  # Allow to work with PDL
    $tmp -> where($tmp<0.0) .= $tmp -> where($tmp<0.0)  + 1.0;
  } elsif ($tmp < 0.0) {
    $tmp += 1.0;
  }

  $left = $tmp;
  $tmp = sqrt($x*$x + $y*$y + $z*$z);