$Astro::Utils::VERSION
=
'0.02'
;
$Astro::Utils::AUTHORITY
=
'cpan:MANWAR'
;
use
vars
qw(@ISA @EXPORT)
;
@ISA
=
qw(Exporter)
;
@EXPORT
=
qw(calculate_equinox calculate_solstice)
;
use
5.006;
my
$PI
= 3.141592653589793;
my
$TYPE
= {
'MAR'
=> 0,
'JUN'
=> 1,
'SEP'
=> 2,
'DEC'
=> 3 };
sub
calculate_equinox {
my
(
$type
,
$timezone
,
$year
) =
@_
;
die
"ERROR: Year should be between -1000 and 3000.\n"
unless
(
defined
$year
&& (
$year
>= -1000) && (
$year
<= 3000));
if
(
$timezone
=~ /^utc$/i) {
return
_calc_utc_equinox(
$TYPE
->{
uc
(
$type
)},
$year
);
}
elsif
(
$timezone
=~ /^tdt$/i) {
return
_calc_tdt_equinox(
$TYPE
->{
uc
(
$type
)},
$year
);
}
else
{
die
"ERROR: Invalid timezone [$timezone] received.\n"
;
}
}
sub
calculate_solstice {
my
(
$type
,
$timezone
,
$year
) =
@_
;
die
"ERROR: Year should be between -1000 and 3000.\n"
unless
(
defined
$year
&& (
$year
>= -1000) && (
$year
<= 3000));
if
(
$timezone
=~ /^utc$/i) {
return
_calc_utc_solstice(
$TYPE
->{
uc
(
$type
)},
$year
);
}
elsif
(
$timezone
=~ /^tdt$/i) {
return
_calc_tdt_solstice(
$TYPE
->{
uc
(
$type
)},
$year
);
}
else
{
die
"ERROR: Invalid timezone [$timezone] received.\n"
;
}
}
sub
_calc_utc_equinox {
my
(
$k
,
$year
) =
@_
;
die
"ERROR: Invalid type for equinox.\n"
unless
(
defined
$k
&& (
$k
=~ /^[0|2]$/));
my
$jd
= _calc_jd(
$k
,
$year
);
return
_jd_to_utc(
$jd
);
}
sub
_calc_utc_solstice {
my
(
$k
,
$year
) =
@_
;
die
"ERROR: Invalid type for solstice.\n"
unless
(
defined
$k
&& (
$k
=~ /^[1|3]$/));
my
$jd
= _calc_jd(
$k
,
$year
);
return
_jd_to_utc(
$jd
);
}
sub
_calc_tdt_equinox {
my
(
$k
,
$year
) =
@_
;
die
"ERROR: Invalid type for equinox.\n"
unless
(
defined
$k
&& (
$k
=~ /^[0|2]$/));
my
$jd
= _calc_jd(
$k
,
$year
);
return
_jd_to_tdt(
$jd
);
}
sub
_calc_tdt_solstice {
my
(
$k
,
$year
) =
@_
;
die
"ERROR: Invalid type for solstice.\n"
unless
(
defined
$k
&& (
$k
=~ /^[1|3]$/));
my
$jd
= _calc_jd(
$k
,
$year
);
return
_jd_to_tdt(
$jd
);
}
sub
_calc_jd {
my
(
$k
,
$year
) =
@_
;
my
$jde0
= _jde0(
$k
,
$year
);
my
$t
= (
$jde0
- 2451545.0) / 36525;
my
$w
= 35999.373 *
$t
- 2.47;
my
$dl
= 1 + 0.0334 * _cos(
$w
) + 0.0007 * _cos(2 *
$w
);
my
$s
= _periodic_term_24(
$t
);
return
(
$jde0
+ ( (0.00001 *
$s
) /
$dl
));
}
sub
_jd_to_utc {
my
(
$jd
) =
@_
;
my
(
$yr
,
$mon
,
$day
,
$hr
,
$min
,
$sec
) = _process(
$jd
);
my
$date
= DateTime->new(
year
=>
$yr
,
month
=>
$mon
,
day
=>
$day
,
hour
=>
$hr
,
minute
=>
$min
,
second
=>
$sec
);
my
$first
= 1620;
my
$last
= 2002;
my
@tbl
= (121,112,103,95,88,82,77,72,68,63,60,56,53,51,48,46,44,42,40,38,
35,33,31,29,26,24,22,20,18,16,14,12,11,10,9,8,7,7,7,7,
7,7,8,8,9,9,9,9,9,10,10,10,10,10,10,10,10,11,11,11,
11,11,12,12,12,12,13,13,13,14,14,14,14,15,15,15,15,15,16,16,
16,16,16,16,16,16,15,15,14,13,
13.1,12.5,12.2,12.0,12.0,12.0,12.0,12.0,12.0,11.9,11.6,11.0,10.2,9.2,8.2,
7.1,6.2,5.6,5.4,5.3,5.4,5.6,5.9,6.2,6.5,6.8,7.1,7.3,7.5,7.6,
7.7,7.3,6.2,5.2,2.7,1.4,-1.2,-2.8,-3.8,-4.8,-5.5,-5.3,-5.6,-5.7,-5.9,
-6.0,-6.3,-6.5,-6.2,-4.7,-2.8,-0.1,2.6,5.3,7.7,10.4,13.3,16.0,18.2,20.2,
21.1,22.4,23.5,23.8,24.3,24.0,23.9,23.9,23.7,24.0,24.3,25.3,26.2,27.3,28.2,
29.1,30.0,30.7,31.4,32.2,33.1,34.0,35.0,36.5,38.3,40.2,42.2,44.5,46.5,48.5,
50.5,52.5,53.8,54.9,55.8,56.9,58.3,60.0,61.6,63.0,63.8,64.3
);
my
$delta_t
= 0;
my
$t
= (
$yr
- 2000) / 100;
if
(
$yr
>=
$first
&&
$yr
<=
$last
) {
if
(
$yr
% 2) {
$delta_t
= (
$tbl
[(
$yr
-
$first
- 1) / 2] +
$tbl
[(
$yr
-
$first
+ 1) / 2] ) / 2;
}
else
{
$delta_t
=
$tbl
[(
$yr
-
$first
) / 2];
}
}
elsif
(
$yr
< 948) {
$delta_t
= 2177 + 497 *
$t
+ (44.1 * _pow(
$t
, 2));
}
elsif
(
$yr
>= 948) {
$delta_t
= 102 + 102 *
$t
+ (25.3 * _pow(
$t
, 2));
if
(
$yr
>= 2000 &&
$yr
<= 2100) {
$delta_t
+= 0.37 * (
$yr
- 2100);
}
}
$date
->subtract(
seconds
=>
$delta_t
);
return
sprintf
(
"%04d-%02d-%02d %02d:%02d:%02d"
,
$date
->year,
$date
->month,
$date
->day,
$date
->hour,
$date
->minute,
$date
->second);
}
sub
_jd_to_tdt {
my
(
$jd
) =
@_
;
my
(
$yr
,
$mon
,
$day
,
$hr
,
$min
,
$sec
) = _process(
$jd
);
return
sprintf
(
"%04d-%02d-%02d %02d:%02d:%02d"
,
$yr
,
$mon
,
$day
,
$hr
,
$min
,
$sec
);
}
sub
_process {
my
(
$jd
) =
@_
;
my
(
$a
,
$alpha
);
my
$z
=
int
(
$jd
+ 0.5);
my
$f
= (
$jd
+ 0.5) -
$z
;
if
(
$z
< 2299161) {
$a
=
$z
;
}
else
{
$alpha
=
int
((
$z
- 1867216.25) / 36524.25);
$a
=
$z
+ 1 +
$alpha
-
int
(
$alpha
/ 4);
}
my
$b
=
$a
+ 1524;
my
$c
=
int
((
$b
- 122.1) / 365.25);
my
$d
=
int
(365.25 *
$c
);
my
$e
=
int
((
$b
-
$d
) / 30.6001);
my
$dt
=
$b
-
$d
-
int
(30.6001 *
$e
) +
$f
;
my
$mon
=
$e
- ((
$e
< 13.5)?(1):(13));
my
$yr
=
$c
- ((
$mon
> 2.5)?(4716):(4715));
my
$day
=
int
(
$dt
);
my
$h
= 24 * (
$dt
-
$day
);
my
$hr
=
int
(
$h
);
my
$m
= 60 * (
$h
-
$hr
);
my
$min
=
int
(
$m
);
my
$sec
=
int
(60 * (
$m
-
$min
));
return
(
$yr
,
$mon
,
$day
,
$hr
,
$min
,
$sec
);
}
sub
_periodic_term_24 {
my
(
$t
) =
@_
;
my
@a
= (485,203,199,182,156,136,77,74,70,58,52,
50,45,44,29,18,17,16,14,12,12,12,9,8);
my
@b
= (324.96,337.23,342.08,27.85,73.14,171.52,
222.54,296.72,243.58,119.81,297.17,21.02,
247.54,325.15,60.93,155.12,288.79,198.04,
199.76,95.39,287.11,320.81,227.73,15.45);
my
@c
= (1934.136,32964.467,20.186,445267.112,45036.886,
22518.443,65928.934,3034.906,9037.513,
33718.147,150.678,2281.226,29929.562,31555.956,
4443.417,67555.328,4562.452,62894.029,
31436.921,14577.848,31931.756,34777.259,
1222.114,16859.074);
my
$s
= 0;
foreach
my
$i
(0..23) {
$s
+= (
$a
[
$i
] * _cos(
$b
[
$i
] + (
$c
[
$i
] *
$t
)));
}
return
$s
;
}
sub
_jde0 {
my
(
$k
,
$year
) =
@_
;
my
(
$jde0
,
$y
);
if
(
$year
>= -1000 &&
$year
<= 1000) {
$y
=
$year
/ 1000;
if
(
$k
== 0) {
$jde0
= 1721139.29189 +
(365242.13740 *
$y
) +
(0.06134 * _pow(
$y
, 2)) +
(0.00111 * _pow(
$y
, 3)) -
(0.00071 * _pow(
$y
, 4));
}
elsif
(
$k
== 1) {
$jde0
= 1721233.25401 +
(365241.72562 *
$y
) -
(0.05323 * _pow(
$y
, 2)) +
(0.00907 * _pow(
$y
, 3)) +
(0.00025 * _pow(
$y
, 4));
}
elsif
(
$k
== 2) {
$jde0
= 1721325.70455 +
(365242.49558 *
$y
) -
(0.11677 * _pow(
$y
, 2)) -
(0.00297 * _pow(
$y
, 3)) +
(0.00074 * _pow(
$y
, 4));
}
elsif
(
$k
== 3) {
$jde0
= 1721414.39987 +
(365242.88257 *
$y
) -
(0.00769 * _pow(
$y
, 2)) -
(0.00933 * _pow(
$y
, 3)) -
(0.00006 * _pow(
$y
, 4));
}
}
elsif
(
$year
> 1000 &&
$year
<= 3000) {
$y
= (
$year
- 2000) / 1000;
if
(
$k
== 0) {
$jde0
= 2451623.80984 +
(365242.37404 *
$y
) +
(0.05169 * _pow(
$y
, 2)) -
(0.00411 * _pow(
$y
, 3)) -
(0.00057 * _pow(
$y
, 4));
}
elsif
(
$k
== 1) {
$jde0
= 2451716.56767 +
(365241.62603 *
$y
) +
(0.00325 * _pow(
$y
, 2)) +
(0.00888 * _pow(
$y
, 3)) -
(0.00030 * _pow(
$y
, 4));
}
elsif
(
$k
== 2) {
$jde0
= 2451810.21715 +
(365242.01767 *
$y
) -
(0.11575 * _pow(
$y
, 2)) +
(0.00337 * _pow(
$y
, 3)) +
(0.00078 * _pow(
$y
, 4));
}
elsif
(
$k
== 3) {
$jde0
= 2451900.05952 +
(365242.74049 *
$y
) -
(0.06223 * _pow(
$y
, 2)) -
(0.00823 * _pow(
$y
, 3)) +
(0.00032 * _pow(
$y
, 4));
}
}
return
$jde0
;
}
sub
_pow {
my
(
$n
,
$m
) =
@_
;
my
$r
= 1;
foreach
(1..
$m
) {
$r
*=
$n
;
}
return
$r
;
}
sub
_cos {
my
(
$degree
) =
@_
;
return
cos
((
$degree
*
$PI
)/180);
}
1;