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
toDegree rad = rad * 180 / pi

deg2rad = pi/180; -- deg2rad
rad2deg = 180 / pi -- rad2deg

--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
out x -- adjust (0,360) by adding/subtracting 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 *

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 Radians = Double
type Degree = Double

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

type Longitude = Double
type Latitude = Double

{-fromDegree :: Degree -> Radians
fromDegree deg = deg * pi / 180-}

toDegree :: Radians -> Degree
toDegree rad = rad * 180 / pi

deg2rad :: Double
deg2rad = pi/180; -- deg2rad
-- rad2deg = 180 / pi -- rad2deg

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
-- | adjust (0,360) by adding/subtracting 360

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

stl :: Time -> Double
stl t = mA + (1.916 sin(2deg2radmA)) + (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

StackExchange Code Review Q#13413, answer score: 7