[Gambas-user] searching sunset/sunrise function

[Ron]

>
> If you have the snippets, I'm interested in converting to gambas
>
>
> ---------- Original Message -----------
> From: Ron <ron at ...1740...>
> To: mailing list for gambas users <gambas-user at lists.sourceforge.net>
> Sent: Mon, 19 May 2008 12:01:33 +0200
> Subject: [Gambas-user] searching sunset/sunrise function
>
>   
>> 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.
>>     

Nando,

I have attached the module made so far.

CalcSunTimes() is returning fair results, as far as I an see, but 
haven't checked it much, you're welcome to make it better. (see link 
above routine)
CalcMoonPhase() is the beginning of converting code (from link in 
module), but the phase it returns is false, puzzle as to why, maybe 
usage of wrong datatype, or round()s? would be nice to get this fixed.

Alas I'm no Gambas nor Math/Astronomy wizard.

Regards,
Ron.


-------------- next part --------------
' Gambas module file
' adapted from source found here http://www.drbeat.li/php/source.php?src=php/sunrise.inc.php
' left out UTC code part
' Time_Sunrise = MSunMoon.CalcSunTimes(ilongitude, ilatitude, itimezone, 1, 0)
' Time_Sunset = MSunMoon.CalcSunTimes(ilongitude, ilatitude, itimezone, 0, 0)
' PRINT "Sunrise at " & Time_Sunrise & ", Sunset at " & Time_Sunset

PUBLIC FUNCTION CalcSunTimes(lon AS Float, lat AS Float, timezone AS Float, isRise AS Boolean, twilight AS Integer) AS String
  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 
  IF (twilight = 0) THEN R = -0.0145439 'effective for sunrise/sunset
  IF (twilight = 1) THEN R = -0.104528 'civil twilight (brightest)
  IF (twilight = 2) THEN R = -0.207912 'nautical twilight
  IF (twilight = 3) THEN R = -0.309017 'astronomical twilight (darkest)

  '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

  'calculate sunrise or sunset?
  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) THEN 
    p += d
  ELSE IF (m > a) THEN 
    p += b
  ENDIF

  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)

  RETURN 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


'calculate moonphase
'tried to convert from javascript source inside this page http://home.att.net/~srschmitt/script_moon_phase.html

PUBLIC FUNCTION CalcMoonPhase()
  DIM yy, mm, k1, k2, k3, jd AS Integer
  DIM ip, dp, np, rp AS Float
  DIM AG AS Float 'moon's age
  DIM DI AS Float 'moon's distance in earth radii
  DIM LA AS Float 'moon's ecliptic latitude
  DIM LO AS Float 'moon 's ecliptic longitude
  DIM Phase, Zodiac AS String
  DIM Y, D, M AS Integer
  
  y = Year(Now)
  m = Month(Now)
  d = Day(Now)
  
    'calculate the Julian Date at 12 h UT
    YY = Y - Round((12 - M)) / 10
    MM = M + 9
    IF (MM >= 12) THEN MM = MM - 12
    
    K1 = Round(365.25 * (YY + 4712))
    K2 = Round(30.6 * MM + 0.5)
    K3 = Round(Round((YY / 100) + 49) * 0.75) - 38
    
    JD = K1 + K2 + D + 59 'for dates in Julian calendar
    IF (JD > 2299160) THEN JD = JD - K3 'for Gregorian calendar

    'calculate moon's age in days

    IP = (JD - 2451550.1) / 29.530588853
    PRINT "ip = " & ip
    
    AG = IP * 29.53

    IF (AG < 1.84566) THEN 
      Phase = "NEW"
    ELSE IF (AG < 5.53699) THEN 
      Phase = "Evening crescent"
    ELSE IF (AG < 9.22831) THEN 
      Phase = "First quarter"
    ELSE IF (AG < 12.91963) THEN 
      Phase = "Waxing gibbous"
    ELSE IF (AG < 16.61096) THEN 
      Phase = "FULL"
    ELSE IF (AG < 20.30228) THEN 
      Phase = "Waning gibbous"
    ELSE IF (AG < 23.99361) THEN 
      Phase = "Last quarter"
    ELSE IF (AG < 27.68493) THEN 
      Phase = "Morning crescent"
    ELSE 
      Phase = "NEW"
    ENDIF 

PRINT "Phase = " & Phase

'giving wrong phase....;-(

END

FUNCTION normalize(v AS Float) AS Float
    v = v - Round(v)
    IF (v < 0) THEN v = v + 1
    RETURN v
END