You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

446 lines
18 KiB

//
// EDSunriseSet.m
//
// Created by Ernesto García on 20/08/11.
// Copyright 2011 Ernesto García. All rights reserved.
//
// C/C++ sun calculations created by Paul Schlyter
// sunriset.c
// http://stjarnhimlen.se/english.html
// SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
// the length of the day at any date and latitude
// Written as DAYLEN.C, 1989-08-16
// Modified to SUNRISET.C, 1992-12-01
// (c) Paul Schlyter, 1989, 1992
// Released to the public domain by Paul Schlyter, December 1992
//
#import "EDSunriseSet.h"
//
// Defines from sunriset.c
//
#define INV360 ( 1.0 / 360.0 )
#define RADEG ( 180.0 / M_PI )
#define DEGRAD ( M_PI / 180.0 )
/* The trigonometric functions in degrees */
#define sind(x) sin((x)*DEGRAD)
#define cosd(x) cos((x)*DEGRAD)
#define tand(x) tan((x)*DEGRAD)
#define atand(x) (RADEG*atan(x))
#define asind(x) (RADEG*asin(x))
#define acosd(x) (RADEG*acos(x))
#define atan2d(y,x) (RADEG*atan2(y,x))
/* A macro to compute the number of days elapsed since 2000 Jan 0.0 */
/* (which is equal to 1999 Dec 31, 0h UT) */
#define days_since_2000_Jan_0(y,m,d) \
(367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
#if defined(__IPHONE_8_0) || defined (__MAC_10_10)
#define EDGregorianCalendar NSCalendarIdentifierGregorian
#else
#define EDGregorianCalendar NSGregorianCalendar
#endif
#pragma mark - Readwrite accessors only private
@interface EDSunriseSet()
@property (nonatomic) double latitude;
@property (nonatomic) double longitude;
@property (nonatomic, strong) NSTimeZone *timezone;
@property (nonatomic, strong) NSCalendar *calendar;
@property (nonatomic, strong) NSTimeZone *utcTimeZone;
@property (readwrite, strong) NSDate *date;
@property (readwrite, strong) NSDate *sunset;
@property (readwrite, strong) NSDate *sunrise;
@property (readwrite, strong) NSDate *civilTwilightStart;
@property (readwrite, strong) NSDate *civilTwilightEnd;
@property (readwrite, strong) NSDate *nauticalTwilightStart;
@property (readwrite, strong) NSDate *nauticalTwilightEnd;
@property (readwrite, strong) NSDate *astronomicalTwilightStart;
@property (readwrite, strong) NSDate *astronomicalTwilightEnd;
@property (readwrite, strong) NSDateComponents* localSunrise;
@property (readwrite, strong) NSDateComponents* localSunset;
@property (readwrite, strong) NSDateComponents* localCivilTwilightStart;
@property (readwrite, strong) NSDateComponents* localCivilTwilightEnd;
@property (readwrite, strong) NSDateComponents* localNauticalTwilightStart;
@property (readwrite, strong) NSDateComponents* localNauticalTwilightEnd;
@property (readwrite, strong) NSDateComponents* localAstronomicalTwilightStart;
@property (readwrite, strong) NSDateComponents* localAstronomicalTwilightEnd;
@end
#pragma mark - Calculations from sunriset.c
@implementation EDSunriseSet(Calculations)
/*****************************************/
/* Reduce angle to within 0..360 degrees */
/*****************************************/
-(double) revolution:(double) x
{
return( x - 360.0 * floor( x * INV360 ) );
}
/*********************************************/
/* Reduce angle to within -180..+180 degrees */
/*********************************************/
-(double) rev180:(double) x
{
return( x - 360.0 * floor( x * INV360 + 0.5 ) );
}
-(double) GMST0:(double) d
{
double sidtim0;
/* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr */
/* L = M + w, as defined in sunpos(). Since I'm too lazy to */
/* add these numbers, I'll let the C compiler do it for me. */
/* Any decent C compiler will add the constants at compile */
/* time, imposing no runtime or code overhead. */
sidtim0 = [self revolution: ( 180.0 + 356.0470 + 282.9404 ) +
( 0.9856002585 + 4.70935E-5 ) * d];
return sidtim0;
}
/******************************************************/
/* Computes the Sun's ecliptic longitude and distance */
/* at an instant given in d, number of days since */
/* 2000 Jan 0.0. The Sun's ecliptic latitude is not */
/* computed, since it's always very near 0. */
/******************************************************/
-(void) sunposAtDay:(double)d longitude:(double*)lon r:(double *)r
{
double M, /* Mean anomaly of the Sun */
w, /* Mean longitude of perihelion */
/* Note: Sun's mean longitude = M + w */
e, /* Eccentricity of Earth's orbit */
E, /* Eccentric anomaly */
x, y, /* x, y coordinates in orbit */
v; /* True anomaly */
/* Compute mean elements */
M = [self revolution:( 356.0470 + 0.9856002585 * d )];
w = 282.9404 + 4.70935E-5 * d;
e = 0.016709 - 1.151E-9 * d;
/* Compute true longitude and radius vector */
E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
x = cosd(E) - e;
y = sqrt( 1.0 - e*e ) * sind(E);
*r = sqrt( x*x + y*y ); /* Solar distance */
v = atan2d( y, x ); /* True anomaly */
*lon = v + w; /* True solar longitude */
if ( *lon >= 360.0 )
*lon -= 360.0; /* Make it 0..360 degrees */
}
-(void) sun_RA_decAtDay:(double)d RA:(double*)RA decl:(double *)dec r:(double *)r
{
double lon, obl_ecl;
double xs, ys;
double xe, ye, ze;
/* Compute Sun's ecliptical coordinates */
//sunpos( d, &lon, r );
[self sunposAtDay:d longitude:&lon r:r];
/* Compute ecliptic rectangular coordinates */
xs = *r * cosd(lon);
ys = *r * sind(lon);
/* Compute obliquity of ecliptic (inclination of Earth's axis) */
obl_ecl = 23.4393 - 3.563E-7 * d;
/* Convert to equatorial rectangular coordinates - x is unchanged */
xe = xs;
ye = ys * cosd(obl_ecl);
ze = ys * sind(obl_ecl);
/* Convert to spherical coordinates */
*RA = atan2d( ye, xe );
*dec = atan2d( ze, sqrt(xe*xe + ye*ye) );
} /* sun_RA_dec */
#define sun_rise_set(year,month,day,lon,lat,rise,set) \
__sunriset__( year, month, day, lon, lat, -35.0/60.0, 1, rise, set )
-(int)sunRiseSetForYear:(int)year month:(int)month day:(int)day longitude:(double)lon latitude:(double)lat
trise:(double *)trise tset:(double *)tset
{
return [self sunRiseSetHelperForYear:year month:month day:day longitude:lon latitude:lat altitude:(-35.0/60.0)
upper_limb:1 trise:trise tset:tset];
}
/*
#define civil_twilight(year,month,day,lon,lat,start,end) \
__sunriset__( year, month, day, lon, lat, -6.0, 0, start, end )
*/
-(int) civilTwilightForYear:(int)year month:(int)month day:(int)day longitude:(double)lon latitude:(double)lat
trise:(double *)trise tset:(double *)tset
{
return [self sunRiseSetHelperForYear:year month:month day:day longitude:lon latitude:lat altitude:-6.0
upper_limb:0 trise:trise tset:tset];
}
/*
#define nautical_twilight(year,month,day,lon,lat,start,end) \
__sunriset__( year, month, day, lon, lat, -12.0, 0, start, end )
*/
-(int) nauticalTwilightForYear:(int)year month:(int)month day:(int)day longitude:(double)lon latitude:(double)lat
trise:(double *)trise tset:(double *)tset
{
return [self sunRiseSetHelperForYear:year month:month day:day longitude:lon latitude:lat altitude:-12.0
upper_limb:0 trise:trise tset:tset];
}
/*
#define astronomical_twilight(year,month,day,lon,lat,start,end) \
__sunriset__( year, month, day, lon, lat, -18.0, 0, start, end )
*/
-(int) astronomicalTwilightForYear:(int)year month:(int)month day:(int)day longitude:(double)lon latitude:(double)lat
trise:(double *)trise tset:(double *)tset
{
return [self sunRiseSetHelperForYear:year month:month day:day longitude:lon latitude:lat altitude:-18.0
upper_limb:0 trise:trise tset:tset];
}
/***************************************************************************/
/* Note: year,month,date = calendar date, 1801-2099 only. */
/* Eastern longitude positive, Western longitude negative */
/* Northern latitude positive, Southern latitude negative */
/* The longitude value IS critical in this function! */
/* altit = the altitude which the Sun should cross */
/* Set to -35/60 degrees for rise/set, -6 degrees */
/* for civil, -12 degrees for nautical and -18 */
/* degrees for astronomical twilight. */
/* upper_limb: non-zero -> upper limb, zero -> center */
/* Set to non-zero (e.g. 1) when computing rise/set */
/* times, and to zero when computing start/end of */
/* twilight. */
/* *rise = where to store the rise time */
/* *set = where to store the set time */
/* Both times are relative to the specified altitude, */
/* and thus this function can be used to comupte */
/* various twilight times, as well as rise/set times */
/* Return value: 0 = sun rises/sets this day, times stored at */
/* *trise and *tset. */
/* +1 = sun above the specified "horizon" 24 hours. */
/* *trise set to time when the sun is at south, */
/* minus 12 hours while *tset is set to the south */
/* time plus 12 hours. "Day" length = 24 hours */
/* -1 = sun is below the specified "horizon" 24 hours */
/* "Day" length = 0 hours, *trise and *tset are */
/* both set to the time when the sun is at south. */
/* */
/**********************************************************************/
-(int)sunRiseSetHelperForYear:(int)year month:(int)month day:(int)day longitude:(double)lon latitude:(double)lat
altitude:(double)altit upper_limb:(int)upper_limb trise:(double *)trise tset:(double *)tset
{
double d, /* Days since 2000 Jan 0.0 (negative before) */
sr, /* Solar distance, astronomical units */
sRA, /* Sun's Right Ascension */
sdec, /* Sun's declination */
sradius, /* Sun's apparent radius */
t, /* Diurnal arc */
tsouth, /* Time when Sun is at south */
sidtime; /* Local sidereal time */
int rc = 0; /* Return cde from function - usually 0 */
/* Compute d of 12h local mean solar time */
d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
/* Compute local sideral time of this moment */
//sidtime = revolution( GMST0(d) + 180.0 + lon );
sidtime = [self revolution:[self GMST0:d] + 180.0 + lon];
/* Compute Sun's RA + Decl at this moment */
//sun_RA_dec( d, &sRA, &sdec, &sr );
[self sun_RA_decAtDay:d RA: &sRA decl:&sdec r:&sr];
/* Compute time when Sun is at south - in hours UT */
//tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
tsouth = 12.0 - [self rev180:sidtime - sRA] / 15.0;
/* Compute the Sun's apparent radius, degrees */
sradius = 0.2666 / sr;
/* Do correction to upper limb, if necessary */
if ( upper_limb )
altit -= sradius;
/* Compute the diurnal arc that the Sun traverses to reach */
/* the specified altitide altit: */
{
double cost;
cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
( cosd(lat) * cosd(sdec) );
if ( cost >= 1.0 )
rc = -1, t = 0.0; /* Sun always below altit */
else if ( cost <= -1.0 )
rc = +1, t = 12.0; /* Sun always above altit */
else
t = acosd(cost)/15.0; /* The diurnal arc, hours */
}
/* Store rise and set times - in hours UT */
*trise = tsouth - t;
*tset = tsouth + t;
return rc;
} /* __sunriset__ */
@end
#pragma mark - Private Implementation
@implementation EDSunriseSet(Private)
static const int kSecondsInHour= 60.0*60.0;
-(NSDate*)utcTime:(NSDateComponents*)dateComponents withOffset:(NSTimeInterval)interval
{
(self.calendar).timeZone = self.utcTimeZone;
return [[self.calendar dateFromComponents:dateComponents] dateByAddingTimeInterval:(NSTimeInterval)(interval)];
}
-(NSDateComponents*)localTime:(NSDate*)refDate
{
(self.calendar).timeZone = self.timezone;
// Return only hour, minute, seconds
NSDateComponents *dc = [self.calendar components:( NSCalendarUnitHour | NSCalendarUnitMinute | NSCalendarUnitSecond) fromDate:refDate] ;
return dc;
}
- (instancetype) init {
[super doesNotRecognizeSelector:_cmd];
return nil;
}
-(NSString *)description
{
return [NSString stringWithFormat:
@"Date: %@\nTimeZone: %@\n"
@"Local Sunrise: %@\n"
@"Local Sunset: %@\n"
@"Local Civil Twilight Start: %@\n"
@"Local Civil Twilight End: %@\n"
@"Local Nautical Twilight Start: %@\n"
@"Local Nautical Twilight End: %@\n"
@"Local Astronomical Twilight Start: %@\n"
@"Local Astronomical Twilight End: %@\n",
self.date.description, self.timezone.name,
self.localSunrise.description, self.localSunset.description,
self.localCivilTwilightStart, self.localCivilTwilightEnd,
self.localNauticalTwilightStart, self.localNauticalTwilightEnd,
self.localAstronomicalTwilightStart, self.localAstronomicalTwilightEnd
];
}
#pragma mark - Calculation methods
-(void)calculateSunriseSunset
{
// Get date components
(self.calendar).timeZone = self.timezone;
NSDateComponents *dateComponents = [self.calendar components:( NSCalendarUnitYear | NSCalendarUnitMonth | NSCalendarUnitDay ) fromDate:self.date];
// Calculate sunrise and sunset
double rise=0.0, set=0.0;
[self sunRiseSetForYear:(int)dateComponents.year month:(int)dateComponents.month day:(int)dateComponents.day longitude:self.longitude latitude:self.latitude
trise:&rise tset:&set ];
NSTimeInterval secondsRise = rise*kSecondsInHour;
NSTimeInterval secondsSet = set*kSecondsInHour;
self.sunrise = [self utcTime:dateComponents withOffset:(NSTimeInterval)secondsRise];
self.sunset = [self utcTime:dateComponents withOffset:(NSTimeInterval)secondsSet];
self.localSunrise = [self localTime:self.sunrise];
self.localSunset = [self localTime:self.sunset];
}
-(void)calculateTwilight
{
// Get date components
(self.calendar).timeZone = self.timezone;
NSDateComponents *dateComponents = [self.calendar components:( NSCalendarUnitYear | NSCalendarUnitMonth | NSCalendarUnitDay ) fromDate:self.date];
double start=0.0, end=0.0;
// Civil twilight
[self civilTwilightForYear:(int)dateComponents.year month:(int)dateComponents.month day:(int)dateComponents.day longitude:self.longitude latitude:self.latitude
trise:&start tset:&end ];
self.civilTwilightStart = [self utcTime:dateComponents withOffset:(NSTimeInterval)(start*kSecondsInHour)];
self.civilTwilightEnd = [self utcTime:dateComponents withOffset:(NSTimeInterval)(end*kSecondsInHour)];
self.localCivilTwilightStart = [self localTime:self.civilTwilightStart];
self.localCivilTwilightEnd = [self localTime:self.civilTwilightEnd];
// Nautical twilight
[self nauticalTwilightForYear:(int)dateComponents.year month:(int)dateComponents.month day:(int)dateComponents.day longitude:self.longitude latitude:self.latitude
trise:&start tset:&end ];
self.nauticalTwilightStart = [self utcTime:dateComponents withOffset:(NSTimeInterval)(start*kSecondsInHour)];
self.nauticalTwilightEnd = [self utcTime:dateComponents withOffset:(NSTimeInterval)(end*kSecondsInHour)];
self.localNauticalTwilightStart = [self localTime:self.nauticalTwilightStart];
self.localNauticalTwilightEnd = [self localTime:self.nauticalTwilightEnd];
// Astronomical twilight
[self astronomicalTwilightForYear:(int)dateComponents.year month:(int)dateComponents.month day:(int)dateComponents.day longitude:self.longitude latitude:self.latitude
trise:&start tset:&end ];
self.astronomicalTwilightStart = [self utcTime:dateComponents withOffset:(NSTimeInterval)(start*kSecondsInHour)];
self.astronomicalTwilightEnd = [self utcTime:dateComponents withOffset:(NSTimeInterval)(end*kSecondsInHour)];
self.localAstronomicalTwilightStart = [self localTime:self.astronomicalTwilightStart];
self.localAstronomicalTwilightEnd = [self localTime:self.astronomicalTwilightEnd];
}
-(void)calculate
{
[self calculateSunriseSunset];
[self calculateTwilight];
}
@end
#pragma mark - Public Implementation
@implementation EDSunriseSet
#pragma mark - Initialization
-(EDSunriseSet*)initWithDate:(NSDate*)date timezone:(NSTimeZone*)tz latitude:(double)latitude longitude:(double)longitude {
self = [super init];
if( self )
{
self.latitude = latitude;
self.longitude = longitude;
self.timezone = tz;
self.date = date;
self.calendar = [[NSCalendar alloc] initWithCalendarIdentifier:EDGregorianCalendar];
self.utcTimeZone = [NSTimeZone timeZoneWithAbbreviation:@"UTC"];
[self calculate];
}
return self;
}
+(EDSunriseSet*)sunrisesetWithDate:(NSDate*)date timezone:(NSTimeZone*)tz latitude:(double)latitude longitude:(double)longitude {
return [[EDSunriseSet alloc] initWithDate:date timezone:tz latitude:latitude longitude:longitude];
}
@end