[Gambas-user] searching sunset/sunrise function

[Ron]

Ron schreef:
> Hi,
>
> Does anybody know if  a sunset/sunrise calculator functions/snippet in 
> gambas exists by any chance?
>
> I have found some vb and vbscript snippets, but porting them is not 
> easy, or they give wrong results.
>
> Regards,
> Ron.
>
>   

Ok replying to my own question.

I found this http://www.drbeat.li/php/source.php?src=php/sunrise.inc.php

And hacked together this code:
----------
' Gambas module file

PUBLIC FUNCTION CalcSunTimes(lon AS Float, lat AS Float, timezone AS 
Float, isRise AS Boolean)
  DIM a, b, c, d, e, f, g, j, k, l, m, p, q, R, s, t, u, v AS Float
  DIM yday, ihour, imin, isec AS Integer

  'twilight setting
  'effective' => -0.0145439 (sunrise/sunset)
  'civil' => -0.104528 (civil twilight)
  'nautical' => -0.207912 (nautical twilight)
  'astronomical' => -0.309017 (astronomical twilight)
  R = -0.0145439

  'multiples of Pi
  a = 0.5 * Pi
  b = Pi
  c = 1.5 * Pi
  d = 2 * Pi

  'convert coordinates and timezone to radians
  e = lat * b / 180
  f = lon * b / 180
  g = timezone * d / 24

  'sunrise
  j = IIf(isRise, a, c)
 
  'calculate day of year
  yday = DateDiff("01/01/" & Year(Now()), Now(), gb.day)

  k = yday + (j - f) / d
  l = k * 0.017202 - 0.0574039 'solar mean anomoly
  m = l + 0.0334405 * Sin(l) 'solar true longitude
  m += 4.93289 + 3.49066E-4 * Sin(2 * l)

  'quadrant determination
  m = norm(m, d)
  IF ((m / a) - CInt(m / a) == 0) THEN
    m += 4.84814E-6
    p = Sin(m) / Cos(m) 'solar right ascension
    p = ATan2(0.91746 * p, 1)
  END IF

  'quadrant adjustment
  IF (m > c)
    p += d
  ELSE IF (m > a)
    p += b
  END IF

  q = 0.39782 * Sin(m) 'solar declination
  q /= Sqr(- q * q + 1)
  q = ATan2(q, 1)

  s = R - Sin(q) * Sin(e)
  s /= Cos(q) * Cos(e)
  IF (Abs(s) > 1) THEN PRINT "(Midnight Sun)"

  s /= Sqr(- s * s + 1)
  s = a - ATan2(s, 1)

  IF (isRise) THEN s = d - s

  t = s + p - 0.0172028 * k - 1.73364 'local apparent time
  u = t - f 'universal time
  v = u + g 'wall clock time

  'quadrant Determination
  v = norm(v, d)

  'scale from radians to hours
  v *= 24 / d

  'local time
  ihour = CInt(v)
  iMin = CInt((v - ihour) * 60)

  IF isRise THEN
    PRINT "Sunrise at ";
  ELSE
    PRINT "Sunset at ";
  END IF
 
  PRINT Format(ihour, "##") & ":" & Format(imin, "0#")
END

FUNCTION norm(a AS Float, b AS Float) AS Float
  WHILE (a < 0)
    a += b
  WEND
  WHILE (a >= b)
    a -= b
  WEND
  RETURN a
END
------------------