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:
What should I do with
sun
arguments? To force all of them to beDouble
?
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?lat
andlon
could beFloat
orDouble
or Fractional. What should I choose and why?zenit
could be Fractional, but could it also could be string taking "civil", or 1-4 and inside 1 = official 90.5?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 nestedwhere
orlet
. Each function defined insun
could be accessed through wholesun
. Is this considered good practice?End result would be
(rise, set)
times. My thoughts were to create as much as possible common functions, and just feed themt_rise
andt_set
. Is it ok? Does anything better come to mind?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.
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 beNum 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/
requireFractional
typeclass. ButIntegral
havequot
andrem
. How do you combine those two and createnn
than can takeNum
? My first thoughts were to create 2 functions:nni
for Integral andnnf
for Fractional. But is there a better way?