package Astro::Coord;
use strict;
=head1 NAME
Astro::Coord - Astronomical coordinate transformations
=head1 SYNOPSIS
use Astro::Coord;
($l, $b) = fk4gal($ra, $dec);
($az, $el) = eqazel($ha, $dec, $latitude);
=head1 DESCRIPTION
Astro::Coord contains an assorted set Perl routines for coordinate
conversions, such as hour angle to elevation and J2000 to B1950.
=head1 AUTHOR
Chris Phillips Chris.Phillips@csiro.au
=head1 FUNCTIONS
=cut
BEGIN {
use Exporter ();
use vars qw( $VERSION @ISA @EXPORT @EXPORT_OK @EXPORT_FAIL
$bepoch );
$VERSION = '1.43';
@ISA = qw(Exporter);
@EXPORT = qw( xy2azel azel2xy eqazel J2000todate
fk4fk5 fk5fk4 fk4gal galfk4 j2gal
coord_convert
haset_ewxy ewxy_tlos haset_azel azel_tlos
antenna_rise pol2r r2pol
);
@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);
if (ref($tmp) =~ /PDL/) { # Allow to work with PDL
$right = &PDL::Math::asin($z/$tmp)/(2.0*$PI);
} else {
$right = asin($z/$tmp)/(2.0*$PI);
}
return ($left, $right);
}
=item B<xy2azel>
($az, $el) = xy2azel($x, $y);
Converts a telescope position in X,Y coordinates into Az,El coordinates
$x, $y The X and Y coordinates (turns)
$az, $el The azimuth and elevation (turns)
=cut
sub xy2azel ($$) {
my ($x, $y) = @_;
# Convert a position in X,Y to Az,El
my @polar = pol2r($x, $y);
my $temp = $polar[0];
$polar[0] = $polar[1];
$polar[1] = $polar[2];
$polar[2] = $temp;
return (r2pol(@polar));
}
=item B<azel2xy>
($x, $y) = azel2xy($az, $el);
Converts a position in Az,El coordinates into X,Y coordinates
$az, $el The azimuth and elevation (turns)
$x, $y The X and Y coordinate (turns)
=cut
sub azel2xy ($$) {
my ($az, $el) = @_;
# Convert a position in Az,El to X,Y
my @polar = pol2r($az, $el);
my $temp = $polar[1];
$polar[1] = $polar[0];
$polar[0] = $polar[2];
$polar[2] = $temp;
my ($x, $y) = r2pol(@polar);
if ($x>0.5) {
$x -= 1.0;
}
if ($y>0.5) {
$y -= 1.0;
}
return ($x, $y);
}
=item B<eqazel>
($ha, $dec) = eqazel($az, $el, $latitude);
($az, $el) = eqazel($ha, $dec, $latitude);
($ha, $dec) = eqazel($az, $el, $latitude, $allownegative);
Converts HA/Dec coordinates to Az/El and vice versa.
$ha, $dec Hour angle and declination of source (turns)
$az, $el Azimuth and elevation of source (turns)
$latitude Latitude of the observatory (turns)
$allownegative If true, allow negative $ha or $az on return (Optional)
Note:
The ha,dec and az,el conversion is symmetrical
=cut
sub eqazel ($$$;$) {
my $sphi = sin(turn2rad($_[2]));
my $cphi = cos(turn2rad($_[2]));
my $sleft = sin(turn2rad($_[0]));
my $cleft = cos(turn2rad($_[0]));
my $sright = sin(turn2rad($_[1]));
my $cright = cos(turn2rad($_[1]));
my $left_out = atan2(-$sleft,-$cleft*$sphi+$sright*$cphi/$cright)/(2.0*$PI);
$left_out = ($left_out < 0.0) ? $left_out + 1.0 : $left_out
if (!(defined $_[3] && $_[3]));
my $right_out= asin($cleft*$cright*$cphi + $sright*$sphi)/(2.0*$PI);
return($left_out, $right_out);
}
=item B<fk4fk5>
($JRA, $JDec) = fk4fk5($BRA, $BDec);
(@fk5) = fk4fk5(@fk4);
Converts an FK4 (B1950) position to the equivalent FK5 (J2000)
position.
$BRA,$BDec fk4/B1950 position (turns)
$JRA,$Dec fk5/J2000 position (turns)
@fk4 fk4/B1950 position (as a 3-vector)
@fk5 fk5/J2000 position (as a 3-vector)
Note:
This code is based on similar routines from the Fortran SLALIB
package, so are quite accurate, but subject to a restrictive
license (see README).
=cut
sub fk4fk5 (@) {
# - - - - - -
# F K 4 5 Z
# - - - - - -
#
# Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
# proper motion in the FK5 frame (double precision)
#
# This routine converts stars from the old, Bessel-Newcomb, FK4
# system to the new, IAU 1976, FK5, Fricke system, in such a
# way that the FK5 proper motion is zero. Because such a star
# has, in general, a non-zero proper motion in the FK4 system,
# the routine requires the epoch at which the position in the
# FK4 system was determined.
#
# The method is from Appendix 2 of Ref 1, but using the constants
# of Ref 4.
#
# Given:
# R1950,D1950 dp B1950.0 FK4 RA,Dec at epoch (rad)
# BEPOCH dp Besselian epoch (e.g. 1979.3D0)
#
# Returned:
# R2000,D2000 dp J2000.0 FK5 RA,Dec (rad)
#
# Notes:
#
# 1) The epoch BEPOCH is strictly speaking Besselian, but
# if a Julian epoch is supplied the result will be
# affected only to a negligible extent.
#
# 2) Conversion from Besselian epoch 1950.0 to Julian epoch
# 2000.0 only is provided for. Conversions involving other
# epochs will require use of the appropriate precession,
# proper motion, and E-terms routines before and/or
# after FK45Z is called.
#
# 3) In the FK4 catalogue the proper motions of stars within
# 10 degrees of the poles do not embody the differential
# E-term effect and should, strictly speaking, be handled
# in a different manner from stars outside these regions.
# However, given the general lack of homogeneity of the star
# data available for routine astrometry, the difficulties of
# handling positions that may have been determined from
# astrometric fields spanning the polar and non-polar regions,
# the likelihood that the differential E-terms effect was not
# taken into account when allowing for proper motion in past
# astrometry, and the undesirability of a discontinuity in
# the algorithm, the decision has been made in this routine to
# include the effect of differential E-terms on the proper
# motions for all stars, whether polar or not. At epoch 2000,
# and measuring on the sky rather than in terms of dRA, the
# errors resulting from this simplification are less than
# 1 milliarcsecond in position and 1 milliarcsecond per
# century in proper motion.
#
# References:
#
# 1 Aoki,S., et al, 1983. Astron.Astrophys., 128, 263.
#
# 2 Smith, C.A. et al, 1989. "The transformation of astrometric
# catalog systems to the equinox J2000.0". Astron.J. 97, 265.
#
# 3 Yallop, B.D. et al, 1989. "Transformation of mean star places
# from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
# Astron.J. 97, 274.
#
# 4 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
# the Astronomical Almanac", ISBN 0-935702-68-7.
#
# Called: sla_DCS2C, sla_EPJ, sla_EPB2D, sla_DCC2S, sla_DRANRM
#
# P.T.Wallace Starlink 21 September 1998
#
# Copyright (C) 1998 Rutherford Appleton Laboratory
my ($rect, $w, $i, $j);
my (@r0, @a1, @v1, @v2); # Position and position+velocity vectors
if (@_==3) { # Rectangular coordinates passed
@r0 = @_;
$rect = 1;
} elsif (@_==2) { # Sperical coordinates
@r0 = pol2r($_[0],$_[1]); # Spherical to Cartesian
$rect = 0;
} elsif (@_>3) {
croak "Too many arguments for Astro::fk4fk5 ";
} else {
croak "Not enough arguments for Astro::fk4fk5 ";
}
# Adjust vector A to give zero proper motion in FK5
$w=($bepoch-1950)/PMF;
for ($i=0; $i<3; $i++) {
$a1[$i]=$a[$i]+$w*$ad[$i];
}
# Remove e-terms
$w=$r0[0]*$a1[0]+$r0[1]*$a1[1]+$r0[2]*$a1[2];
for ($i=0; $i<3; $i++) {
$v1[$i]=$r0[$i]-$a1[$i]+$w*$r0[$i];
}
# Convert position vector to Fricke system
for ($i=0; $i<6; $i++) {
$w=0;
for ($j=0; $j<3; $j++) {
#warn "DEBUG: [$i,$j]\n";
$w=$w+$em[$i][$j]*$v1[$j];
$v2[$i]=$w
}
}
# Allow for fictitious proper motion in FK4
$w=(epj(epb2d($bepoch))-2000)/PMF;
for ($i=0; $i<3; $i++) {
$v2[$i]=$v2[$i]+$w*$v2[$i+3];
}
if ($rect) {
return @v2[0..2];
} else {
# Revert to spherical coordinates
return r2pol(@v2[0..2]);
}
}
=item B<fk4fk5r>
@fk5 = fk4fk5r(@fk4);
Converts an FK4 (B1950) position to the equivalent FK5 (J2000) position.
Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
@fk4 fk4 position (as a 3-vector, turns)
@fk5 fk5 position (as a 3-vector, turns)
Note:
Just a wrapper to fk4fk5 which now handler polar and rectangular
arguments
=cut
sub fk4fk5r (@) {
return fk4fk5(@_);
}
#sub fk4fk5r (@) {
# # First check that we have 3 arguments
# if (scalar @_ < 3) {
# croak 'Not enough arguments for Astro::Coord::fk4fk5r at ';
# } elsif (scalar @_ > 3) {
# croak 'Too many arguments for Astro::Coord::fk4fk5r at ';
# }
#
# my ($i, $j, @temp, @fk5);
# my $w = 0.0;
#
# # Add the eterms
# for ($i=0 ; $i<3 ; $i++) {
# $w += $_[$i] * $eterm[$i];
# }
# for ($i=0 ; $i<3 ; $i++) {
# $temp[$i] = $_[$i] - $eterm[$i] + $w * $_[$i];
# }
#
# # Precess from FK4 to FK5
# for ($i=0 ; $i<3 ; $i++) {
# $fk5[$i] = 0.0;
# for ($j=0 ; $j<3 ; $j++) {
# $fk5[$i] += $btoj[$i][$j] * $temp[$j];
# }
# }
#
# return(@fk5);
#}
=item B<fk5fk4>
($JRA, $JDec) = fk4fk5($BRA, $BDec);
($@fk5) = fk4fk5(@fk4);
Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
$JRA,$Dec fk5/J2000 position (turns)
$BRA,$BDec fk4/B1950 position (turns)
@fk5 fk5/J2000 position (as a 3-vector)
@fk4 fk4/B1950 position (as a 3-vector)
Note:
This code is based on similar routines from the Fortran SLALIB
package, so are quite accurate, but subject to a restrictive
license (see README).
=cut
sub fk5fk4 (@) {
#+
# - - - - - -
# F K 5 2 4
# - - - - - -
#
# Convert J2000.0 FK5 star data to B1950.0 FK4 (double precision)
#
# This routine converts stars from the new, IAU 1976, FK5, Fricke
# system, to the old, Bessel-Newcomb, FK4 system. The precepts
# of Smith et al (Ref 1) are followed, using the implementation
# by Yallop et al (Ref 2) of a matrix method due to Standish.
# Kinoshita's development of Andoyer's post-Newcomb precession is
# used. The numerical constants from Seidelmann et al (Ref 3) are
# used canonically.
#
# Given: (all J2000.0,FK5)
# R2000,D2000 dp J2000.0 RA,Dec (rad)
# DR2000,DD2000 dp J2000.0 proper motions (rad/Jul.yr)
# P2000 dp parallax (arcsec)
# V2000 dp radial velocity (km/s, +ve = moving away)
#
# Returned: (all B1950.0,FK4)
# R1950,D1950 dp B1950.0 RA,Dec (rad)
# DR1950,DD1950 dp B1950.0 proper motions (rad/trop.yr)
# P1950 dp parallax (arcsec)
# V1950 dp radial velocity (km/s, +ve = moving away)
#
# Notes:
#
# 1) The proper motions in RA are dRA/dt rather than
# cos(Dec)#dRA/dt, and are per year rather than per century.
#
# 2) Note that conversion from Julian epoch 2000.0 to Besselian
# epoch 1950.0 only is provided for. Conversions involving
# other epochs will require use of the appropriate precession,
# proper motion, and E-terms routines before and/or after
# FK524 is called.
#
# 3) In the FK4 catalogue the proper motions of stars within
# 10 degrees of the poles do not embody the differential
# E-term effect and should, strictly speaking, be handled
# in a different manner from stars outside these regions.
# However, given the general lack of homogeneity of the star
# data available for routine astrometry, the difficulties of
# handling positions that may have been determined from
# astrometric fields spanning the polar and non-polar regions,
# the likelihood that the differential E-terms effect was not
# taken into account when allowing for proper motion in past
# astrometry, and the undesirability of a discontinuity in
# the algorithm, the decision has been made in this routine to
# include the effect of differential E-terms on the proper
# motions for all stars, whether polar or not. At epoch 2000,
# and measuring on the sky rather than in terms of dRA, the
# errors resulting from this simplification are less than
# 1 milliarcsecond in position and 1 milliarcsecond per
# century in proper motion.
#
# References:
#
# 1 Smith, C.A. et al, 1989. "The transformation of astrometric
# catalog systems to the equinox J2000.0". Astron.J. 97, 265.
#
# 2 Yallop, B.D. et al, 1989. "Transformation of mean star places
# from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
# Astron.J. 97, 274.
#
# 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
# the Astronomical Almanac", ISBN 0-935702-68-7.
#
# P.T.Wallace Starlink 19 December 1993
#
# Copyright (C) 1995 Rutherford Appleton Laboratory
#-
my ($rect, @v1, @v2);
if (@_==3) { # Rectangular coordinates passed
@v1 = @_;
$rect = 1;
} elsif (@_==2) { # Sperical coordinates
@v1 = pol2r($_[0],$_[1]); # Spherical to Cartesian
$rect = 0;
} elsif (@_>2) {
croak "Too many arguments for Astro::fk5fk4 ";
} else {
croak "Not enough arguments for Astro::fk5fk4 ";
}
# Miscellaneous
my ($w, $x, $y, $z, $wd, $rxyz);
my ($ur, $ud, $xd, $yd, $zd);
my ($i,$j);
# Convert position+velocity vector to BN system
for ($i=0; $i<6; $i++) {
$w=0.0;
##for ($j=0; $j<6; $j++) {
for ($j=0; $j<3; $j++) {
$w=$w+$emi[$i][$j]*$v1[$j];
}
$v2[$i]=$w;
}
# Position vector components and magnitude
$x=$v2[0];
$y=$v2[1];
$z=$v2[2];
$rxyz=sqrt($x*$x+$y*$y+$z*$z);
# Apply E-terms to position
$w=$x*$a[0]+$y*$a[1]+$z*$a[2];
$x=$x+$a[0]*$rxyz-$w*$x;
$y=$y+$a[1]*$rxyz-$w*$y;
$z=$z+$a[2]*$rxyz-$w*$z;
# Recompute magnitude
$rxyz=sqrt($x*$x+$y*$y+$z*$z);
# Apply E-terms to both position and velocity
$x=$v2[0];
$y=$v2[1];
$z=$v2[2];
$w=$x*$a[0]+$y*$a[1]+$z*$a[2];
$wd=$x*$a[3]+$y*$a[4]+$z*$a[5];
$x=$x+$a[0]*$rxyz-$w*$x;
$y=$y+$a[1]*$rxyz-$w*$y;
$z=$z+$a[2]*$rxyz-$w*$z;
$xd=$v2[3]+$a[3]*$rxyz-$wd*$x;
$yd=$v2[4]+$a[4]*$rxyz-$wd*$y;
$zd=$v2[5]+$a[5]*$rxyz-$wd*$z;
my @r;
if ($rect) {
@r = ($x, $y, $z);
} else {
@r= r2pol($x, $y, $z);
}
# my $rxysq =$x*$x+$y*$y;
# my $rxy = sqrt($rxysq);
# if ($rxy>TINY) {
# $ur=($x*$yd-$y*$xd)/$rxysq;
# $ud=($zd*$rxysq-$z*($x*$xd+$y*$yd))/(($rxysq+$z*$z)*$rxy);
# }
#
## Return results
# my $dr1950=$ur/PMF;
# my $dd1950=$ud/PMF;
return(@r);
}
=item B<fk5fk4r>
@fk4 = fk5fk4r(@fk5);
Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
@fk4 fk4 position (as a 3-vector, turns)
@fk5 fk5 position (as a 3-vector, turns)
Note:
Just a wrapper to fk5fk4 which now handler polar and rectangular
arguments
=cut
sub fk5fk4r (@) {
return fk5fk4(@_);
}
#sub fk5fk4 (@) {
## - - - - - -
## F K 5 4 Z
## - - - - - -
##
## Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming
## zero proper motion and parallax (double precision)
##
## This routine converts star positions from the new, IAU 1976,
## FK5, Fricke system to the old, Bessel-Newcomb, FK4 system.
##
## Given:
## R2000,D2000 dp J2000.0 FK5 RA,Dec (rad)
## BEPOCH dp Besselian epoch (e.g. 1950D0)
##
## Returned:
## R1950,D1950 dp B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH
## DR1950,DD1950 dp B1950.0 FK4 proper motions (rad/trop.yr)
##
## Notes:
##
## 1) The proper motion in RA is dRA/dt rather than cos(Dec)#dRA/dt.
##
## 2) Conversion from Julian epoch 2000.0 to Besselian epoch 1950.0
## only is provided for. Conversions involving other epochs will
## require use of the appropriate precession routines before and
## after this routine is called.
##
## 3) Unlike in the sla_FK524 routine, the FK5 proper motions, the
## parallax and the radial velocity are presumed zero.
##
## 4) It is the intention that FK5 should be a close approximation
## to an inertial frame, so that distant objects have zero proper
## motion; such objects have (in general) non-zero proper motion
## in FK4, and this routine returns those fictitious proper
## motions.
##
## 5) The position returned by this routine is in the B1950
## reference frame but at Besselian epoch BEPOCH. For
## comparison with catalogues the BEPOCH argument will
## frequently be 1950D0.
##
## Called: sla_FK524, sla_PM
##
## P.T.Wallace Starlink 10 April 1990
##
## Copyright (C) 1995 Rutherford Appleton Laboratory
#
# my $bepoch = 1950.0;
#
# my $rect;
# if (@_>3) {
# croak "Too many arguments for Astro::fk5fk4 ";
# } elsif (@_<2) {
# croak "Not enough arguments for Astro::fk5fk4 ";
# }
# my @r2000 = @_;
#
# # fk5 equinox j2000 (any epoch) to fk4 equinox b1950 epoch b1950
# my (@r1950) = fk524(@r2000);
# my $dd1950 = pop @r1950;
# my $dr1950 = pop @r1950;
#
# ## fictitious proper motion to epoch bepoch
# #my ($r1950, $d1950) = pm($r,$d,$dr1950,$dd1950,0.0,0.0,1950,$bepoch);
# return @r1950;
#}
#=item B<fk5fk4r>
#
# @fk4 = fk5fk4r(@fk5);
#
# Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
# Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
# @fk5 fk5 position (as a 3-vector, turns)
# @fk4 fk4 position (as a 3-vector, turns)
#
#=cut
#
#sub fk5fk4r(@) {
#
# # First check that we have 3 arguments
# if (scalar @_ < 3) {
# croak 'Not enough arguments for Astro::Coord::fk5fk4r at ';
# } elsif (scalar @_ > 3) {
# croak 'Too many arguments for Astro::Coord::fk5fk4r at ';
# }
#
# my ($i, $j, @fk4);
# my $w = 0.0;
#
# # Precess. Note : the same matrix is used as for the FK4 -> FK5
# # transformation, but we have transposed it within the
# # for loop
#
# for ($i=0 ; $i<3 ; $i++) {
# $fk4[$i] = 0.0;
# for ($j=0 ; $j<3 ; $j++) {
# $fk4[$i] += $btoj[$j][$i] * $_[$j];
# }
# }
#
# # Allow for e-terms
# for ($i=0 ; $i<3 ; $i++) {
# $w += $_[$i] * $eterm[$i];
# }
# $w += 1.0;
# for ($i=0 ; $i<3 ; $i++) {
# $fk4[$i] = ($fk4[$i] + $eterm[$i])/$w;
# }
#
# return(@fk4);
#}
=item B<fk4galr>
@gal = fk4galr(@fk4)
Converts an FK4 position (B1950.0) to the IAU 1958 Galactic
coordinate system
Note: convert equitoral positions to/from 3-vectors using pol2r and r2pol.
@fk4 fk4 position to convert (as a 3-vector, turns)
@gal Galactic position (as a 3-vector, turns)
Returns undef if too few or two many arguments are passed.
Reference : Blaauw et al., 1960, MNRAS, 121, 123.
=cut
# Within 1e-7 arcsec of SLALIB slaEg50
sub fk4galr(@) {
# First check that we have 3 arguments
if (scalar @_ < 3) {
croak 'Not enough arguments for Astro::Coord::fk4galr at ';
} elsif (scalar @_ > 3) {
croak 'Too many arguments for Astro::Coord::fk4galr at ';
}
my ($i, $j, @temp, @gal);
my $w = 0.0;
# Allow for e-terms
for ($i=0 ; $i<3 ; $i++) {
$w += $_[$i] * $eterm[$i];
}
for ($i=0 ; $i<3 ; $i++) {
$temp[$i] = $_[$i] - $eterm[$i] + $w * $_[$i];
}
# Precess
for ($i=0 ; $i<3 ; $i++) {
$gal[$i] = 0.0;
for ($j=0 ; $j<3 ; $j++) {
$gal[$i] += $etog[$i][$j] * $temp[$j];
}
}
return(@gal);
}
=item B<galfk4>
($bRA, $bDec) = galfk4($l, $b);
@fk4 = galfk4(@gal);
Converts an IAU 1958 Galactic position to the FK4 coordinate system (B1950)
Notes: Converts equitoral positions to/from 3-vectors using pol2r and r2pol.
$BRA,$BDec fk4/B1950 position (turns)
$l, $b Galactic longitude and latitude
@gal Galactic position (as a 3-vector, turns)
@fk4 fk4 position (as a 3-vector, turns)
Reference : Blaauw et al., 1960, MNRAS, 121, 123.
=cut
# Within 1e-7 arcsec of SLALIB slaGe50
sub galfk4(@) {
my (@r, $rect);
if (@_==3) { # Rectangular coordinates passed
@r = @_;
$rect = 1;
} elsif (@_==2) { # Sperical coordinates
@r = pol2r($_[0],$_[1]); # Spherical to Cartesian
$rect = 0;
} elsif (@_>3) {
croak "Too many arguments for Astro::galfk4 at";
} else {
croak "Not enough arguments for Astro::galfk4 at";
}
my ($i, $j, @fk4);
my $w = 0.0;
# Precess. Note : the same matrix is used as for the FK4 -> Galactic
# transformation, but we have transposed it within the
# for loop
for ($i=0 ; $i<3 ; $i++) {
$fk4[$i] = 0.0;
for ($j=0 ; $j<3 ; $j++) {
$fk4[$i] += $etog[$j][$i] * $r[$j];
}
}
# Allow for e-terms */
for ($i=0 ; $i<3 ; $i++) {
$w += $r[$i] * $eterm[$i];
}
$w += 1.0;
for ($i=0 ; $i<3 ; $i++) {
$fk4[$i] = ($fk4[$i] + $eterm[$i])/$w;
}
if ($rect) {
return @fk4;
} else {
return r2pol(@fk4);
}
}
sub galfk4r(@) {galfk4(@_)};
#=item B<fk4fk5>
#
# ($JRA, $JDec) = fk4fk5($BRA, $BDec);
#
# Converts an FK4 (B1950) position to the equivalent FK5 (J2000) position.
# **LOW PRECISION**
# $BRA,$BDec fk4/B1950 position (turns)
# $JRA,$Dec fk5/J2000 position (turns)
#
#=cut
#
#sub fk4fk5 ($$) {
# return r2pol(fk4fk5r(pol2r(shift,shift)));
#}
#=item B<fk5fk4>
#
# ($BRA, $BDec) = fk5fk4($JRA, $JDec);
#
# Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
# **LOW PRECISION**
# $JRA,$Dec fk5/J2000 position (turns)
# $BRA,$BDec fk4/B1950 position (turns)
#
#=cut
#
#sub fk5fk4 ($$) {
# return r2pol(fk5fk4r(pol2r(shift,shift)));
#}
=item B<fk4gal>
($l, $b) = fk4gal($ra, $dec);
Converts an FK4 position (B1950) to the IAU 1958 Galactic
coordinate system
($ra, $dec) fk4 position to convert (turns)
($l, $b) Galactic position (turns)
Reference : Blaauw et al., 1960, MNRAS, 121, 123.
=cut
sub fk4gal ($$) {
return r2pol(fk4galr(pol2r(shift,shift)));
}
#=item B<galfk4>
#
# ($ra, $dec) = galfk4($l, $b);
#
# Converts an IAU 1958 Galactic coordinate system position
# to FK4 (B1950).
# ($l, $b) Galactic position (turns)
# ($ra, $dec) fk4 position to convert (turns)
# Reference : Blaauw et al., 1960, MNRAS, 121, 123.
#
#=cut
#
#sub galfk4 ($$) {
# return r2pol(galfk4r(pol2r(shift,shift)));
#}
=item B<ephem_vars>
($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($jd)
Given the Julian day ($jd) this routine calculates the ephemeris
values required by the prcmat and nutate routines
The returned values are :
$omega - Longitude of the ascending node of the Moons mean orbit on
the ecliptic, measured from the mean equinox of date.
$rma - Mean anomaly of the Sun.
$mlanom - Mean anomaly of the Moon.
$F - L - omega, where L is the mean longitude of the Moon.
$D - Mean elongation of the Moon from the Sun.
$eps0 - Mean obilquity of the ecliptic.
=cut
=item B<J2000todate>
($DRA, $DDec) = J2000todate($JRA, $JDec, $mjd);
@date = J2000todate(@J2000, $mjd);
Converts an J2000 position date coordinate
$DRA,$DDec Date coordinate (turns)
$JRA,$Dec J2000 position (turns)
@date Date coordinate (as a 3-vector)
@J2000 J2000 position (as a 3-vector)
=cut
# Untested
sub J2000todate(@) {
my ($rect);
my (@J2000, @date); # Position vectors
my $mjd = pop @_;
if (@_==3) { # Rectangular coordinates passed
@J2000 = @_;
$rect = 1;
} elsif (@_==2) { # Sperical coordinates
@J2000 = pol2r($_[0],$_[1]); # Spherical to Cartesian
$rect = 0;
} elsif (@_>3) {
croak "Too many arguments for Astro::Coord::J2000todate ";
} else {
croak "Not enough arguments for Astro::Coord::J2000todate ";
}
# compute the general precession matrix.
my @gp = precsn(JULIAN_DAY_J2000, $mjd+2400000.5);
# Determine ephemeris quantities
my ($deps, $dpsi);
my @nu = ();
my ($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($mjd+2400000.5);
($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
my @prcmat = ();
for (my $i=0 ; $i<3 ; $i++) {
for (my $j=0 ; $j<3 ; $j++) {
my $xx = 0.0;
for (my $k=0 ; $k<3 ; $k++) {
$xx = $xx + $gp[$i][$k] * $nu[$k][$j];
}
$prcmat[$i][$j] = $xx;
}
}
for (my $i=0 ; $i<3 ; $i++) {
$date[$i] = 0.0;
for (my $j=0 ; $j<3 ; $j++) {
$date[$i] += $prcmat[$i][$j] * $J2000[$j];
}
}
if ($rect) {
return @date;
} else {
# Revert to spherical coordinates
return r2pol(@date);
}
}
sub ephem_vars ($) {
my $epoch = shift;
# Calculates values required internally by prcmat and for nutate from
# the passed Julian Day
# Calculate the interval to/from J2000 in Julian Centuries
my $jcents = ($epoch-(JULIAN_DAY_J2000))/JULIAN_DAYS_IN_CENTURY;
# Calculate the longitude of the mean ascending node of the lunar
# orbit on the ecliptic [A.A. Suppl. 1984, p S26]
my $omega = (((0.000000039 * $jcents + 0.000036143) *
$jcents - 33.757045934) *
$jcents + 2.182438624)/(2.0*$PI);
$omega = fmod($omega,1.0);
if ($omega < 0.0) {
$omega += 1.0;
}
# Calculate the mean anomaly. [A.A suppl. 1984, p S26]
my $manom = (6.240035939 - ((5.818e-8 * $jcents +2.797e-6) * $jcents -
628.301956024) * $jcents)/(2.0*$PI);
$manom = fmod($manom,1.0);
if ($manom < 0.0) {
$manom += 1.0;
}
# Calculate the mean anomaly of the Moon. [A.A. Suppl, 1984, p S26]
my $mlanom = (((0.000000310 * $jcents + 0.000151795) * $jcents
+8328.691422884) * $jcents + 2.355548394)/(2.0*$PI);
$mlanom = fmod($mlanom,1.0);
if ($mlanom < 0.0) {
$mlanom += 1.0;
}
# Calculate the longitude of the moon from ascending node.
# [A.A. Suppl, 1984, p S26]
my $F = (((0.000000053 * $jcents - 0.000064272) * $jcents + 8433.466158318)
* $jcents + 1.627901934)/(2.0*$PI);
$F = fmod($F,1.0);
if ($F < 0.0) {
$F += 1.0;
}
# Calculate the mean elongation of the moon from the sun.
# [A.A suppl. 1984, p S26]
my $D = (((0.000000092 * $jcents + 0.000033409) * $jcents + 7771.377146171)
* $jcents + 5.198469514)/(2.0*$PI);
$D = fmod($D,1.0);
if ($D < 0.0) {
$D += 1.0;
}
# Calculate the mean obliquity of the ecliptic = mean obliquity.
# [A.A suppl. 1984, p S26]
my $eps0 = (((0.000000009 * $jcents - 0.000000003) * $jcents - 0.000226966)
* $jcents + 0.409092804)/(2.0*$PI);
return($omega, $manom, $mlanom, $F, $D, $eps0)
}
=item B<nutate>
($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
To calculate the nutation in longitude and obliquity according to
the 1980 IAU Theory of Nutation including terms with amplitudes
greater than 0.01 arcsecond. The nutation matrix is used to
compute true place from mean place: true vector = N x mean place
vector where the three components of each vector are the direction
cosines wrt the mean equinox and equator.
/ 1 -dpsi.cos(eps) -dpsi.sin(eps) \
| |
N = | dpsi.cos(eps) 1 -deps |
| |
\ dpsi.sin(eps) deps 1 /
The required inputs are : (NOTE: these are the values returned by ephem_vars)
$omega - Longitude of the ascending node of the Moons mean orbit on
the ecliptic, measured from the mean equinox of date.
$rma - Mean anomaly of the Sun.
$mlanom - Mean anomaly of the Moon.
$F - L - omega, where L is the mean longitude of the Moon.
$D - Mean elongation of the Moon from the Sun.
$eps0 - Mean obilquity of the ecliptic.
The returned values are :
$deps - nutation in obliquity
$dpsi - nutation in longitude (scalar)
@nu - nutation matrix (array [3][3])
=cut
sub nutate ($$$$$$) {
my ($omega, $F, $D, $manom, $mlanom, $eps0) = @_;
my $arg1 = $omega;
my $arg2 = 2.0 * $omega;
my $arg9 = 2.0 * ($F-$D+$omega);
my $arg10 = $manom;
my $arg11 = $arg9 + $arg10;
my $arg12 = $arg9 - $arg10;
my $arg13 = $arg9 - $arg1;
my $arg31 = 2.0 * ($F+$omega);
my $arg32 = $mlanom;
my $arg33 = $arg31 - $arg1;
my $arg34 = $arg31 + $arg32;
my $arg35 = $mlanom - 2.0 * $D;
my $arg36 = $arg31 - $arg32;
my $dpsi = (-0.000083386 * sin($arg1*2.0*$PI)
+0.000001000 * sin($arg2*2.0*$PI)
-0.000006393 * sin($arg9*2.0*$PI)
+0.000000691 * sin($arg10*2.0*$PI)
-0.000000251 * sin($arg11*2.0*$PI)
+0.000000105 * sin($arg12*2.0*$PI)
+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
and $ref0) are not always required for the coordinate conversion.
In particular if the conversion is between two coordinate systems
which are fixed with respect to the celestial sphere (RA/Dec J2000,
B1950 or Galactic), or two coordinate systems which are fixed with
respect to the antenna (X/Y and Az/El) then these parameters are not
used (NOTE: they must always be passed, even if they only hold 0 or
undef as the routine will return undef if it is not passed 8
parameters). The RA/Dec date system is defined with respect to the
celestial sphere, but varies with time. The table below shows which
of $mjd, $longitude, $latitude and $ref0 are used for a given
conversion. If in doubt you should determing the correct values for
all input parameters, no checking is done in the routine that the
passed values are sensible.
Conversion $mjd $longitude $latitude $ref0
------------------------------------------------------------------------
Galactic, Galactic,
RA/Dec J2000,B1950 <->RA/Dec J2000, B1950 N N N N
Galactic,
RA/Dec J2000,B1950 <->RA/Dec date Y N N N
Galactic,
RA/Dec J2000,B1950,<->HA/Dec Y Y N N
date
Galactic,
RA/Dec J2000,B1950,<->X/Y, Az/El Y Y Y Y
date
X/Y, Az/El <->X/Y, Az/El N N N N
X/Y, Az/El <->HA/Dec N N Y Y
NOTE : The method used for refraction correction is asymptotic at
an elevation of 0 degrees.
The required inputs are :
$input_left - The left/longitude input coordinate (turns)
$input_right - The right/latitude input coordinate (turns)
$input_mode - The mode of the input coordinates (0-6)
$output_mode - The mode to convert the coordinates to.
$mjd - The time as modified Julian day (if necessary) at
which to perform the conversion
$longitude - The longitude of the location/observatory (if necessary)
at which to perform the conversion (turns)
$latitude - The latitude of the location/observatory (if necessary)
at which to perform the conversion (turns)
$ref0 - The refraction constant (if in doubt use 0.00005).
The returned values are :
$output_left - The left/longitude output coordinate (turns)
$output_right - The right/latitude output coordinate (turns)
=cut
sub coord_convert ($$$$;$$$$) {
my ($input_left, $input_right, $input_mode, $output_mode, $mjd, $longitude,
$latitude, $ref0) = @_;
# Some required constants
my ($EWXY, $AZEL, $HADEC, $DATE, $J2000, $B1950, $GALACTIC) = 0..6;
# First check what the input and output modes are.
if (($input_mode < $EWXY) || ($input_mode > $GALACTIC)) {
carp "Invalid input coordinate mode : $input_mode\n".
"Valid inputs are numbers in the range 0-6, which corrspond to X/Y, ".
"Az/El,\n HA/Dec, RA/Dec (date), RA/Dec (J2000), RA/Dec (B1950), ".
"Galactic.";
return undef;
}
if (($output_mode < $EWXY) || ($output_mode > $GALACTIC)) {
carp "Invalid output coordinate mode : $output_mode\n".
"Valid outputs are numbers in the range 0-6, which corrspond to X/Y, ".
"Az/El,\n HA/Dec, RA/Dec (date), RA/Dec (J2000), RA/Dec (B1950), ".
"Galactic.";
return undef;
}
# Check we have the correct parameters passed
# Need mjd
if ((($input_mode>=$DATE && $output_mode<=$DATE) ||
($input_mode<=$DATE && $output_mode>=$DATE)) &&
!(defined($mjd))) {
carp '$mjd parametr missing';
return undef;
}
# Need longitude
if ((($input_mode>=$HADEC && $output_mode<=$AZEL) ||
($input_mode<=$HADEC && $output_mode>=$HADEC)) &&
!(defined($longitude))) {
carp '$longitude parametr missing';
return undef;
}
# Need latitude
if ((($input_mode>=$HADEC && $output_mode<$HADEC) ||
($input_mode<=$AZEL && $output_mode>$AZEL)) &&
!(defined($latitude))) {
carp '$latitude parameter missing';
return undef;
}
# Need ref0
if ((($input_mode>=$HADEC && $output_mode<$HADEC) ||
($input_mode<=$AZEL && $output_mode>$AZEL)) &&
!(defined($ref0))) {
carp '$ref0 parameter missing';
return undef;
}
# If necessary determine ephemeris quantities (if either of the modes are
# date, HA/Dec, AzEl or EWXY).
my ($omega, $rma, $mlanom, $F, $D, $eps0, $deps, $dpsi);
my @nu = ();
if (($input_mode<=$DATE && $output_mode>=$DATE) ||
($input_mode>=$DATE && $output_mode<=$DATE)) {
($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($mjd+2400000.5);
($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
}
my @vonc = ();
if (($input_mode<=$HADEC && $output_mode>=$DATE) ||
($input_mode>=$DATE && $output_mode<=$HADEC)) {
# Calculate the interval to/from J2000 in Julian Centuries
my $jcents = ($mjd+2400000.5-(JULIAN_DAY_J2000))/JULIAN_DAYS_IN_CENTURY;
# Compute the eccentricity of the Earth's orbit (in radians)
# [Explanatory supplement to the Astronomical Ephemeris 1961, p 98]
my $e = (-0.000000126 * $jcents - 0.00004205) * $jcents + 0.016709114;
# Compute the eccentric anomaly, by iteratively solving :
# ea = e*sin(ea) - rma
my $ea = $rma;
my $xx;
do {
$xx = $ea;
$ea = $xx + ($rma - $xx + $e*sin($xx)) / (1.0 - $e*cos($xx));
} while (abs($ea -$xx) > 1.0e-9);
# Compute the mean longitude of perihelion, in radians
# (reference as for `e').
my $perihl = ((0.00000005817764*$jcents + 0.000008077) * $jcents
+ 0.030010190) * $jcents + 1.796613066;
# Calculate the equation of the equinoxes
#my $eqenx = $dpsi * cos(($eps0+$deps)*2.0*$PI);
# Compute the abberation vector
my $eps = $eps0 + $deps;
$xx = 0.00009936508 / (1.0 - $e*cos($ea));
my $efac = sqrt(1.0 - $e*$e);
$vonc[0] = $xx * (-cos($perihl)*sin($ea) - $efac*sin($perihl)*cos($ea));
$vonc[1] = $xx * (-sin($perihl)*cos($eps)*sin($ea) +
$efac*cos($perihl)*cos($eps)*cos($ea));
$vonc[2] = $xx * (-sin($perihl)*sin($eps)*sin($ea) +
$efac*cos($perihl)*sin($eps)*cos($ea));
}
my @prcmat = ();
if (($input_mode<=$DATE && $output_mode>=$J2000) ||
($input_mode>=$J2000 && $output_mode<=$DATE)) {
# compute the general precession matrix. */
my @gp = precsn(JULIAN_DAY_J2000, $mjd+2400000.5);
# The matrices returned from nutate (nu) and precsn (gp) can be used
# to convert J2000 coordinates to date by :
# (coords at date) = gp * nu * (coords at J2000)
# gp and nu can be combined to give the required precession matrix
for (my $i=0 ; $i<3 ; $i++) {
for (my $j=0 ; $j<3 ; $j++) {
my $xx = 0.0;
for (my $k=0 ; $k<3 ; $k++) {
$xx = $xx + $gp[$i][$k] * $nu[$k][$j];
}
$prcmat[$i][$j] = $xx;
}
}
}
my $lmst;
if (($input_mode<=$HADEC && $output_mode>=$DATE) ||
($output_mode<=$HADEC && $input_mode>=$DATE)) {
$lmst = mjd2lst($mjd, $longitude);
}
# Perform the conversion
my (@lb, @b1950, @j2000, @date, $ra, $ha, $dec, $az, $el, $x, $y);
if ($input_mode == $GALACTIC) {
@lb = pol2r($input_left, $input_right);
} elsif ($input_mode == $B1950) {
@b1950 = pol2r($input_left, $input_right);
} elsif ($input_mode == $J2000) {
@j2000 = pol2r($input_left, $input_right);
} elsif ($input_mode == $DATE) {
@date = pol2r($input_left, $input_right);
} elsif ($input_mode == $HADEC) {
$ha = $input_left;
$dec = $input_right;
} elsif ($input_mode == $AZEL) {
$az = $input_left;
$el = $input_right;
} else {
$x = $input_left;
$y = $input_right;
}
# Conversion is to a "lower" mode
if ($output_mode < $input_mode) {
# Convert from Galactic to B1950
if ($input_mode == $GALACTIC) {
@b1950 = galfk4r(@lb);
}
# Convert from B1950 to J2000
if (($input_mode >= $B1950) && ($output_mode < $B1950)) {
@j2000 = fk4fk5r(@b1950);
}
# Precess from J2000 to date
if (($input_mode >= $J2000) && ($output_mode < $J2000)) {
for (my $i=0 ; $i<3 ; $i++) {
$date[$i] = 0.0;
for (my $j=0 ; $j<3 ; $j++) {
$date[$i] += $prcmat[$i][$j] * $j2000[$j];
}
}
}
# Convert from date to HA/Dec
if (($input_mode >= $DATE) && ($output_mode < $DATE)) {
# Convert to geometrical equitorial coordinates
for (my $i=0 ; $i<3 ; $i++) {
$date[$i] += $vonc[$i];
}
# Convert from retangular back to polar coordinates
($ra, $dec) = r2pol(@date);
# Convert to hour angle
$ha = $lmst - $ra;
if ($ha < 0.0) {
$ha += 1.0;
}
}
# Convert from HA/Dec to Az/El
if (($input_mode >= $HADEC) && ($output_mode < $HADEC)) {
($az, $el) = eqazel($ha, $dec, $latitude);
# Correct for refraction
$el += $ref0/tan($el*2.0*$PI);
}
# Convert from Az/El to X/Y
if (($input_mode >= $AZEL) && ($output_mode < $AZEL)) {
($x, $y) = azel2xy($az, $el);
}
} else {
# Convert from X/Y to Az/El
if (($input_mode == $EWXY) && ($output_mode > $EWXY)) {
($az, $el) = xy2azel($x, $y);
}
# Convert from Az/El to HA/Dec
if (($input_mode <= $AZEL) && ($output_mode > $AZEL)) {
# First numerically invert the refraction correction
my $upper = $el - $ref0/tan($el*2.0*$PI);
my $lower = $el - 1.5*$ref0/tan($el*2.0*$PI);
my $root = ($lower+$upper)/2.0;
my $niter = 0;
do {
if ($root + $ref0/tan($root*2.0*$PI) - $el > 0.0) {
$upper = $root;
} else {
$lower = $root;
}
$root = ($lower+$upper)/2.0;
$niter++;
} while (($niter <= 10) && (($upper-$root) > 7.0e-8));
$el = $root;
# Now do the conversion
($ha, $dec) = eqazel($az, $el, $latitude);
}
# Convert from HA/Dec to date
if (($input_mode <= $HADEC) && ($output_mode > $HADEC)) {
$ra = $lmst - $ha;
if ($ra < 0.0) {
$ra += 1.0;
}
@date = pol2r($ra, $dec);
# Remove the abberation vector
for (my $i=0 ; $i<3 ; $i++) {
$date[$i] -= $vonc[$i];
}
}
# precess from date to J2000
if (($input_mode <= $DATE) && ($output_mode > $DATE)) {
for (my $i=0 ; $i<3 ; $i++) {
$j2000[$i] = 0.0;
for (my $j=0 ; $j<3 ; $j++) {
$j2000[$i] += $prcmat[$j][$i] * $date[$j];
}
}
}
# Convert from J2000 to B1950
if (($input_mode <= $J2000) && ($output_mode > $J2000)) {
@b1950 = fk5fk4(@j2000);
}
# Convert from B1950 to Galactic
if (($input_mode <= $B1950) && ($output_mode >= $B1950)) {
@lb = fk4galr(@b1950);
}
}
if ($output_mode == $EWXY) {
return($x, $y);
} elsif ($output_mode == $AZEL) {
return($az, $el);
} elsif ($output_mode == $HADEC) {
return($ha, $dec);
} elsif ($output_mode == $DATE) {
return(r2pol(@date));
} elsif ($output_mode == $J2000) {
return(r2pol(@j2000));
} elsif ($output_mode == $B1950) {
return(r2pol(@b1950));
} elsif ($output_mode == $GALACTIC) {
return(r2pol(@lb));
}
}
=item B<haset_ewxy>
$haset = haset_ewxy($declination, $latitude, %limits);
This routine takes the $declination of the source, and the $latitude of the
EWXY mounted antenna and calculates the hour angle at which the source
will set. It is then trivial to calculate the time until the source
sets, simply by subtracting the current hour angle of the source from
the hour angle at which it sets.
The required inputs are :
$declination - The declination of the source (turns)
$latitude - The latitude of the observatory (turns)
%limits - A reference to a hash holding the EWXY antenna limits
The following keys must be defined XLOW, XLOW_KEYHOLE,
XHIGH, XHIGH_KEYHOLE, YLOW, YLOW_KEYHOLE, YHIGH,
YHIGH_KEYHOLE (all values shoule be in turns)
The returned value is :
$haset - The hour angle (turns) at which a source at this
declination sets for an EWXY mounted antenna with the
given limits at the given latitude
NOTE: returns undef if %limits hash is missing any of the required keys
=cut
sub haset_ewxy($$\%) {
my ($declination, $latitude, $limitsref) = @_;
# Check that all the required keys are present
if ((!exists $limitsref->{XLOW}) || (!exists $limitsref->{XLOW_KEYHOLE}) ||
(!exists $limitsref->{XHIGH}) || (!exists $limitsref->{XHIGH_KEYHOLE}) ||
(!exists $limitsref->{YLOW}) || (!exists $limitsref->{YLOW_KEYHOLE}) ||
(!exists $limitsref->{YHIGH}) || (!exists $limitsref->{YHIGH_KEYHOLE})) {
carp 'Missing key in %limits';
return undef;
}
# Local variables
my ($pole, $pxlim, $exlim, $hix, $hixk, $lowx, $lowxk);
if ($latitude < 0.0) {
$pole = -90.0/360.0;
$pxlim = $limitsref->{XLOW};
$exlim = $limitsref->{XHIGH};
} else {
$pole = 90.0/360.0;
$pxlim = $limitsref->{XHIGH};
$exlim = $limitsref->{XLOW};
}
my $dec_never = $latitude + $exlim;
my $dec_always = $pole - ($latitude + $pxlim - $pole);
if ((($latitude < 0.0) && ($declination > $dec_never)) ||
(($latitude > 0.0) && ($declination < $dec_never))) {
# Source is never up
return(0.0);
} elsif ((($latitude < 0.0) && ($declination < $dec_always)) ||
(($latitude > 0.0) && ($declination > $dec_always))) {
# Source is always up
return(1.0);
} else {
# Up some of the time - calculate the ghastly constants and
# do everything in radians from here on.
$declination = 2.0 * $PI * $declination;
$latitude = 2.0 * $PI * $latitude;
my $k0 = -cos($declination);
my $k1 = sin($declination)*cos($latitude);
my $k2 = sin($declination)*sin($latitude);
my $k3 = cos($declination)*sin($latitude);
my $k4 = cos($declination)*cos($latitude);
my $k5 = $k4 * $k1 + $k2 * $k3;
my $x = 2.0 * $PI * $limitsref->{XLOW_KEYHOLE};
my $dec_split = asin(cos(2.0 * $PI * $limitsref->{YLOW}) *
(cos($x) * sin($latitude) + sin($x) *
cos($latitude)));
if ($latitude > 0.0) {
# Set up for northern antenna
$hix = $limitsref->{XLOW};
$hixk = $limitsref->{XLOW_KEYHOLE};
$lowx = $limitsref->{XHIGH};
$lowxk = $limitsref->{XHIGH_KEYHOLE};
} else {
# Set up for southern antenna
$hix = $limitsref->{XHIGH};
$hixk = $limitsref->{XHIGH_KEYHOLE};
$lowx = $limitsref->{XLOW};
$lowxk = $limitsref->{XLOW_KEYHOLE};
}
if ((($declination > $dec_split) && ($latitude < 0.0)) ||
(($declination < $dec_split) && ($latitude > 0.0))) {
# We are on the equatorial side
my $x = 2.0 * $PI * $hix;
my $y = -1.0 * abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y))/
($k3 + $k4))/(2.0 * $PI));
} else {
my $x = 2.0 * $PI * $hixk;
my $y = -1.0 * abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
return(asin(sin(2.0 * $PI * $limitsref->{YLOW_KEYHOLE}) /
$k0)/(2.0 * $PI));
} elsif (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW})) {
return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
($k3 + $k4))/(2.0 * $PI));
} else {
return(asin(sin(2.0 * $PI*$limitsref->{YLOW}) / $k0) /
(2.0 * $PI));
}
}
} else {
# We are on the polar side
my $x = 2.0 * $PI * $lowx;
my $y = abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
($k3 + $k4))/(2.0 * $PI));
} else {
my $x = 2.0 * $PI * $lowxk;
my $y = -1.0 * abs(acos($k5 /($k4 * sin($x) + $k3 * cos($x))));
if (abs($y) < abs(2.0 * $PI* $limitsref->{YLOW_KEYHOLE})) {
return(asin(sin(2.0 * $PI * $limitsref->{YLOW_KEYHOLE}) /
$k0)/(2.0 * $PI));
} elsif (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW})) {
return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
($k3 + $k4))/(2.0 * $PI));
} else {
return(asin(sin(2.0 * $PI * $limitsref->{YLOW}) / $k0)/
(2.0 * $PI));
}
}
}
}
}
=item B<ewxy_tlos>
$tlos = ewxy_tlos($hour_angle, $declination, $latitude, %limits);
This routine calculates the time left on-source (tlos) for a source
at $hour_angle, $declination for an EWXY mount antenna at $latitude.
The required inputs are :
$hour_angle - The current hour angle of the source (turns)
$declination - The declination of the source (turns)
$latitude - The latitude of the observatory (turns)
\%limits - A reference to a hash holding the EWXY antenna limits
The following keys must be defined XLOW, XLOW_KEYHOLE,
XHIGH, XHIGH_KEYHOLE, YLOW, YLOW_KEYHOLE, YHIGH,
YHIGH_KEYHOLE (all values should be in turns)
The returned value is :
$tlos - The time left on-source (turns)
=cut
sub ewxy_tlos($$$\%) {
my ($hour_angle, $declination, $latitude, $limitsref) = @_;
my $haset = haset_ewxy($declination, $latitude, %$limitsref);
return(undef) if (!defined $haset);
$haset -= $hour_angle if (($haset > 0.0) && ($haset < 1.0));
$haset += 1.0 if ($haset < 0.0);
return $haset;
}
=item B<haset_azel>
$haset = haset_azel($declination, $latitude, %limits);
This routine takes the $declination of the source, and the
$latitude of the Az/El mounted antenna and calculates the hour
angle at which the source will set. It is then trivial to
calculate the time until the source sets, simply by subtracting the
current hour angle of the source from the hour angle at which it
sets. This routine assumes that the antenna is able to rotate
through 360 degrees in azimuth.
The required inputs are :
$declination - The declination of the source (turns)
$latitude - The latitude of the observatory (turns)
\%limits - A reference to a hash holding the Az/El antenna limits
The following keys must be defined ELLOW (all values should
be in turns)
The returned value is :
$haset - The hour angle (turns) at which a source at this
declination sets for an Az/El mounted antenna with the
given limits at the given latitude
NOTE: returns undef if the %limits hash is missing any of the required keys
=cut
sub haset_azel($$\%) {
my ($declination, $latitude, $limitsref) = @_;
# Check that all the required keys are present
if (!exists $limitsref->{ELLOW}) {
carp 'Missing key in %limits';
return undef ;
}
my $cos_haset = (cos($PI / 2.0 - $limitsref->{ELLOW} * 2.0 *
$PI) - sin($latitude * 2.0 * $PI) *
sin($declination * 2.0 * $PI))/
(cos($declination * 2.0 * $PI)
*cos($latitude * 2.0 * $PI));
if ($cos_haset > 1.0) {
# The source never rises
return(0.0);
} elsif ($cos_haset < -1.0) {
# The source never sets
return(1.0);
} else {
return(acos($cos_haset)/(2.0*$PI));
}
}
=item B<azel_tlos>
$tlos = azel_tlos($hour_angle, $declination, $latitude, \%limits);
This routine calculates the time left on-source (tlos) for a source
at $hour_angle, $declination for an Az/El mount antenna at $latitude.
The required inputs are :
$hour_angle - The current hour angle of the source (turns)
$declination - The declination of the source (turns)
$latitude - The latitude of the observatory (turns)
%limits - A reference to a hash holding the Az/El antenna limits
The following keys must be defined ELLOW (all values
should be in turns)
The returned value is :
$tlos - The time left on-source (turns)
=cut
sub azel_tlos($$$\%) {
my ($hour_angle, $declination, $latitude, $limitsref) = @_;
# Calculate the time left onsource
my $haset = haset_azel($declination, $latitude, %$limitsref);
if (!defined $haset) {return(undef)};
if (($haset > 0.0) && ($haset < 1.0)) { $haset -= $hour_angle; }
if ($haset < 0.0) { $haset += 1.0; }
return($haset);
}
=item B<antenna_rise>
$ha_set = antenna_rise($declination, $latitude, $mount, \%limits);
Given the $declination of the source, the $latitude of the antenna,
the type of the antenna $mount and a reference to a hash holding
information on the antenna limits, this routine calculates the hour
angle at which the source sets for the antenna. The hour angle at
which it rises is simply the negative of that at which it sets.
These values in turn can be used to calculate the LMST at which the
source rises/sets and from that the UT at which the source
rises/sets on a given day, or to calculate the native coordinates
at which the source rises/sets.
If you want to calculate source rise times above arbitrary elevation,
use the routine rise.
The required inputs are :
$declination - The declination of the source (turns)
$latitude - The latitude of the observatory (turns)
$mount - The type of antenna mount, 0 => EWXY mount, 1 => Az/El,
any other number will cause the routine to return
undef
%limits - A reference to a hash holding the antenna limits
For an EWXY antenna there must be keys for all the
limits (i.e. XLOW, XLOW_KEYHOLE, XHIGH, XHIGH_KEYHOLE,
YLOW, YLOW_KEYHOLE, YHIGH, YHIGH_KEYHOLE). For an Az/El
antenna there must be a key for ELLOW (all values should
be in turns).
The returned values are :
$ha_set - The hour angle at which the source sets (turns). The hour
angle at which the source rises is simply the negative of this
value.
=cut
sub antenna_rise($$$$) {
my ($declination, $latitude, $mount, $limitsref) = @_;
# Check that the mount type is either EWXY (0) or AZEL (1)
if (($mount != 0) && ($mount != 1)) {
carp 'mount must equal 0 or 1';
return undef;
}
if ($mount == 0) {
return(haset_ewxy($declination, $latitude, %$limitsref));
} elsif ($mount == 1) {
return(haset_azel($declination, $latitude, %$limitsref));
}
}
my @b2g = ([-0.054875539726, 0.494109453312, -0.867666135858],
[-0.873437108010, -0.444829589425, -0.198076386122],
[-0.483834985808, 0.746982251810, 0.455983795705]);
#my @b2g = ([ -0.0548777621, +0.4941083214, -0.8676666398],
# [ -0.8734369591, -0.4448308610, -0.1980741871],
# [ -0.4838350026, +0.7469822433, +0.4559837919]);
sub j2gal($$) {
my ($ra,$dec) = @_;
my @r = pol2r($ra,$dec);
my @g = (0,0,0);
for (my $i=0; $i<3; $i++) {
for (my $j=0; $j<3; $j++) {
$g[$i]+= $b2g[$j][$i] * $r[$j];
}
}
return r2pol(@g);
}
# SLALIB support routines
sub epb2d ($) {
# - - - - - -
# E P B 2 D
# - - - - - -
#
# Conversion of Besselian Epoch to Modified Julian Date
# (double precision)
#
# Given:
# EPB dp Besselian Epoch
#
# The result is the Modified Julian Date (JD - 2400000.5).
#
# Reference:
# Lieske,J.H., 1979. Astron.Astrophys.,73,282.
#
# P.T.Wallace Starlink February 1984
#
# Copyright (C) 1995 Rutherford Appleton Laboratory
my $epb = shift;
return 15019.81352 + ($epb-1900)*365.242198781;
}
sub epj ($) {
# - - - -
# E P J
# - - - -
#
# Conversion of Modified Julian Date to Julian Epoch (double precision)
#
# Given:
# DATE dp Modified Julian Date (JD - 2400000.5)
#
# The result is the Julian Epoch.
#
# Reference:
# Lieske,J.H., 1979. Astron.Astrophys.,73,282.
#
# P.T.Wallace Starlink February 1984
#
# Copyright (C) 1995 Rutherford Appleton Laboratory
my $date = shift;
return 2000 + ($date-51544.5)/365.25;
}
sub pm ($$$$$$$$$$) {
# - - -
# P M
# - - -
#
# Apply corrections for proper motion to a star RA,Dec
# (double precision)
#
# References:
# 1984 Astronomical Almanac, pp B39-B41.
# (also Lederle & Schwan, Astron. Astrophys. 134,
# 1-6, 1984)
#
# Given:
# R0,D0 dp RA,Dec at epoch EP0 (rad)
# PR,PD dp proper motions: RA,Dec changes per year of epoch
# EP0 dp start epoch in years (e.g. Julian epoch)
# EP1 dp end epoch in years (same system as EP0)
#
# Returned:
# R1,D1 dp RA,Dec at epoch EP1 (rad)
#
# Called:
# sla_DCS2C spherical to Cartesian
# sla_DCC2S Cartesian to spherical
# sla_DRANRM normalize angle 0-2Pi
#
# Note:
# The proper motions in RA are dRA/dt rather than
# cos(Dec)*dRA/dt, and are in the same coordinate
# system as R0,D0.
#
# P.T.Wallace Starlink 23 August 1996
#
# Copyright (C) 1996 Rutherford Appleton Laboratory
my ($r0, $d0, $pr, $pd, $ep0, $ep1) = @_;
# Km/s to AU/year multiplied by arc seconds to radians
use constant VFR => 0.21094502*0.484813681109535994e-5;
my (@em, $t);
# Spherical to Cartesian
my @p = pol2r($r0,$d0);
# Space motion (radians per year)
$em[0]=-$pr*$p[1]-$pd*cos($r0)*sin($d0);
$em[1]= $pr*$p[0]-$pd*sin($r0)*sin($d0);
$em[2]= $pd*cos($d0);
# Apply the motion
$t=$ep1-$ep0;
for (my $i = 0; $i<3; $i++) {
$p[$i]=$p[$i]+$t*$em[$i];
}
# Cartesian to spherical
return r2pol(@p);
}
1;
__END__