# How to improve these Haskell functions?

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.
(if someone figure out better title, please edit)

This is what I have done so far. Ver 0.3. Questions are below the code.

``````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)

]
``````

This code looks ugly. Really ugly :(

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 or... What to choose and why?

3. `zenit` could be Fractional, but 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? Anything better comes to your 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 there has to be 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))
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
``````
-
Thanks :) You really helped me. –  CoR Jul 11 '12 at 1:19

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 Data.Fixed   -- mod'

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

-- 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 =
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
``````
-