view release on metacpan or  search on metacpan

lib/Astro/Sunrise.pm  view on Meta::CPAN

    my $t;
    if ( $cost >= 1.0 ) {
      if ($polar eq 'retval') {
        return ('night', 'night');
      }
      carp "Sun never rises!!\n";
      $t = 0.0;    # Sun always below altit
    }
    elsif ( $cost <= -1.0 ) {
      if ($polar eq 'retval') {
        return ('day', 'day');
      }
      carp "Sun never sets!!\n";
      $t = 12.0;    # Sun always above altit
    }
    else {
      my $arc = acosd($cost);    # The diurnal arc
      $t = $arc / $ang_spd;      # Time to traverse the diurnal arc, hours
      if ($trace) {
        printf $trace "Diurnal arc $arc -> $t hours (%s)\n", _fmt_dur($t);
      }
    }

    # Store rise and set times - in hours UT

    my $hour_rise_ut = $tsouth - $t;
    my $hour_set_ut  = $tsouth + $t;
    if ($trace) {
      printf $trace "For day $d (%s), sunrise at $hour_rise_ut (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon),
                   _fmt_hr($hour_rise_ut, $lon);
      printf $trace "For day $d (%s), sunset  at $hour_set_ut (%s)\n",  _fmt_hr(24 * ($d - int($d)), $lon),
                   _fmt_hr($hour_set_ut , $lon);
    }
    return($hour_rise_ut, $hour_set_ut);
}

#########################################################################################################
#
#
# FUNCTIONAL SEQUENCE for GMST0
#
# _GIVEN
# Day number
#
# _THEN
#
# computes GMST0, the Greenwich Mean Sidereal Time
# at 0h UT (i.e. the sidereal time at the Greenwich meridian at
# 0h UT).  GMST is then the sidereal time at Greenwich at any
# time of the day..
#
#
# _RETURN
#
# Sidtime
#
sub GMST0 {
    my ($d) = @_;

    my $sidtim0 = revolution(   ( 180.0 + 356.0470 + 282.9404 )
                              + ( 0.9856002585 + 4.70935E-5 ) * $d );
    return $sidtim0;

}

#
#
# FUNCTIONAL SEQUENCE for sun_RA_dec
#
# _GIVEN
# day number, $r and $lon (from sunpos)
#
# _THEN
#
# compute RA and dec
#
#
# _RETURN
#
# Sun's Right Ascension (RA), Declination (dec) and distance (r)
#
#
sub sun_RA_dec {
    my ($d, $lon_noon, $trace) = @_;

    # Compute Sun's ecliptical coordinates
    my ( $r, $lon ) = sunpos($d);
    if ($trace) {
      printf $trace "For day $d (%s), solar noon at ecliptic longitude $lon %s\n", _fmt_hr(24 * ($d - int($d)), $lon_noon), _fmt_angle($lon);
    }

    # Compute ecliptic rectangular coordinates (z=0)
    my $x = $r * cosd($lon);
    my $y = $r * sind($lon);

    # Compute obliquity of ecliptic (inclination of Earth's axis)
    my $obl_ecl = 23.4393 - 3.563E-7 * $d;

    # Convert to equatorial rectangular coordinates - x is unchanged
    my $z = $y * sind($obl_ecl);
    $y    = $y * cosd($obl_ecl);

    # Convert to spherical coordinates
    my $RA  = atan2d( $y, $x );
    my $dec = atan2d( $z, sqrt( $x * $x + $y * $y ) );

    return ( $RA, $dec, $r );

}    # sun_RA_dec


#
#
# FUNCTIONAL SEQUENCE for sunpos
#
# _GIVEN
#  day number
#
# _THEN
#
# Computes the Sun's ecliptic longitude and distance
# at an instant given in d, number of days since
# 2000 Jan 0.0.
#
#
# _RETURN
#
# ecliptic longitude and distance
# ie. $True_solar_longitude, $Solar_distance
#
sub sunpos {
    my ($d) = @_;

    #                       Mean anomaly of the Sun
    #                       Mean longitude of perihelion
    #                         Note: Sun's mean longitude = M + w
    #                       Eccentricity of Earth's orbit
    #                       Eccentric anomaly
    #                       x, y coordinates in orbit
    #                       True anomaly

    # Compute mean elements
    my $Mean_anomaly_of_sun = revolution( 356.0470 + 0.9856002585 * $d );
    my $Mean_longitude_of_perihelion = 282.9404 + 4.70935E-5 * $d;
    my $Eccentricity_of_Earth_orbit  = 0.016709 - 1.151E-9 * $d;

    # Compute true longitude and radius vector
    my $Eccentric_anomaly =   $Mean_anomaly_of_sun
                            + $Eccentricity_of_Earth_orbit * $RADEG
                               * sind($Mean_anomaly_of_sun)
                               * ( 1.0 + $Eccentricity_of_Earth_orbit * cosd($Mean_anomaly_of_sun) );

    my $x = cosd($Eccentric_anomaly) - $Eccentricity_of_Earth_orbit;

    my $y = sqrt( 1.0 - $Eccentricity_of_Earth_orbit * $Eccentricity_of_Earth_orbit )
            * sind($Eccentric_anomaly);

    my $Solar_distance = sqrt( $x * $x + $y * $y );    # Solar distance
    my $True_anomaly = atan2d( $y, $x );               # True anomaly

    my $True_solar_longitude = $True_anomaly + $Mean_longitude_of_perihelion;    # True solar longitude

    if ( $True_solar_longitude >= 360.0 ) {
      $True_solar_longitude -= 360.0;    # Make it 0..360 degrees
    }

    return ( $Solar_distance, $True_solar_longitude );
}



sub sind {
    sin( ( $_[0] ) * $DEGRAD );
}

sub cosd {
    cos( ( $_[0] ) * $DEGRAD );
}

sub tand {
    tan( ( $_[0] ) * $DEGRAD );
}

sub atand {
    ( $RADEG * atan( $_[0] ) );
}

sub asind {
    ( $RADEG * asin( $_[0] ) );
}

sub acosd {
    ( $RADEG * acos( $_[0] ) );
}

sub atan2d {
    ( $RADEG * atan2( $_[0], $_[1] ) );
}

#
#
# FUNCTIONAL SEQUENCE for revolution
#
# _GIVEN