← 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:45 2010

File/wise/base/static/lib/perl5/site_perl/5.10.0/Astro/Time.pm
Statements Executed277411
Total Time0.397828000000279 seconds

Subroutines — ordered by exclusive time
Calls P F Exclusive
Time		 Inclusive
Time		 Subroutine
69345310.563330.56333Astro::Time::turn2rad
00000Astro::Time::BEGIN
00000Astro::Time::cal2dayno
00000Astro::Time::cal2lst
00000Astro::Time::cal2mjd
00000Astro::Time::dayOK
00000Astro::Time::dayno2cal
00000Astro::Time::dayno2lst
00000Astro::Time::dayno2mjd
00000Astro::Time::daynoOK
00000Astro::Time::deg2rad
00000Astro::Time::deg2str
00000Astro::Time::deg2turn
00000Astro::Time::gst
00000Astro::Time::hms2time
00000Astro::Time::jd2mjd
00000Astro::Time::leap
00000Astro::Time::lst2mjd
00000Astro::Time::mjd2cal
00000Astro::Time::mjd2dayno
00000Astro::Time::mjd2jd
00000Astro::Time::mjd2lst
00000Astro::Time::mjd2time
00000Astro::Time::monthOK
00000Astro::Time::nint
00000Astro::Time::now2mjd
00000Astro::Time::rad2deg
00000Astro::Time::rad2str
00000Astro::Time::rad2turn
00000Astro::Time::rise
00000Astro::Time::str2deg
00000Astro::Time::str2rad
00000Astro::Time::str2turn
00000Astro::Time::time2hms
00000Astro::Time::tomorrow
00000Astro::Time::turn2deg
00000Astro::Time::turn2str
00000Astro::Time::utOK
00000Astro::Time::yesterday

LineStmts.Exclusive
Time		Avg.Code
1package Astro::Time;
238.5e-52.8e-5use strict;
# spent 10µs making 1 call to strict::import
3
4=head1 NAME
5
6Astro::Time - Time based astronomical routines
7
8=head1 SYNOPSIS
9
10 use Astro::Time;
11
12 $dayno = cal2dayno($day, $month, $year);
13 print "It's a leap year!\n" if (leap($year));
14 $lmst = mjd2lst($mjd, $longitude, $dUT1);
15 $turns = str2turn($string, 'H');
16 $str = turn2str($turn, 'D', $sig);
17
18=head1 DESCRIPTION
19
20Astro::Time contains an assorted set Perl routines for time based
21conversions, such as conversion between calendar dates and Modified
22Julian day and conversion of UT to local sidereal time. Include are
23routines for conversion between numerical and string representation of
24angles.
25
26=head1 AUTHOR
27
28Chris Phillips Chris.Phillips@csiro.au
29
30=head1 FUNCTIONS
31
32=cut
33
34
35BEGIN {
3632.6e-58.7e-6 use Exporter ();
3718.0e-68.0e-6 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK @EXPORT_FAIL
# spent 102µs making 1 call to vars::import
3828.0e-54.0e-5 $PI $StrSep $StrZero $Quiet );
3911.0e-61.0e-6 $VERSION = '1.15';
4017.0e-67.0e-6 @ISA = qw(Exporter);
41
4211.3e-51.3e-5 @EXPORT = qw( cal2dayno dayno2cal leap yesterday tomorrow
43 mjd2cal cal2mjd mjd2dayno dayno2mjd now2mjd
44 jd2mjd mjd2jd mjd2time
45 gst mjd2lst cal2lst dayno2lst rise lst2mjd
46 turn2str deg2str rad2str str2turn str2deg str2rad
47 hms2time time2hms
48 deg2rad rad2deg turn2rad rad2turn turn2deg deg2turn
49 $PI );
5012.0e-62.0e-6 @EXPORT_OK = qw ( daynoOK monthOK dayOK utOK nint $StrSep $StrZero );
5112.0e-62.0e-6 @EXPORT_FAIL = qw ( @days );
52
5333.2e-51.1e-5 use Carp;
# spent 51µs making 1 call to Exporter::import
5432.9e-59.7e-6 use POSIX qw( fmod floor ceil acos );
# spent 82µs making 1 call to POSIX::import
5510.002640.00264}
56
5711.0e-61.0e-6$PI = 3.1415926535897932384626433832795028841971693993751;
58
5911.0e-61.0e-6$StrZero = 0;
6011.0e-61.0e-6$StrSep = ':';
61
62100my $debug = 0; # Used for debugging str2turn
63
6411.0e-61.0e-6$Quiet = 0;
65
6612.0e-62.0e-6my @days = (31,28,31,30,31,30,31,31,30,31,30,31);
67
68# Is the dayno valid?
69sub daynoOK ($$) {
70 my ($dayno, $year) = @_;
71
72 if ($dayno<1 || $dayno>366 || ($dayno>365 && !leap($year))) {
73 carp '$dayno out of range' if (!$Quiet);
74 return 0;
75 } else {
76 return 1;
77 }
78}
79
80# Is the month valid?
81sub monthOK ($) {
82 my $month = shift;
83
84 if ($month > 12 || $month < 1) {
85 carp '$month out of range' if (!$Quiet);
86 return 0;
87 } else {
88 return 1
89 }
90}
91
92# IS the day of month OK? (assumes month IS ok - should be checked first)
93sub dayOK ($$$) {
94 my ($day, $month, $year) = @_;
95
96 $month--; # For array indexing
97 if (leap($year)) {
98 $days[1] = 29;
99 } else {
100 $days[1] = 28;
101 }
102
103 if ($day < 1 || $day > $days[$month]) {
104 carp '$day out of range' if (!$Quiet);
105 return 0;
106 } else {
107 return 1;
108 }
109}
110
111# Is the day fraction OK?
112sub utOK ($) {
113 my $ut = shift;
114
115 if ($ut < 0.0 || $ut >= 1.0) {
116 carp '$ut out of range' if (!$Quiet);
117 return 0;
118 } else {
119 return 1;
120 }
121}
122
123# Return the nearest integer (ie round)
124sub nint ($) {
125 my ($x) = @_;
126 ($x<0.0) ? return(ceil($x-0.5)) : return(floor($x+0.5))
127}
128
129=item B<turn2str>
130
131 $str = turn2str($turn, $mode, $sig);
132 $str = turn2str($turn, $mode, $sig, $strsep);
133
134 Convert fraction of a turn into string representation
135 $turn Angle in turns
136 $mode Mode of string to convert to:
137 'H' for hours
138 'D' for degrees
139 $sig number of significant figures
140 $strsep String separator (override for default $Astro::Time::StrSep)
141Note:
142 The behavior can be modified by the following two variables:
143 $Astro::Time::StrZero Minimum number of leading digits (zero padded
144 if needed)
145 $Astro::Time::StrSep (Overridden by optional fourth argument)
146 Deliminator used in string (Default ':')
147 This may also equal one of a number of special values:
148 'HMS' 12H45M12.3S or 170D34M56.2S
149 'hms' 12h45m12.3s or 170d34m56.2s
150 'deg' 170d34'56.2"
151
152=cut
153
154sub turn2str ($$$;$) {
155 my($turn, $mode, $sig, $strsep) = @_;
156 $mode = uc $mode;
157 if (($mode ne 'H') && ($mode ne 'D')) {
158 carp 'turn2str: $mode must equal \'H\' or \'D\'';
159 return undef;
160 }
161 $strsep = $StrSep if (!defined $strsep);
162
163 my ($angle, $str, $sign, $wholesec, $secfract, $min);
164
165 if ($mode eq 'H') {
166 $angle = $turn * 24;
167 } else {
168 $angle = $turn * 360;
169 }
170
171 if ($angle < 0.0) {
172 $sign = -1;
173 $angle = -$angle;
174 } else {
175 $sign = 1;
176 }
177
178 my $wholeangle = (int $angle);
179
180 $angle -= $wholeangle;
181 $angle *= 3600;
182
183 # Get second fraction
184 $wholesec = int $angle;
185 $secfract = $angle - $wholesec;
186
187 $wholesec %= 60;
188 $min = ($angle-$wholesec - $secfract)/60.0;
189 $secfract = int ($secfract * 10**$sig + 0.5); # Add 0.5 to ensure rounding
190
191 # Check we have not rounded too far
192 if ($secfract >= 10**$sig) {
193 $secfract -= 10**$sig;
194 $wholesec++;
195 if ($wholesec >= 60.0) {
196 $wholesec -= 60;
197 $min++;
198 if ($min >= 60.0) {
199 $min -= 60;
200 $wholeangle++;
201 }
202 }
203 }
204
205 my $angleform;
206 if ($StrZero > 0) {
207 $angleform = "%0$StrZero";
208 } else {
209 $angleform = '%';
210 }
211
212 my ($sep1, $sep2, $sep3);
213 if ($strsep eq 'HMS') {
214 if ($mode eq 'H') {
215 $sep1 = 'H';
216 } else {
217 $sep1 = 'D';
218 }
219 $sep2 = 'M';
220 $sep3 = 'S';
221 } elsif ($strsep eq 'hms') {
222 if ($mode eq 'H') {
223 $sep1 = 'h';
224 } else {
225 $sep1 = 'd';
226 }
227 $sep2 = 'm';
228 $sep3 = 's';
229 } elsif ($strsep eq 'deg') { # What if $mode eq 'H'??
230 $sep1 = 'd';
231 $sep2 = "'";
232 $sep3 = '"';
233 } else {
234 $sep1 = $sep2 = $strsep;
235 $sep3 = '';
236 }
237
238 if ($sig > 0) {
239 $str = sprintf("${angleform}d$sep1%02d".
240 "$sep2%02d.%0${sig}d$sep3",
241 $wholeangle, $min, $wholesec, $secfract);
242 } else {
243 $str = sprintf("${angleform}d$sep1%02d".
244 "$sep2%02d$sep3",
245 $wholeangle, $min, $wholesec);
246 }
247
248 if ($sign == -1) {
249 $str = '-'.$str;
250 }
251 return $str;
252}
253
254=item B<deg2str>
255
256 $str=deg2str($deg, $mode, $sig);
257
258 Convert degrees into string representation
259 $deg angle in degrees
260 $mode mode of string to convert to:
261 'H' for hours
262 'D' for degrees
263 $sig number of significant figures
264 See note for turn2str
265
266=cut
267
268sub deg2str ($$$;$) {
269
270 my($deg, $mode, $sig, $strsep) = @_;
271 return turn2str($deg/(360), $mode, $sig, $strsep);
272}
273
274=item B<rad2str>
275
276 $str=rad2str($rad, $mode, $sig);
277
278 Convert radians into string representation
279 $rad angle in radians
280 $mode mode of string to convert to:
281 'H' for hours
282 'D' for degrees
283 $sig is number of significant figures
284 See note for turn2str
285
286=cut
287
288sub rad2str ($$$;$) {
289
290 my($rad, $mode, $sig, $strsep) = @_;
291 return turn2str($rad/(2*$PI), $mode, $sig, $strsep);
292}
293
294=item B<str2turn>
295
296 $turns = str2turn($string,$mode);
297
298 Convert angle from string representation into fraction of a turn
299 $string a : or space delimited angle
300 $mode type of angle
301 'H' if $string is in hours,min,sec
302 'D' if $string is in deg,arcmin,arcsec
303 The format of $string can be fairly flexible e.g.:
304 12.2
305 12:34
306 12:34:45.1
307 -23 34 12.3
308 -34 34.3
309
310 Note: You cannot mix spaces and :
311
312=cut
313
314sub str2turn ($$) {
315 my($str,$mode) = @_;
316
317 if (! defined $str) {
318 carp 'Use of uninitialized value at';
319 return undef;
320 }
321 $mode = uc $mode;
322 if (($mode ne "H") && ($mode ne "D")) {
323 carp 'str2turn: $mode must equal "H" or "D"';
324 return undef;
325 }
326
327 my $sign = 1.0;
328 my $angle = 0.0;
329 my $min = 0.0;
330 my $sec = 0.0;
331
332 # Does it match dd:dd:dd.d form
333 if ($str =~ /^\s*(?:([+-])\s*)? # Sign (optional)
334 (\d*): # Hours degrees
335 (\d{0,2})(?:: # Minutes
336 (\d{0,2}(?:\.\d*)?))? # Seconds and fractions (both optional)
337 /x) {
338 print STDERR "Matches dd:dd:dd.d\n" if $debug;
339 $sign = -1 if (defined($1) && $1 eq "-");
340 $angle = $2 if ($2);
341 $min = $3 if ($3);
342 $sec = $4 if ($4);
343
344 # Does it match hms form 12h33m34.6s
345 } elsif ($str =~ /^\s*(?:([+-])\s*)? # Sign (optional)
346 (\d+)\s*(h)\s* # Hours
347 (?:(\d{0,2})\s*m\s* # Minutes optional
348 (?:(\d{0,2} # Seconds and fractions (optional)
349 (?:\.\d*)?)\s*s)?)?
350
351 /x) {
352 print STDERR "Matches hms\n" if $debug;
353 $sign = -1 if (defined($1) && $1 eq "-");
354 $angle = $2 if ($2);
355 $mode = 'H';
356 $min = $4 if ($4);
357 $sec = $5 if ($5);
358
359 # Does it match dms form 12d33m34.6s or 12d33'34.6"
360 } elsif ($str =~ /^\s*(?:([+-])\s*)? # Sign (optional)
361 (\d+)\s*([d])\s* # Degrees
362 (?:(\d{0,2})\s*[m']\s* # Minutes optional
363 (?:(\d{0,2} # Seconds and fractions (optional)
364 (?:\.\d*)?)\s*[s"])?)?
365 /x) {
366 print STDERR "Matches dms\n" if $debug;
367 $sign = -1 if (defined($1) && $1 eq "-");
368 $angle = $2 if ($2);
369 #$mode = uc($3);
370 $mode = 'D';
371 $min = $4 if ($4);
372 $sec = $5 if ($5);
373
374 # Does is match dd dd dd.d form
375 } elsif ($str =~ /^\s*(?:([+-])\s*)? # Sign (optional)
376 (\d+)\s+ # Hours degrees
377 (\d{0,2})(?:\s+ # Minutes
378 (\d{0,2}(?:\.\d*)?))? # Seconds and fractions
379 /x) {
380 print STDERR "Matches dd dd dd.d\n" if $debug;
381 $sign = -1 if (defined($1) && $1 eq "-");
382 $angle = $2 if ($2);
383 $min = $3 if ($3);
384 $sec = $4 if ($4);
385
386 # Does it match dd.d form
387 } elsif ($str =~ /^\s*(?:([+-])\s*)?(\d+(?:\.\d*)?)/) {
388 print STDERR "Matches dd.d\n" if $debug;
389 $sign = -1 if (defined($1) && $1 eq "-");
390 $angle = $2 if ($2);
391 } else {
392 return undef;
393 }
394
395 my $factor;
396 if ($mode eq "H") {
397 $factor = 24.0;
398 } else {
399 $factor = 360.0;
400 }
401
402 print "Got ($sign) $angle/$min/$sec [$mode]\n" if $debug;
403
404 return $sign * (($angle + ($min + $sec/60.0)/60.0)/ $factor);
405}
406
407=item B<str2deg>
408
409 $degrees=str2deg($string,$mode);
410
411 Convert angle from string representation into degrees
412 $string a : or space delimited angle
413 $mode 'H' if $string is in hours,min,sec
414 'D' if $string is in deg,arcmin,arcsec
415 See note for str2turn
416
417=cut
418
419sub str2deg ($$) {
420 my($str, $mode) = @_;
421 return str2turn($str, $mode) * 360;
422}
423
424=item B<str2rad>
425
426 $radians=str2rad($string,$mode);
427
428 Convert angle from string representation into radians
429 $string a : or space delimited angle
430 $mode 'H' if $string is in hours,min,sec
431 'D' if $string is in deg,arcmin,arcsec
432 See note for str2turn
433
434=cut
435
436sub str2rad ($$) {
437 my($str, $mode) = @_;
438 return str2turn($str, $mode) * 2*$PI;
439}
440
441=item B<hms2time>
442
443 ($time) = hms2time($hour, $minute, $second);
444 ($time) = hms2time($hour, $minute, $second, $mode);
445
446 Returns the day fraction given hours minutes and seconds (or degrees)
447 $time Day fraction
448 $hour Hours
449 $minutes Minutes
450 $second Seconds
451 $mode 'H' or 'D' to interpret as hours or degrees (default
452 hours)
453
454=cut
455
456sub hms2time ($$$;$) {
457
458 my($hour, $minute, $second, $mode) = @_;
459 $mode = 'H' if (!defined $mode);
460
461 my $factor;
462 if ($mode eq 'H' || $mode eq 'h') {
463 $factor = 24.0;
464 } elsif ($mode eq 'D' || $mode eq 'd') {
465 $factor = 360.0;
466 } else {
467 carp 'Illegal $mode value';
468 return undef;
469 }
470
471 return (($second/60 + $minute)/60 + $hour)/$factor;
472}
473
474=item B<time2hms>
475
476 ($sign, $hour, $minute, $second) = time2hms($time, $mode, $sig);
477
478 Returns hours (or degrees), minutes and seconds given the day fraction
479 $sign Sign of angle ('+' or '-')
480 $hour Hours
481 $minutes Minutes
482 $second Seconds
483 $time Day fraction
484 $mode Return degrees or Hours?
485 'H' for hours
486 'D' for degrees
487 $sig Number of significant digits for $second
488
489=cut
490
491sub time2hms ($$$) {
492
493 my($turn, $mode, $sig) = @_;
494 $mode = uc $mode;
495 if (($mode ne 'H') && ($mode ne 'D')) {
496 carp 'time2hms: $mode must equal \'H\' or \'D\'';
497 return undef;
498 }
499
500 my ($angle, $str, $sign, $wholesec, $secfract, $min);
501
502 if ($mode eq 'H') {
503 $angle = $turn * 24;
504 } else {
505 $angle = $turn * 360;
506 }
507
508 if ($angle < 0.0) {
509 $sign = '-';
510 $angle = -$angle;
511 } else {
512 $sign = '+';
513 }
514
515 my $wholeangle = (int $angle);
516
517 $angle -= $wholeangle;
518 $angle *= 3600;
519
520 # Get second fraction
521 $wholesec = int $angle;
522 $secfract = $angle - $wholesec;
523
524 $wholesec %= 60;
525 $min = ($angle-$wholesec - $secfract)/60.0;
526 $secfract = int ($secfract * 10**$sig + 0.5); # Add 0.5 to ensure rounding
527
528 # Check we have not rounded too far
529 if ($secfract >= 10**$sig) {
530 $secfract -= 10**$sig;
531 $wholesec++;
532 if ($wholesec >= 60.0) {
533 $wholesec -= 60;
534 $min++;
535 if ($min >= 60.0) {
536 $min -= 60;
537 $wholeangle++;
538 }
539 }
540 }
541 my $format = sprintf '%%0%dd', $sig;
542 $secfract = sprintf($format, $secfract);
543
544 if ($sig > 0) {
545 return($sign, $wholeangle, $min, "$wholesec.$secfract");
546 } else {
547 return($sign, $wholeangle, $min, $wholesec);
548 }
549}
550
551=item B<deg2rad>
552
553 $rad=deg2rad($deg);
554 Convert degrees to radians
555
556=cut
557
558sub deg2rad ($) {
559 return $_[0] / 180*$PI;
560}
561
562=item B<rad2deg>
563
564 $deg=rad2deg($rad);
565 Convert radians to degrees
566
567=cut
568
569sub rad2deg ($) {
570 return $_[0] * 180/$PI;
571}
572
573=item B<turn2rad>
574
575 $rad=turn2rad($turn);
576 Convert turns to radians
577
578=cut
579
580#sub turn2rad ($) {
581# return $_[0] * 2*$PI;
582#}
583
584
# spent 563ms within Astro::Time::turn2rad which was called 69345 times, avg 8µs/call: # 27738 times (241ms+0) by Astro::Coord::pol2r at line 129 of /wise/base/static/lib/perl5/site_perl/5.10.0/Astro/Coord.pm, avg 9µs/call # 27738 times (215ms+0) by Astro::Coord::pol2r at line 130 of /wise/base/static/lib/perl5/site_perl/5.10.0/Astro/Coord.pm, avg 8µs/call # 13869 times (107ms+0) by Astro::Coord::pol2r at line 131 of /wise/base/static/lib/perl5/site_perl/5.10.0/Astro/Coord.pm, avg 8µs/call
sub turn2rad {
585693450.031734.6e-7 my @ret;
586693450.101961.5e-6 foreach (@_) {
587693450.140792.0e-6 push @ret, $_ * 2*$PI;
588 }
589693450.119221.7e-6 if (@ret==1) {
590 return $ret[0];
591 } elsif (@ret==0) {
592 return undef;
593 } else {
594 return @ret;
595 }
596}
597
598=item B<rad2turn>
599
600 $turn=rad2turn($rad);
601 Convert radians to turns
602
603=cut
604
605#sub rad2turn ($) {
606# return $_[0] / (2*$PI);
607#}
608
609sub rad2turn {
610 my @ret;
611 foreach (@_) {
612 push @ret, $_/(2*$PI);
613 }
614 if (@ret==1) {
615 return $ret[0];
616 } elsif (@ret==0) {
617 return undef;
618 } else {
619 return @ret;
620 }
621}
622
623=item B<turn2deg>
624
625 $deg=turn2deg($turn);
626 Convert turns to radians
627
628=cut
629
630#sub turn2deg ($) {
631# return $_[0] * 360.0;
632#}
633
634sub turn2deg {
635 my @ret;
636 foreach (@_) {
637 push @ret, $_*360.0;
638 }
639 if (@ret==1) {
640 return $ret[0];
641 } elsif (@ret==0) {
642 return undef;
643 } else {
644 return @ret;
645 }
646}
647
648=item B<deg2turn>
649
650 $turn=deg2turn($deg);
651 Convert degrees to turns
652
653=cut
654
655#sub deg2turn ($) {
656# return $_[0] / 360.0;
657#}
658
659sub deg2turn {
660 my @ret;
661 foreach (@_) {
662 push @ret, $_/360.0;
663 }
664 if (@ret==1) {
665 return $ret[0];
666 } elsif (@ret==0) {
667 return undef;
668 } else {
669 return @ret;
670 }
671}
672
673=item B<cal2dayno>
674
675 $dayno = cal2dayno($day, $month, $year);
676
677 Returns the day number corresponding to $day of $month in $year.
678
679=cut
680
681# VERIFYED
682sub cal2dayno ($$$) {
683
684 my ($day, $month, $year) = @_;
685 return undef if (!monthOK($month));
686 return undef if (!dayOK($day, $month, $year));
687
688 $month--; # For array indexing
689
690 if (leap($year)) {
691 $days[1] = 29;
692 } else {
693 $days[1] = 28;
694 }
695
696 my $mon;
697 my $dayno = $day;
698 for ($mon=0; $mon<$month; $mon++) {
699 $dayno += $days[$mon];
700 }
701
702 return($dayno);
703}
704
705=item B<dayno2cal>
706
707 ($day, $month) = dayno2cal($dayno, $year);
708
709 Return the $day and $month corresponding to $dayno of $year.
710
711=cut
712
713# Verified
714sub dayno2cal ($$) {
715
716 my($dayno, $year) = @_;
717 return undef if (!daynoOK($dayno, $year));
718
719 if (leap($year)) {
720 $days[1] = 29;
721 } else {
722 $days[1] = 28;
723 }
724
725 my $month = 0;
726 my $end = $days[$month];
727 while ($dayno>$end) {
728 $month++;
729 $end+= $days[$month];
730 }
731 $end -= $days[$month];
732 my $day = $dayno - $end;
733 $month++;
734
735 return($day, $month);
736}
737
738=item B<leap>
739
740 $isleapyear = leap($year);
741
742 Returns true if $year is a leap year.
743 $year year in full
744
745=cut
746
747# NOT Verified
748sub leap ($) {
749 my $year = shift;
750 return (((!($year%4))&&($year%100))||(!($year%400)));
751}
752
753=item B<yesterday>
754
755 ($dayno, $year) = yesterday($dayno, $year);
756 ($day, $month, $year) = yesterday($day, $month, $year);
757
758 Winds back the day number by one, taking account of year wraps.
759 $dayno Day number of year
760 $year Year
761 $month Month
762 $day Day of month
763
764=cut
765
766# Verified
767sub yesterday($$;$) {
768
769 my ($day, $month, $dayno, $year);
770
771 if (scalar(@_)==2) {
772 ($dayno, $year) = @_;
773 return undef if (!daynoOK($dayno, $year));
774 } else {
775 ($day, $month, $year) = @_;
776 return undef if (!monthOK($month));
777 return undef if (!dayOK($day, $month, $year));
778 $dayno = cal2dayno($day, $month, $year);
779 }
780
781 $dayno--;
782 if ($dayno==0) {
783 $year--;
784 if (leap($year)) {
785 $dayno = 366;
786 } else {
787 $dayno = 365;
788 }
789 }
790
791 if (scalar(@_)==2) {
792 return ($dayno, $year);
793 } else {
794 ($day, $month) = dayno2cal($dayno, $year);
795 return($day, $month, $year);
796 }
797}
798
799=item B<tomorrow>
800
801 ($dayno, $year) = tomorrow($dayno, $year);
802 ($day, $month, $year) = tomorrow($day, $month, $year);
803
804 Advances the day number by one, taking account of year wraps.
805 $dayno Day number of year
806 $year Year
807 $month Month
808 $day Day of month
809
810=cut
811
812# Verified
813sub tomorrow($$;$) {
814
815 my ($day, $month, $dayno, $year);
816
817 if (scalar(@_)==2) {
818 ($dayno, $year) = @_;
819 return undef if (!daynoOK($dayno, $year));
820 } else {
821 ($day, $month, $year) = @_;
822 return undef if (!monthOK($month));
823 return undef if (!dayOK($day, $month, $year));
824 $dayno = cal2dayno($day, $month, $year);
825 }
826
827 $dayno++;
828 if (($dayno==366 && !leap($year)) || $dayno==367) {
829 $year++;
830 $dayno = 1;
831 }
832
833 if (scalar(@_)==2) {
834 return ($dayno, $year);
835 } else {
836 ($day, $month) = dayno2cal($dayno, $year);
837 return($day, $month, $year);
838 }
839}
840
841=item B<mjd2cal>
842
843 ($day, $month, $year, $ut) = mjd2cal($mjd);
844
845 Converts a modified Julian day number into calendar date (universal
846 time). (based on the slalib routine sla_djcl).
847 $mjd Modified Julian day (JD-2400000.5)
848 $day Day of the month.
849 $month Month of the year.
850 $year Year
851 $ut UT day fraction
852
853=cut
854
855# VERIFIED
856sub mjd2cal($) {
857
858 my $mjd = shift;
859
860 my $ut = fmod($mjd,1.0);
861 if ($ut<0.0) {
862 $ut += 1.0;
863 $mjd -= 1;
864 }
865
86630.001180.00039 use integer; # Calculations require integer division and modulation
# spent 10µs making 1 call to integer::import
867 # Get the integral Julian Day number
868 my $jd = nint($mjd + 2400001);
869
870 # Do some rather cryptic calculations
871
872 my $temp1 = 4*($jd+((6*(((4*$jd-17918)/146097)))/4+1)/2-37);
873 my $temp2 = 10*((($temp1-237)%1461)/4)+5;
874
875 my $year = $temp1/1461-4712;
876 my $month =(($temp2/306+2)%12)+1;
877 my $day = ($temp2%306)/10+1;
878
879 return($day, $month, $year, $ut);
880}
881
882=item B<cal2mjd>
883
884 $mjd = cal2mjd($day, $month, $year, $ut);
885
886 Converts a calendar date (universal time) into modified Julian day
887 number.
888 $day Day of the month.
889 $month Month of the year.
890 $year Year
891 $ut UT dayfraction
892 $mjd Modified Julian day (JD-2400000.5)
893
894=cut
895
896# Verified
897sub cal2mjd($$$;$) {
898 my($day, $month, $year, $ut) = @_;
899 $ut = 0.0 if (!defined $ut);
900
901 return undef if (!monthOK($month));
902 return undef if (!dayOK($day, $month, $year));
903 return undef if (!utOK($ut));
904
905 my ($m, $y);
906 if ($month <= 2) {
907 $m = int($month+9);
908 $y = int($year-1);
909 } else {
910 $m = int($month-3);
911 $y = int($year);
912 }
913 my $c = int($y/100);
914 $y = $y-$c*100;
915 my $x1 = int(146097.0*$c/4.0);
916 my $x2 = int(1461.0*$y/4.0);
917 my $x3 = int((153.0*$m+2.0)/5.0);
918 return($x1+$x2+$x3+$day-678882+$ut);
919}
920
921=item B<mjd2dayno>
922
923 ($dayno, $year, $ut) = mjd2dayno($mjd);
924
925 Converts a modified Julian day number into year and dayno (universal
926 time).
927 $mjd Modified Julian day (JD-2400000.5)
928 $year Year
929 $dayno Dayno of year
930
931=cut
932
933# NOT Verified
934sub mjd2dayno($) {
935
936 my $mjd = shift;
937 my ($day, $month, $year, $ut) = mjd2cal($mjd);
938
939 return (cal2dayno($day,$month,$year), $year, $ut);
940}
941
942=item B<dayno2mjd>
943
944 $mjd = dayno2mjd($dayno, $year, $ut);
945
946 Converts a dayno and year to modified Julian day
947 $mjd Modified Julian day (JD-2400000.5)
948 $year Year
949 $dayno Dayno of year
950
951=cut
952
953# Not verified - wrapper to cal2mjd
954sub dayno2mjd($$;$) {
955
956 my ($dayno, $year, $ut) = @_;
957 $ut = 0.0 if (!defined $ut);
958
959 return undef if (!daynoOK($dayno,$year));
960 return undef if (!utOK($ut));
961
962 my ($day, $month) = dayno2cal($dayno, $year);
963
964 return cal2mjd($day, $month, $year, $ut);
965}
966
967=item B<now2mjd>
968
969 $mjd = now2mjd()
970
971=cut
972
973# Not verified - just wrapper
974sub now2mjd() {
975 my ($s, $m, $h, $day, $month, $year) = gmtime();
976 $month++;
977 $year += 1900;
978 return(cal2mjd($day, $month, $year, ((($s/60+$m)/60+$h)/24)));
979}
980
981=item B<jd2mjd>
982
983 $mjd = jd2mjd($jd);
984
985 Converts a Julian day to Modified Julian day
986 $jd Julian day
987 $mjd Modified Julian day
988
989=cut
990
991sub jd2mjd($) {
992 return (shift)-2400000.5;
993}
994
995=item B<mjd2jd>
996
997 $jd = mjd2jd($mjd);
998
999 Converts a Modified Julian day to Julian day
1000 $mjd Modified Julian day
1001 $jd Julian day
1002
1003=cut
1004
1005sub mjd2jd($) {
1006 return (shift)+2400000.5;
1007}
1008
1009=item B<mjd2ime>
1010
1011 $str = mjd2time($mjd);
1012
1013 Converts a Modified Julian day to a formatted string
1014 $mjd Modified Julian day
1015 $str Formatted time
1016
1017=cut
1018
1019sub mjd2time($) {
1020 my ($dayno, $year, $ut) = mjd2dayno($_[0]);
1021 return sprintf("$year %03d/%s", $dayno, turn2str($ut, 'H', 0));
1022}
1023
1024=item B<gst>
1025
1026 $gst = gst($mjd);
1027 $gmst = gst($mjd, $dUT1);
1028 $gtst = gst($mjd, $dUT1, $eqenx);
1029
1030 Converts a modified Julian day number to Greenwich sidereal time
1031 $mjd modified Julian day (JD-2400000.5)
1032 $dUT1 difference between UTC and UT1 (UT1 = UTC + dUT1) (seconds)
1033 $eqenx Equation of the equinoxes (not yet supported)
1034 $gst Greenwich sidereal time (turns)
1035 $gmst Greenwich mean sidereal time (turns)
1036 $gtst Greenwich true sidereal time (turns)
1037
1038=cut
1039
1040# Verified
1041sub gst($;$$) {
1042 my ($mjd, $dUT1, $eqenx) = @_;
1043 $dUT1 = 0.0 if (! defined $dUT1);
1044 if ($dUT1 > 0.5 || $dUT1 < -0.5) {
1045 carp '$dUT1 out of range';
1046 return undef;
1047 }
1048 if (defined $eqenx) {
1049 croak '$eqenx is not supported yet';
1050 }
1051
1052 my $JULIAN_DAY_J2000 = 2451545.0;
1053 my $JULIAN_DAYS_IN_CENTURY = 36525.0;
1054 my $SOLAR_TO_SIDEREAL = 1.002737909350795;
1055
1056 my $a=101.0 + 24110.54841/86400.0;
1057 my $b=8640184.812866/86400.0;
1058 my $e=0.093104/86400.0;
1059 my $d=0.0000062/86400.0;
1060 my $tu = (int($mjd)-($JULIAN_DAY_J2000-2400000.5))/$JULIAN_DAYS_IN_CENTURY;
1061 my $sidtim = $a + $tu*($b + $tu*($e - $tu*$d));
1062 $sidtim -= int($sidtim);
1063 if ($sidtim < 0.0) {$sidtim += 1.0};
1064 my $gmst = $sidtim + ($mjd - int($mjd) + $dUT1/86400.0)*$SOLAR_TO_SIDEREAL;
1065 while ($gmst<0.0) {
1066 $gmst += 1.0;
1067 }
1068 while ($gmst>1.0) {
1069 $gmst -= 1.0;
1070 }
1071
1072 return $gmst;
1073}
1074
1075# Not verified - wrapper to gmst
1076
1077=item B<mjd2lst>
1078
1079 $lst = mjd2lst($mjd, $longitude);
1080 $lmst = mjd2lst($mjd, $longitude, $dUT1);
1081 $ltst = mjd2lst($mjd, $longitude, $dUT1, $eqenx);
1082
1083 Converts a modified Julian day number into local sidereal time (lst),
1084 local mean sidereal time (lmst) or local true sidereal time (ltst).
1085 Unless high precisions is required dUT1 can be omitted (it will always
1086 be in the range -0.5 to 0.5 seconds).
1087 $mjd Modified Julian day (JD-2400000.5)
1088 $longitude Longitude for which the LST is required (turns)
1089 $dUT1 Difference between UTC and UT1 (UT1 = UTC + dUT1)(seconds)
1090 $eqenx Equation of the equinoxes (not yet supported)
1091 $lst Local sidereal time (turns)
1092 $lmst Local mean sidereal time (turns)
1093 $ltst Local true sidereal time (turns)
1094
1095=cut
1096
1097sub mjd2lst($$;$$) {
1098 my ($mjd, $longitude, $dUT1, $eqenx) = @_;
1099
1100 my $lst = gst($mjd, $dUT1, $eqenx);
1101 return undef if (!defined $lst);
1102 $lst += $longitude;
1103
1104 while ($lst>1.0) {
1105 $lst-= 1;
1106 }
1107 while ($lst < 0.0) {
1108 $lst += 1;
1109 }
1110 return $lst;
1111}
1112
1113=item B<cal2lst>
1114
1115 $lst = cal2lst($day, $month, $year, $ut, $longitude);
1116 $lmst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1);
1117 $ltst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1, $eqenx);
1118
1119 Wrapper to mjd2lst using calendar date rather than mjd
1120
1121=cut
1122
1123sub cal2lst($$$$$;$$) {
1124 my ($day, $month, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
1125 my $mjd = cal2mjd($day, $month, $year, $ut);
1126 return undef if (!defined $mjd);
1127
1128 return mjd2lst($mjd, $longitude, $dUT1, $eqenx);
1129}
1130
1131=item B<dayno2lst>
1132
1133 $lst = dayno2lst($dayno, $year, $ut, $longitude);
1134 $lmst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1);
1135 $ltst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1, $eqenx);
1136
1137 Wrapper to mjd2lst using calendar date rather than mjd
1138
1139=cut
1140
1141sub dayno2lst($$$$;$$) {
1142 my ($dayno, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
1143 my $mjd = dayno2mjd($dayno, $year, $ut);
1144 return undef if (!defined $mjd);
1145
1146 return mjd2lst($mjd, $longitude, $dUT1, $eqenx);
1147}
1148
1149# Not verified
1150
1151=item B<rise>
1152
1153 ($lst_rise, $lst_set) = rise($ra, $dec, $obslat, $el_limit);
1154
1155 Return the lst rise and set time of the given source
1156 $lst_rise, $lst_set Rise and set time (turns)
1157 $ra, $dec RA and Dec of source (turns)
1158 $obslat Latitude of observatory (turns)
1159 $el_limit Elevation limit of observatory
1160 (turns, 0 horizontal)
1161 Returns 'Circumpolar' if source circumpolar
1162 Returns undef if source never rises
1163
1164 Uses the formula:
1165 cos $z_limit = sin $obslat * sin $dec + cos $obslat * cos $dec
1166 * cos $HA
1167 where:
1168 $z_limit is the zenith angle limit corresponding to $el_limit
1169 $HA is the Hour Angle of the source
1170NOTE: For maximum accuracy source coordinated should be precessed to
1171 the current date.
1172
1173=cut
1174
1175sub rise ($$$$) {
1176 #print "rise: Got @_\n";
1177 my ($ra, $dec, $obslat, $el_limit) = @_;
1178 $ra = turn2rad($ra);
1179 $dec = turn2rad($dec);
1180 $obslat = turn2rad($obslat);
1181 $el_limit = turn2rad($el_limit);
1182
1183 my $z_limit = $PI/2-$el_limit;
1184
1185 #print "Check it\n";
1186
1187 # Check whether the source ever rises or is circumpolar
1188 my $z = acos(sin($obslat)*sin($dec) + cos($obslat)*cos($dec)); # Highest point
1189 return (undef) if ($z>$z_limit);
1190
1191 #print "Got here\n";
1192
1193 $z = acos(sin($obslat)*sin($dec) - cos($obslat)*cos($dec)); # Lowest point
1194
1195 return ('Circumpolar') if ($z<$z_limit);
1196
1197 my $cos_ha = (cos($z_limit) - sin($obslat)*sin($dec))
1198 /(cos($obslat)*cos($dec));
1199 my $ha = acos($cos_ha);
1200
1201 my $lst_rise = $ra - $ha;
1202 my $lst_set = $ra + $ha;
1203
1204 $lst_rise += 2*$PI if ($lst_rise < 0.0);
1205 $lst_set -= 2*$PI if ($lst_set >= 2*$PI);
1206
1207 return rad2turn($lst_rise), rad2turn($lst_set);
1208}
1209
1210=item B<lst2mjd>
1211
1212 $mjd = lst2mjd($lmst, $dayno, $year, $longitude);
1213 $mjd = lst2mjd($lmst, $dayno, $year, $longitude, $dUT1);
1214
1215 This routine calculates the modified Julian day number corresponding
1216 to the local mean sidereal time $lmst at $longitude, on a given UT
1217 day number ($dayno). Unless high precision is required dUT1 can be
1218 omitted.
1219
1220 The required inputs are :
1221 $lmst - The local mean sidereal time (turns)
1222 $dayno - The UT day of year for which to do the conversion
1223 $year - The year for which to do the conversion
1224 $longitude - The longitude of the observatory (turns)
1225 $dUT1 - Difference between UTC and UT1 (UT1 = UTC + dUT1)
1226 (seconds)
1227 $mjd The modified Julian day corresponding to $lmst on $dayno
1228
1229=cut
1230
1231sub lst2mjd($$$$;$) {
1232 my ($lmst, $dayno, $year, $longitude, $dUT1) = @_;
1233 $dUT1 = 0.0 if (!defined $dUT1);
1234
1235 my $SOLAR_TO_SIDEREAL = 1.002737909350795;
1236
1237 my $mjd = dayno2mjd($dayno, $year, $dUT1);
1238
1239 # Time in turns from passed lmst to lmst at the start of $dayno
1240 my $delay = $lmst-mjd2lst($mjd, $longitude);
1241
1242 if ($delay < 0.0) {
1243 $delay += 1.0;
1244 }
1245
1246 return($mjd + $delay/$SOLAR_TO_SIDEREAL);
1247}
1248
124917.0e-67.0e-61;
1250
1251__END__
1252

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