← Index
Performance Profile   « block view • line view • sub view »
For /wise/base/deliv/dev/bin/getfix
  Run on Thu May 20 15:30:03 2010
Reported on Thu May 20 16:25:28 2010

File/wise/base/static/lib/perl5/site_perl/5.10.0/Astro/Coord.pm
Statements Executed1012491
Total Time1.76655099999729 seconds

Subroutines — ordered by exclusive time
Calls P F Exclusive
Time		 Inclusive
Time		 Subroutine
13869111.079862.18973Astro::Coord::fk5fk4
13869110.249080.81241Astro::Coord::pol2r
13869110.234670.29746Astro::Coord::r2pol
00000Astro::Coord::BEGIN
00000Astro::Coord::J2000todate
00000Astro::Coord::antenna_rise
00000Astro::Coord::azel2xy
00000Astro::Coord::azel_tlos
00000Astro::Coord::coord_convert
00000Astro::Coord::epb2d
00000Astro::Coord::ephem_vars
00000Astro::Coord::epj
00000Astro::Coord::eqazel
00000Astro::Coord::ewxy_tlos
00000Astro::Coord::fk4fk5
00000Astro::Coord::fk4fk5r
00000Astro::Coord::fk4gal
00000Astro::Coord::fk4galr
00000Astro::Coord::fk5fk4r
00000Astro::Coord::galfk4
00000Astro::Coord::galfk4r
00000Astro::Coord::haset_azel
00000Astro::Coord::haset_ewxy
00000Astro::Coord::j2gal
00000Astro::Coord::nutate
00000Astro::Coord::pm
00000Astro::Coord::precsn
00000Astro::Coord::xy2azel

LineStmts.Exclusive
Time		Avg.Code
1package Astro::Coord;
236.0e-52.0e-5use strict;
# spent 24µs making 1 call to strict::import
3
4=head1 NAME
5
6Astro::Coord - Astronomical coordinate transformations
7
8=head1 SYNOPSIS
9
10 use Astro::Coord;
11
12 ($l, $b) = fk4gal($ra, $dec);
13 ($az, $el) = eqazel($ha, $dec, $latitude);
14
15=head1 DESCRIPTION
16
17Astro::Coord contains an assorted set Perl routines for coordinate
18conversions, such as hour angle to elevation and J2000 to B1950.
19
20=head1 AUTHOR
21
22Chris Phillips Chris.Phillips@csiro.au
23
24=head1 FUNCTIONS
25
26=cut
27
28
29BEGIN {
3032.9e-59.7e-6 use Exporter ();
31 use vars qw( $VERSION @ISA @EXPORT @EXPORT_OK @EXPORT_FAIL
# spent 74µs making 1 call to vars::import
3237.4e-52.5e-5 $bepoch );
3351.4e-52.8e-6 $VERSION = '1.42';
34
35 @ISA = qw(Exporter);
36
37 @EXPORT = qw( xy2azel azel2xy eqazel J2000todate
38 fk4fk5 fk5fk4 fk4gal galfk4 j2gal
39 coord_convert
40 haset_ewxy ewxy_tlos haset_azel azel_tlos
41 antenna_rise pol2r r2pol
42 );
43 @EXPORT_OK = qw ( fk4fk5r fk5fk4r fk4galr galfk4r
44 ephem_vars nutate precsn $bepoch );
45 @EXPORT_FAIL = qw ( );
46
4733.3e-51.1e-5 use Carp;
# spent 57µs making 1 call to Exporter::import
4833.1e-51.0e-5 use POSIX qw( asin acos fmod tan );
# spent 80µs making 1 call to POSIX::import
49
5030.000790.00026 use Astro::Time qw( $PI rad2turn turn2rad mjd2lst );
# spent 275µs making 1 call to Exporter::import
5113.4e-53.4e-5}
52
5311.0e-61.0e-6$bepoch = 1950.0;
54
5533.4e-51.1e-5use constant JULIAN_DAY_J2000 => 2451545.0;
# spent 66µs making 1 call to constant::import
5636.3e-52.1e-5use constant JULIAN_DAYS_IN_CENTURY => 36525.0;
# spent 42µs making 1 call to constant::import
57
58# The E-terms vector for FK4 <--> other coordinate system transforms
59# (used in fk4fk5 fk5fk4 fk4gal galfk4)
6011.0e-61.0e-6my @eterm = (-1.62557E-06, -0.31919E-06, -0.13843E-06);
61
62## The precession matrix for FK4 <--> FK5 conversions (used in
63## fk4fk5 and fk5fk4)
64#my @btoj = ([+0.999925678186902,-0.011182059642247,-0.004857946558960],
65# [+0.011182059571766,+0.999937478448132,-0.000027176441185],
66# [+0.004857946721186,-0.000027147426498,+0.999988199738770]);
67
68# The precession matrix for FK4 <--> Galactic conversions (used in
69# fk4gal and galfk4)
7013.0e-63.0e-6my @etog = ([-0.066988739415,-0.872755765852,-0.483538914632],
71 [+0.492728466075,-0.450346958020,+0.744584633283],
72 [-0.867600811151,-0.188374601723,+0.460199784784]);
73
74# Values used in SLALIB routines
75
7634.3e-51.4e-5use constant D2PI => 6.283185307179586476925287;
# spent 51µs making 1 call to constant::import
77
78# Radians per year to arcsec per century
7933.0e-51.0e-5use constant PMF => 100*60*60*360/D2PI;
# spent 45µs making 1 call to constant::import
80
81# Small number to avoid arithmetic problems
8232.8e-59.3e-6use constant TINY => 1e-30;
# spent 40µs making 1 call to constant::import
83
84# Km per sec to AU per tropical century
85# = 86400 * 36524.2198782 / 149597870
8630.005620.00187use constant VF => 21.095;
# spent 40µs making 1 call to constant::import
87
88# Vectors A and Adot, and matrix M
8911.0e-61.0e-6my @a = ( -1.62557e-6, -0.31919e-6, -0.13843e-6,
90 +1.245e-3, -1.580e-3, -0.659e-3);
91
9211.0e-61.0e-6my @ad =(+1.245e-3, -1.580e-3, -0.659e-3);
93
9414.0e-64.0e-6my @em = ([+0.9999256782, -0.0111820611, -0.0048579477],
95 [+0.0111820610, +0.9999374784, -0.0000271765],
96 [+0.0048579479, -0.0000271474, +0.9999881997],
97 [-0.000551, -0.238565, +0.435739],
98 [+0.238514, -0.002667, -0.008541],
99 [-0.435623, +0.012254, +0.002117]);
100
10111.2e-51.2e-5my @emi = ([+0.9999256795, +0.0111814828, +0.0048590039,
102 -0.00000242389840, -0.00000002710544, -0.00000001177742],
103 [-0.0111814828, +0.9999374849, -0.0000271771,
104 +0.00000002710544, -0.00000242392702, +0.00000000006585],
105 [-0.0048590040, -0.0000271557, +0.9999881946,
106 +0.00000001177742, +0.00000000006585, -0.00000242404995],
107 [-0.000551, +0.238509, -0.435614,
108 +0.99990432, +0.01118145, +0.00485852],
109 [-0.238560, -0.002667, +0.012254,
110 -0.01118145, +0.99991613, -0.00002717],
111 [+0.435730, -0.008541, +0.002117,
112 -0.00485852, -0.00002716, +0.99996684]);
113
114=item B<pol2r>
115
116 ($x, $y, $z) = pol2r($polar1, $polar2);
117
118 Converts a position in polar coordinates into rectangular coordinates
119 $polar1, $polar2 The polar coordinates to convert (turns)
120 $x, $y, $z The rectangular coordinates
121
122=cut
123
124
# spent 812ms (249+563) within Astro::Coord::pol2r which was called 13869 times, avg 59µs/call: # 13869 times (249ms+563ms) by Astro::Coord::fk5fk4 at line 545, avg 59µs/call
sub pol2r ($$) {
125832140.386984.7e-6 my ($p1, $p2) = @_;
126
127 # Converts polar coordinates into rectangluar
128 my @rect;
129 $rect[0] = cos(turn2rad($p1))*cos(turn2rad($p2));
# spent 241ms making 27738 calls to Astro::Time::turn2rad, avg 9µs/call
130 $rect[1] = sin(turn2rad($p1))*cos(turn2rad($p2));
# spent 215ms making 27738 calls to Astro::Time::turn2rad, avg 8µs/call
131 $rect[2] = sin(turn2rad($p2));
# spent 107ms making 13869 calls to Astro::Time::turn2rad, avg 8µs/call
132 return(@rect);
133}
134
135=item B<r2pol>
136
137 ($polar1, $polar2) = r2pol($x, $y, $z);
138
139 Converts a position in rectangular coordinates into polar coordinates
140 $x, $y, $y The rectangular coordinates to convert
141 $polar1, $polar2 The polar coordinates (turns);
142 Returns undef if too few or too many arguments are passed.
143
144=cut
145
146
# spent 297ms (235+62.8) within Astro::Coord::r2pol which was called 13869 times, avg 21µs/call: # 13869 times (235ms+62.8ms) by Astro::Coord::fk5fk4 at line 600, avg 21µs/call
sub r2pol (@) {
147 # First check that we have 3 arguments
1481386900.251211.8e-6 if (scalar @_ != 3) {
149 carp 'Inconsistent arguments';
150 return undef ;
151 }
152 my ($x, $y, $z) = @_;
153
154 # Converts rectangular coordinates to polar
155 my ($tmp, $left, $right);
156 $tmp = atan2($y, $x)/(2.0*$PI);
157
158 if (ref($tmp) =~ /PDL/ ) { # Allow to work with PDL
159 $tmp -> where($tmp<0.0) .= $tmp -> where($tmp<0.0) + 1.0;
160 } elsif ($tmp < 0.0) {
161 $tmp += 1.0;
162 }
163
164 $left = $tmp;
165 $tmp = sqrt($x*$x + $y*$y + $z*$z);
166
167 if (ref($tmp) =~ /PDL/) { # Allow to work with PDL
168 $right = &PDL::Math::asin($z/$tmp)/(2.0*$PI);
169 } else {
170 $right = asin($z/$tmp)/(2.0*$PI);
# spent 62.8ms making 13869 calls to POSIX::asin, avg 5µs/call
171 }
172
173 return ($left, $right);
174}
175
176=item B<xy2azel>
177
178 ($az, $el) = xy2azel($x, $y);
179
180 Converts a telescope position in X,Y coordinates into Az,El coordinates
181 $x, $y The X and Y coordinates (turns)
182 $az, $el The azimuth and elevation (turns)
183
184=cut
185
186sub xy2azel ($$) {
187 my ($x, $y) = @_;
188
189 # Convert a position in X,Y to Az,El
190 my @polar = pol2r($x, $y);
191 my $temp = $polar[0];
192 $polar[0] = $polar[1];
193 $polar[1] = $polar[2];
194 $polar[2] = $temp;
195 return (r2pol(@polar));
196}
197
198=item B<azel2xy>
199
200 ($x, $y) = azel2xy($az, $el);
201
202 Converts a position in Az,El coordinates into X,Y coordinates
203 $az, $el The azimuth and elevation (turns)
204 $x, $y The X and Y coordinate (turns)
205
206=cut
207
208sub azel2xy ($$) {
209 my ($az, $el) = @_;
210
211 # Convert a position in Az,El to X,Y
212 my @polar = pol2r($az, $el);
213 my $temp = $polar[1];
214 $polar[1] = $polar[0];
215 $polar[0] = $polar[2];
216 $polar[2] = $temp;
217 my ($x, $y) = r2pol(@polar);
218 if ($x>0.5) {
219 $x -= 1.0;
220 }
221 if ($y>0.5) {
222 $y -= 1.0;
223 }
224 return ($x, $y);
225}
226
227=item B<eqazel>
228
229 ($ha, $dec) = eqazel($az, $el, $latitude);
230 ($az, $el) = eqazel($ha, $dec, $latitude);
231
232 Converts HA/Dec coordinates to Az/El and vice versa.
233 $ha, $dec Hour angle and declination of source (turns)
234 $az, $el Azimuth and elevation of source (turns)
235 $latitude Latitude of the observatory (turns)
236 Note:
237 The ha,dec and az,el conversion is symmetrical
238
239=cut
240
241sub eqazel ($$$) {
242 my $sphi = sin(turn2rad($_[2]));
243 my $cphi = cos(turn2rad($_[2]));
244 my $sleft = sin(turn2rad($_[0]));
245 my $cleft = cos(turn2rad($_[0]));
246 my $sright = sin(turn2rad($_[1]));
247 my $cright = cos(turn2rad($_[1]));
248 my $left_out = atan2(-$sleft,-$cleft*$sphi+$sright*$cphi/$cright)/(2.0*$PI);
249 $left_out = ($left_out < 0.0) ? $left_out + 1.0 : $left_out;
250 my $right_out= asin($cleft*$cright*$cphi + $sright*$sphi)/(2.0*$PI);
251
252 return($left_out, $right_out);
253
254}
255
256=item B<fk4fk5>
257
258 ($JRA, $JDec) = fk4fk5($BRA, $BDec);
259 (@fk5) = fk4fk5(@fk4);
260
261 Converts an FK4 (B1950) position to the equivalent FK5 (J2000)
262 position.
263 $BRA,$BDec fk4/B1950 position (turns)
264 $JRA,$Dec fk5/J2000 position (turns)
265 @fk4 fk4/B1950 position (as a 3-vector)
266 @fk5 fk5/J2000 position (as a 3-vector)
267 Note:
268 This code is based on similar routines from the Fortran SLALIB
269 package, so are quite accurate, but subject to a restrictive
270 license (see README).
271
272=cut
273
274sub fk4fk5 (@) {
275# - - - - - -
276# F K 4 5 Z
277# - - - - - -
278#
279# Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
280# proper motion in the FK5 frame (double precision)
281#
282# This routine converts stars from the old, Bessel-Newcomb, FK4
283# system to the new, IAU 1976, FK5, Fricke system, in such a
284# way that the FK5 proper motion is zero. Because such a star
285# has, in general, a non-zero proper motion in the FK4 system,
286# the routine requires the epoch at which the position in the
287# FK4 system was determined.
288#
289# The method is from Appendix 2 of Ref 1, but using the constants
290# of Ref 4.
291#
292# Given:
293# R1950,D1950 dp B1950.0 FK4 RA,Dec at epoch (rad)
294# BEPOCH dp Besselian epoch (e.g. 1979.3D0)
295#
296# Returned:
297# R2000,D2000 dp J2000.0 FK5 RA,Dec (rad)
298#
299# Notes:
300#
301# 1) The epoch BEPOCH is strictly speaking Besselian, but
302# if a Julian epoch is supplied the result will be
303# affected only to a negligible extent.
304#
305# 2) Conversion from Besselian epoch 1950.0 to Julian epoch
306# 2000.0 only is provided for. Conversions involving other
307# epochs will require use of the appropriate precession,
308# proper motion, and E-terms routines before and/or
309# after FK45Z is called.
310#
311# 3) In the FK4 catalogue the proper motions of stars within
312# 10 degrees of the poles do not embody the differential
313# E-term effect and should, strictly speaking, be handled
314# in a different manner from stars outside these regions.
315# However, given the general lack of homogeneity of the star
316# data available for routine astrometry, the difficulties of
317# handling positions that may have been determined from
318# astrometric fields spanning the polar and non-polar regions,
319# the likelihood that the differential E-terms effect was not
320# taken into account when allowing for proper motion in past
321# astrometry, and the undesirability of a discontinuity in
322# the algorithm, the decision has been made in this routine to
323# include the effect of differential E-terms on the proper
324# motions for all stars, whether polar or not. At epoch 2000,
325# and measuring on the sky rather than in terms of dRA, the
326# errors resulting from this simplification are less than
327# 1 milliarcsecond in position and 1 milliarcsecond per
328# century in proper motion.
329#
330# References:
331#
332# 1 Aoki,S., et al, 1983. Astron.Astrophys., 128, 263.
333#
334# 2 Smith, C.A. et al, 1989. "The transformation of astrometric
335# catalog systems to the equinox J2000.0". Astron.J. 97, 265.
336#
337# 3 Yallop, B.D. et al, 1989. "Transformation of mean star places
338# from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
339# Astron.J. 97, 274.
340#
341# 4 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
342# the Astronomical Almanac", ISBN 0-935702-68-7.
343#
344# Called: sla_DCS2C, sla_EPJ, sla_EPB2D, sla_DCC2S, sla_DRANRM
345#
346# P.T.Wallace Starlink 21 September 1998
347#
348# Copyright (C) 1998 Rutherford Appleton Laboratory
349
350
351 my ($rect, $w, $i, $j);
352 my (@r0, @a1, @v1, @v2); # Position and position+velocity vectors
353
354 if (@_==3) { # Rectangular coordinates passed
355 @r0 = @_;
356 $rect = 1;
357 } elsif (@_==2) { # Sperical coordinates
358 @r0 = pol2r($_[0],$_[1]); # Spherical to Cartesian
359 $rect = 0;
360 } elsif (@_>3) {
361 croak "Too many arguments for Astro::fk4fk5 ";
362 } else {
363 croak "Not enough arguments for Astro::fk4fk5 ";
364 }
365
366 # Adjust vector A to give zero proper motion in FK5
367 $w=($bepoch-1950)/PMF;
368 for ($i=0; $i<3; $i++) {
369 $a1[$i]=$a[$i]+$w*$ad[$i];
370 }
371 # Remove e-terms
372 $w=$r0[0]*$a1[0]+$r0[1]*$a1[1]+$r0[2]*$a1[2];
373 for ($i=0; $i<3; $i++) {
374 $v1[$i]=$r0[$i]-$a1[$i]+$w*$r0[$i];
375 }
376
377 # Convert position vector to Fricke system
378 for ($i=0; $i<6; $i++) {
379 $w=0;
380 for ($j=0; $j<3; $j++) {
381 #warn "DEBUG: [$i,$j]\n";
382 $w=$w+$em[$i][$j]*$v1[$j];
383 $v2[$i]=$w
384 }
385 }
386
387 # Allow for fictitious proper motion in FK4
388 $w=(epj(epb2d($bepoch))-2000)/PMF;
389 for ($i=0; $i<3; $i++) {
390 $v2[$i]=$v2[$i]+$w*$v2[$i+3];
391 }
392
393 if ($rect) {
394 return @v2[0..2];
395 } else {
396 # Revert to spherical coordinates
397 return r2pol(@v2[0..2]);
398 }
399}
400
401=item B<fk4fk5r>
402
403 @fk5 = fk4fk5r(@fk4);
404
405 Converts an FK4 (B1950) position to the equivalent FK5 (J2000) position.
406 Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
407 @fk4 fk4 position (as a 3-vector, turns)
408 @fk5 fk5 position (as a 3-vector, turns)
409 Note:
410 Just a wrapper to fk4fk5 which now handler polar and rectangular
411 arguments
412
413=cut
414
415sub fk4fk5r (@) {
416 return fk4fk5(@_);
417}
418
419#sub fk4fk5r (@) {
420# # First check that we have 3 arguments
421# if (scalar @_ < 3) {
422# croak 'Not enough arguments for Astro::Coord::fk4fk5r at ';
423# } elsif (scalar @_ > 3) {
424# croak 'Too many arguments for Astro::Coord::fk4fk5r at ';
425# }
426#
427# my ($i, $j, @temp, @fk5);
428# my $w = 0.0;
429#
430# # Add the eterms
431# for ($i=0 ; $i<3 ; $i++) {
432# $w += $_[$i] * $eterm[$i];
433# }
434# for ($i=0 ; $i<3 ; $i++) {
435# $temp[$i] = $_[$i] - $eterm[$i] + $w * $_[$i];
436# }
437#
438# # Precess from FK4 to FK5
439# for ($i=0 ; $i<3 ; $i++) {
440# $fk5[$i] = 0.0;
441# for ($j=0 ; $j<3 ; $j++) {
442# $fk5[$i] += $btoj[$i][$j] * $temp[$j];
443# }
444# }
445#
446# return(@fk5);
447#}
448
449=item B<fk5fk4>
450
451 ($JRA, $JDec) = fk4fk5($BRA, $BDec);
452 ($@fk5) = fk4fk5(@fk4);
453
454 Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
455 $JRA,$Dec fk5/J2000 position (turns)
456 $BRA,$BDec fk4/B1950 position (turns)
457 @fk5 fk5/J2000 position (as a 3-vector)
458 @fk4 fk4/B1950 position (as a 3-vector)
459 Note:
460 This code is based on similar routines from the Fortran SLALIB
461 package, so are quite accurate, but subject to a restrictive
462 license (see README).
463
464=cut
465
466
# spent 2.19s (1.08+1.11) within Astro::Coord::fk5fk4 which was called 13869 times, avg 158µs/call: # 13869 times (1.08s+1.11s) by WISE::CoUtils::_fk5fk4 at line 44 of /wise/base/deliv/dev/lib/perl/WISE/CoUtils.pm, avg 158µs/call
sub fk5fk4 (@) {
467#+
468# - - - - - -
469# F K 5 2 4
470# - - - - - -
471#
472# Convert J2000.0 FK5 star data to B1950.0 FK4 (double precision)
473#
474# This routine converts stars from the new, IAU 1976, FK5, Fricke
475# system, to the old, Bessel-Newcomb, FK4 system. The precepts
476# of Smith et al (Ref 1) are followed, using the implementation
477# by Yallop et al (Ref 2) of a matrix method due to Standish.
478# Kinoshita's development of Andoyer's post-Newcomb precession is
479# used. The numerical constants from Seidelmann et al (Ref 3) are
480# used canonically.
481#
482# Given: (all J2000.0,FK5)
483# R2000,D2000 dp J2000.0 RA,Dec (rad)
484# DR2000,DD2000 dp J2000.0 proper motions (rad/Jul.yr)
485# P2000 dp parallax (arcsec)
486# V2000 dp radial velocity (km/s, +ve = moving away)
487#
488# Returned: (all B1950.0,FK4)
489# R1950,D1950 dp B1950.0 RA,Dec (rad)
490# DR1950,DD1950 dp B1950.0 proper motions (rad/trop.yr)
491# P1950 dp parallax (arcsec)
492# V1950 dp radial velocity (km/s, +ve = moving away)
493#
494# Notes:
495#
496# 1) The proper motions in RA are dRA/dt rather than
497# cos(Dec)#dRA/dt, and are per year rather than per century.
498#
499# 2) Note that conversion from Julian epoch 2000.0 to Besselian
500# epoch 1950.0 only is provided for. Conversions involving
501# other epochs will require use of the appropriate precession,
502# proper motion, and E-terms routines before and/or after
503# FK524 is called.
504#
505# 3) In the FK4 catalogue the proper motions of stars within
506# 10 degrees of the poles do not embody the differential
507# E-term effect and should, strictly speaking, be handled
508# in a different manner from stars outside these regions.
509# However, given the general lack of homogeneity of the star
510# data available for routine astrometry, the difficulties of
511# handling positions that may have been determined from
512# astrometric fields spanning the polar and non-polar regions,
513# the likelihood that the differential E-terms effect was not
514# taken into account when allowing for proper motion in past
515# astrometry, and the undesirability of a discontinuity in
516# the algorithm, the decision has been made in this routine to
517# include the effect of differential E-terms on the proper
518# motions for all stars, whether polar or not. At epoch 2000,
519# and measuring on the sky rather than in terms of dRA, the
520# errors resulting from this simplification are less than
521# 1 milliarcsecond in position and 1 milliarcsecond per
522# century in proper motion.
523#
524# References:
525#
526# 1 Smith, C.A. et al, 1989. "The transformation of astrometric
527# catalog systems to the equinox J2000.0". Astron.J. 97, 265.
528#
529# 2 Yallop, B.D. et al, 1989. "Transformation of mean star places
530# from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
531# Astron.J. 97, 274.
532#
533# 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to
534# the Astronomical Almanac", ISBN 0-935702-68-7.
535#
536# P.T.Wallace Starlink 19 December 1993
537#
538# Copyright (C) 1995 Rutherford Appleton Laboratory
539#-
5407905331.121311.4e-6 my ($rect, @v1, @v2);
541 if (@_==3) { # Rectangular coordinates passed
542 @v1 = @_;
543 $rect = 1;
544 } elsif (@_==2) { # Sperical coordinates
545 @v1 = pol2r($_[0],$_[1]); # Spherical to Cartesian
# spent 812ms making 13869 calls to Astro::Coord::pol2r, avg 59µs/call
546 $rect = 0;
547 } elsif (@_>2) {
548 croak "Too many arguments for Astro::fk5fk4 ";
549 } else {
550 croak "Not enough arguments for Astro::fk5fk4 ";
551 }
552
553# Miscellaneous
554 my ($w, $x, $y, $z, $wd, $rxyz);
555 my ($ur, $ud, $xd, $yd, $zd);
556 my ($i,$j);
557
558 # Convert position+velocity vector to BN system
559 for ($i=0; $i<6; $i++) {
560 $w=0.0;
561 ##for ($j=0; $j<6; $j++) {
562 for ($j=0; $j<3; $j++) {
563 $w=$w+$emi[$i][$j]*$v1[$j];
564 }
565 $v2[$i]=$w;
566 }
567
568# Position vector components and magnitude
569 $x=$v2[0];
570 $y=$v2[1];
571 $z=$v2[2];
572 $rxyz=sqrt($x*$x+$y*$y+$z*$z);
573
574# Apply E-terms to position
575 $w=$x*$a[0]+$y*$a[1]+$z*$a[2];
576 $x=$x+$a[0]*$rxyz-$w*$x;
577 $y=$y+$a[1]*$rxyz-$w*$y;
578 $z=$z+$a[2]*$rxyz-$w*$z;
579
580# Recompute magnitude
581 $rxyz=sqrt($x*$x+$y*$y+$z*$z);
582
583# Apply E-terms to both position and velocity
584 $x=$v2[0];
585 $y=$v2[1];
586 $z=$v2[2];
587 $w=$x*$a[0]+$y*$a[1]+$z*$a[2];
588 $wd=$x*$a[3]+$y*$a[4]+$z*$a[5];
589 $x=$x+$a[0]*$rxyz-$w*$x;
590 $y=$y+$a[1]*$rxyz-$w*$y;
591 $z=$z+$a[2]*$rxyz-$w*$z;
592 $xd=$v2[3]+$a[3]*$rxyz-$wd*$x;
593 $yd=$v2[4]+$a[4]*$rxyz-$wd*$y;
594 $zd=$v2[5]+$a[5]*$rxyz-$wd*$z;
595
596 my @r;
597 if ($rect) {
598 @r = ($x, $y, $z);
599 } else {
600 @r= r2pol($x, $y, $z);
# spent 297ms making 13869 calls to Astro::Coord::r2pol, avg 21µs/call
601 }
602
603# my $rxysq =$x*$x+$y*$y;
604# my $rxy = sqrt($rxysq);
605# if ($rxy>TINY) {
606# $ur=($x*$yd-$y*$xd)/$rxysq;
607# $ud=($zd*$rxysq-$z*($x*$xd+$y*$yd))/(($rxysq+$z*$z)*$rxy);
608# }
609#
610## Return results
611# my $dr1950=$ur/PMF;
612# my $dd1950=$ud/PMF;
613
614 return(@r);
615}
616
617
618=item B<fk5fk4r>
619
620 @fk4 = fk5fk4r(@fk5);
621
622 Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
623 Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
624 @fk4 fk4 position (as a 3-vector, turns)
625 @fk5 fk5 position (as a 3-vector, turns)
626 Note:
627 Just a wrapper to fk5fk4 which now handler polar and rectangular
628 arguments
629
630=cut
631
632sub fk5fk4r (@) {
633 return fk5fk4(@_);
634}
635
636
637#sub fk5fk4 (@) {
638## - - - - - -
639## F K 5 4 Z
640## - - - - - -
641##
642## Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming
643## zero proper motion and parallax (double precision)
644##
645## This routine converts star positions from the new, IAU 1976,
646## FK5, Fricke system to the old, Bessel-Newcomb, FK4 system.
647##
648## Given:
649## R2000,D2000 dp J2000.0 FK5 RA,Dec (rad)
650## BEPOCH dp Besselian epoch (e.g. 1950D0)
651##
652## Returned:
653## R1950,D1950 dp B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH
654## DR1950,DD1950 dp B1950.0 FK4 proper motions (rad/trop.yr)
655##
656## Notes:
657##
658## 1) The proper motion in RA is dRA/dt rather than cos(Dec)#dRA/dt.
659##
660## 2) Conversion from Julian epoch 2000.0 to Besselian epoch 1950.0
661## only is provided for. Conversions involving other epochs will
662## require use of the appropriate precession routines before and
663## after this routine is called.
664##
665## 3) Unlike in the sla_FK524 routine, the FK5 proper motions, the
666## parallax and the radial velocity are presumed zero.
667##
668## 4) It is the intention that FK5 should be a close approximation
669## to an inertial frame, so that distant objects have zero proper
670## motion; such objects have (in general) non-zero proper motion
671## in FK4, and this routine returns those fictitious proper
672## motions.
673##
674## 5) The position returned by this routine is in the B1950
675## reference frame but at Besselian epoch BEPOCH. For
676## comparison with catalogues the BEPOCH argument will
677## frequently be 1950D0.
678##
679## Called: sla_FK524, sla_PM
680##
681## P.T.Wallace Starlink 10 April 1990
682##
683## Copyright (C) 1995 Rutherford Appleton Laboratory
684#
685# my $bepoch = 1950.0;
686#
687# my $rect;
688# if (@_>3) {
689# croak "Too many arguments for Astro::fk5fk4 ";
690# } elsif (@_<2) {
691# croak "Not enough arguments for Astro::fk5fk4 ";
692# }
693# my @r2000 = @_;
694#
695# # fk5 equinox j2000 (any epoch) to fk4 equinox b1950 epoch b1950
696# my (@r1950) = fk524(@r2000);
697# my $dd1950 = pop @r1950;
698# my $dr1950 = pop @r1950;
699#
700# ## fictitious proper motion to epoch bepoch
701# #my ($r1950, $d1950) = pm($r,$d,$dr1950,$dd1950,0.0,0.0,1950,$bepoch);
702# return @r1950;
703#}
704
705#=item B<fk5fk4r>
706#
707# @fk4 = fk5fk4r(@fk5);
708#
709# Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
710# Note: Convert equitoral positions to/from 3-vectors using pol2r and r2pol.
711# @fk5 fk5 position (as a 3-vector, turns)
712# @fk4 fk4 position (as a 3-vector, turns)
713#
714#=cut
715#
716#sub fk5fk4r(@) {
717#
718# # First check that we have 3 arguments
719# if (scalar @_ < 3) {
720# croak 'Not enough arguments for Astro::Coord::fk5fk4r at ';
721# } elsif (scalar @_ > 3) {
722# croak 'Too many arguments for Astro::Coord::fk5fk4r at ';
723# }
724#
725# my ($i, $j, @fk4);
726# my $w = 0.0;
727#
728# # Precess. Note : the same matrix is used as for the FK4 -> FK5
729# # transformation, but we have transposed it within the
730# # for loop
731#
732# for ($i=0 ; $i<3 ; $i++) {
733# $fk4[$i] = 0.0;
734# for ($j=0 ; $j<3 ; $j++) {
735# $fk4[$i] += $btoj[$j][$i] * $_[$j];
736# }
737# }
738#
739# # Allow for e-terms
740# for ($i=0 ; $i<3 ; $i++) {
741# $w += $_[$i] * $eterm[$i];
742# }
743# $w += 1.0;
744# for ($i=0 ; $i<3 ; $i++) {
745# $fk4[$i] = ($fk4[$i] + $eterm[$i])/$w;
746# }
747#
748# return(@fk4);
749#}
750
751=item B<fk4galr>
752
753 @gal = fk4galr(@fk4)
754
755 Converts an FK4 position (B1950.0) to the IAU 1958 Galactic
756 coordinate system
757 Note: convert equitoral positions to/from 3-vectors using pol2r and r2pol.
758 @fk4 fk4 position to convert (as a 3-vector, turns)
759 @gal Galactic position (as a 3-vector, turns)
760 Returns undef if too few or two many arguments are passed.
761 Reference : Blaauw et al., 1960, MNRAS, 121, 123.
762
763=cut
764
765# Within 1e-7 arcsec of SLALIB slaEg50
766sub fk4galr(@) {
767 # First check that we have 3 arguments
768 if (scalar @_ < 3) {
769 croak 'Not enough arguments for Astro::Coord::fk4galr at ';
770 } elsif (scalar @_ > 3) {
771 croak 'Too many arguments for Astro::Coord::fk4galr at ';
772 }
773
774 my ($i, $j, @temp, @gal);
775 my $w = 0.0;
776
777 # Allow for e-terms
778 for ($i=0 ; $i<3 ; $i++) {
779 $w += $_[$i] * $eterm[$i];
780 }
781 for ($i=0 ; $i<3 ; $i++) {
782 $temp[$i] = $_[$i] - $eterm[$i] + $w * $_[$i];
783 }
784
785
786 # Precess
787 for ($i=0 ; $i<3 ; $i++) {
788 $gal[$i] = 0.0;
789 for ($j=0 ; $j<3 ; $j++) {
790 $gal[$i] += $etog[$i][$j] * $temp[$j];
791 }
792 }
793
794 return(@gal);
795}
796
797=item B<galfk4>
798
799 ($bRA, $bDec) = galfk4($l, $b);
800 @fk4 = galfk4(@gal);
801
802 Converts an IAU 1958 Galactic position to the FK4 coordinate system (B1950)
803 Notes: Converts equitoral positions to/from 3-vectors using pol2r and r2pol.
804 $BRA,$BDec fk4/B1950 position (turns)
805 $l, $b Galactic longitude and latitude
806 @gal Galactic position (as a 3-vector, turns)
807 @fk4 fk4 position (as a 3-vector, turns)
808 Reference : Blaauw et al., 1960, MNRAS, 121, 123.
809
810=cut
811
812# Within 1e-7 arcsec of SLALIB slaGe50
813sub galfk4(@) {
814 my (@r, $rect);
815
816 if (@_==3) { # Rectangular coordinates passed
817 @r = @_;
818 $rect = 1;
819 } elsif (@_==2) { # Sperical coordinates
820 @r = pol2r($_[0],$_[1]); # Spherical to Cartesian
821 $rect = 0;
822 } elsif (@_>3) {
823 croak "Too many arguments for Astro::galfk4 at";
824 } else {
825 croak "Not enough arguments for Astro::galfk4 at";
826 }
827
828 my ($i, $j, @fk4);
829 my $w = 0.0;
830
831 # Precess. Note : the same matrix is used as for the FK4 -> Galactic
832 # transformation, but we have transposed it within the
833 # for loop
834 for ($i=0 ; $i<3 ; $i++) {
835 $fk4[$i] = 0.0;
836 for ($j=0 ; $j<3 ; $j++) {
837 $fk4[$i] += $etog[$j][$i] * $r[$j];
838 }
839 }
840
841 # Allow for e-terms */
842 for ($i=0 ; $i<3 ; $i++) {
843 $w += $r[$i] * $eterm[$i];
844 }
845 $w += 1.0;
846 for ($i=0 ; $i<3 ; $i++) {
847 $fk4[$i] = ($fk4[$i] + $eterm[$i])/$w;
848 }
849
850 if ($rect) {
851 return @fk4;
852 } else {
853 return r2pol(@fk4);
854 }
855}
856
857sub galfk4r(@) {galfk4(@_)};
858
859#=item B<fk4fk5>
860#
861# ($JRA, $JDec) = fk4fk5($BRA, $BDec);
862#
863# Converts an FK4 (B1950) position to the equivalent FK5 (J2000) position.
864# **LOW PRECISION**
865# $BRA,$BDec fk4/B1950 position (turns)
866# $JRA,$Dec fk5/J2000 position (turns)
867#
868#=cut
869#
870#sub fk4fk5 ($$) {
871# return r2pol(fk4fk5r(pol2r(shift,shift)));
872#}
873
874#=item B<fk5fk4>
875#
876# ($BRA, $BDec) = fk5fk4($JRA, $JDec);
877#
878# Converts an FK5 (J2000) position to the equivalent FK4 (B1950) position.
879# **LOW PRECISION**
880# $JRA,$Dec fk5/J2000 position (turns)
881# $BRA,$BDec fk4/B1950 position (turns)
882#
883#=cut
884#
885#sub fk5fk4 ($$) {
886# return r2pol(fk5fk4r(pol2r(shift,shift)));
887#}
888
889=item B<fk4gal>
890
891 ($l, $b) = fk4gal($ra, $dec);
892
893 Converts an FK4 position (B1950) to the IAU 1958 Galactic
894 coordinate system
895 ($ra, $dec) fk4 position to convert (turns)
896 ($l, $b) Galactic position (turns)
897 Reference : Blaauw et al., 1960, MNRAS, 121, 123.
898
899=cut
900
901sub fk4gal ($$) {
902 return r2pol(fk4galr(pol2r(shift,shift)));
903}
904
905#=item B<galfk4>
906#
907# ($ra, $dec) = galfk4($l, $b);
908#
909# Converts an IAU 1958 Galactic coordinate system position
910# to FK4 (B1950).
911# ($l, $b) Galactic position (turns)
912# ($ra, $dec) fk4 position to convert (turns)
913# Reference : Blaauw et al., 1960, MNRAS, 121, 123.
914#
915#=cut
916#
917#sub galfk4 ($$) {
918# return r2pol(galfk4r(pol2r(shift,shift)));
919#}
920
921=item B<ephem_vars>
922
923 ($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($jd)
924
925 Given the Julian day ($jd) this routine calculates the ephemeris
926 values required by the prcmat and nutate routines
927
928 The returned values are :
929 $omega - Longitude of the ascending node of the Moons mean orbit on
930 the ecliptic, measured from the mean equinox of date.
931 $rma - Mean anomaly of the Sun.
932 $mlanom - Mean anomaly of the Moon.
933 $F - L - omega, where L is the mean longitude of the Moon.
934 $D - Mean elongation of the Moon from the Sun.
935 $eps0 - Mean obilquity of the ecliptic.
936
937=cut
938
939=item B<J2000todate>
940
941
942 ($DRA, $DDec) = J2000todate($JRA, $JDec, $mjd);
943 @date = J2000todate(@J2000, $mjd);
944
945 Converts an J2000 position date coordinate
946
947 $DRA,$DDec Date coordinate (turns)
948 $JRA,$Dec J2000 position (turns)
949 @date Date coordinate (as a 3-vector)
950 @J2000 J2000 position (as a 3-vector)
951
952=cut
953
954# Untested
955sub J2000todate(@) {
956
957 my ($rect);
958 my (@J2000, @date); # Position vectors
959
960 my $mjd = pop @_;
961 if (@_==3) { # Rectangular coordinates passed
962 @J2000 = @_;
963 $rect = 1;
964 } elsif (@_==2) { # Sperical coordinates
965 @J2000 = pol2r($_[0],$_[1]); # Spherical to Cartesian
966 $rect = 0;
967 } elsif (@_>3) {
968 croak "Too many arguments for Astro::Coord::J2000todate ";
969 } else {
970 croak "Not enough arguments for Astro::Coord::J2000todate ";
971 }
972
973 # compute the general precession matrix.
974 my @gp = precsn(JULIAN_DAY_J2000, $mjd+2400000.5);
975
976 # Determine ephemeris quantities
977 my ($deps, $dpsi);
978 my @nu = ();
979 my ($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($mjd+2400000.5);
980 ($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
981
982 my @prcmat = ();
983 for (my $i=0 ; $i<3 ; $i++) {
984 for (my $j=0 ; $j<3 ; $j++) {
985 my $xx = 0.0;
986 for (my $k=0 ; $k<3 ; $k++) {
987 $xx = $xx + $gp[$i][$k] * $nu[$k][$j];
988 }
989 $prcmat[$i][$j] = $xx;
990 }
991 }
992
993 for (my $i=0 ; $i<3 ; $i++) {
994 $date[$i] = 0.0;
995 for (my $j=0 ; $j<3 ; $j++) {
996 $date[$i] += $prcmat[$i][$j] * $J2000[$j];
997 }
998 }
999
1000 if ($rect) {
1001 return @date;
1002 } else {
1003 # Revert to spherical coordinates
1004 return r2pol(@date);
1005 }
1006}
1007
1008sub ephem_vars ($) {
1009 my $epoch = shift;
1010
1011 # Calculates values required internally by prcmat and for nutate from
1012 # the passed Julian Day
1013
1014 # Calculate the interval to/from J2000 in Julian Centuries
1015 my $jcents = ($epoch-(JULIAN_DAY_J2000))/JULIAN_DAYS_IN_CENTURY;
1016
1017 # Calculate the longitude of the mean ascending node of the lunar
1018 # orbit on the ecliptic [A.A. Suppl. 1984, p S26]
1019 my $omega = (((0.000000039 * $jcents + 0.000036143) *
1020 $jcents - 33.757045934) *
1021 $jcents + 2.182438624)/(2.0*$PI);
1022 $omega = fmod($omega,1.0);
1023 if ($omega < 0.0) {
1024 $omega += 1.0;
1025 }
1026
1027 # Calculate the mean anomaly. [A.A suppl. 1984, p S26]
1028 my $manom = (6.240035939 - ((5.818e-8 * $jcents +2.797e-6) * $jcents -
1029 628.301956024) * $jcents)/(2.0*$PI);
1030 $manom = fmod($manom,1.0);
1031 if ($manom < 0.0) {
1032 $manom += 1.0;
1033 }
1034
1035 # Calculate the mean anomaly of the Moon. [A.A. Suppl, 1984, p S26]
1036 my $mlanom = (((0.000000310 * $jcents + 0.000151795) * $jcents
1037 +8328.691422884) * $jcents + 2.355548394)/(2.0*$PI);
1038 $mlanom = fmod($mlanom,1.0);
1039 if ($mlanom < 0.0) {
1040 $mlanom += 1.0;
1041 }
1042
1043 # Calculate the longitude of the moon from ascending node.
1044 # [A.A. Suppl, 1984, p S26]
1045 my $F = (((0.000000053 * $jcents - 0.000064272) * $jcents + 8433.466158318)
1046 * $jcents + 1.627901934)/(2.0*$PI);
1047 $F = fmod($F,1.0);
1048 if ($F < 0.0) {
1049 $F += 1.0;
1050 }
1051
1052 # Calculate the mean elongation of the moon from the sun.
1053 # [A.A suppl. 1984, p S26]
1054 my $D = (((0.000000092 * $jcents + 0.000033409) * $jcents + 7771.377146171)
1055 * $jcents + 5.198469514)/(2.0*$PI);
1056 $D = fmod($D,1.0);
1057 if ($D < 0.0) {
1058 $D += 1.0;
1059 }
1060
1061 # Calculate the mean obliquity of the ecliptic = mean obliquity.
1062 # [A.A suppl. 1984, p S26]
1063 my $eps0 = (((0.000000009 * $jcents - 0.000000003) * $jcents - 0.000226966)
1064 * $jcents + 0.409092804)/(2.0*$PI);
1065 return($omega, $manom, $mlanom, $F, $D, $eps0)
1066}
1067
1068=item B<nutate>
1069
1070 ($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
1071
1072 To calculate the nutation in longitude and obliquity according to
1073 the 1980 IAU Theory of Nutation including terms with amplitudes
1074 greater than 0.01 arcsecond. The nutation matrix is used to
1075 compute true place from mean place: true vector = N x mean place
1076 vector where the three components of each vector are the direction
1077 cosines wrt the mean equinox and equator.
1078
1079 / 1 -dpsi.cos(eps) -dpsi.sin(eps) \
1080 | |
1081 N = | dpsi.cos(eps) 1 -deps |
1082 | |
1083 \ dpsi.sin(eps) deps 1 /
1084
1085 The required inputs are : (NOTE: these are the values returned by ephem_vars)
1086 $omega - Longitude of the ascending node of the Moons mean orbit on
1087 the ecliptic, measured from the mean equinox of date.
1088 $rma - Mean anomaly of the Sun.
1089 $mlanom - Mean anomaly of the Moon.
1090 $F - L - omega, where L is the mean longitude of the Moon.
1091 $D - Mean elongation of the Moon from the Sun.
1092 $eps0 - Mean obilquity of the ecliptic.
1093
1094 The returned values are :
1095 $deps - nutation in obliquity
1096 $dpsi - nutation in longitude (scalar)
1097 @nu - nutation matrix (array [3][3])
1098
1099=cut
1100
1101sub nutate ($$$$$$) {
1102 my ($omega, $F, $D, $manom, $mlanom, $eps0) = @_;
1103
1104 my $arg1 = $omega;
1105 my $arg2 = 2.0 * $omega;
1106 my $arg9 = 2.0 * ($F-$D+$omega);
1107 my $arg10 = $manom;
1108 my $arg11 = $arg9 + $arg10;
1109 my $arg12 = $arg9 - $arg10;
1110 my $arg13 = $arg9 - $arg1;
1111 my $arg31 = 2.0 * ($F+$omega);
1112 my $arg32 = $mlanom;
1113 my $arg33 = $arg31 - $arg1;
1114 my $arg34 = $arg31 + $arg32;
1115 my $arg35 = $mlanom - 2.0 * $D;
1116 my $arg36 = $arg31 - $arg32;
1117
1118 my $dpsi = (-0.000083386 * sin($arg1*2.0*$PI)
1119 +0.000001000 * sin($arg2*2.0*$PI)
1120 -0.000006393 * sin($arg9*2.0*$PI)
1121 +0.000000691 * sin($arg10*2.0*$PI)
1122 -0.000000251 * sin($arg11*2.0*$PI)
1123 +0.000000105 * sin($arg12*2.0*$PI)
1124 +0.000000063 * sin($arg13*2.0*$PI)
1125 -0.000001102 * sin($arg31*2.0*$PI)
1126 +0.000000345 * sin($arg32*2.0*$PI)
1127 -0.000000187 * sin($arg33*2.0*$PI)
1128 -0.000000146 * sin($arg34*2.0*$PI)
1129 -0.000000077 * sin($arg35*2.0*$PI)
1130 +0.000000060 * sin($arg36*2.0*$PI))/(2.0*$PI);
1131
1132 my $deps = ( 0.000044615 * cos($arg1*2.0*$PI)
1133 -0.000000434 * cos($arg2*2.0*$PI)
1134 +0.000002781 * cos($arg9*2.0*$PI)
1135 +0.000000109 * cos($arg11*2.0*$PI)
1136 +0.000000474 * cos($arg31*2.0*$PI)
1137 +0.000000097 * cos($arg33*2.0*$PI)
1138 +0.000000063 * cos($arg34*2.0*$PI))/(2.0*$PI);
1139 my $eps = $eps0 + $deps;
1140
1141 my @N = ([1.0, -($dpsi)*(2.0*$PI)*cos($eps*2.0*$PI),
1142 -($dpsi)*(2.0*$PI)*sin($eps*2.0*$PI)],
1143 [0.0, 1.0, -($deps)*(2.0*$PI)],
1144 [0.0, ($deps)*(2.0*$PI), 1.0]);
1145 $N[1][0] = -1.0*$N[0][1];
1146 $N[2][0] = -1.0*$N[0][2];
1147 return($deps, $dpsi, @N);
1148}
1149
1150=item B<precsn>
1151
1152 @gp = precsn($jd_start, $jd_stop);
1153
1154 To calculate the precession matrix P for dates AFTER 1984.0 (JD =
1155 2445700.5) Given the position of an object referred to the equator
1156 and equinox of the epoch $jd_start its position referred to the
1157 equator and equinox of epoch $jd_stop can be calculated as follows :
1158
1159 1) Express the position as a direction cosine 3-vector (V1)
1160 (use pol2r to do this).
1161 2) The corresponding vector V2 for epoch jd_end is V2 = P.V1
1162
1163 The required inputs are :
1164 $jd_start - The Julian day of the current epoch of the coordinates.
1165 $jd_end - The Julian day at the required epoch for the conversion.
1166
1167 The returned values are :
1168 @gp - The required precession matrix (array [3][3])
1169
1170=cut
1171
1172sub precsn ($$) {
1173 my ($jd_start, $jd_end) = @_;
1174
1175 my @a = (0.011180860865024,
1176 0.000006770713945,
1177 -0.000000000673891,
1178 0.000001463555541,
1179 -0.000000001667759,
1180 0.000000087256766);
1181 my @b = (0.011180860865024,
1182 0.000006770713945,
1183 -0.000000000673891,
1184 0.000005307158404,
1185 0.000000000319977,
1186 0.000000088250634);
1187 my @d = (0.009717173455170,
1188 -0.000004136915141,
1189 -0.000000001052046,
1190 0.000002068457570,
1191 0.000000001052046,
1192 -0.000000202812107);
1193
1194 my $t = ($jd_start - JULIAN_DAY_J2000)/JULIAN_DAYS_IN_CENTURY;
1195 my $st = ($jd_end - $jd_start)/JULIAN_DAYS_IN_CENTURY;
1196 my $t2 = $t * $t;
1197 my $st2 = $st * $st;
1198 my $st3 = $st2 * $st;
1199
1200 # Calculate the Equatorial precession parameters
1201 # (ref. USNO Circular no. 163 1981,
1202 # Lieske et al., Astron. & Astrophys., 58, 1 1977)
1203
1204 my $zeta = ($a[0] + $a[1]*$t + $a[2]*$t2) * $st +
1205 ($a[3] + $a[4]*$t) * $st2 + $a[5] * $st3;
1206 my $z = ($b[0] + $b[1]*$t + $b[2]*$t2) * $st +
1207 ($b[3] + $b[4]*$t) * $st2 + $b[5] * $st3;
1208 my $theta = ($d[0] + $d[1]*$t + $d[2]*$t2) * $st -
1209 ($d[3] + $d[4]*$t) * $st2 + $d[5] * $st3;
1210
1211 # Calculate the P matrix
1212 my @precession = ([0.0, 0.0, 0.0],
1213 [0.0, 0.0, 0.0],
1214 [0.0, 0.0, 0.0]);
1215 $precession[0][0] = cos($zeta)*cos($z)*cos($theta) - sin($zeta)*sin($z);
1216 $precession[0][1] = -sin($zeta)*cos($z)*cos($theta) - cos($zeta)*sin($z);
1217 $precession[0][2] = -cos($z)*sin($theta);
1218 $precession[1][0] = cos($zeta)*sin($z)*cos($theta) + sin($zeta)*cos($z);
1219 $precession[1][1] = -sin($zeta)*sin($z)*cos($theta) + cos($zeta)*cos($z);
1220 $precession[1][2] = -sin($z)*sin($theta);
1221 $precession[2][0] = cos($zeta)*sin($theta);
1222 $precession[2][1] = -sin($zeta)*sin($theta);
1223 $precession[2][2] = cos($theta);
1224
1225 return(@precession);
1226}
1227
1228=item B<coord_convert>
1229
1230 ($output_left, $output_right) = coord_convert($input_left, $input_right,
1231 $input_mode, $output_mode,
1232 $mjd, $longitude, $latitude,
1233 $ref0);
1234
1235 A routine for converting between any of the following coordinate systems :
1236 Coordinate system input/output mode
1237 ----------------- -----------------
1238 X, Y (East-West mounted) 0
1239 Azimuth, Elevation 1
1240 Hour Angle, Declination 2
1241 Right Ascension, Declination (date, J2000 or B1950) 3,4,5
1242 Galactic (B1950) 6
1243
1244 The last four parameters in the call ($mjd, $longitude, $latitude
1245 and $ref0) are not always required for the coordinate conversion.
1246 In particular if the conversion is between two coordinate systems
1247 which are fixed with respect to the celestial sphere (RA/Dec J2000,
1248 B1950 or Galactic), or two coordinate systems which are fixed with
1249 respect to the antenna (X/Y and Az/El) then these parameters are not
1250 used (NOTE: they must always be passed, even if they only hold 0 or
1251 undef as the routine will return undef if it is not passed 8
1252 parameters). The RA/Dec date system is defined with respect to the
1253 celestial sphere, but varies with time. The table below shows which
1254 of $mjd, $longitude, $latitude and $ref0 are used for a given
1255 conversion. If in doubt you should determing the correct values for
1256 all input parameters, no checking is done in the routine that the
1257 passed values are sensible.
1258
1259 Conversion $mjd $longitude $latitude $ref0
1260 ------------------------------------------------------------------------
1261 Galactic, Galactic,
1262 RA/Dec J2000,B1950 <->RA/Dec J2000, B1950 N N N N
1263
1264 Galactic,
1265 RA/Dec J2000,B1950 <->RA/Dec date Y N N N
1266
1267 Galactic,
1268 RA/Dec J2000,B1950,<->HA/Dec Y Y N N
1269 date
1270
1271 Galactic,
1272 RA/Dec J2000,B1950,<->X/Y, Az/El Y Y Y Y
1273 date
1274
1275 X/Y, Az/El <->X/Y, Az/El N N N N
1276
1277 X/Y, Az/El <->HA/Dec N N Y Y
1278
1279
1280 NOTE : The method used for refraction correction is asymptotic at
1281 an elevation of 0 degrees.
1282
1283 The required inputs are :
1284 $input_left - The left/longitude input coordinate (turns)
1285 $input_right - The right/latitude input coordinate (turns)
1286 $input_mode - The mode of the input coordinates (0-6)
1287 $output_mode - The mode to convert the coordinates to.
1288 $mjd - The time as modified Julian day (if necessary) at
1289 which to perform the conversion
1290 $longitude - The longitude of the location/observatory (if necessary)
1291 at which to perform the conversion (turns)
1292 $latitude - The latitude of the location/observatory (if necessary)
1293 at which to perform the conversion (turns)
1294 $ref0 - The refraction constant (if in doubt use 0.00005).
1295
1296 The returned values are :
1297 $output_left - The left/longitude output coordinate (turns)
1298 $output_right - The right/latitude output coordinate (turns)
1299
1300=cut
1301
1302sub coord_convert ($$$$;$$$$) {
1303 my ($input_left, $input_right, $input_mode, $output_mode, $mjd, $longitude,
1304 $latitude, $ref0) = @_;
1305
1306 # Some required constants
1307 my ($EWXY, $AZEL, $HADEC, $DATE, $J2000, $B1950, $GALACTIC) = 0..6;
1308
1309 # First check what the input and output modes are.
1310 if (($input_mode < $EWXY) || ($input_mode > $GALACTIC)) {
1311 carp "Invalid input coordinate mode : $input_mode\n".
1312 "Valid inputs are numbers in the range 0-6, which corrspond to X/Y, ".
1313 "Az/El,\n HA/Dec, RA/Dec (date), RA/Dec (J2000), RA/Dec (B1950), ".
1314 "Galactic.";
1315 return undef;
1316 }
1317 if (($output_mode < $EWXY) || ($output_mode > $GALACTIC)) {
1318 carp "Invalid output coordinate mode : $output_mode\n".
1319 "Valid outputs are numbers in the range 0-6, which corrspond to X/Y, ".
1320 "Az/El,\n HA/Dec, RA/Dec (date), RA/Dec (J2000), RA/Dec (B1950), ".
1321 "Galactic.";
1322 return undef;
1323 }
1324
1325 # Check we have the correct parameters passed
1326
1327 # Need mjd
1328 if ((($input_mode>=$DATE && $output_mode<=$DATE) ||
1329 ($input_mode<=$DATE && $output_mode>=$DATE)) &&
1330 !(defined($mjd))) {
1331 carp '$mjd parametr missing';
1332 return undef;
1333 }
1334 # Need longitude
1335 if ((($input_mode>=$HADEC && $output_mode<=$AZEL) ||
1336 ($input_mode<=$HADEC && $output_mode>=$HADEC)) &&
1337 !(defined($longitude))) {
1338 carp '$longitude parametr missing';
1339 return undef;
1340 }
1341 # Need latitude
1342 if ((($input_mode>=$HADEC && $output_mode<$HADEC) ||
1343 ($input_mode<=$AZEL && $output_mode>$AZEL)) &&
1344 !(defined($latitude))) {
1345 carp '$latitude parameter missing';
1346 return undef;
1347 }
1348 # Need ref0
1349 if ((($input_mode>=$HADEC && $output_mode<$HADEC) ||
1350 ($input_mode<=$AZEL && $output_mode>$AZEL)) &&
1351 !(defined($ref0))) {
1352 carp '$ref0 parameter missing';
1353 return undef;
1354 }
1355
1356 # If necessary determine ephemeris quantities (if either of the modes are
1357 # date, HA/Dec, AzEl or EWXY).
1358 my ($omega, $rma, $mlanom, $F, $D, $eps0, $deps, $dpsi);
1359 my @nu = ();
1360
1361 if (($input_mode<=$DATE && $output_mode>=$DATE) ||
1362 ($input_mode>=$DATE && $output_mode<=$DATE)) {
1363 ($omega, $rma, $mlanom, $F, $D, $eps0) = ephem_vars($mjd+2400000.5);
1364 ($deps, $dpsi, @nu) = nutate($omega, $F, $D, $rma, $mlanom, $eps0);
1365 }
1366
1367 my @vonc = ();
1368 if (($input_mode<=$HADEC && $output_mode>=$DATE) ||
1369 ($input_mode>=$DATE && $output_mode<=$HADEC)) {
1370
1371 # Calculate the interval to/from J2000 in Julian Centuries
1372 my $jcents = ($mjd+2400000.5-(JULIAN_DAY_J2000))/JULIAN_DAYS_IN_CENTURY;
1373
1374 # Compute the eccentricity of the Earth's orbit (in radians)
1375 # [Explanatory supplement to the Astronomical Ephemeris 1961, p 98]
1376 my $e = (-0.000000126 * $jcents - 0.00004205) * $jcents + 0.016709114;
1377
1378 # Compute the eccentric anomaly, by iteratively solving :
1379 # ea = e*sin(ea) - rma
1380 my $ea = $rma;
1381 my $xx;
1382 do {
1383 $xx = $ea;
1384 $ea = $xx + ($rma - $xx + $e*sin($xx)) / (1.0 - $e*cos($xx));
1385 } while (abs($ea -$xx) > 1.0e-9);
1386
1387 # Compute the mean longitude of perihelion, in radians
1388 # (reference as for `e').
1389 my $perihl = ((0.00000005817764*$jcents + 0.000008077) * $jcents
1390 + 0.030010190) * $jcents + 1.796613066;
1391
1392 # Calculate the equation of the equinoxes
1393 #my $eqenx = $dpsi * cos(($eps0+$deps)*2.0*$PI);
1394
1395 # Compute the abberation vector
1396 my $eps = $eps0 + $deps;
1397 $xx = 0.00009936508 / (1.0 - $e*cos($ea));
1398 my $efac = sqrt(1.0 - $e*$e);
1399 $vonc[0] = $xx * (-cos($perihl)*sin($ea) - $efac*sin($perihl)*cos($ea));
1400 $vonc[1] = $xx * (-sin($perihl)*cos($eps)*sin($ea) +
1401 $efac*cos($perihl)*cos($eps)*cos($ea));
1402 $vonc[2] = $xx * (-sin($perihl)*sin($eps)*sin($ea) +
1403 $efac*cos($perihl)*sin($eps)*cos($ea));
1404
1405 }
1406
1407 my @prcmat = ();
1408 if (($input_mode<=$DATE && $output_mode>=$J2000) ||
1409 ($input_mode>=$J2000 && $output_mode<=$DATE)) {
1410
1411 # compute the general precession matrix. */
1412 my @gp = precsn(JULIAN_DAY_J2000, $mjd+2400000.5);
1413
1414 # The matrices returned from nutate (nu) and precsn (gp) can be used
1415 # to convert J2000 coordinates to date by :
1416 # (coords at date) = gp * nu * (coords at J2000)
1417 # gp and nu can be combined to give the required precession matrix
1418
1419 for (my $i=0 ; $i<3 ; $i++) {
1420 for (my $j=0 ; $j<3 ; $j++) {
1421 my $xx = 0.0;
1422 for (my $k=0 ; $k<3 ; $k++) {
1423 $xx = $xx + $gp[$i][$k] * $nu[$k][$j];
1424 }
1425 $prcmat[$i][$j] = $xx;
1426 }
1427 }
1428 }
1429
1430 my $lmst;
1431 if (($input_mode<=$HADEC && $output_mode>=$DATE) ||
1432 ($output_mode<=$HADEC && $input_mode>=$DATE)) {
1433 $lmst = mjd2lst($mjd, $longitude);
1434 }
1435
1436 # Perform the conversion
1437 my (@lb, @b1950, @j2000, @date, $ra, $ha, $dec, $az, $el, $x, $y);
1438 if ($input_mode == $GALACTIC) {
1439 @lb = pol2r($input_left, $input_right);
1440 } elsif ($input_mode == $B1950) {
1441 @b1950 = pol2r($input_left, $input_right);
1442 } elsif ($input_mode == $J2000) {
1443 @j2000 = pol2r($input_left, $input_right);
1444 } elsif ($input_mode == $DATE) {
1445 @date = pol2r($input_left, $input_right);
1446 } elsif ($input_mode == $HADEC) {
1447 $ha = $input_left;
1448 $dec = $input_right;
1449 } elsif ($input_mode == $AZEL) {
1450 $az = $input_left;
1451 $el = $input_right;
1452 } else {
1453 $x = $input_left;
1454 $y = $input_right;
1455 }
1456
1457 # Conversion is to a "lower" mode
1458 if ($output_mode < $input_mode) {
1459
1460 # Convert from Galactic to B1950
1461 if ($input_mode == $GALACTIC) {
1462 @b1950 = galfk4r(@lb);
1463 }
1464
1465 # Convert from B1950 to J2000
1466 if (($input_mode >= $B1950) && ($output_mode < $B1950)) {
1467 @j2000 = fk4fk5r(@b1950);
1468 }
1469
1470 # Precess from J2000 to date
1471 if (($input_mode >= $J2000) && ($output_mode < $J2000)) {
1472 for (my $i=0 ; $i<3 ; $i++) {
1473 $date[$i] = 0.0;
1474 for (my $j=0 ; $j<3 ; $j++) {
1475 $date[$i] += $prcmat[$i][$j] * $j2000[$j];
1476 }
1477 }
1478 }
1479
1480 # Convert from date to HA/Dec
1481 if (($input_mode >= $DATE) && ($output_mode < $DATE)) {
1482
1483 # Convert to geometrical equitorial coordinates
1484 for (my $i=0 ; $i<3 ; $i++) {
1485 $date[$i] += $vonc[$i];
1486 }
1487
1488 # Convert from retangular back to polar coordinates
1489 ($ra, $dec) = r2pol(@date);
1490
1491 # Convert to hour angle
1492 $ha = $lmst - $ra;
1493 if ($ha < 0.0) {
1494 $ha += 1.0;
1495 }
1496 }
1497
1498 # Convert from HA/Dec to Az/El
1499 if (($input_mode >= $HADEC) && ($output_mode < $HADEC)) {
1500 ($az, $el) = eqazel($ha, $dec, $latitude);
1501
1502 # Correct for refraction
1503 $el += $ref0/tan($el*2.0*$PI);
1504 }
1505
1506 # Convert from Az/El to X/Y
1507 if (($input_mode >= $AZEL) && ($output_mode < $AZEL)) {
1508 ($x, $y) = azel2xy($az, $el);
1509 }
1510 } else {
1511 # Convert from X/Y to Az/El
1512 if (($input_mode == $EWXY) && ($output_mode > $EWXY)) {
1513 ($az, $el) = xy2azel($x, $y);
1514 }
1515
1516 # Convert from Az/El to HA/Dec
1517 if (($input_mode <= $AZEL) && ($output_mode > $AZEL)) {
1518
1519 # First numerically invert the refraction correction
1520 my $upper = $el - $ref0/tan($el*2.0*$PI);
1521 my $lower = $el - 1.5*$ref0/tan($el*2.0*$PI);
1522 my $root = ($lower+$upper)/2.0;
1523 my $niter = 0;
1524 do {
1525 if ($root + $ref0/tan($root*2.0*$PI) - $el > 0.0) {
1526 $upper = $root;
1527 } else {
1528 $lower = $root;
1529 }
1530 $root = ($lower+$upper)/2.0;
1531 $niter++;
1532 } while (($niter <= 10) && (($upper-$root) > 7.0e-8));
1533 $el = $root;
1534
1535 # Now do the conversion
1536 ($ha, $dec) = eqazel($az, $el, $latitude);
1537 }
1538
1539 # Convert from HA/Dec to date
1540 if (($input_mode <= $HADEC) && ($output_mode > $HADEC)) {
1541 $ra = $lmst - $ha;
1542 if ($ra < 0.0) {
1543 $ra += 1.0;
1544 }
1545 @date = pol2r($ra, $dec);
1546
1547 # Remove the abberation vector
1548 for (my $i=0 ; $i<3 ; $i++) {
1549 $date[$i] -= $vonc[$i];
1550 }
1551 }
1552
1553 # precess from date to J2000
1554 if (($input_mode <= $DATE) && ($output_mode > $DATE)) {
1555 for (my $i=0 ; $i<3 ; $i++) {
1556 $j2000[$i] = 0.0;
1557 for (my $j=0 ; $j<3 ; $j++) {
1558 $j2000[$i] += $prcmat[$j][$i] * $date[$j];
1559 }
1560 }
1561 }
1562
1563 # Convert from J2000 to B1950
1564 if (($input_mode <= $J2000) && ($output_mode > $J2000)) {
1565 @b1950 = fk5fk4(@j2000);
1566 }
1567
1568 # Convert from B1950 to Galactic
1569 if (($input_mode <= $B1950) && ($output_mode >= $B1950)) {
1570 @lb = fk4galr(@b1950);
1571 }
1572 }
1573
1574 if ($output_mode == $EWXY) {
1575 return($x, $y);
1576 } elsif ($output_mode == $AZEL) {
1577 return($az, $el);
1578 } elsif ($output_mode == $HADEC) {
1579 return($ha, $dec);
1580 } elsif ($output_mode == $DATE) {
1581 return(r2pol(@date));
1582 } elsif ($output_mode == $J2000) {
1583 return(r2pol(@j2000));
1584 } elsif ($output_mode == $B1950) {
1585 return(r2pol(@b1950));
1586 } elsif ($output_mode == $GALACTIC) {
1587 return(r2pol(@lb));
1588 }
1589}
1590
1591=item B<haset_ewxy>
1592
1593 $haset = haset_ewxy($declination, $latitude, %limits);
1594
1595 This routine takes the $declination of the source, and the $latitude of the
1596 EWXY mounted antenna and calculates the hour angle at which the source
1597 will set. It is then trivial to calculate the time until the source
1598 sets, simply by subtracting the current hour angle of the source from
1599 the hour angle at which it sets.
1600
1601 The required inputs are :
1602 $declination - The declination of the source (turns)
1603 $latitude - The latitude of the observatory (turns)
1604 %limits - A reference to a hash holding the EWXY antenna limits
1605 The following keys must be defined XLOW, XLOW_KEYHOLE,
1606 XHIGH, XHIGH_KEYHOLE, YLOW, YLOW_KEYHOLE, YHIGH,
1607 YHIGH_KEYHOLE (all values shoule be in turns)
1608
1609 The returned value is :
1610 $haset - The hour angle (turns) at which a source at this
1611 declination sets for an EWXY mounted antenna with the
1612 given limits at the given latitude
1613
1614 NOTE: returns undef if %limits hash is missing any of the required keys
1615
1616=cut
1617
1618sub haset_ewxy($$\%) {
1619
1620 my ($declination, $latitude, $limitsref) = @_;
1621
1622 # Check that all the required keys are present
1623 if ((!exists $limitsref->{XLOW}) || (!exists $limitsref->{XLOW_KEYHOLE}) ||
1624 (!exists $limitsref->{XHIGH}) || (!exists $limitsref->{XHIGH_KEYHOLE}) ||
1625 (!exists $limitsref->{YLOW}) || (!exists $limitsref->{YLOW_KEYHOLE}) ||
1626 (!exists $limitsref->{YHIGH}) || (!exists $limitsref->{YHIGH_KEYHOLE})) {
1627 carp 'Missing key in %limits';
1628 return undef;
1629 }
1630
1631 # Local variables
1632 my ($pole, $pxlim, $exlim, $hix, $hixk, $lowx, $lowxk);
1633
1634 if ($latitude < 0.0) {
1635 $pole = -90.0/360.0;
1636 $pxlim = $limitsref->{XLOW};
1637 $exlim = $limitsref->{XHIGH};
1638 } else {
1639 $pole = 90.0/360.0;
1640 $pxlim = $limitsref->{XHIGH};
1641 $exlim = $limitsref->{XLOW};
1642 }
1643 my $dec_never = $latitude + $exlim;
1644 my $dec_always = $pole - ($latitude + $pxlim - $pole);
1645
1646 if ((($latitude < 0.0) && ($declination > $dec_never)) ||
1647 (($latitude > 0.0) && ($declination < $dec_never))) {
1648
1649 # Source is never up
1650 return(0.0);
1651 } elsif ((($latitude < 0.0) && ($declination < $dec_always)) ||
1652 (($latitude > 0.0) && ($declination > $dec_always))) {
1653
1654 # Source is always up
1655 return(1.0);
1656 } else {
1657
1658 # Up some of the time - calculate the ghastly constants and
1659 # do everything in radians from here on.
1660 $declination = 2.0 * $PI * $declination;
1661 $latitude = 2.0 * $PI * $latitude;
1662 my $k0 = -cos($declination);
1663 my $k1 = sin($declination)*cos($latitude);
1664 my $k2 = sin($declination)*sin($latitude);
1665 my $k3 = cos($declination)*sin($latitude);
1666 my $k4 = cos($declination)*cos($latitude);
1667 my $k5 = $k4 * $k1 + $k2 * $k3;
1668 my $x = 2.0 * $PI * $limitsref->{XLOW_KEYHOLE};
1669 my $dec_split = asin(cos(2.0 * $PI * $limitsref->{YLOW}) *
1670 (cos($x) * sin($latitude) + sin($x) *
1671 cos($latitude)));
1672 if ($latitude > 0.0) {
1673
1674 # Set up for northern antenna
1675 $hix = $limitsref->{XLOW};
1676 $hixk = $limitsref->{XLOW_KEYHOLE};
1677 $lowx = $limitsref->{XHIGH};
1678 $lowxk = $limitsref->{XHIGH_KEYHOLE};
1679
1680 } else {
1681
1682 # Set up for southern antenna
1683 $hix = $limitsref->{XHIGH};
1684 $hixk = $limitsref->{XHIGH_KEYHOLE};
1685 $lowx = $limitsref->{XLOW};
1686 $lowxk = $limitsref->{XLOW_KEYHOLE};
1687 }
1688
1689 if ((($declination > $dec_split) && ($latitude < 0.0)) ||
1690 (($declination < $dec_split) && ($latitude > 0.0))) {
1691
1692 # We are on the equatorial side
1693 my $x = 2.0 * $PI * $hix;
1694 my $y = -1.0 * abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
1695 if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
1696 return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y))/
1697 ($k3 + $k4))/(2.0 * $PI));
1698 } else {
1699 my $x = 2.0 * $PI * $hixk;
1700 my $y = -1.0 * abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
1701 if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
1702 return(asin(sin(2.0 * $PI * $limitsref->{YLOW_KEYHOLE}) /
1703 $k0)/(2.0 * $PI));
1704 } elsif (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW})) {
1705 return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
1706 ($k3 + $k4))/(2.0 * $PI));
1707 } else {
1708 return(asin(sin(2.0 * $PI*$limitsref->{YLOW}) / $k0) /
1709 (2.0 * $PI));
1710 }
1711 }
1712 } else {
1713
1714 # We are on the polar side
1715 my $x = 2.0 * $PI * $lowx;
1716 my $y = abs(acos($k5 / ($k4 * sin($x) + $k3 * cos($x))));
1717 if (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW_KEYHOLE})) {
1718 return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
1719 ($k3 + $k4))/(2.0 * $PI));
1720 } else {
1721 my $x = 2.0 * $PI * $lowxk;
1722 my $y = -1.0 * abs(acos($k5 /($k4 * sin($x) + $k3 * cos($x))));
1723 if (abs($y) < abs(2.0 * $PI* $limitsref->{YLOW_KEYHOLE})) {
1724 return(asin(sin(2.0 * $PI * $limitsref->{YLOW_KEYHOLE}) /
1725 $k0)/(2.0 * $PI));
1726 } elsif (abs($y) < abs(2.0 * $PI * $limitsref->{YLOW})) {
1727 return(acos(($k1 - $k2 + cos($x) * cos($y) - sin($x) * cos($y)) /
1728 ($k3 + $k4))/(2.0 * $PI));
1729 } else {
1730 return(asin(sin(2.0 * $PI * $limitsref->{YLOW}) / $k0)/
1731 (2.0 * $PI));
1732 }
1733 }
1734 }
1735 }
1736}
1737
1738=item B<ewxy_tlos>
1739
1740 $tlos = ewxy_tlos($hour_angle, $declination, $latitude, %limits);
1741
1742 This routine calculates the time left on-source (tlos) for a source
1743 at $hour_angle, $declination for an EWXY mount antenna at $latitude.
1744
1745 The required inputs are :
1746 $hour_angle - The current hour angle of the source (turns)
1747 $declination - The declination of the source (turns)
1748 $latitude - The latitude of the observatory (turns)
1749 \%limits - A reference to a hash holding the EWXY antenna limits
1750 The following keys must be defined XLOW, XLOW_KEYHOLE,
1751 XHIGH, XHIGH_KEYHOLE, YLOW, YLOW_KEYHOLE, YHIGH,
1752 YHIGH_KEYHOLE (all values should be in turns)
1753
1754 The returned value is :
1755 $tlos - The time left on-source (turns)
1756
1757=cut
1758
1759sub ewxy_tlos($$$\%) {
1760
1761 my ($hour_angle, $declination, $latitude, $limitsref) = @_;
1762
1763 my $haset = haset_ewxy($declination, $latitude, %$limitsref);
1764 return(undef) if (!defined $haset);
1765 $haset -= $hour_angle if (($haset > 0.0) && ($haset < 1.0));
1766 $haset += 1.0 if ($haset < 0.0);
1767
1768 return $haset;
1769}
1770
1771=item B<haset_azel>
1772
1773 $haset = haset_azel($declination, $latitude, %limits);
1774
1775 This routine takes the $declination of the source, and the
1776 $latitude of the Az/El mounted antenna and calculates the hour
1777 angle at which the source will set. It is then trivial to
1778 calculate the time until the source sets, simply by subtracting the
1779 current hour angle of the source from the hour angle at which it
1780 sets. This routine assumes that the antenna is able to rotate
1781 through 360 degrees in azimuth.
1782
1783 The required inputs are :
1784 $declination - The declination of the source (turns)
1785 $latitude - The latitude of the observatory (turns)
1786 \%limits - A reference to a hash holding the Az/El antenna limits
1787 The following keys must be defined ELLOW (all values should
1788 be in turns)
1789
1790 The returned value is :
1791 $haset - The hour angle (turns) at which a source at this
1792 declination sets for an Az/El mounted antenna with the
1793 given limits at the given latitude
1794
1795 NOTE: returns undef if the %limits hash is missing any of the required keys
1796
1797=cut
1798
1799sub haset_azel($$\%) {
1800
1801 my ($declination, $latitude, $limitsref) = @_;
1802
1803 # Check that all the required keys are present
1804 if (!exists $limitsref->{ELLOW}) {
1805 carp 'Missing key in %limits';
1806 return undef ;
1807 }
1808
1809 my $cos_haset = (cos($PI / 2.0 - $limitsref->{ELLOW} * 2.0 *
1810 $PI) - sin($latitude * 2.0 * $PI) *
1811 sin($declination * 2.0 * $PI))/
1812 (cos($declination * 2.0 * $PI)
1813 *cos($latitude * 2.0 * $PI));
1814 if ($cos_haset > 1.0) {
1815
1816 # The source never rises
1817 return(0.0);
1818 } elsif ($cos_haset < -1.0) {
1819
1820 # The source never sets
1821 return(1.0);
1822 } else {
1823
1824 return(acos($cos_haset)/(2.0*$PI));
1825 }
1826}
1827
1828=item B<azel_tlos>
1829
1830 $tlos = azel_tlos($hour_angle, $declination, $latitude, \%limits);
1831
1832 This routine calculates the time left on-source (tlos) for a source
1833 at $hour_angle, $declination for an Az/El mount antenna at $latitude.
1834
1835 The required inputs are :
1836 $hour_angle - The current hour angle of the source (turns)
1837 $declination - The declination of the source (turns)
1838 $latitude - The latitude of the observatory (turns)
1839 %limits - A reference to a hash holding the Az/El antenna limits
1840 The following keys must be defined ELLOW (all values
1841 should be in turns)
1842
1843 The returned value is :
1844 $tlos - The time left on-source (turns)
1845
1846=cut
1847
1848sub azel_tlos($$$\%) {
1849 my ($hour_angle, $declination, $latitude, $limitsref) = @_;
1850
1851 # Calculate the time left onsource
1852 my $haset = haset_azel($declination, $latitude, %$limitsref);
1853 if (!defined $haset) {return(undef)};
1854 if (($haset > 0.0) && ($haset < 1.0)) { $haset -= $hour_angle; }
1855 if ($haset < 0.0) { $haset += 1.0; }
1856
1857 return($haset);
1858}
1859
1860=item B<antenna_rise>
1861
1862 $ha_set = antenna_rise($declination, $latitude, $mount, \%limits);
1863
1864 Given the $declination of the source, the $latitude of the antenna,
1865 the type of the antenna $mount and a reference to a hash holding
1866 information on the antenna limits, this routine calculates the hour
1867 angle at which the source sets for the antenna. The hour angle at
1868 which it rises is simply the negative of that at which it sets.
1869 These values in turn can be used to calculate the LMST at which the
1870 source rises/sets and from that the UT at which the source
1871 rises/sets on a given day, or to calculate the native coordinates
1872 at which the source rises/sets.
1873
1874 If you want to calculate source rise times above arbitrary elevation,
1875 use the routine rise.
1876
1877 The required inputs are :
1878 $declination - The declination of the source (turns)
1879 $latitude - The latitude of the observatory (turns)
1880 $mount - The type of antenna mount, 0 => EWXY mount, 1 => Az/El,
1881 any other number will cause the routine to return
1882 undef
1883 %limits - A reference to a hash holding the antenna limits
1884 For an EWXY antenna there must be keys for all the
1885 limits (i.e. XLOW, XLOW_KEYHOLE, XHIGH, XHIGH_KEYHOLE,
1886 YLOW, YLOW_KEYHOLE, YHIGH, YHIGH_KEYHOLE). For an Az/El
1887 antenna there must be a key for ELLOW (all values should
1888 be in turns).
1889
1890 The returned values are :
1891 $ha_set - The hour angle at which the source sets (turns). The hour
1892 angle at which the source rises is simply the negative of this
1893 value.
1894
1895=cut
1896
1897sub antenna_rise($$$$) {
1898
1899 my ($declination, $latitude, $mount, $limitsref) = @_;
1900
1901 # Check that the mount type is either EWXY (0) or AZEL (1)
1902 if (($mount != 0) && ($mount != 1)) {
1903 carp 'mount must equal 0 or 1';
1904 return undef;
1905 }
1906
1907 if ($mount == 0) {
1908 return(haset_ewxy($declination, $latitude, %$limitsref));
1909 } elsif ($mount == 1) {
1910 return(haset_azel($declination, $latitude, %$limitsref));
1911 }
1912}
1913
191412.0e-62.0e-6my @b2g = ([-0.054875539726, 0.494109453312, -0.867666135858],
1915 [-0.873437108010, -0.444829589425, -0.198076386122],
1916 [-0.483834985808, 0.746982251810, 0.455983795705]);
1917
1918#my @b2g = ([ -0.0548777621, +0.4941083214, -0.8676666398],
1919# [ -0.8734369591, -0.4448308610, -0.1980741871],
1920# [ -0.4838350026, +0.7469822433, +0.4559837919]);
1921
1922sub j2gal($$) {
1923 my ($ra,$dec) = @_;
1924 my @r = pol2r($ra,$dec);
1925 my @g = (0,0,0);
1926 for (my $i=0; $i<3; $i++) {
1927 for (my $j=0; $j<3; $j++) {
1928 $g[$i]+= $b2g[$j][$i] * $r[$j];
1929 }
1930 }
1931 return r2pol(@g);
1932}
1933
1934# SLALIB support routines
1935
1936sub epb2d ($) {
1937# - - - - - -
1938# E P B 2 D
1939# - - - - - -
1940#
1941# Conversion of Besselian Epoch to Modified Julian Date
1942# (double precision)
1943#
1944# Given:
1945# EPB dp Besselian Epoch
1946#
1947# The result is the Modified Julian Date (JD - 2400000.5).
1948#
1949# Reference:
1950# Lieske,J.H., 1979. Astron.Astrophys.,73,282.
1951#
1952# P.T.Wallace Starlink February 1984
1953#
1954# Copyright (C) 1995 Rutherford Appleton Laboratory
1955
1956 my $epb = shift;
1957
1958 return 15019.81352 + ($epb-1900)*365.242198781;
1959}
1960
1961sub epj ($) {
1962# - - - -
1963# E P J
1964# - - - -
1965#
1966# Conversion of Modified Julian Date to Julian Epoch (double precision)
1967#
1968# Given:
1969# DATE dp Modified Julian Date (JD - 2400000.5)
1970#
1971# The result is the Julian Epoch.
1972#
1973# Reference:
1974# Lieske,J.H., 1979. Astron.Astrophys.,73,282.
1975#
1976# P.T.Wallace Starlink February 1984
1977#
1978# Copyright (C) 1995 Rutherford Appleton Laboratory
1979 my $date = shift;
1980
1981 return 2000 + ($date-51544.5)/365.25;
1982}
1983
1984sub pm ($$$$$$$$$$) {
1985# - - -
1986# P M
1987# - - -
1988#
1989# Apply corrections for proper motion to a star RA,Dec
1990# (double precision)
1991#
1992# References:
1993# 1984 Astronomical Almanac, pp B39-B41.
1994# (also Lederle & Schwan, Astron. Astrophys. 134,
1995# 1-6, 1984)
1996#
1997# Given:
1998# R0,D0 dp RA,Dec at epoch EP0 (rad)
1999# PR,PD dp proper motions: RA,Dec changes per year of epoch
2000# EP0 dp start epoch in years (e.g. Julian epoch)
2001# EP1 dp end epoch in years (same system as EP0)
2002#
2003# Returned:
2004# R1,D1 dp RA,Dec at epoch EP1 (rad)
2005#
2006# Called:
2007# sla_DCS2C spherical to Cartesian
2008# sla_DCC2S Cartesian to spherical
2009# sla_DRANRM normalize angle 0-2Pi
2010#
2011# Note:
2012# The proper motions in RA are dRA/dt rather than
2013# cos(Dec)*dRA/dt, and are in the same coordinate
2014# system as R0,D0.
2015#
2016# P.T.Wallace Starlink 23 August 1996
2017#
2018# Copyright (C) 1996 Rutherford Appleton Laboratory
2019
2020 my ($r0, $d0, $pr, $pd, $ep0, $ep1) = @_;
2021
2022 # Km/s to AU/year multiplied by arc seconds to radians
202330.000144.5e-5 use constant VFR => 0.21094502*0.484813681109535994e-5;
# spent 60µs making 1 call to constant::import
2024
2025 my (@em, $t);
2026
2027 # Spherical to Cartesian
2028 my @p = pol2r($r0,$d0);
2029
2030 # Space motion (radians per year)
2031 $em[0]=-$pr*$p[1]-$pd*cos($r0)*sin($d0);
2032 $em[1]= $pr*$p[0]-$pd*sin($r0)*sin($d0);
2033 $em[2]= $pd*cos($d0);
2034
2035 # Apply the motion
2036 $t=$ep1-$ep0;
2037 for (my $i = 0; $i<3; $i++) {
2038 $p[$i]=$p[$i]+$t*$em[$i];
2039 }
2040
2041 # Cartesian to spherical
2042 return r2pol(@p);
2043}
2044
2045
204611.6e-51.6e-51;
2047
2048__END__
2049

Report produced by the NYTProf 2.07 Perl profiler, developed by Tim Bunce and Adam Kaplan.