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
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)
]
This code looks ugly. Really ugly :(
Here are my questions:
What should I do with
sunarguments? 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?latandloncould be Float or Double or Fractional or... What to choose and why?zenitcould be Fractional, but 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 insun. After a while I figured out that I do not need multiple nestedwhereorlet.
Each function defined insuncould 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_riseandt_set
Is it ok? Anything better comes to your 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.
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?
