# Calculating sun rise/set time at various places on Earth

This is the algorithm for calculating sun rise/set time at various places on Earth. I took it as an example of multiple functions inside of one top function.

This is Ver 0.3:

fromDegree deg = deg * pi / 180

--nn :: RealFrac a => a -> a -> a -> a
nn value start end = let
width = end - start
offsetValue = value - start   -- value relative to 0
in (offsetValue - (fromIntegral (floor(offsetValue / width)) * width)) + start
-- + start to reset back to start of original range

{-
zenith:                Sun's zenith for sunrise/sunset
offical      = 90 degrees 50'
civil        = 96 degrees
nautical     = 102 degrees
astronomical = 108 degrees

longitude is positive for East and negative for West
-}

-- sun :: Double  -> Double -> Double -> Double -> Double -> Double -> Double
sun year month day lat lon zenit local = let

-- OK
-- 1. day of year
doy = n1 - (n2 * n3) + day - 30
where
n1 = fromIntegral (floor(275 * month / 9)) :: Double
n2 = fromIntegral (floor((month + 9) / 12)) :: Double
n3 = fromIntegral (1 + floor ( (year - fromIntegral(4 * floor(year / 4)) + 2) / 3)) :: Double

-- OK
-- 2. convert the longitude to hour value and calculate an approximate time
lngHour = lon / 15
t_rise = doy + ((6 - lngHour) / 24)
t_set = doy + ((18 - lngHour) / 24)

-- OK
-- 3. calculate the Sun's mean anomaly
m t = (0.9856 * t) - 3.289

-- OK
--4. calculate the Sun's true longitude, UGLY!
l t = out res
where
| x < 0 = x + 360
| x > 360 = x - 360
| otherwise = x
res = m t + (1.916 * sin(m t * deg2rad)) + (0.020 * sin(2 * m t * deg2rad)) + 282.634

l' t = nn stl 0 360     -- little bit better
where
stl = m t + (1.916 * sin(m t * deg2rad)) + (0.020 * sin(2 * m t * deg2rad)) + 282.634

-- test OK
m_rise = m t_rise
m_set = m t_set
stl1 = m_rise + (1.916 * sin(m_rise * deg2rad)) + (0.020 * sin(2 * m_rise * deg2rad)) + 282.634
stl = stl1 - 360
stl' = nn stl1 0 360
-- end test

-- OK
--5a. calculate the Sun's right ascension
ra' t = out res    -- ugly, not used any more
where
out x
| x < 0 = x + 360
| x > 360 = x - 360
| otherwise = x
res = toDegree $atan(0.91764 * tan(l t * deg2rad)) -- ??? --5b. right ascension value needs to be in the same quadrant as l ra t = (ra' t + (lQuadrant t - raQuadrant t)) /15 -- 5c... /15 where lQuadrant t = fromIntegral ((floor(l t / 90)) * 90) :: Double raQuadrant t = fromIntegral((floor(ra' t / 90)) * 90) :: Double ra' t = nn raas 0 360 -- better where raas = toDegree$ atan(0.91764 * tan(l t * deg2rad))

-- OK, done in 5b
--5c. right ascension value needs to be converted into hours
--  ra'' t = ra2 t / 15

-- ??? deg/rad conversion not needed?
--6. calculate the Sun's declination
sinDec t = 0.39782 * sin(l t * deg2rad)
cosDec t = cos(asin(sinDec t))

-- ??? what values are apropriate for this?
--7a. calculate the Sun's local hour angle
--  if (cosH >  1) sun never rises
--  if (cosH < -1) sun never sets

-- OK and 7b
--8. calculate local mean time of rising/setting

lmtr t = hr t + ra t - (0.06571 * t) - 6.622
where
hr t = (360 - toDegree (acos(cosH t))) / 15

lmts t = hs t + ra t - (0.06571 * t) - 6.622
where
hs t = (toDegree (acos(cosH t))) /15

--OK
--NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
utr t = out res
where
out x
| x < 0 = x + 24
| x > 24 = x - 24
| otherwise = x
res =  lmtr t - lngHour

uts t = out res
where
out x
| x < 0 = x + 24
| x > 24 = x - 24
| otherwise = x
res = lmts t - lngHour

utr' t = (nn (lmtr t - lngHour) 0 24)  + local
uts' t = (nn (lmts t - lngHour) 0 24) + local

-- OK done in 9.
--10. convert UT value to local time zone of latitude/longitude
--  localT = UT + localOffset

in mapM_ putStrLn [   -- or print
" day: " ++ show doy,
" lngHour: " ++ show lngHour,
" t rise: " ++ show t_rise,
" tset: " ++ show t_set,
" m rise: " ++ show m_rise,
" m set: " ++ show m_set,
" stl rise: " ++ show stl,
" stl' rise: " ++ show stl',
" l rise: " ++ show (l t_rise),
" l' rise: " ++ show (l' t_rise),
" l set: " ++ show (l t_set),
" l' set: " ++ show (l' t_set),
" ra rise: " ++ show (ra t_rise),
" ra set: " ++ show (ra t_set),
" sinDec set: " ++ show (sinDec t_set),
" cosDec set: " ++ show (cosDec t_set),
" cosH set: " ++ show (cosH t_set),
" cosH rise: " ++ show (cosH t_rise),
" lmtr r : " ++ show (lmtr t_rise),
" lmts s: " ++ show (lmts t_set),
" ut rise : " ++ show (utr t_rise),
" ut set: " ++ show (uts t_set),
" utr' rise : " ++ show (utr' t_rise),
" uts' set: " ++ show (uts' t_set)

]


Here are my questions:

1. What should I do with sun arguments? To force all of them to be Double?
First three arguments "natural" type would be Integral. But is it ok to have function with mixed type arguments? What are Haskell's convention on that?

2. lat and lon could be Float or Double or Fractional. What should I choose and why?

3. zenit could be Fractional, but could it also could be string taking "civil", or 1-4 and inside 1 = official 90.5?

4. What is the best or easiest way to put together a lot of complex functions like in this example? From where do you start? I didn't want to pollute global name space, so I put them all in sun. After a while I figured out that I do not need multiple nested where or let. Each function defined in sun could be accessed through whole sun. Is this considered good practice?

5. End result would be (rise, set) times. My thoughts were to create as much as possible common functions, and just feed them t_rise and t_set. Is it ok? Does anything better come to mind?

6. Yes, I know it gives erroneous result. But that's not my main concern. My biggest issue is to find out "proper way" of doing complex things like this. "OK" are just my markers for what I think is correct.

7. nn normalizes any number to an arbitrary range by assuming the range wraps around when going below min or above max. start, end: -180,180; 0,360; 0,24; -Math.PI, Math.PI. It was a good function in JavaScript. How can I make it be Num a => a -> a -> a -> a?

--nn :: RealFrac a => a -> a -> a -> a
nn value start end = let
width = end - start
offsetValue = value - start   -- value relative to 0
in (offsetValue - (fromIntegral (floor(offsetValue / width)) * width)) + start


floor and / require Fractional typeclass. But Integral have quot and rem. How do you combine those two and create nn than can take Num? My first thoughts were to create 2 functions: nni for Integral and nnf for Fractional. But is there a better way?

first thing good question and good observation - the code is ugly ;-), but here it starts getting better; You have a lot duplicate functionality and plain code duplication.

Here are my questions:

What should I do with sun arguments? To force all of them to be Double? First three arguments "natural" type would be Integral. But is it ok to have function with mixed type arguments? What are Haskell's convention

Best way I can think about it is to give the natural types and adapt the functions, as good as possible.

type Radians = Double
type Degree  = Double

type Year  = Int
type Month = Int
type Day   = Int
type Time  = Double

type Longitude = Double
type Latitude  = Double


lat and lon could be Float or Double or Fractional or... What to choose and why?

I chose Double just for simplicity.

zenit could be Fractional, but it also could be string taking "civil", or 1-4 and inside 1 = official 90.5...

At first try to do it statically and provide functionality later on.

What is the best or easiest way to put together a lot of complex functions like in this example? From where do you start? I didn't want to pollute global name space, so I put them all in sun. After a while I figured out that I do not need multiple nested where or let. Each function defined in sun could be accessed through whole sun. Is this considered good practice?

At first you do not litter namespace with global functions !!! Not at all, you just make a module and snap - no namespace problem - export only your function of choice and noone will care if you have a thousand helper functions. You do this with:

module SunPos (sun) where


And secondly, please use type annotations for top level functions - it is so much more easy to read and reason about your code, and the type checker kicks in. With any decent editor you have a syntax checking plugin to help you prevent errors, in vim it is Syntastic I know, others like emacs should have one too.

Now for the matter of complexity, just start simple - all simple functions. toDegree is a great example, doy too. almost everything you put into your big function sun can be extracted. This makes it also easier to write tests for your function and not be bugged with print "foo", but have automated tests with HUnit or Quickcheck, so if you change any funcitonality it can be tested and you know where to look for bugs.

It is really hard to find a bug in a 200 line of code function, but in 2 lines of code… . And try to restrict yourself to 80 characters of textwidth, as for example here on stackexchange you have to scroll, code with more than 80 characters per line, which is annoying, and easily solved.

Next thing in Haskell it is common practise to use CamelCase instead of underscores. And please use names in your functions that are a bit more self explaining. Today editors help you with autocompletion and your harddrive has no problem with a few extra characters to memorize. m -> meanAnomaly for example, l -> trueLong or nn -> ????, I changed a few but not nearly enough.

so here is what i have corrected by now - it is not a complete but a compiling code i think.

module SunPos (sun) where

type Degree  = Double

type Year  = Int
type Month = Int
type Day   = Int
type Time  = Double

type Longitude = Double
type Latitude  = Double

fromDegree deg = deg * pi / 180-}

nn :: Double -> Double -> Double -> Double
nn value start end = offsetValue - (windingNum * width) + start
where width       = end - start
offsetValue = value - start   -- value relative to 0
windingNum  = fromIntegral (floor (offsetValue / width)::Int)
-- + start to reset back to start of original range

{-
zenith:                Sun's zenith for sunrise/sunset
offical      = 90 degrees 50'
civil        = 96 degrees
nautical     = 102 degrees
astronomical = 108 degrees

longitude is positive for East and negative for West
-}

-- | 1. day of year
doy :: Year -> Month -> Day -> Double
doy year month day = fromIntegral (n1 - (n2 * n3) + day - 30)
where n1 =  275 * month div  9
n2 = (month + 9) div 12
n3 = 1 + ((year - (4 * (year div 4)) + 2) div 3)

meanAnomaly :: Time -> Time
meanAnomaly t = (0.9856 * t) - 3.289

trueLong :: Time -> Double
trueLong t = out res
where res = stl t

out :: Double -> Double
out x = nn x 0.0 360.0

stl :: Time -> Double
stl t = mA + (1.916 * sin(2*deg2rad*mA)) + (0.020 * sin(deg2rad*mA)) + 282.634
where  mA  = meanAnomaly t

l' :: Time -> Double
l' t          = nn (stl t) 0 360     -- little bit better
--
-- | calculate the Sun's right ascension ra [0,360) by adding/subtracting 360
ra :: Time -> Double
ra t = (ra' + (lQuadrant - raQuadrant)) /15    -- 5c... /15
where lQuadrant  = fromIntegral (floor(trueLong t / 90) * 90 ::Int)
raQuadrant = fromIntegral(floor(ra' / 90) * 90 ::Int)
ra' = out raas
raas = toDegree $atan(0.91764 * tan(trueLong t * deg2rad)) -- | local mean time of rising/setting lmtr :: Time -> Double -> Latitude -> Double lmtr t zenit lat = hr + ra t - (0.06571 * t) - 6.622 where hr = (360 - toDegree (acos cosHour')) / 15 cosHour' = cosHour t zenit lat lmts :: Time -> Double -> Latitude -> Double lmts t zenit lat = hs + ra t - (0.06571 * t) - 6.622 where hs = toDegree (acos cosHour') /15 cosHour' = cosHour t zenit lat -- calculate the cosine of Sun's local hour angle cosHour :: Time -> Double -> Latitude -> Double cosHour t zenit lat = (cZen - sLat ) / (cosDec t * cos(lat * deg2rad)) where cZen = cos(zenit*deg2rad) sLat = sinDec t * sin(lat*deg2rad) sinDec :: Time -> Double sinDec t = 0.39782 * sin(trueLong t * deg2rad) cosDec :: Time -> Double cosDec t = cos(asin(sinDec t)) longHour :: Longitude -> Double longHour lon = lon / 15 sun :: Year -> Month -> Day -> Latitude -> Longitude -> Double -> Double -> String sun year month day lat lon zenit local = unlines [ " day: " ++ show doy', " longHour: " ++ show (longHour lon), " t rise: " ++ show timeOfRise, " tset: " ++ show timeOfSet, " m rise: " ++ show meanAnomalyOfRise, " m set: " ++ show meanAnomalyOfSet, " stl' rise: " ++ show stl', " l rise: " ++ show (trueLong timeOfRise), " l' rise: " ++ show (l' timeOfRise), " l set: " ++ show (trueLong timeOfSet), " l' set: " ++ show (l' timeOfSet), " ra rise: " ++ show (ra timeOfRise), " ra set: " ++ show (ra timeOfSet), " sinDec set: " ++ show (sinDec timeOfSet), " cosDec set: " ++ show (cosDec timeOfSet), " lmtr: " ++ show lR, " lmts: " ++ show lS, " ut rise : " ++ show utr, " ut set: " ++ show uts] where doy' = doy year month day timeOfRise = doy' + (( 6 - longHour lon) / 24) timeOfSet = doy' + ((18 - longHour lon) / 24) meanAnomalyOfRise = meanAnomaly timeOfRise meanAnomalyOfSet = meanAnomaly timeOfSet stl1 = meanAnomalyOfRise + stl meanAnomalyOfRise stl' = nn stl1 0 360 utr = nn lR 0 24 + local lR = lmtr timeOfRise zenit lat - longHour lon uts = nn lS 0 24 + local lS = lmts timeOfSet zenit lat - longHour lon  Haskell is really beautiful and really ugly language at the same time. It's pretty useless in current state (everything is a function) simply because not everything is a function... What I did here was just finding common patterns and making functions for that patterns. Like a tree. Leafs are almost universal functions like normalize number or degree transformations. Branches become functions used just in sun. And at the end there is sun function which combines branches and leafs. It took me some real time to refactor and debug sun function. Fully functional. Full module. Use it for what ever you want :) module Sky (sun, toHourMin) where --import MyLib.Functions -- nn nni, toDegree, deg2rad, rad2deg, doy import Data.Fixed -- mod' type Radians = Double type Degree = Double {-type Year = Integral type Month = Integral type Day = Integral -} type Longitude = Degree type Latitude = Degree type Zenith = Degree type Time = Double type LocalTime = Double type GMTtime = Double ------------------ common helper functions------------------------------------- -- normalise number x nn start end x = let width = end - start offsetValue = x - start -- x relative to 0 in (offsetValue - (fromIntegral (floor(offsetValue / width)) * width)) + start norm360 :: Degree -> Degree norm360 x = nn 0 360 x norm24 :: Time -> Time norm24 x = nn 0 24 x toRadian deg = deg * pi / 180 toDegree rad = rad * 180 / pi deg2rad = pi/180 rad2deg = 180 / pi -- upgrade it to Real a arcAngleRight max start end | (diff <= 0) = -diff | otherwise = max - diff where m2 = max/2 diff = (start + m2 - end) mod' max - m2; doy year month day = n1 - (n2 * n3) + day - 30 where n1 = quot (275 * month) 9 -- infixl 7: n1 = 275 * month quot 9 n2 = quot (month + 9) 12 n3 = 1 + quot (year - 4 * quot year 4 + 2) 3 toHourMin time = hour ++ ":" ++ min where (h, m) = properFraction time -- (Integral, Realfrac) 1.99 => (1, 0.99) hour = show h m' = floor (m * 60) min | m' >= 10 = show m' | otherwise = "0" ++ show m' ------------------ common sun functions --------------------------------------- longitudeHour :: Longitude -> Time longitudeHour l = l / 15 sunMeanAnomaly :: Time -> Degree -- 3. sunMeanAnomaly time = (0.9856 * time) - 3.289 sunTrueLongitude :: Time -> Degree -- 4. sunTrueLongitude time = norm360 stl where m = sunMeanAnomaly time stl = m + (1.916 * sin (m*deg2rad)) + (0.02 * sin((2*m)*deg2rad)) + 282.634 -- stl = m + (1.916 * sin m) + (0.02 * sin(2*m)) + 282.634 sunRightAscention :: Time -> Degree sunRightAscention time = ra / 15 -- 5c where stl = sunTrueLongitude time nra = norm360$ toDegree $atan(0.91764 * tan(stl*deg2rad)) -- 5a lQ = fromIntegral (floor(stl/90) * 90) -- 5b raQ = fromIntegral (floor(nra/90) * 90) -- 5b ra = nra + (lQ - raQ) -- 5b -------------------- require common sun f ------------------------------------- -- if (cosH > 1) sun never rises -- if (cosH < -1) sun never sets --7a. calculate the Sun's local hour angle; aka: cosH sunLocalHourAngle :: Time -> Latitude -> Zenith -> Degree sunLocalHourAngle time lat zenith = (cos(zenith*deg2rad) - (sinDec * sin(lat*deg2rad))) / (cosDec * cos(lat*deg2rad)) where stl = sunTrueLongitude time -- sun declination sinDec = 0.39782 * sin(stl*deg2rad) -- 6 cosDec = cos(asin sinDec) -- 6 {- Local mean rising / setting time are real local times of dusk / sunset but many cities and countries incorporate manu time zones in just one. othervise there would be 4 min time difference for every longitude degree. 360/24H = 15 15 degrees of longitude separation between each of the 24 primary time zones. 15 / 60 = 0.25 or 1 = 4min realLocalRiseTime = toHourMin . norm24$ lmtRise aproxTimeRise

Because we live in politicly corrected timezones, this is not usefull.
-}
--8. local mean Rising time
lmtRise time lat zen = hRise + ra - (0.06571 * time) - 6.622
where
ra = sunRightAscention time
cosH = sunLocalHourAngle time lat zen
hRise = (360 - toDegree (acos cosH)) / 15   -- 7b.

-- local mean Setting time
lmtSet time lat zen = hSet + ra - (0.06571 * time) - 6.622
where
ra = sunRightAscention time
cosH = sunLocalHourAngle time lat zen
hSet = (toDegree (acos cosH )) /15      -- 7b.

{-
zenithh:                Sun's zenithh for sunrise/sunset
offical      = 90 degrees 50'
civil        = 96 degrees
nautical     = 102 degrees
astronomical = 108 degrees

longitude is positive for East and negative for West

return double or  (hh,mm)
-}

--sun ::  Integral a => a ->  a ->  a -> Latitude -> Longitude -> Zenith -> Time -> Time -> ()
sun year month dayy lat lon zenith timeZone dST = let

day :: Time
day = fromIntegral $doy year month dayy lonHour :: Time lonHour = longitudeHour lon aproxTimeRise :: Time aproxTimeRise = day + ((6.0 - lonHour) / 24.0) aproxTimeSet :: Time aproxTimeSet = day + ((18.0 - lonHour) / 24.0) -- if (sunLocalHourAngle aproxTimeRise > 1) sun never rises -- if (sunLocalHourAngle aproxTimeSet < -1) sun never sets lRise :: Time lRise | noRise > 1 = 0 -- sun never rises | noRise < -1 = 24 -- sun never sets | otherwise = lmtRise aproxTimeRise lat zenith where noRise = sunLocalHourAngle aproxTimeRise lat zenith lSet :: Time lSet | noRise > 1 = 0 -- sun never rises | noRise < -1 = 24 -- sun never sets | otherwise = lmtSet aproxTimeSet lat zenith where noRise = sunLocalHourAngle aproxTimeSet lat zenith --9. adjust back to UTC/GMT 0-24 gmtRise :: GMTtime gmtRise = norm24 (lRise - lonHour) gmtSet :: GMTtime gmtSet = norm24$ lSet - lonHour

localRise :: LocalTime
localRise = norm24 $gmtRise + timeZone + dST localSet :: LocalTime localSet = norm24$ gmtSet + timeZone + dST

outT :: Time -> Time -> Time
outT _ 0 = 0           -- sun never rises
outT _ 24 = 24        -- sun never sets
outT time _ = time

dayLength :: Time
dayLength
| lRise ==0 && lSet == 0 = 0             -- sun never rises
| lRise == 24 && lSet == 24 = 24       -- sun never sets
| lRise /= 0 && lRise /= 24 && lSet /= 0 && lSet /= 24 = arcAngleRight 24 gmtRise gmtSet
| otherwise = 0

in ((outT localRise lRise, outT localSet lSet, dayLength),
(outT gmtRise lRise, outT gmtSet lSet, dayLength))


How to call it:

ny =  ((r,s,l),(ur, us, p))
-- or: toHourMin r ++ "   set: " ++ toHourMin s ++ "  l: " ++ toHourMin l
where
((r,s,l),(ur, us, p)) = sun 2012 7 8 40.7 (-74) 90.5 (-5) 1

• @200_success: If you've asked me that 3 years ago, my answer would be more eloquent. But now... Haskel is really beautiful and really ugly language at the same time. It's pretty useless in current state (everything is a function) simply because not everything is a function... What I did here was just finding common patterns and making functions for that patterns. Like a tree. Leafs are almost universal functions like normalize number or degree transformations. Branches become functions used just in sun. And at the end there is sun function which combines branches and leafs. – CoR Jun 9 '15 at 8:39
• If there were no modules, this approach would be very wrong. Haskel didn't have modules from beginning. Modules were added later. epsilonhalbe answer helped a lot. He pointed me to avoid global pollution, to re-factor the code. Doing that I created something like tree structure. After a while almost all leafs ended in my functions.hs file because they are universal. But helper functions stayed in sun module because they will probably be used only here and nowhere else. – CoR Jun 9 '15 at 8:45