Calculating sun rise/set time at various places on Earth

Question

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 < 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 [0,360) by adding/subtracting 360
  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
  cosH t = (cos(zenit*deg2rad) - (sinDec t * sin(lat*deg2rad))) / (cosDec t * cos(lat * deg2rad))
--  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
--9. adjust back to UTC
--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?

Answer 1

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(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

Answer 2

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