8 กันยายน 2558

โค้ดฟังชั่นการหาเวลาพระอาทิตย์ขึ้น-ตกจริง (Visual Basic)

การเรียกใช้งานฟังก์ชั่น  SunTimes (ค.ศ., เดือน, วันที่, ลองจิจูด, ละติจูด, Time Zone, Daylight Savings)
ค่าคงที่และค่าอื่นๆที่เกี่ยวข้อง
Const pi = 3.1415926535897932384626433832795
RADEG=180/pi
DEGRAD=pi/180

Function SunTimes (strYear, strMonth, strDay, longitude, latitude, TZ, isdst)
Dim d, n, i, w, m, l, e, e1, a, xv, yv, v, xs, ys, xe, ecl, lonsun, ye, ze, ra, dec, h
Dim GMST0, UT_Sun_in_south, LHA, hour_rise, hour_set, min_rise, min_set

'calculate days since 2000 jan 1
  d = ( 367 * (strYear) - int((7 * ((strYear) + (((strMonth) + 9) / 12))) / 4) + int((275 * (strMonth)) / 9) + (strDay) - 730530)

' Orbital elements of the Sun:
  N = 0.0
  i = 0.0
  w = 282.9404 + 4.70935E-5 * d

  a = 1.000000
  e = 0.016709 - 1.151E-9 * d

  M = 356.0470 + 0.9856002585 * d
  M = rev(M)

  ecl = 23.4393 - 3.563E-7 * d
  L = w + M
  if (L < 0 OR L > 360) then
    L = rev(L)
  end if

' position of the Sun
  E1 = M + e*(180/pi) * sind(M) * ( 1.0 + e * cosd(M) )
  xv = cosd(E1) - e
  yv = sqrt(1.0 - e * e) * sind(E1)

  v = atan2d(yv, xv)
  r = sqrt(xv * xv + yv * yv)
  lonsun = v + w
  if (lonsun < 0 OR lonsun > 360) then
    lonsun = rev(lonsun)
  end if
  xs = r * cosd(lonsun)
  ys = r * sind(lonsun)
  xe = xs
  ye = ys * cosd(ecl)
  ze = ys * sind(ecl)
  RA = atan2d(ye, xe)
  Dec = atan2d(ze, (sqrt((xe * xe)+(ye * ye))))
  h=-0.833

  GMST0 = L + 180
  if (GMST0 < 0 OR GMST0 > 360) then
    GMST0 = rev(GMST0)
  end if

  UT_Sun_in_south = ( RA - GMST0 - longitude ) / 15.0
  if (UT_Sun_in_south < 0) then
    UT_Sun_in_south=UT_Sun_in_south + 24
  end if

  LHA= (sind(h) - (sind(latitude) * sind(Dec)))/(cosd(latitude) * cosd(Dec))
  if (LHA > -1 AND LHA < 1) then
    LHA=acosd(LHA)/15
  else
    SunTimes = "No sunrise,No sunset"
    exit function
  end if
  hour_rise=UT_Sun_in_south - LHA
  hour_set=UT_Sun_in_south + LHA
  min_rise=int((hour_rise-int(hour_rise)) * 60)
  min_set=int((hour_set-int(hour_set)) * 60)

  hour_rise=(int(hour_rise) + (TZ + isdst))
  hour_set=(int(hour_set) + (TZ + isdst))
  if (min_rise < 10) then
    min_rise = right("0000" & min_rise, 2)
  end if
  if (min_set < 10) then
    min_set = right("0000" & min_set, 2)
  end if

  SunTimes = hour_rise & ":" & min_rise & "," & hour_set & ":" & min_set
End Function

' Support Functions
Function sind(qqq)
  sind = sin((qqq) * DEGRAD)
End Function

Function cosd(qqq)
 cosd = cos((qqq) * DEGRAD)
End Function

Function tand(qqq)
 tand = tan((qqq) * DEGRAD)
End Function

Function atand(qqq)
 atand = (RADEG * atan(qqq))
End Function

Function asind(qqq)
 asind = (RADEG * asin(qqq))
End Function

Function acosd(qqq)
 acosd = (RADEG * acos(qqq))
End Function

Function atan2d (qqq, qqq1)
 atan2d = (RADEG * atan2(qqq, qqq1))
End Function

Function rev(qqq)
Dim x
 x = (qqq - int(qqq/360.0) * 360.0)
 if (x <= 0) then
  x = x + 360
 end if
 rev = x
End Function

Function atan2(ys,xs)
' from the Draw Arrows tip
' @ http://www.devx.com/upload/free/features/VBPJ/techtips/techtips11.pdf
' by —Jim Deutch, Syracuse, New York

Dim theta
 If xs <> 0 Then
  theta = Atn(ys / xs)
  If xs < 0 Then
   theta = theta + pi
  End If
 Else
  If ys < 0 Then
   theta = 3 * pi / 2 '90
  Else
   theta = pi / 2 '270
  End If
 End If
 atan2 = theta
End Function

function acos(x)
 acos = Atn(-X / Sqrt(-X * X + 1)) + 2 * Atn(1)
end function

function sqrt(x)
 if x > 0 then
  sqrt = Sqr(x)
 else
  sqrt = 0
 end if

end function
เขียนโดย