This forum uses cookies
This forum makes use of cookies to store your login information if you are registered, and your last visit if you are not. Cookies are small text documents stored on your computer; the cookies set by this forum can only be used on this website and pose no security risk. Cookies on this forum also track the specific topics you have read and when you last read them. Please confirm that you accept these cookies being set.

How to calculate sun position?
#11
Try this one instead:
Code:
function suncalc(lat, lng, time)
  -- sun calculations are based on http://aa.quae.nl/en/reken/zonpositie.html formulas
  local PI   = math.pi
  local sin  = math.sin
  local cos  = math.cos
  local tan  = math.tan
  local asin = math.asin
  local atan = math.atan2
  local acos = math.acos
  local deg  = math.deg
  local rad  = PI / 180
  local e = rad * 23.4397 -- obliquity of the Earth
  local daysec = 60 * 60 * 24
  local J1970 = 2440588
  local J2000 = 2451545

  local function toDays(time)
    return time / daysec - 0.5 + J1970 - J2000
  end

  local function rightAscension(l, b)
    return atan(sin(l) * cos(e) - tan(b) * sin(e), cos(l))
  end

  local function declination(l, b)
    return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l))
  end

  local function azimuth(H, phi, dec)
    return atan(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi))
  end

  local function altitude(H, phi, dec)
    return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H))
  end

  local function siderealTime(d, lw)
    return rad * (280.16 + 360.9856235 * d) - lw
  end

  local function astroRefraction(h)
    if h < 0 then -- the following formula works for positive altitudes only.
      h = 0 -- if h = -0.08901179 a div/0 would occur.
    end

    -- formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
    -- 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
    return 0.0002967 / math.tan(h + 0.00312536 / (h + 0.08901179))
  end

  -- general sun calculations
  local function solarMeanAnomaly(d)
    return rad * (357.5291 + 0.98560028 * d)
  end

  local function eclipticLongitude(M)
    local C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)) -- equation of center
    local P = rad * 102.9372 -- perihelion of the Earth
    return M + C + P + PI
  end

  local function sunCoords(d)
    local M = solarMeanAnomaly(d)
    local L = eclipticLongitude(M)
    return declination(L, 0), rightAscension(L, 0)
  end

  local lw      = rad * -lng
  local phi     = rad * lat
  local d       = toDays(time or os.time())
  local dec, ra = sunCoords(d)
  local H       = siderealTime(d, lw) - ra

  local alt, az = altitude(H, phi, dec), azimuth(H, phi, dec)
  return deg(alt), 180 + deg(az)
end
Reply


Messages In This Thread
How to calculate sun position? - by buuuudzik - 13.12.2015, 16:23
RE: How to calculate sun position? - by admin - 14.12.2015, 08:07
RE: How to calculate sun position? - by admin - 14.12.2015, 10:20
RE: How to calculate sun position? - by admin - 15.12.2015, 13:24
RE: How to calculate sun position? - by admin - 21.12.2015, 07:08
RE: How to calculate sun position? - by admin - 11.05.2017, 07:49
RE: How to calculate sun position? - by admin - 18.04.2018, 07:26
RE: How to calculate sun position? - by admin - 02.07.2020, 12:43

Forum Jump: