How to calculate sunset?

[emosenkis]

I'm trying to use this library to calculate sunset times for locations on Earth on a given date but I don't have a background in astronomy and I can't figure out how to calculate four of the values needed as input to the astro::transit::time function. Here's what I have so far:

extern crate astro;
extern crate chrono;

use astro::coords::GeographPoint;
use astro::time;
use astro::transit::{self, TransitBody, TransitType};
use chrono::{Datelike, Timelike, UTC};

fn main() {
    let coords: (f64, f64) = (40.7128, -74.0059);  // GPS coords of New York City
    let now = UTC::now();
    let day_of_month = time::DayOfMonth {
        day: now.day() as u8,
        hr: now.hour() as u8,
        min: now.minute() as u8,
        sec: now.second() as f64,
        time_zone: 0.0,
    };
    let date = time::Date {
        year: now.year() as i16,                                                                                                     
        month: now.month() as u8,                                                                                                    
        decimal_day: time::decimal_day(&day_of_month),                                                                               
        cal_type: time::CalType::Gregorian,                                                                                          
    };                                                                                                                               
    let julian_day = time::julian_day(&date);                                                                                        
    let delta_t = time::delta_t(date.year as i32, date.month);                                                                       
    let julian_ephm_day = time::julian_ephemeris_day(julian_day, delta_t);                                                           
    let sunset_time = astro::transit::time(&TransitType::Set,
                                           &TransitBody::Sun,
                                           &GeographPoint {
                                               lat: coords.0.to_radians(),
                                               long: coords.1.to_radians(),
                                           },
                                           // How do I calculate the next four args?
                                           eq_point1,
                                           eq_point2,
                                           eq_point3,
                                           apprnt_greenwhich_sidr,
                                           // End of missing args
                                           delta_t,
                                           0f64); // moon_eq_hz_parallax not needed
    println!("{:?}", sunset_time);                                                                                                   
}

Any pointers would be greatly appreciated. Thanks!

[saurvs]

Hi Eitan,

For obtaining the eq_points, the basic idea is to calculate the heliocentric coordinates of the Earth from astro::planet::heliocent_coords, and then apply a series of transformations and corrections to get the equatorial coordinates of the Sun. The whole process is fairly lengthy, but is pretty well documented here. Every step in that document is essentially a function in this library. By step 3.7, you'll have the Sun's apparent ecliptic longitude and latitude, which you can then convert to equatorial coordinates using the eq_frm_ecl! macro. I've been wanting to write a simple macro that calculates the geocentric coordinates of the Sun by just taking the JD; I just haven't gotten around to it yet.

Regarding apprnt_greenwhich_sidr, you can use the apprnt_sidr! macro, but you should probably manually calculate the quantities needed for astro::time::apprnt_sidr so you can reuse the nutation values for the calculation of the Sun's equatorial coordinates.

[emosenkis]

Thanks. I've tried implementing that, but I'm getting results that don't seem at all related to the actual sunset time. Can you spot my mistake?

#[macro_use]
extern crate astro;
extern crate chrono;

use astro::coords::{EqPoint, GeographPoint};
use astro::{aberr, nutation, planet, sun, time};
use astro::transit::{TransitBody, TransitType};
use chrono::{Datelike, Timelike, UTC};

fn main() {
    let coords: (f64, f64) = (40.7128, -74.0059);  // GPS coords of New York City
    let now = UTC::now();
    let day_of_month = time::DayOfMonth {
        day: now.day() as u8,
        hr: 0,
        min: 0,
        sec: 0.0,  // Should be TT 0, which is +/- 32.184 secs, I think
        time_zone: 0.0,
    };
    println!("{}", time::decimal_day(&day_of_month));
    let date = time::Date {
        year: now.year() as i16,                                                                                                     
        month: time::Month::from_int(now.month() as u8),                                                                             
        decimal_day: time::decimal_day(&day_of_month),                                                                               
        cal_type: time::CalType::Gregorian,                                                                                          
    };                                                                                                                               
    let julian_day = time::julian_day(&date);                                                                                        
    println!("Julian day: {}", julian_day);                                                                                          
    let delta_t = time::delta_t(date.year as i32, date.month as u8);                                                                 
    println!("Delta-T: {}", delta_t);                                                                                                
    let julian_ephm_day = time::julian_ephemeris_day(julian_day, delta_t);                                                           
    println!("Julian_ephm_day: {}", julian_ephm_day);                                                                                
    let (nut_in_long, nut_in_oblq) = nutation::nutation(julian_ephm_day);                                                            
    let mean_eclip_oblq = astro::ecliptic::mn_oblq_laskar(julian_ephm_day);                                                          
    let true_eclip_oblq = mean_eclip_oblq + nut_in_oblq;                                                                             
    let mean_sidr = time::mn_sidr(julian_ephm_day);                                                                                  
    let apprnt_greenwhich_sidr = time::apprnt_sidr(mean_sidr, nut_in_long, true_eclip_oblq);                                         
    let eq_point1 = sun_eq_coords(julian_ephm_day - 1f64);                                                                           
    let eq_point2 = sun_eq_coords(julian_ephm_day);                                                                                  
    let eq_point3 = sun_eq_coords(julian_ephm_day + 1f64);                                                                           
    let sunset_time = astro::transit::time(&TransitType::Set,                                                                        
                                           &TransitBody::Sun,                                                                        
                                           &GeographPoint {                                                                          
                                               lat: coords.0.to_radians(),                                                           
                                               long: coords.1.to_radians(),                                                          
                                           },                                                                                        
                                           &eq_point1,                                                                               
                                           &eq_point2,                                                                               
                                           &eq_point3,                                                                               
                                           apprnt_greenwhich_sidr,                                                                   
                                           delta_t,                                                                                  
                                           0f64); // moon_eq_hz_parallax not needed                                                  
    println!("{:?}", sunset_time);                                                                                                   
}                                                                                                                                    
                                                                                                                                     
fn sun_eq_coords(julian_ephm_day: f64) -> EqPoint {                                                                                  
    let (earth_long, earth_lat, _sun_earth_dist) = planet::heliocent_coords(&planet::Planet::Earth, julian_ephm_day);                
    let geocentric_long = (earth_long + 180f64) % 360f64;                                                                            
    let geocentric_lat = -earth_lat;                                                                                                 
    let (nut_in_long, nut_in_oblq) = nutation::nutation(julian_ephm_day);                                                            
    let (_sun_ecl_point, sun_earth_dist) = sun::geocent_ecl_pos(julian_ephm_day);
    let mean_eclip_oblq = astro::ecliptic::mn_oblq_laskar(julian_ephm_day);
    let true_eclip_oblq = mean_eclip_oblq + nut_in_oblq;
    let aberration_correction = aberr::sol_aberr(sun_earth_dist);
    let apprnt_sun_long = geocentric_long + nut_in_long + aberration_correction;
    let (asc, dec) = eq_frm_ecl!(apprnt_sun_long, geocentric_lat, true_eclip_oblq); 
    EqPoint { asc: asc, dec: dec }
}