Calculating Equation of Center in Swift

Question

I am working on a pet project that involves calculating sunrise / sunset data. I am struggling with implementing the following formula in Swift:

Equation of Center:

C = (1.9148 * sin(meanSolarAnomaly)) + (0.0200 * sin(2 * meanSolarAnomaly)) + (0.0003 * sin(3 * meanSolarAnomaly))

Here is the answer I should be getting for my given Lat/Lon:

C = 1.9148 * sin(18.30143135945) + 0.0200 * sin(2 * 18.30143135945) + 0.0003 * sin(3 * 18.30143135945) = 0.61344892821988

Here is my code, which is not giving the correct value as the final value:

// meanSolarAnomaly has a value of 18.30143135945036 at this point
let center = (1.9148 * sin(meanSolarAnomaly)) + (0.0200 * sin(2*meanSolarAnomaly)) + (0.0003 * sin(3*meanSolarAnomaly))

My code says center = -1.015867439183884 while the correct center = 0.61344892821988

I have double and triple checked the equation but I can't seem to spot my error. I am hoping it is a simple syntax mistake but will be embarrassed if it is...

I am working from the equation and answers provided here.

EDIT Here is the complete code:

//: Playground - noun: a place where people can play

import UIKit
//calculator to determine what time of day Sunset / Sunrise will occur
func jdFromDate(date : NSDate) -> Double {
    let JD_JAN_1_1970_0000GMT = 2440587.5
    return JD_JAN_1_1970_0000GMT + date.timeIntervalSince1970 / 86400
}

func dateFromJd(jd : Double) -> NSDate {
    let JD_JAN_1_1970_0000GMT = 2440587.5
    return  NSDate(timeIntervalSince1970: (jd - JD_JAN_1_1970_0000GMT) * 86400)
}
let julianOffset = 2451545 as Double
let refDateFormatter = NSDateFormatter()
refDateFormatter.dateFormat = "MM-dd-yyyy"

let today = NSDate()

let julianDaysToToday = round(jdFromDate(today))

//get the lat/lon variables set (Tampa in example)
let lat = 27.9681
let lon = 82.4764

//now we need to calculate julian cycle
let nRaw = (julianDaysToToday - julianOffset - 0.0009) - (lon / 360)
let n = round(nRaw)

//n now contains the julian cycle
//next we must calculate the julian date of solar noon (approximately)
//J* = 2451545 + 0.0009 + (lw/360) + n

let jSolarNoon = julianOffset + 0.0009 + (lon/360) + n

//next calculate the mean solar anomaly
//M = [357.5291 + 0.98560028 * (J* - 2451545)] mod 360
let meanSolarAnomaly = (357.5291 + 0.98560028 * (jSolarNoon - julianOffset)) % 360

//next calculate the equation of center

let center = (1.9148 * sin(meanSolarAnomaly)) + (0.0200 * sin(2*meanSolarAnomaly)) + (0.0003 * sin(3*meanSolarAnomaly))

//Now, using Center  and Mean, calculate the ecliptical longitude of the sun.
//λ = (M + 102.9372 + C + 180) mod 360
let eclLonOfSun = (meanSolarAnomaly + 102.9372 + center + 180) % 360

//now we can finally get an accurate julian date for solar noon
let jTransit = jSolarNoon + (0.0053 * sin(meanSolarAnomaly)) - (0.0069 * sin(2 * eclLonOfSun))

//To calculate the hour angle we need to find the declination of the sun
//δ = arcsin( sin(λ) * sin(23.45) )
let declinationOfSun = asin(sin(eclLonOfSun) * sin(23.45))

//now calculate the hour angle
//H = arccos( [sin(-0.83) - sin(ln) * sin(δ)] / [cos(ln) * cos(δ)] )
let hourCosNum = sin(-0.83) - sin(lat) * sin(declinationOfSun)
let hourDenom = cos(lat)*cos(declinationOfSun)
let hourAngle = acos(hourCosNum)/hourDenom

//time to go back through the approximation again using the hour angle
let finalJulianApproximation = 2451545 + 0.0009 + ((hourAngle + lon)/360) + n

//The values of M and λ from above don't really change from solar noon to sunset, so there is no need to recalculate them before calculating sunset.
let jSet = finalJulianApproximation + (0.0053 * sin(meanSolarAnomaly)) - (0.0069 * sin(2*eclLonOfSun))

let sunset = dateFromJd(jSet)

Answer 1

As @kennytm has suggested, mean anomaly (of the sun or anything else) is an angle. Angles in Swift (and C, from which the maths libraries come) are all radians, whilst astronomers talk in degrees. Here is your code in Playground:

var meanSolarAnomaly = 18.30143135945036
var c = (1.9148 * sin(meanSolarAnomaly)) + (0.0200 * sin(2 * meanSolarAnomaly)) + (0.0003 * sin(3 * meanSolarAnomaly)) 
// = -1.01586743918389 - wrong answer

meanSolarAnomaly = meanSolarAnomaly * M_PI / 180.0 
// Convert it to radians

c = (1.9148 * sin(meanSolarAnomaly)) + (0.0200 * sin(2 * meanSolarAnomaly)) + (0.0003 * sin(3 * meanSolarAnomaly)) 
// = 0.6134489282198807 - right answer