Astro
view release on metacpan or search on metacpan
Astro/Coord.pm view on Meta::CPAN
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);
Astro/Coord.pm view on Meta::CPAN
+0.000000063 * sin($arg13*2.0*$PI)
-0.000001102 * sin($arg31*2.0*$PI)
+0.000000345 * sin($arg32*2.0*$PI)
-0.000000187 * sin($arg33*2.0*$PI)
-0.000000146 * sin($arg34*2.0*$PI)
-0.000000077 * sin($arg35*2.0*$PI)
+0.000000060 * sin($arg36*2.0*$PI))/(2.0*$PI);
my $deps = ( 0.000044615 * cos($arg1*2.0*$PI)
-0.000000434 * cos($arg2*2.0*$PI)
+0.000002781 * cos($arg9*2.0*$PI)
+0.000000109 * cos($arg11*2.0*$PI)
+0.000000474 * cos($arg31*2.0*$PI)
+0.000000097 * cos($arg33*2.0*$PI)
+0.000000063 * cos($arg34*2.0*$PI))/(2.0*$PI);
my $eps = $eps0 + $deps;
my @N = ([1.0, -($dpsi)*(2.0*$PI)*cos($eps*2.0*$PI),
-($dpsi)*(2.0*$PI)*sin($eps*2.0*$PI)],
[0.0, 1.0, -($deps)*(2.0*$PI)],
[0.0, ($deps)*(2.0*$PI), 1.0]);
$N[1][0] = -1.0*$N[0][1];
$N[2][0] = -1.0*$N[0][2];
return($deps, $dpsi, @N);
}
=item B<precsn>
@gp = precsn($jd_start, $jd_stop);
To calculate the precession matrix P for dates AFTER 1984.0 (JD =
2445700.5) Given the position of an object referred to the equator
and equinox of the epoch $jd_start its position referred to the
equator and equinox of epoch $jd_stop can be calculated as follows :
1) Express the position as a direction cosine 3-vector (V1)
(use pol2r to do this).
2) The corresponding vector V2 for epoch jd_end is V2 = P.V1
The required inputs are :
$jd_start - The Julian day of the current epoch of the coordinates.
$jd_end - The Julian day at the required epoch for the conversion.
The returned values are :
@gp - The required precession matrix (array [3][3])
=cut
sub precsn ($$) {
my ($jd_start, $jd_end) = @_;
my @a = (0.011180860865024,
0.000006770713945,
-0.000000000673891,
0.000001463555541,
-0.000000001667759,
0.000000087256766);
my @b = (0.011180860865024,
0.000006770713945,
-0.000000000673891,
0.000005307158404,
0.000000000319977,
0.000000088250634);
my @d = (0.009717173455170,
-0.000004136915141,
-0.000000001052046,
0.000002068457570,
0.000000001052046,
-0.000000202812107);
my $t = ($jd_start - JULIAN_DAY_J2000)/JULIAN_DAYS_IN_CENTURY;
my $st = ($jd_end - $jd_start)/JULIAN_DAYS_IN_CENTURY;
my $t2 = $t * $t;
my $st2 = $st * $st;
my $st3 = $st2 * $st;
# Calculate the Equatorial precession parameters
# (ref. USNO Circular no. 163 1981,
# Lieske et al., Astron. & Astrophys., 58, 1 1977)
my $zeta = ($a[0] + $a[1]*$t + $a[2]*$t2) * $st +
($a[3] + $a[4]*$t) * $st2 + $a[5] * $st3;
my $z = ($b[0] + $b[1]*$t + $b[2]*$t2) * $st +
($b[3] + $b[4]*$t) * $st2 + $b[5] * $st3;
my $theta = ($d[0] + $d[1]*$t + $d[2]*$t2) * $st -
($d[3] + $d[4]*$t) * $st2 + $d[5] * $st3;
# Calculate the P matrix
my @precession = ([0.0, 0.0, 0.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0]);
$precession[0][0] = cos($zeta)*cos($z)*cos($theta) - sin($zeta)*sin($z);
$precession[0][1] = -sin($zeta)*cos($z)*cos($theta) - cos($zeta)*sin($z);
$precession[0][2] = -cos($z)*sin($theta);
$precession[1][0] = cos($zeta)*sin($z)*cos($theta) + sin($zeta)*cos($z);
$precession[1][1] = -sin($zeta)*sin($z)*cos($theta) + cos($zeta)*cos($z);
$precession[1][2] = -sin($z)*sin($theta);
$precession[2][0] = cos($zeta)*sin($theta);
$precession[2][1] = -sin($zeta)*sin($theta);
$precession[2][2] = cos($theta);
return(@precession);
}
=item B<coord_convert>
($output_left, $output_right) = coord_convert($input_left, $input_right,
$input_mode, $output_mode,
$mjd, $longitude, $latitude,
$ref0);
A routine for converting between any of the following coordinate systems :
Coordinate system input/output mode
----------------- -----------------
X, Y (East-West mounted) 0
Azimuth, Elevation 1
Hour Angle, Declination 2
Right Ascension, Declination (date, J2000 or B1950) 3,4,5
Galactic (B1950) 6
The last four parameters in the call ($mjd, $longitude, $latitude
( run in 1.718 second using v1.01-cache-2.11-cpan-5837b0d9d2c )