Astro

 view release on metacpan or  search on metacpan

Astro/Coord.pm  view on Meta::CPAN

 coordinate system
 Note: convert equitoral positions to/from 3-vectors using pol2r and r2pol.
   @fk4     fk4 position to convert (as a 3-vector, turns)
   @gal     Galactic position (as a 3-vector, turns)
 Returns undef if too few or two many arguments are passed.
 Reference : Blaauw et al., 1960, MNRAS, 121, 123.

=cut

# Within 1e-7 arcsec of SLALIB slaEg50
sub fk4galr(@) {
  # First check that we have 3 arguments
  if (scalar @_ < 3) {
    croak 'Not enough arguments for Astro::Coord::fk4galr at ';
  } elsif (scalar @_ > 3) {
    croak 'Too many arguments for Astro::Coord::fk4galr at ';
  }

  my ($i, $j, @temp, @gal);
  my $w = 0.0;

Astro/Coord.pm  view on Meta::CPAN

 Notes: Converts equitoral positions to/from 3-vectors using pol2r and r2pol.
   $BRA,$BDec  fk4/B1950 position (turns)
   $l, $b      Galactic longitude and latitude
   @gal        Galactic position (as a 3-vector, turns)
   @fk4        fk4 position (as a  3-vector, turns)
 Reference : Blaauw et al., 1960, MNRAS, 121, 123.

=cut

# Within 1e-7 arcsec of SLALIB slaGe50
sub galfk4(@) {
  my (@r, $rect);

  if (@_==3) { # Rectangular coordinates passed
    @r = @_;
    $rect = 1;
  } elsif (@_==2) { # Sperical coordinates
    @r = pol2r($_[0],$_[1]); #  Spherical to Cartesian
    $rect = 0;
  } elsif (@_>3) {
    croak "Too many arguments for Astro::galfk4 at";

Astro/Coord.pm  view on Meta::CPAN

    $fk4[$i] = ($fk4[$i] + $eterm[$i])/$w;
  }

  if ($rect) {
    return @fk4;
  } else {
    return r2pol(@fk4);
  }
}

sub galfk4r(@) {galfk4(@_)};

#=item B<fk4fk5>
#
# ($JRA, $JDec) = fk4fk5($BRA, $BDec);
#
# Converts an FK4 (B1950) position to the equivalent FK5 (J2000) position.
#   **LOW PRECISION**
#   $BRA,$BDec     fk4/B1950 position (turns)
#   $JRA,$Dec      fk5/J2000 position (turns)
#

Astro/Coord.pm  view on Meta::CPAN

 Converts an J2000 position date coordinate

   $DRA,$DDec     Date coordinate (turns)
   $JRA,$Dec      J2000 position (turns)
   @date          Date coordinate (as a 3-vector)
   @J2000         J2000 position (as a 3-vector)

=cut

# Untested
sub J2000todate(@) {

  my ($rect);
  my (@J2000, @date); #  Position  vectors

  my $mjd = pop @_;
  if (@_==3) { # Rectangular coordinates passed
    @J2000 = @_;
    $rect = 1;
  } elsif (@_==2) { # Sperical coordinates
    @J2000 = pol2r($_[0],$_[1]); #  Spherical to Cartesian

Astro/Coord.pm  view on Meta::CPAN


  The returned value is :
    $haset       - The hour angle (turns) at which a source at this 
                   declination sets for an EWXY mounted antenna with the 
                   given limits at the given latitude

  NOTE: returns undef if %limits hash is missing any of the required keys

=cut

sub haset_ewxy($$\%) {

  my ($declination, $latitude, $limitsref) = @_;

  # Check that all the required keys are present
  if ((!exists $limitsref->{XLOW}) || (!exists $limitsref->{XLOW_KEYHOLE}) ||
      (!exists $limitsref->{XHIGH}) || (!exists $limitsref->{XHIGH_KEYHOLE}) ||
      (!exists $limitsref->{YLOW}) || (!exists $limitsref->{YLOW_KEYHOLE}) ||
      (!exists $limitsref->{YHIGH}) || (!exists $limitsref->{YHIGH_KEYHOLE})) {
    carp 'Missing key in %limits';
   return undef;

Astro/Coord.pm  view on Meta::CPAN

    \%limits     - A reference to a hash holding the EWXY antenna limits
                   The following keys must be defined XLOW, XLOW_KEYHOLE,
		   XHIGH, XHIGH_KEYHOLE, YLOW, YLOW_KEYHOLE, YHIGH,
		   YHIGH_KEYHOLE (all values should be in turns)

  The returned value is :
    $tlos        - The time left on-source (turns)

=cut

sub ewxy_tlos($$$\%) {

  my ($hour_angle, $declination, $latitude, $limitsref) = @_;

  my $haset = haset_ewxy($declination, $latitude, %$limitsref);
  return(undef) if (!defined $haset);
  $haset -= $hour_angle if (($haset > 0.0) && ($haset < 1.0));
  $haset += 1.0 if ($haset < 0.0);

  return $haset;
}

Astro/Coord.pm  view on Meta::CPAN


  The returned value is :
    $haset       - The hour angle (turns) at which a source at this
                   declination sets for an Az/El mounted antenna with the
                   given limits at the given latitude

  NOTE: returns undef if the %limits hash is missing any of the required keys

=cut

sub haset_azel($$\%) {

  my ($declination,  $latitude, $limitsref) = @_;

  # Check that all the required keys are present
  if (!exists $limitsref->{ELLOW}) {
    carp 'Missing key in %limits';
    return undef ;
  }

  my $cos_haset = (cos($PI / 2.0 - $limitsref->{ELLOW} * 2.0 *

Astro/Coord.pm  view on Meta::CPAN

    $latitude    - The latitude of the observatory (turns)
    %limits     - A reference to a hash holding the Az/El antenna limits
                   The following keys must be defined ELLOW (all values
                   should be in turns)

  The returned value is :
    $tlos        - The time left on-source (turns)

=cut

sub azel_tlos($$$\%) {
  my ($hour_angle, $declination, $latitude, $limitsref) = @_;

  # Calculate the time left onsource
  my $haset = haset_azel($declination, $latitude, %$limitsref);
  if (!defined $haset) {return(undef)};
  if (($haset > 0.0) && ($haset < 1.0)) { $haset -= $hour_angle; }
  if ($haset < 0.0) { $haset += 1.0; }

  return($haset);
}

Astro/Coord.pm  view on Meta::CPAN

	           antenna there must be a key for ELLOW (all values should
                   be in turns).

  The returned values are :
    $ha_set  - The hour angle at which the source sets (turns).  The hour
               angle at which the source rises is simply the negative of this
               value.

=cut

sub antenna_rise($$$$) {

  my ($declination, $latitude, $mount, $limitsref) = @_;

  # Check that the mount type is either EWXY (0) or AZEL (1)
  if (($mount != 0) && ($mount != 1)) {
    carp 'mount must equal 0 or 1';
    return undef;
  }

  if ($mount == 0) {

Astro/Coord.pm  view on Meta::CPAN

}

my @b2g = ([-0.054875539726,  0.494109453312, -0.867666135858],
	   [-0.873437108010, -0.444829589425, -0.198076386122],
	   [-0.483834985808,  0.746982251810,  0.455983795705]);

#my @b2g = ([ -0.0548777621, +0.4941083214, -0.8676666398],
#	   [ -0.8734369591, -0.4448308610, -0.1980741871],
#	   [ -0.4838350026, +0.7469822433, +0.4559837919]);

sub j2gal($$) {
  my ($ra,$dec) = @_;
  my @r = pol2r($ra,$dec);
  my @g = (0,0,0);
  for (my $i=0; $i<3; $i++) {
    for (my $j=0; $j<3; $j++) {
      $g[$i]+= $b2g[$j][$i] * $r[$j];
    }
  }
  return r2pol(@g);
}

Astro/Misc.pm  view on Meta::CPAN

=item B<lum2spectral>

  $spectral_type = lum2spectral($luminosity);

 Calculate the spectral type of a ZAMS star from its luminosity
 Based on Thompson 1984 ApJ 283 165 Table 1
   $luminosity   Star luminosity (normalised to Lsun)

=cut

sub lum2spectral($) {
  my $lum = log10(shift);

  my $n = scalar (@ThompsonData);
  if ($lum < $ThompsonData[0][LUM]) {
    return "<$ThompsonData[0][SPEC]";
  } elsif ($lum > $ThompsonData[$n-1][LUM]) {
    return ">$ThompsonData[$n-1][SPEC]";
  };

  $n = 1;

Astro/Time.pm  view on Meta::CPAN


 Winds back the day number by one, taking account of year wraps.
   $dayno       Day number of year
   $year        Year
   $month       Month
   $day         Day of month

=cut

# Verified
sub yesterday($$;$) {

  my ($day, $month, $dayno, $year);
  
  if (scalar(@_)==2) {
    ($dayno, $year) = @_;
    return undef if (!daynoOK($dayno, $year));
  } else {
    ($day, $month, $year) = @_;
    return undef if (!monthOK($month));
    return undef if (!dayOK($day, $month, $year));

Astro/Time.pm  view on Meta::CPAN


 Advances the day number by one, taking account of year wraps.
   $dayno       Day number of year
   $year        Year
   $month       Month
   $day         Day of month

=cut

# Verified
sub tomorrow($$;$) {

  my ($day, $month, $dayno, $year);

  if (scalar(@_)==2) {
    ($dayno, $year) = @_;
    return undef if (!daynoOK($dayno, $year));
  } else {
    ($day, $month, $year) = @_;
    return undef if (!monthOK($month));
    return undef if (!dayOK($day, $month, $year));

Astro/Time.pm  view on Meta::CPAN

 time). (based on the slalib routine sla_djcl).
    $mjd     Modified Julian day (JD-2400000.5)
    $day     Day of the month.
    $month   Month of the year.
    $year    Year
    $ut      UT day fraction

=cut

# VERIFIED
sub mjd2cal($) {

  my $mjd = shift;

  my $ut = fmod($mjd,1.0);
  if ($ut<0.0) {
    $ut += 1.0;
    $mjd -= 1;
  }

  use integer;  # Calculations require integer division and modulation

Astro/Time.pm  view on Meta::CPAN

 number.
    $day     Day of the month.
    $month   Month of the year.
    $year    Year
    $ut      UT dayfraction
    $mjd     Modified Julian day (JD-2400000.5)

=cut

# Verified
sub cal2mjd($$$;$) {
  my($day, $month, $year, $ut) = @_;
  $ut = 0.0 if (!defined $ut);

  return undef if (!monthOK($month));
  return undef if (!dayOK($day, $month, $year));
  return undef if (!utOK($ut));

  my ($m, $y);
  if ($month <= 2) {
    $m = int($month+9);

Astro/Time.pm  view on Meta::CPAN


 Converts a modified Julian day number into year and dayno (universal 
 time).
    $mjd     Modified Julian day (JD-2400000.5)
    $year    Year
    $dayno   Dayno of year

=cut

# NOT Verified
sub mjd2dayno($) {

  my $mjd = shift;
  my ($day, $month, $year, $ut) = mjd2cal($mjd);

  return  (cal2dayno($day,$month,$year), $year, $ut);
}

=item B<dayno2mjd>

  $mjd = dayno2mjd($dayno, $year, $ut);

 Converts a dayno and year to modified Julian day
    $mjd     Modified Julian day (JD-2400000.5)
    $year    Year
    $dayno   Dayno of year

=cut

# Not verified - wrapper to cal2mjd
sub dayno2mjd($$;$) {

  my ($dayno, $year, $ut) = @_;
  $ut = 0.0 if (!defined $ut);

  return undef if (!daynoOK($dayno,$year));
  return undef if (!utOK($ut));

  my ($day, $month) = dayno2cal($dayno, $year);

  return cal2mjd($day, $month, $year, $ut);
}

=item B<now2mjd>

  $mjd = now2mjd()

=cut

# Not verified - just wrapper
sub now2mjd() {
  my ($s, $m, $h, $day, $month, $year) = gmtime();
  $month++;
  $year += 1900;
  return(cal2mjd($day, $month, $year, ((($s/60+$m)/60+$h)/24)));
}

=item B<jd2mjd>

  $mjd = jd2mjd($jd);

 Converts a Julian day to Modified Julian day
    $jd      Julian day
    $mjd     Modified Julian day

=cut

sub jd2mjd($) {
  return (shift)-2400000.5;
}

=item B<mjd2jd>

  $jd = mjd2jd($mjd);

 Converts a Modified Julian day to Julian day
    $mjd     Modified Julian day
    $jd      Julian day

=cut

sub mjd2jd($) {
  return (shift)+2400000.5;
}

=item B<mjd2time>

  $str = mjd2time($mjd);
  $str = mjd2time($mjd, $np);

 Converts a Modified Julian day to a formatted string
    $mjd     Modified Julian day
    $str     Formatted time
    $np      Number of significant digits for fraction of a sec. Default 0

=cut

sub mjd2time($;$) {
  my ($dayno, $year, $ut) = mjd2dayno(shift);
  my $np = shift;
  $np = 0 if (! defined $np);
  return sprintf("$year %03d/%s", $dayno, turn2str($ut, 'H', $np));
}

=item B<mjd2vextime>

  $str = mjd2vextime($mjd);
  $str = mjd2vextime($mjd, $np);

 Converts a Modified Julian day to a vex formatted string
    $mjd     Modified Julian day
    $str     Formatted time
    $np      Number of significant digits for fraction of a sec. Default 0

=cut

sub mjd2vextime($;$) {
  my ($dayno, $year, $ut) = mjd2dayno(shift);
  my $np = shift;
  $np = 0 if (! defined $np);
  return sprintf("%dy%03dd%s", $year, $dayno, turn2str($ut, 'H', $np, 'hms'));
}

=item B<mjd2epoch>

  $time = mjd2epoch($mjd);

 Converts a Modified Julian day to unix Epoch (seconds sinve 1 Jan 1970)
 Rounded to the nearest second
    $mjd     Modified Julian day
    $tie     Seconds since 1 Jan 1970

=cut

sub mjd2epoch($) {
  my $mjd = shift;
  my $epoch = ($mjd - 40587)*24*60*60;
  return int($epoch + $epoch/abs($epoch*2)); # Work even if epoch is negative
}

=item B<gst>

  $gst  = gst($mjd);
  $gmst = gst($mjd, $dUT1);
  $gtst = gst($mjd, $dUT1, $eqenx);

Astro/Time.pm  view on Meta::CPAN

   $mjd     modified Julian day (JD-2400000.5)
   $dUT1    difference between UTC and UT1 (UT1 = UTC + dUT1) (seconds)
   $eqenx   Equation of the equinoxes (not yet supported)
   $gst     Greenwich sidereal time (turns)
   $gmst    Greenwich mean sidereal time (turns)
   $gtst    Greenwich true sidereal time (turns)

=cut

# Verified
sub gst($;$$) {
  my ($mjd, $dUT1, $eqenx) = @_;
  $dUT1 = 0.0 if (! defined $dUT1);
  if ($dUT1 > 0.5 || $dUT1 < -0.5) {
    carp '$dUT1 out of range';
    return undef;
  }
  if (defined $eqenx) {
    croak '$eqenx is not supported yet';
  }

Astro/Time.pm  view on Meta::CPAN

   $mjd         Modified Julian day (JD-2400000.5)
   $longitude   Longitude for which the LST is required (turns)
   $dUT1        Difference between UTC and UT1 (UT1 = UTC + dUT1)(seconds)
   $eqenx       Equation of the equinoxes (not yet supported)
   $lst         Local sidereal time (turns)
   $lmst        Local mean sidereal time (turns)
   $ltst        Local true sidereal time (turns)

=cut

sub mjd2lst($$;$$) {
  my ($mjd, $longitude, $dUT1, $eqenx) = @_;

  my $lst = gst($mjd, $dUT1, $eqenx);
  return undef if (!defined $lst);
  $lst += $longitude;

  while ($lst>1.0) {
    $lst-= 1;
  }
  while ($lst < 0.0) {

Astro/Time.pm  view on Meta::CPAN

=item B<cal2lst>

  $lst = cal2lst($day, $month, $year, $ut, $longitude);
  $lmst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1);
  $ltst = cal2lst($day, $month, $year, $ut, $longitude, $dUT1, $eqenx);

 Wrapper to mjd2lst using calendar date rather than mjd

=cut

sub cal2lst($$$$$;$$) {
  my ($day, $month, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
  my $mjd = cal2mjd($day, $month, $year, $ut);
  return undef if (!defined $mjd);

  return mjd2lst($mjd, $longitude, $dUT1, $eqenx);
}

=item B<dayno2lst>

  $lst = dayno2lst($dayno, $year, $ut, $longitude);
  $lmst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1);
  $ltst = dayno2lst($dayno, $year, $ut, $longitude, $dUT1, $eqenx);

 Wrapper to mjd2lst using calendar date rather than mjd

=cut

sub dayno2lst($$$$;$$) {
  my ($dayno, $year, $ut, $longitude, $dUT1, $eqenx) = @_;
  my $mjd = dayno2mjd($dayno, $year, $ut);
  return undef if (!defined $mjd);

  return mjd2lst($mjd, $longitude,  $dUT1, $eqenx);
}

# Not verified

=item B<rise>

Astro/Time.pm  view on Meta::CPAN

    $lmst      - The local mean sidereal time (turns)
    $dayno     - The UT day of year for which to do the conversion
    $year      - The year for which to do the conversion
    $longitude - The longitude of the observatory (turns)
    $dUT1      - Difference between UTC and UT1 (UT1 = UTC + dUT1) 
                                                                (seconds)
    $mjd         The modified Julian day corresponding to $lmst on $dayno

=cut

sub lst2mjd($$$$;$) {
  my ($lmst, $dayno, $year, $longitude, $dUT1) = @_;
  $dUT1 = 0.0 if (!defined $dUT1);

  my $SOLAR_TO_SIDEREAL = 1.002737909350795;

  my $mjd = dayno2mjd($dayno, $year, $dUT1);

  # Time in turns from passed lmst to lmst at the start of $dayno
  my $delay = $lmst-mjd2lst($mjd, $longitude);

Astro/Time.pm  view on Meta::CPAN

  This routine returns the name of the given month (as a number 1..12), 
  where 1 is January. The default is a 3 character version of the month
  ('Jan', 'Feb', etc) in the second form the full month is returned


  The required inputs are :
    $month      - The month in question with 1 == January.

=cut

sub month2str($;$) {
  my ($mon, $long) = @_;

  return undef if (!monthOK($mon));

  if ($long) {
    return $MonthStr[$mon-1];
  } else {
    return $MonthShortStr[$mon-1];
  }
}

Astro/Time.pm  view on Meta::CPAN

  $weekdaystr = mjd2weekdaystr($mjd);

 Returns the name of the weekday correspondig to the given MJD.
 May not work for historical dates.

    $mjd     Modified Julian day (JD-2400000.5)

=cut


sub mjd2weekdaystr($;$) {
  my ($mjd, $long) = @_;
  my $dow = mjd2weekday($mjd);
  if ($long) {
    return $WeekStr[$dow];
  } else {
    return $WeekShortStr[$dow];
  }
}

=item B<str2month>

Astro/Time.pm  view on Meta::CPAN

  Given the name of a month (in English), this routine returns the
  an integer between 1 and 12, where 1 is January. Full month names of
  3 character abbreviations are acceptable. Minumum matching (e.g. "Marc")
  is not supported.

  The required inputs are :
    $month      - Name of the month ('Jan', 'January', 'Feb', 'February' etc)

=cut

sub str2month($) {
  my $month = uc(shift);

  for (my $i=0; $i<12; $i++) {
    if ($month eq uc($MonthStr[$i]) || $month eq uc($MonthShortStr[$i])) {
      return $i+1;
    }
  }
  return undef;
}



( run in 0.425 second using v1.01-cache-2.11-cpan-65fba6d93b7 )