Sun and Moon Position in Lat/Lon from Date/time VB6

I am looking for code in visual basic 6 to get sun and moon position in latitude and longitude.
I have converted the code from this link: http://www.experts-exchange.com/Programming/Languages/CPP/Q_23501542.html?sfQueryTermInfo=1+10+posit+sun 
but it is not very acurate.

Does any one have  source code for both sun and moon position in visual basic 6?
LVL 6
iscodeAsked:
Who is Participating?
 
mdouganCommented:
OK, so, the formulas for calculating Alt and AZ from one of our earlier links here:
http://bodmas.org/kepler/sun.html

Works out for us... first, these formulas use a different "days" value... they want the number of days before 2000 which is the -877 number that you got before.  Whereas our previous calculations use the number of days since the epoch which is 2000 days or something like that.  So, I've added a couple of variables up near the top... one is  TheDate which you put the entire yyyy-MM-dd hh:mm:ss, we then use that to calculate DaysBefore

Once we have that, it is relatively easy to plug values we've already calculated into the formulas for creating LST and HA, which then helps us calculate ALT and AZ.

Down at the bottom, I've hardcoded the  observer_lat and the observer_lon but those are values that you'd probably get from a textbox on the screen, or from a GPS... if you get them from the screen, though, you have to insure that you end up with the lat/lon formatted as decimals with West longitude expressed as a negative number or South Latitude expressed as a negative number.

Note, if their result values are slightly different than our results, it is because they round the DaysBefore  to -877.04167   whereas when we calculate it is something like -877.041666666666679

If you want to get their exact numbers, then go to where we are calculating LST and uncomment the line with the rounded value and comment out the line that uses DaysBefore.
        Dim DAY0 As Double
        Dim M0 As Double
        Dim L0 As Double
        Dim EQC1 As Double
        Dim EQC2 As Double
        Dim OBLQ As Double
 
        Dim TheDate As Date
        Dim TheYear As Double
        Dim TheDay As Double
        Dim TheHour As Double
        Dim TheMinute As Double
        Dim TheSecond As Double
        Dim DaysBefore As Double
 
        Dim YRDAY As Double
        Dim FRDAY As Double
        Dim TIM As Double
 
        DAY0 = 726482  ' Day number of Epoch (origin approx. year BC 1)
        M0 = 356.6349  ' Mean anomaly at Epoch
        L0 = 279.4033  ' Mean longitude along ecliptic
        EQC1 = 1.9151  ' Equation of centre harmonic coefficient
        EQC2 = 0.02    ' Equation of centre harmonic coefficient
        OBLQ = 23.4406 ' Obliquity of ecliptic (inclination)
 
        'Position of the Sun at 11:00 UT on 1997 August 7th
        TheDate = "1997-08-07 11:00:00"
        TheYear = 1997
        TheDay = DatePart("y", "1997-08-7")
        TheHour = 11
        TheMinute = 0
        TheSecond = 0
 
        'LAT  15.0302 North
        'LON 204.9828 West
        'TheDate = "1985-08-12 01:45:00"
        'TheYear = 1985
        'TheDay = DatePart("y", "1985-08-12")
        'TheHour = 1
        'TheMinute = 45
        'TheSecond = 0
 
 
        'LAT  15.8453 North
        'LON 057.7893 West
        'TheDate = "2009-05-03 15:48:00"
        'TheYear = 2009
        'TheDay = DatePart("y", "2009-05-03")
        'TheHour = 15
        'TheMinute = 48
        'TheSecond = 0
 
 
        ' Step 1. Calculate number of days since epoch
        'Enter Year (e.g.1985), Day number (Jan 1 = 1), hours, minutes, seconds:
        YRDAY = Int((TheYear - 1) * 365.25) - DAY0                     ' Days from year difference
        FRDAY = (TheHour + ((TheMinute + (TheSecond / 60)) / 60)) / 24   ' Fraction of day
        TIM = YRDAY + TheDay + FRDAY                                     ' Time in days since epoch
 
        DaysBefore = 367 * Year(TheDate) - 7 * (Year(TheDate) + (Month(TheDate) + 9) \ 12) \ 4 + 275 * Month(TheDate) \ 9 + Microsoft.VisualBasic.Day(TheDate) - 730531.5 + Hour(TheDate) / 24
 
        'Note: These formulae only valid 1900 Mar 01 - 2100 Feb 28.
 
        'Step 2. Compute GHAA and True Longitude
        Dim MA As Double
        Dim ML As Double
        Dim GHAA As Double
        Dim TL As Double
 
        MA = M0 + 360 * TIM / 365.2596413                          ' Sun's Mean Anomaly
        While MA < 0
            MA = MA + 360
        End While
        While MA > 360
            MA = MA - 360
        End While
 
        ML = L0 + 360 * TIM / 365.2421988                              ' Mean Longitude
        While ML < 0
            ML = ML + 360
        End While
        While ML > 360
            ML = ML - 360
        End While
 
 
        GHAA = ML - 180 + 360 * FRDAY                                  ' Greenwich Hour Angle Aries
        While GHAA < 0
            GHAA = GHAA + 360
        End While
        While GHAA > 360
            GHAA = GHAA - 360
        End While
 
        ' True Longitude
        TL = ML + EQC1 * Sin(DegreesToRadians(MA)) + EQC2 * Sin(DegreesToRadians(2 * MA))
 
        'Note:   1. 'Longitude' means angle from ascending node up the ecliptic.
        '        2. The(Sun) 's argument of perigee (AP) is given by  AP = ML - MA
 
        'Step 3a. Compute Unit Vector, Geocentric Rectangular Coordinates
        Dim Xsun As Double
        Dim Ysun As Double
        Dim Zsun As Double
 
        Xsun = Cos(DegreesToRadians(TL))
        Ysun = Sin(DegreesToRadians(TL)) * Cos(DegreesToRadians(OBLQ))
        Zsun = Sin(DegreesToRadians(TL)) * Sin(DegreesToRadians(OBLQ))
 
        'Step 3b. Unit Vector, Geocentric, Greenwich Coordinates
        Dim C As Double
        Dim S As Double
        Dim Xg As Double
        Dim Yg As Double
        Dim Zg As Double
 
        C = Cos(DegreesToRadians(GHAA))
        S = Sin(DegreesToRadians(GHAA))
        Xg = (C * Xsun) + (S * Ysun)
        Yg = (-S * Xsun) + (C * Ysun)
        Zg = Zsun
 
        'Step 4a.  Astronomic Coordinates
        Dim DEC As Double
        Dim RA As Double
 
        DEC = RadiansToDegrees(ArcSin(Zsun))         ' Declination
        RA = ArcTan(DegreesToRadians(Ysun) / DegreesToRadians(Xsun))   ' Right Ascension
 
        If Xsun < 0 Then RA = RA + 180
        If RA < 0 Then RA = RA + 360
 
        'ATN is assumed to yield a value in the range +- 90 degrees.  RA is in the
        'same quadrant as true longitude, TL.
 
        'Step 4b. Geographic Coordinates
        Dim LAT As Double
        Dim LON As Double
        LAT = DEC                                                      ' Sun's Latitude +North
        LON = GHAA - RA                                                ' Sun's Longitude, +West
 
        If LAT >= 0 Then
            'text15.Text = Format(LAT, "##0.0000") & "N"
            text15.Text = LAT.ToString("##0.0000") & "N"
        Else
            'text15.Text = Format(Abs(LAT),"##0.0000") & "S"
            text15.Text = Abs(LAT).ToString("##0.0000") & "S"
        End If
        If LON >= 0 Then
            'text16.Text = Format(LON,"000.0000") & "W"
            text16.Text = LON.ToString("000.0000") & "W"
        Else
            'text16.Text = Format(Abs(LON),"000.0000") & "E"
            text16.Text = Abs(LON).ToString("000.0000") & "E"
        End If
 
        Dim LST As Double
        Dim HA As Double
        Dim observer_lat As Double
        Dim observer_lon As Double
 
        ' Birmingham UK, you might get this from user entered text fields or GPS
        observer_lat = 52.5
        observer_lon = -1.91667
 
        'LST = 280.46061837 + 360.98564736629 * -877.04167 + observer_lon
        LST = 280.46061837 + 360.98564736629 * DaysBefore + observer_lon
        While LST < 0
            LST = LST + 360
        End While
        While LST > 360
            LST = LST - 360
        End While
 
        HA = LST - RA
        ' Instructions say to do this, but then they do not do it
        'While HA < 0
        '    HA = HA + 360
        'End While
        'While HA > 360
        '    HA = HA - 360
        'End While
 
        Dim Alt As Double
 
        Alt = Sin(DegreesToRadians(DEC)) * Sin(DegreesToRadians(observer_lat)) + Cos(DegreesToRadians(DEC)) * Cos(DegreesToRadians(observer_lat)) * Cos(DegreesToRadians(HA))
        Alt = RadiansToDegrees(ArcSin(Alt))
 
        Dim X1 As Double
        Dim Y1 As Double
        Dim AZ As Double
 
        Y1 = -Cos(DegreesToRadians(DEC)) * Cos(DegreesToRadians(observer_lat)) * Sin(DegreesToRadians(HA))
        X1 = Sin(DegreesToRadians(DEC)) - Sin(DegreesToRadians(observer_lat)) * Sin(DegreesToRadians(Alt))
        AZ = ArcTan(Y1 / X1)
 
        If X1 < 0 Then
            AZ = AZ + 180
        ElseIf Y1 < 0 And X1 > 0 Then
            AZ = AZ + 360
        Else
            AZ = AZ
        End If

Open in new window

0
 
d-glitchCommented:
Check out the US Naval Observatory:

     http://www.usno.navy.mil/USNO/astronomical-applications

     http://www.usno.navy.mil/USNO/astronomical-applications/software-products/mica

Software algorithms (in C and Fortran):

     http://www.usno.navy.mil/USNO/astronomical-applications/software-products/novas

And check out the Standards Of Fundamental Astronomy (SOFA) Libraries:

     http://www.iau-sofa.rl.ac.uk/2009_0201.html
0
 
mdouganCommented:
Your question is not really specific enough for anyone to be able to comment on.  The sun and moon don't have a specific latitude/longitude.  At any given moment, in the course of the day, the sun or moon will be over a different lat/lon position on earth.  So, if you want to know what lat/lon, then what date and time are we talking about?

Longitude is a function of time, usually expressed in hours + or - from Greenwich Mean Time.

As far as latitude, the earth axis (in relation to the sun) changes during the course of the year.  At the equinoxes (approx March 21, Sept 21) the declination of the sun is 0 degrees, meaning that the equator is the closest part of earth to the sun.  As you get to the summer solstice (approx June 21st) the sun's declination is 23.4 degrees North, meaning that N23.4 latitude is the closest part of earth to the sun.  And at the winter solstice, it would be 23.4 degrees South.

You could get pretty close to the sun's position by interpolating the changes in the earth's declination using the current date and by calculating the longitude by having an accurate accounting of GMT.  That's one of the techniques that people have been using for navigation for hundreds of years. But, calculating the moon's orbit is much more complex.

If you want to know the exact lat/lon of the sun or moon for a given date/time, then you would normally consult a nautical almanac:

http://en.wikipedia.org/wiki/Nautical_almanac

The math for calculating this is pretty significant, though, and I don't know of a VB example that will produce it.  Though, I'll be curious to learn if anyone else does.

Good Luck!
0
Free Tool: IP Lookup

Get more info about an IP address or domain name, such as organization, abuse contacts and geolocation.

One of a set of tools we are providing to everyone as a way of saying thank you for being a part of the community.

 
iscodeAuthor Commented:
I want the Lat/Lon position that the sun is directly above, as shown on this site : http://www.fourmilab.ch/cgi-bin/Earth , so I can plot it on a worldmap, like I said I have formula in c++ that I have converted to visual basic 6 but its not acurate enough for me.

0
 
mdouganCommented:
Well, it is possible that they are looking up the data in a nautical almanac :-)  But, it is possible to calculate it to, using a more precise formula than the CPP example you provided.  This site has a calculation that should be pretty easy for you to convert to VB.

http://www.stargazing.net/kepler/sun.html

0
 
iscodeAuthor Commented:
pretty easy you say :-)
I have solved the first part I think:
but I am not sure how I do the second part can you give me a little help on this one?
'Position of the Sun at 11:00 UT on 1997 August 7th
'1. Find the days before J2000.0 (d) from the table
Dim test_d As Double
'11 = hour
'24 = day
test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
Text1.Text = test_d
Dim y, m, d, hours, minutes, seconde As Integer
Dim h As Double
y = 1997 'Year(Now)
m = 8 'Month(Now)
d = 7 'Day(Now)
hours = 11 'Hour(Now)
minutes = 0 'Minute(Now)
seconds = 0 'Second(Now)
h = CDbl(hours + (minutes / 60#) + (seconds / 3600#))
dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
Text2.Text = dd
'2. Find the Mean Longitude (L) of the Sun
   L = 280.461 + 0.9856474 * dd
     = -583.99284 + 720  
    (add multiples of 360 to bring in range 0 to 360)
     = 136.00716

Open in new window

0
 
mdouganCommented:
hahaha, ok, maybe not "pretty easy", but doable, I think!

By second part, you mean item 2?  The numbers that they use in their example don't exactly match up for me, but they are close, so, I think I know what they are trying to do:

'2. Find the Mean Longitude (L) of the Sun
Dim test_L as double

   ' according to the order of operations, * or / before + or -, this is the order it would calculate
   test_L = ((0.9856474 * dd) + 280.461)

   ' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
   ' but, you are trying to get L in the range of 0 to 360
   while test_L < 0
       test_L = test_L + 360
   End While
0
 
iscodeAuthor Commented:
Number 2 and 3 gives me correct answer but somehow 4 gives me wrong even though g is correct
Any idea?

'2. Find the Mean Longitude (L) of the Sun
Dim L As Double

' according to the order of operations, * or / before + or -, this is the order it would calculate
L = ((0.9856474 * dd) + 280.461)

' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
' but, you are trying to get L in the range of 0 to 360
While L < 0
 L = L + 360
Wend
Text3.Text = L

'3. Find the Mean anomaly (g) of the Sun
Dim g As Double
g = 357.528 + 0.9856003 * dd
'= -506.88453 + 720
g = g + 720
'= 213.11547
Text4.Text = g

'4. Find the ecliptic longitude (lambda) of the sun
Dim lambda As Double
lambda = L + 1.915 * Sin(g) + 0.02 * Sin(2 * g)
'= 134.97925
Text5.Text = lambda ' = 135,05047 should be 134.97925 !!!
0
 
mdouganCommented:
Yea, that is confusing.  I tried grouping the operations of the formula in different ways but that really shifted the result out of range.  So, then I tried to back into the number.  If you change the value of the 1.915 constant, you can get the correct result...

        lambda = L + 2.060165 * Sin(g) + 0.02 * Sin(2 * g)

But, that's a cludge...  it could be that they calculated incorrectly, and lambda should be 135.05047.

If you have an answer that you know is correct that you can check both formulas against, that is what I'd try.
0
 
iscodeAuthor Commented:
Yes sure is confusing, l find it strange there would be incorect math here even though it looks like it.
Anyway then number 5. is OK no confusion there:

5. Find the obliquity of the ecliptic plane (epsilon)

   epsilon = 23.439 - 0.0000004 * d
           = 23.439351

But in 6. I dont get it to match the calculation, any idea what I am doing wrong here?
6. Find the Right Ascension (alpha) and Declination (delta) of the Sun

   Y = cos(epsilon) * sin(lambda)
   X = cos(lambda)

   a = arctan(Y/X)

   If X < 0 then alpha = a + 180
   If Y < 0 and X > 0 then alpha = a + 360
   else alpha = a

   Y =  0.6489924
   X = -0.7068507

   a = -42.556485
   alpha = -42.556485 + 180 = 137.44352 (degrees)
   
   delta = arcsin(sin(epsilon)*sin(lambda))
         = 16.342193 degrees
 
Dim y, x, a, alpha, delta As Double
y = Cos(epsilon) * Sin(lambda)
x = Cos(lambda)
a = ArcTan(y / x)
 
If x < 0 Then
 alpha = a + 180
End If
 
If y < 0 And x > 0 Then
 alpha = a + 360
Else
 alpha = a
End If
 
'y = 0.6489924
Text7.Text = y
 
'x = -0.7068507
Text8.Text = x
 
'a = -42.556485
Text9.Text = a
 
'alpha = -42.556485 + 180
'alpha = a + 180
'= 137.44352 (degrees)
Text10.Text = alpha
   
delta = ArcSin(Sin(epsilon) * Sin(lambda))
'= 16.342193 degrees
Text11.Text = delta
 
Public Function ArcSin(x As Double) As Double
  'Inverse Sine
    On Error Resume Next
        ArcSin = Atn(x / Sqr(-x * x + 1))
    On Error GoTo 0
End Function
Public Function ArcTan(x As Double) As Double
  'Inverse Tangent
    On Error Resume Next
        ArcTan = Atn(x) * (180 / pi)
    On Error GoTo 0
End Function

Open in new window

0
 
mdouganCommented:
I've been expecting this next question :-)

I was trying to plug through the rest myself.

I think I have most of it figured out, but there were some things I wanted to point out... first, I see you're referencing ArcTan and ArcSin, but those aren't defined in VB, so, are you testing this in VB?  If so, what version?  Or, have you imported a library I don't have?  I switched the references to Atan and Asin, which I think are the same thing (could be wrong).

Second, you had previously declared y as a variable up when defining some date fields, for Year... did you rename that year "y" variable?  I just changed all my references in this to YY, XX, AA etc.

Third, their logic in this block doesn't make sense:

If x < 0 Then
 alpha = a + 180
End If
 
If y < 0 And x > 0 Then
 alpha = a + 360
Else
 alpha = a
End If

Because, if X < 0 then alpha will be set to a + 180, but then in the next if statement it will be overlayed with alpha = a.  So, I tried changing that to:

        If X < 0 Then
            alpha = A + 180
        ElseIf (Y < 0 And X > 0) Then
            alpha = A + 360
        Else
            alpha = A
        End If

Now, all that said, none of the variables for X, Y, A, alpha or delta come out the way you would expect.  The reason for this is that when most trig functions like cos, sin etc are referenced in literature  they are assuming that you are passing an angle in Degrees.  However, all of the corresponding math functions in VB expect the angle in Radians.  There are simple formulas (provided) for converting to and from Degrees/Radians, but the trick is figuring out where... so, after some trial and error, it looks like all of the values the formulas on this webpage are working with, are in Degrees, therefore, we have to convert them to Radians before passing to the Trig function, then convert back to Degrees after the result is returned.

That works great for everything in item 6.  So, I went back and tried to apply the same change to the lambda calculation in item 4 to see if we could get the formula to work right with their constant... but, I couldn't get it to work :-)  But, maybe if you play around with it, you will!
        'Position of the Sun at 11:00 UT on 1997 August 7th
        '1. Find the days before J2000.0 (d) from the table
        Dim test_d As Double
        '11 = hour
        '24 = day
        test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
        text1.Text = test_d
 
        Dim y, m, d, hours, minutes, seconds As Integer
        Dim h As Double
        Dim dd As Double
 
        y = 1997 'Year(Now)
        m = 8 'Month(Now)
        d = 7 'Day(Now)
        hours = 11 'Hour(Now)
        minutes = 0 'Minute(Now)
        seconds = 0 'Second(Now)
        h = CDbl(hours + (minutes / 60.0#) + (seconds / 3600.0#))
        dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
        text2.Text = dd
 
        '2. Find the Mean Longitude (L) of the Sun
        Dim L As Double
 
        ' according to the order of operations, * or / before + or -, this is the order it would calculate
        L = ((0.9856474 * dd) + 280.461)
 
        ' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
        ' but, you are trying to get L in the range of 0 to 360
        While L < 0
            L = L + 360
        End While
        text3.Text = L
 
        '3. Find the Mean anomaly (g) of the Sun
        Dim g As Double
        g = 357.528 + 0.9856003 * dd
        '= -506.88453 + 720
        g = g + 720
        '= 213.11547
        text4.Text = g
 
        '4. Find the ecliptic longitude (lambda) of the sun
        Dim lambda As Double
        'lambda = L + 1.915 * Sin(g) + 0.02 * Sin(2 * g)
        'lambda = L + 1.915 * RadiansToDegrees(Sin(DegreesToRadians(g))) + 0.02 * RadiansToDegrees(Sin(DegreesToRadians(2 * g)))
        lambda = L + 2.060165 * Sin(g) + 0.02 * Sin(2 * g)
 
        '= 134.97925
        text5.Text = lambda ' = 135.05047 should be 134.97925 !!!
 
        '5. Find the obliquity of the ecliptic plane (epsilon)
 
        Dim epsilon As Double
        epsilon = 23.439 - (0.0000004 * dd)
        '= 23.439351
        text6.Text = epsilon
 
        '6. Find the Right Ascension (alpha) and Declination (delta) of the(Sun)
        Dim XX As Double
        Dim YY As Double
        Dim AA As Double
        Dim alpha As Double
        Dim delta As Double
 
        YY = Cos(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))
        XX = Cos(DegreesToRadians(lambda))
        'YY = 0.6489924
        'XX = -0.7068507
 
        'AA = arctan(YY / XX)
        AA = RadiansToDegrees(Atan(YY / XX))
        'AA = -42.556485
 
        If XX < 0 Then
            alpha = AA + 180
        ElseIf (YY < 0 And XX > 0) Then
            alpha = AA + 360
        Else
            alpha = AA
        End If
 
        'alpha = -42.556485 + 180
        '= 137.44352 '(degrees)
        text7.Text = alpha
 
        'delta = arcsin(Sin(epsilon) * Sin(lambda))
        'delta = Asin(Sin(epsilon) * Sin(lambda))
        delta = RadiansToDegrees(Asin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
        '= 16.342193 degrees
        text8.Text = delta
 
    End Sub
 
    Private Function DegreesToRadians(ByVal X As Double) As Double
        Const pi = 3.14159265358979
        'To convert degrees to radians, multiply degrees by pi/180.
        DegreesToRadians = X * pi / 180
    End Function
 
    Private Function RadiansToDegrees(ByVal X As Double) As Double
        Const pi = 3.14159265358979
        'To convert radians to degrees, multiply radians by 180/pi.
        RadiansToDegrees = X * (180 / pi)
    End Function

Open in new window

0
 
iscodeAuthor Commented:
I was using these function, I did post them in last post:

Public Function ArcSin(x As Double) As Double
  'Inverse Sine
    On Error Resume Next
        ArcSin = Atn(x / Sqr(-x * x + 1))
    On Error GoTo 0
End Function
Public Function ArcTan(x As Double) As Double
  'Inverse Tangent
    On Error Resume Next
        ArcTan = Atn(x) * (180 / pi)
    On Error GoTo 0
End Function

>>
Second, you had previously declared y as a variable up when defining some date fields, for Year... did you rename that year "y" variable?  I just changed all my references in this to YY, XX, AA etc.
<<
I changed that to as they use in exsample its confusig enough :-)

Ahh angle in Degrees!!!
you are smart, let me test your code and see how far I can go
0
 
mdouganCommented:
maybe not so smart, just persistent.  I've done a lot of charting where it is necessary to work with trig functions in VB as well as several navigation problems which also use the same functions, so, I've have to wrestle with radians/degrees before.

First time working out the celestial positions though, this is interesting.

I missed your ArcTan and ArcSin functions, sorry.  I googled some link which suggested ArcTan and Atan and ArcSin and Asin were equivalent.  Your functions reference Atn, which I might have thought was the same as Atan... hehe, pretty confusing.  But, it appears using Atan and Asin work in place of ArcTan and ArcSin.
0
 
iscodeAuthor Commented:
yes very interesting subject.
It seams to work now,
Final result

   Right ascension is usually given in hours of time, and both
   figures need to be rounded to a sensible number of decimal
   places.

   alpha =   9.163 hrs      or   9h 09m 46s
   delta = +16.34 degrees   or +16d 20' 32"

   The Interactive Computer Ephemeris gives

   alpha =   9h 09m 45.347s and
   delta = +16d 20' 30.89"
How do they get alpha =  137.44352 = 9.163 hrs?
and is  alpha the longitude? and delta the latitude?

I found another link http://bodmas.org/kepler/sun.html
where they take this to steps 7. and 8. to find the Aziumuth (Z) and altitude (A)
I think we must take this all the way as we are started .-)
0
 
iscodeAuthor Commented:
This is the working code so far
'Position of the Sun at 11:00 UT on 1997 August 7th
        '1. Find the days before J2000.0 (d) from the table
        Dim test_d As Double
        '11 = hour
        '24 = day
        test_d = 11 / 24 + 212 + 7 - 1096.5 '= -877.04167
        Text1.Text = test_d
 
        Dim y, m, d, hours, minutes, seconds As Integer
        Dim h As Double
        Dim dd As Double
 
        y = 1997 'Year(Now)
        m = 8 'Month(Now)
        d = 7 'Day(Now)
        hours = 11 'Hour(Now)
        minutes = 0 'Minute(Now)
        seconds = 0 'Second(Now)
        
'        y = Year(Now)
'        m = Month(Now)
'        d = Day(Now)
'        hours = Hour(Now)
'        minutes = Minute(Now)
'        seconds = Second(Now)
        
        h = CDbl(hours + (minutes / 60#) + (seconds / 3600#))
        dd = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
        Text2.Text = dd
 
        '2. Find the Mean Longitude (L) of the Sun
        Dim L As Double
 
        ' according to the order of operations, * or / before + or -, this is the order it would calculate
        L = ((0.9856474 * dd) + 280.461)
 
        ' assuming L is a negative number... if not, I'm not sure if you'd subtract if L > 360
        ' but, you are trying to get L in the range of 0 to 360
        While L < 0
            L = L + 360
        Wend
        Text3.Text = L
 
        '3. Find the Mean anomaly (g) of the Sun
        Dim g As Double
        g = 357.528 + 0.9856003 * dd
        '= -506.88453 + 720
        g = g + 720
        '= 213.11547
        Text4.Text = g
 
        '4. Find the ecliptic longitude (lambda) of the sun
        Dim lambda As Double
        'lambda = L + 1.915 * Sin(g) + 0.02 * Sin(2 * g)
        'lambda = L + 1.915 * RadiansToDegrees(Sin(DegreesToRadians(g))) + 0.02 * RadiansToDegrees(Sin(DegreesToRadians(2 * g)))
        lambda = L + 2.060165 * Sin(g) + 0.02 * Sin(2 * g)
 
        '= 134.97925
        Text5.Text = lambda ' = 135.05047 should be 134.97925 !!!
 
        '5. Find the obliquity of the ecliptic plane (epsilon)
 
        Dim epsilon As Double
        epsilon = 23.439 - (0.0000004 * dd)
        '= 23.439351
        Text6.Text = epsilon
 
        '6. Find the Right Ascension (alpha) and Declination (delta) of the(Sun)
        Dim XX As Double
        Dim YY As Double
        Dim AA As Double
        Dim alpha As Double
        Dim delta As Double
 
        YY = Cos(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))
        XX = Cos(DegreesToRadians(lambda))
        'YY = 0.6489924
        Text7.Text = YY
        'XX = -0.7068507
        Text8.Text = XX
        AA = ArcTan(YY / XX)
        'AA = RadiansToDegrees(Atan(YY / XX))
        'AA = RadiansToDegrees(ArcTan(YY / XX))
        'AA = -42.556485
        Text9.Text = AA
        If XX < 0 Then
            alpha = AA + 180
        ElseIf (YY < 0 And XX > 0) Then
            alpha = AA + 360
        Else
            alpha = AA
        End If
 
        'alpha = -42.556485 + 180
        '= 137.44352 '(degrees)
        Text10.Text = alpha
 
        'delta = arcsin(Sin(epsilon) * Sin(lambda))
        'delta = Asin(Sin(epsilon) * Sin(lambda))
        'delta = RadiansToDegrees(Asin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
        delta = RadiansToDegrees(ArcSin(Sin(DegreesToRadians(epsilon)) * Sin(DegreesToRadians(lambda))))
        '= 16.342193 degrees
        Text11.Text = delta

Open in new window

0
 
mdouganCommented:
To convert degrees to hours, divide by 15.  360/24 = 15  (so, there are 15 degrees per hour)

Plugging that into the alpha above returns 9.1629

I'm game, let's take it all the way... but I'll follow you ;-)
0
 
iscodeAuthor Commented:
haha OK I am in
I think this solves the lambda problem
lambda = L + 1.915 * Sin(DegreesToRadians(g)) + 0.02 * Sin(DegreesToRadians(2 * g))
0
 
mdouganCommented:
Nice!  I thought I tried that first (without converting back to Degrees), but I might have put the second DegreesToRadians just around the g rather than the 2*g... too many combinations ;-)
0
 
iscodeAuthor Commented:
So is that correct that
alpha = RA = longitude, and delta= latitude?
0
 
mdouganCommented:
Yes, RA stands for Right Ascension which in our example was alpha and is longitude, however, in what units?  Not degrees, as that doesn't compute correctly.  Possibly radians... but, I backed into the number by taking their answer then solving all the other known stuff and the value they want for "long" is:

-1.9166702079166957

Which looks suspiciously like that constant, we had a while back, in step 4, of   1.915

delta should be the same as latitude, but then, when I was looking at the formulas for calculating altitude and azimuth, they were multiplying by the sin of both delta and lat... which made me stop and think for a minute... then it hit me...

Wait, we'd already calculated at what geometric position (GP) the sun would be directly under on the face of the earth at that given time... so, the altitude, if expressed in degrees, has to be exactly 90, and the azimuth, which would be the bearing away from that point would have to be 0 or 360 or perhaps you can't calculate it because you're standing on the GP.  So, I'm thinking this is why the previous example stopped at step 6.

But, to be sure, you should calc out the other formulas... this was as far as I got:
        '7. Find the Aziumuth (Z) and altitude (A) of the Sun, we first
        'find the local sidereal time, then the hour angle of the Sun.
        'Long is the longitude (west negative) of the place, and lat is
        'the latitude. LST and RA are taken in degrees. Note that the
        'calculation of LST needs about 12 places of precision.
 
        ' I tried with two different variable types Decimal and Double
        ' I ended up with the same result, however, I think it might
        ' need to be this decimal type, and the other Double values
        ' typecast to decimals to get greater accuracy, as the value
        ' rounded slightly
        Dim LST_DEC As Decimal
        LST_DEC = 280.46061837D + (360.98564736629D * dd) + (-1.9166702079166957)
        While LST_DEC < 0
            LST_DEC += 360
        End While
        text12.Text = LST_DEC.ToString("###,##0.000000000000")
 
        Dim LST As Double
        LST = CDbl(280.46061837D + (360.98564736629D * dd) + (-1.9166702079166957))
        While LST < 0
            LST += 360
        End While
        text13.Text = LST.ToString("###,##0.000000000000")
        'LST = -316320.911064 degs
        '= 119.088936 degs
 
        Dim ha As Double
        ha = LST - alpha
        While ha < 0
            ha += 360
        End While
        text14.Text = LST.ToString("###,##0.000000000000")
        'ha = 119.088936 - 137.44352degs
        '= -18.354584 degs
        '= 341.64662711177084  (they forgot to 0-360 the value)

Open in new window

0
 
iscodeAuthor Commented:
Oh my this is even more confusing,
I dont have this decimal type in vb6 so I cant try  it, and LST.ToString neither.

first I thought you would do it like this
Dim LST As Double
       LST = 280.46061837 + 360.98564736629 * dd + DegreesToRadians(alpha) '= -316316,594348
      that is close to -316320.911064 :) but  no?


0
 
iscodeAuthor Commented:
Here is something I found
'Right Ascension- a (alpha-- small Greek a, or RA) is like longitude,
'and "meridian lines" of constant RA connect the north and south poles of the sky
'and are everywhere perpendicular to lines of constant declination.

'RA and longitude differ in two ways. First, while the zero meridian on Earth is the one passing Greenwich, England,
'the celestial meridian on which a = 0 is the one passing one of the two points
'where the ecliptic crosses the celestial equator, in this case,
'the "first point in Aries" or "vernal equinox." And second, for historical reasons,
'RA is measured not in degrees but in hours, each hour equal to 15 degrees
'(that is, 24 hours = 360 degrees). Hours are subdivided into minutes and seconds,
'but these are not the same as minutes and seconds into which degrees are divided.
0
 
mdouganCommented:
Hi,

OK, I've got my sample in .NET, I've got VB6 available, but will probably just keep going with what I have.

Rather than Decimal, you should try declaring LST to be a float, but, you might want to look up the limits of the float data type.  In the problem comments, they were saying that you might need as many as 12 decimal positions to accurately calculate LST, and I think that float might handle that.  So, try something like this and let me know how it goes:

        Dim LST As float
        Dim dd_float as float
     
        dd_float = dd

        LST = 280.46061837 + (360.98564736629 * dd_float) + (-1.9166702079166957)

        While LST < 0.0
            LST = LST + 360.0
        End While

        text12.Text = FORMAT(LST,"###,##0.000000000000")

The reason I'm bothering to convert dd to a float is that VB decides what data type to cast the result of a calculation based on the lowest common denominator of the arguments of the calculation.  So, if all your values are floats except for one double, then it will make the result of each intermediate calculation a double... and you'll suffer a loss of precision.  By explicitly converting dd to a float, then it forces VB to make all intermediate results floats and you should keep your precision (... should)

If that doesn't work, then we might have to create a bunch of float variables and assign all of the constants to the float variables to ensure VB uses all floats.

I've also given you the VB6 equivalent to ToString with the Format statement.

On the stuff you found.  Yes, I'd forgot that in navigation, their zero meridan is called Aries, and so often, rather than longitude, it is called Aries Angle, which in our example is called Right Ascension.  I don't know if we're supposed to correct for that or not, or if it was included in our calcs... It might be what that 1.915 number is all about.

Everything else matches what we've already talked about.

(the minutes I'm talking about below are the degrees minutes, not the hours minutes)
One note, at the equator, the distance between a degree of longitude is about 60 nautical miles or about one nautical mile per minute of longitude.  However, even though there are also 60 minutes of longitude at, say, N45 degrees latitude, the distance from one minute to the next is only maybe a half nautical mile, making one degree only 30 nautical miles apart.  At the pole, all the longitude lines converge.

Probably not important for this problem, but just so you know that though longitude lines are evenly spaced around the globe, the actual distance between them has to be calculated based on the latitude.

0
 
iscodeAuthor Commented:
float is not in vb6, here are the types:

" Double
The Double data type is a 64-bit floating point number used when high accuracy is needed. These variables can range from -1.79769313486232e308 to -4.94065645841247e-324 for negative values and from 4.94065645841247e-324 to 1.79769313486232e308 for positive values.

" Long
The Long data type is a 32-bit number which can range from -2,147,483,648 to 2,147,483,647. Long variables can only contain non-fractional integer values. I myself use Long variables over Integers for increased performance. Most Win32 functions use this data type for this reason.

" Single
The Single data type is a 32-bit number ranging from -3.402823e38 to -1.401298e-45 for negative values and from 1.401298e-45 to 3.402823e38 for positive values. When you need fractional numbers within this range, this is the data type to use.

0
 
iscodeAuthor Commented:
I have found what is making us crazy here in the sample they take
" Position of the Sun at 11:00 UT on 1997 August 7th in
Birmingham UK (52.5 latitude, -1.91667 longitude. Longitude is
always taken as negative to the West)"
So the -1.91667 is longitude in Birmingham :-)

I was wondering about that number and why we need it but I saw this from another c++ sample I have
to calculate sunrise and set and we are not doing that right :)
 
0
 
mdouganCommented:
I guess Double is going to be as high as you're going to get in terms of precision, so, stick with that, but perhaps either move all constants into a double variable before the calculation, or cast them as doubles in the calculation to insure that all intermediate variables are typed as double as well.

Yea, this second example added a new twist that I didn't read carefully enough... I thought it was the same exercise as the first one.  In the first one, we were calculating the lat/lon of the geographical position (GP) of the sun for a given date/time (ie, where the sun is directly over head at that date/time).

This second exercise is saying "where in the sky is the sun, as seen from Birmingham England on this date/time".  Assuming that Birmingham isn't in the exact GP of the sun at that moment, then the sun will appear lower than 90 degrees in the sky.  And, rather than being directly South (if we were standing in the Northern Hemisphere and we were at the same longitude as the GP of the sun), or North (if we were standing in the Southern Hemisphere), the sun will have some bearing from our position in Birmingham.

So, that is what this extra calculation does.  It tells you the angle of the sun above the horizon (the altitude), and the aziumuth (bearing) from your position.  So, the LONG and the LAT in that equation are supposed to be the lat/long of Birmingham, not the alpha and delta we calculated.  

This is pretty much the same process that they go through when taking a "sun sight" with a sextant in navigation.  The general idea is that you take the lat/long of your assumed position (where you guess that you are), and you run this same exact calculation to see what the altitude and aziumuth of the sun are for that position.  Then, you use a sextant to measure the angle of the sun off the horizon (altitude), and you take a compass bearing (aziumuth) to see what the actual altitude and aziumuth are.  Then, there are a series of other equations that will tell you the distance and bearing your actual position is away from the assumed position you started with.  Pretty cool.

But, they make this pretty cool calculator called a Celisticomp that will do all this stuff for you....
0
 
iscodeAuthor Commented:
I am trying to visual the picture,
so in the first sample we get alpha and alpha/15 =  RA and delta = DEC
alpha was 137.44352 degrees ´longitude
delta was 16.342193 degrees 'latitude
how do you convert that to get location lat/long where the sun is directly over earth on map
as  right now it is 15.2 N and 32.1 W?

I wanted to use real date and time and test if this is correct now
0
 
mdouganCommented:
I think you only need to calculate RA if you are going to calculate the Altitude and Aziumuth.  I thought that alpha was the longitude and delta was the latitude of the location where the sun is directly over the earth on the map.

To do it for right now, you would input the date and time that it is right now, however the time has to be UTC - universal time.  For me, I think I'm 4 hours behind UTC, so, the code at the top of your routine would be this...

        y = Year(Now)
        m = Month(Now)
        d = Microsoft.VisualBasic.Day(Now)
        hours = (Hour(Now) + 4)
        minutes = Minute(Now)
        seconds = Second(Now)

However, when I run it, I get different results than the results on this page:

http://www.fourmilab.ch/cgi-bin/Earth

Is this what you're checking it against?
0
 
iscodeAuthor Commented:
yes or http://www.usno.navy.mil/USNO/astronomical-applications/data-services/earthview
I dont get corect result either what is the problem?
0
 
mdouganCommented:
Well, I checked with the fourmilab site, and our calculation for delta as the latitude is correct.  It is just that ours is expressed in degrees and tenths of a degree whereas they display degrees and minutes.  If you divide their minute value by 60 and add it to their degree value, we get the same number.

But, alpha isn't giving us our longitude.  It is being used in calculating Local Hour Angle, but I'm at a loss to remember how that helps us get our longitude.

Whenever you want to calculate the longitude of the sun, you should simply have to find out what the UTC time is, then calculate how much time is before or after 12:00 noon in hours and tenths of hours and then multiply by 15.  If the time is before noon, then the longitude is East, else the longitude is West, unless you are over 180 either way, then it switches around.

I just did the calculation for the UTC as appears on that site, and it is close, within about a half a degree... not sure why it wouldn't be more exact.
0
 
mdouganCommented:
By the way, I found this cool Java Applet that does all the calcs we're trying to do... now, if we could just get their code :-)

http://www.jgiesen.de/SME/
0
 
iscodeAuthor Commented:
I am very bad with this conversation degree and hours minute do you have simple conversation for that?

Yes this is just what we need :-) and I wonder how you do the shadow on the map.
0
 
mdouganCommented:
Conversion actually works pretty much the same as Hours-Minutes-Seconds of time.

Each degree of arc is subdivided into 60 minutes of arc, so, a minute is 1/60th of a degree.  
Each minute of arc is subdivided into 60 seconds of arc, so, a second is 1/60th of a minute

If you have a longitude of this for example:

 74° 2'29.15"W

That reads as 74 degrees
2 minutes
29.15 seconds of arc

So, you convert the seconds to minutes by  29.15 / 60 = 0.48583 minutes
Then add that to the minutes = 2.48583
Then convert the minutes to degrees by   2.48583 / 60 = 0.04143 degrees
Then add that to the degrees to get the complete degrees =  74.04143

Looking at the Java applet, the value they are giving as GHA which stands for Greenwich Hour Angle appears to be the Longitude of the sun.  I remember calculating GHA in navigation, so, looked it up and found this page:

http://www.irbs.com/lists/navigation/9712/0002.html

I'm working through the worksheet to see how that matches what we've done in our code... so far, I'm mystified by the formula for calculating the Greenwich Hour Angle of Aries GHAA, which you have to do first... sigh.  I got an A in celestial navigation in school... but that was some time ago.
0
 
mdouganCommented:
projecting the shadow on the map is a whole other story (and probably question!).  The typical map you see is called a Mercator Projection, where they flatten out the earth.  Was it you that asked that whole Mercator thread a while back?  :-D
0
 
iscodeAuthor Commented:
thank you for the conversation lesson I think I have to set that up acording to the website so comparation is right

so we need to do this hour angle to get right lat/on I dont know where I was if not for your clever mind and A is not bad and even if some time ago its like bike you never forget how to :-)
I found this website about this subject too maybe it helps us. http://www.amsat.org/amsat/articles/g3ruh/112.html

Step 2. Compute GHAA and True Longitude
---------------------------------------
MA   = M0   +   360 * TIM / 365.2596413         Sun's Mean Anomaly
ML   = L0   +   360 * TIM / 365.2421988         Mean Longitude
GHAA = ML-180 + 360 * FRDAY                     Greenwich Hour Angle Aries

yes we did the mercator together did you forget already? haha

shadowing the map is a must after we solve the position dont you think? ofcourse in a new question
0
 
iscodeAuthor Commented:
Calculating Local Hour Angle is really tricky business I have really hard time figure this out. I found some interesting formula but they use numbers I dont know how to use in visual basic 6.
http://www.esrl.noaa.gov/gmd/ozwv/dobson/papers/report6/appf.html

POSITION OF THE SUN TO 1' OF ARC PRECISION
i Definitions

Julian day number = number of days elapsed since 12h
Julian date = Julian day number followed by the fraction of a day elapsed since the preceding noon.

ii Variables
          D = number of days from 1900.0 (Julian date 2415020.0)
          T = number of Julian centuries from 1900.0 = D/36525
          L = mean longitude of sun in degrees
          E = equation of time (in seconds of time)
    epsilon = obliquity of the ecliptic
      alpha = right ascension of sun (apparent) in degrees
      delta = declination of sun (apparent) in degrees
        GHA = Greenwich hour angle of sun (in degrees)
         UT = Universal time which is equivalent to GMT in this context (in degrees)

iii Formulae

          D = Julian date - 2415020
          T = D/36525
          L = 279.697 + 36000.769T
          E = - (93.0 + 14.23T - 0.0144T2)sin L - (432.5 - 3.71T - 0.2063T2)cos L
              + (596.9 - 0.81T - 0.0096T2)sin2L - (1.4 + 0.28T)cos 2L
              + (3.8 + 0.60T)sin3L + (19.5 - 0.21T - 0.0103T2)cos 3L - (12.8 - 0.03T)sin 4L
  tan epsilon = 0.43382 - 0.00027T
      alpha = L - E/240
  tan delta = tan(epsilon) sin(alpha)
        GHA = UT + E/240 + 180

iv Worked Example
          Find GHA and declination of the sun on 1976 August 8 at 6h UT
          D = 27978.75
          T = 0.7660164
          L = 136.877ý
          E = -335.4 in seconds of time (= -1.397ý)
  tan epsilon = 0.433613
      alpha = 138.274ý = 138ý 16.4'
      delta = 16.098ý = 16ý 5.9'
        GHA = 268.603ý = 268ý 36.2'


Have you seen this before or can you figure how we can use this numbers?
0
 
mdouganCommented:
I looked over the link you gave in the previous example, and it looks really promissing.  I just haven't had the time to try to plug it into our code sample yet.  But, it looks like it will give us all that we need.

Not sure what you are referring to as difficult numbers in the last example... I see some weird characters such as  ý  but I think that is maybe just a problem with the website or something.

Hopefully, I'll have time to work on it tomorrow!
0
 
iscodeAuthor Commented:
Thats good news if we can use that and mdougan its no rush take your time for this.

About that pages yes it looks like some strange char maybe its the degree ° sign but there are some
math formula setup like  0.0144T2, 0.60T, sin2L I am not familar with.
0
 
mdouganCommented:
I saw the sin2L and my guess was that meant  SIN(2 * L) or something like that.  

But, for the ones that mention 0.60T or the like, I've seen that notation before, but I'm not sure what it means either!  I know that when certain values, such as date/times are translated to XML, they automatically get a T at the end of the value.

I think that the link you gave will give us the longitude and that is what we really need.
0
 
mdouganCommented:
Alright!  I think that link of yours worked!  This is how I converted the formulas from this link:

http://www.amsat.org/amsat/articles/g3ruh/112.html

It should mostly be compatible with VB6... you might have problems with the DatePart function, not sure if VB6 has that same function... also, be sure to use the Format statements instead of the ToString statements.

I tried it with their test code, then with the original date we were using to test with, then with today's date, and matched it to one of those other sites that calculates the Sun's position.  Seems to work!


        Dim DAY0 As Double
        Dim M0 As Double
        Dim L0 As Double
        Dim EQC1 As Double
        Dim EQC2 As Double
        Dim OBLQ As Double
 
        Dim TheYear As Double
        Dim TheDay As Double
        Dim TheHour As Double
        Dim TheMinute As Double
        Dim TheSecond As Double
 
        Dim YRDAY As Double
        Dim FRDAY As Double
        Dim TIM As Double
 
        DAY0 = 726482  ' Day number of Epoch (origin approx. year BC 1)
        M0 = 356.6349  ' Mean anomaly at Epoch
        L0 = 279.4033  ' Mean longitude along ecliptic
        EQC1 = 1.9151  ' Equation of centre harmonic coefficient
        EQC2 = 0.02    ' Equation of centre harmonic coefficient
        OBLQ = 23.4406 ' Obliquity of ecliptic (inclination)
 
        'Position of the Sun at 11:00 UT on 1997 August 7th
        'TheYear = 1997
        'TheDay = DatePart("y", "1997-08-7")
        'TheHour = 11
        'TheMinute = 0
        'TheSecond = 0
 
        'LAT  15.0302 North
        'LON 204.9828 West
        'TheYear = 1985
        'TheDay = DatePart("y", "1985-08-12")
        'TheHour = 1
        'TheMinute = 45
        'TheSecond = 0
 
 
        'LAT  15.8453 North
        'LON 057.7893 West
        TheYear = 2009
        TheDay = DatePart("y", "2009-05-03")
        TheHour = 15
        TheMinute = 48
        TheSecond = 0
 
 
        ' Step 1. Calculate number of days since epoch
        'Enter Year (e.g.1985), Day number (Jan 1 = 1), hours, minutes, seconds:
        YRDAY = Int((TheYear - 1) * 365.25) - DAY0                     ' Days from year difference
        FRDAY = (TheHour + ((TheMinute + (TheSecond / 60)) / 60)) / 24   ' Fraction of day
        TIM = YRDAY + TheDay + FRDAY                                     ' Time in days since epoch
 
        'Note: These formulae only valid 1900 Mar 01 - 2100 Feb 28.
 
        'Step 2. Compute GHAA and True Longitude
        Dim MA As Double
        Dim ML As Double
        Dim GHAA As Double
        Dim TL As Double
 
        MA = M0 + 360 * TIM / 365.2596413                          ' Sun's Mean Anomaly
        While MA < 0
            MA = MA + 360
        End While
        While MA > 360
            MA = MA - 360
        End While
 
        ML = L0 + 360 * TIM / 365.2421988                              ' Mean Longitude
        While ML < 0
            ML = ML + 360
        End While
        While ML > 360
            ML = ML - 360
        End While
 
 
        GHAA = ML - 180 + 360 * FRDAY                                  ' Greenwich Hour Angle Aries
        While GHAA < 0
            GHAA = GHAA + 360
        End While
        While GHAA > 360
            GHAA = GHAA - 360
        End While
 
        ' True Longitude
        TL = ML + EQC1 * Sin(DegreesToRadians(MA)) + EQC2 * Sin(DegreesToRadians(2 * MA))
 
        'Note:   1. 'Longitude' means angle from ascending node up the ecliptic.
        '        2. The(Sun) 's argument of perigee (AP) is given by  AP = ML - MA
 
        'Step 3a. Compute Unit Vector, Geocentric Rectangular Coordinates
        Dim Xsun As Double
        Dim Ysun As Double
        Dim Zsun As Double
 
        Xsun = Cos(DegreesToRadians(TL))
        Ysun = Sin(DegreesToRadians(TL)) * Cos(DegreesToRadians(OBLQ))
        Zsun = Sin(DegreesToRadians(TL)) * Sin(DegreesToRadians(OBLQ))
 
        'Step 3b. Unit Vector, Geocentric, Greenwich Coordinates
        Dim C As Double
        Dim S As Double
        Dim Xg As Double
        Dim Yg As Double
        Dim Zg As Double
 
        C = Cos(DegreesToRadians(GHAA))
        S = Sin(DegreesToRadians(GHAA))
        Xg = (C * Xsun) + (S * Ysun)
        Yg = (-S * Xsun) + (C * Ysun)
        Zg = Zsun
 
        'Step 4a.  Astronomic Coordinates
        Dim DEC As Double
        Dim RA As Double
 
        DEC = RadiansToDegrees(ArcSin(Zsun))         ' Declination
        RA = ArcTan(DegreesToRadians(Ysun) / DegreesToRadians(Xsun))   ' Right Ascension
 
        If Xsun < 0 Then RA = RA + 180
        If RA < 0 Then RA = RA + 360
 
        'ATN is assumed to yield a value in the range +- 90 degrees.  RA is in the
        'same quadrant as true longitude, TL.
 
        'Step 4b. Geographic Coordinates
        Dim LAT As Double
        Dim LON As Double
        LAT = DEC                                                      ' Sun's Latitude +North
        LON = GHAA - RA                                                ' Sun's Longitude, +West
 
        If LAT >= 0 Then
            'text15.Text = Format(LAT, "##0.0000") & "N"
            text15.Text = LAT.ToString("##0.0000") & "N"
        Else
            'text15.Text = Format(Abs(LAT),"##0.0000") & "S"
            text15.Text = Abs(LAT).ToString("##0.0000") & "S"
        End If
        If LON >= 0 Then
            'text16.Text = Format(LON,"000.0000") & "W"
            text16.Text = LON.ToString("000.0000") & "W"
        Else
            'text16.Text = Format(Abs(LON),"000.0000") & "E"
            text16.Text = Abs(LON).ToString("000.0000") & "E"
        End If

Open in new window

0
 
iscodeAuthor Commented:
I dont get it to work with live date, datepart is working though gives me 123 days
RA is always 0
GHAA is not matching the java applet on that website site
mine is now 139,88766 and on the website is 99°03'2
the dates and time are:
        TheYear = Year(Now) '2009
        TheDay = DatePart("y", Date) 'DatePart("y", "2009-05-03")
        TheHour = Hour(Now) '15
        TheMinute = Minute(Now) '48
        TheSecond = Second(Now) '0

 And yours worked, how is your Arctan function?
0
 
iscodeAuthor Commented:
ok I see my problem I get FRDAY = 7,29166666666667E-02 when theyr  test is FRDAY = 0.072917
any idea how I fix that?
0
 
mdouganCommented:
I think your FRDAY value is correct... if I remember correctly, the E-02 at the end of your number means "move the decimal position two positions to the left.  So, your value would be:

.07916666666

If you were to use the format statement:

msgbox(FORMAT(FRDAY, "##0.000000"))  

Then I bet you see the same value.
0
 
mdouganCommented:
123 days sounds about right... what that TheDay variable is supposed to hold is what number of day of the year is it... so, May 3rd is probably about the 123rd day of the year.

What I did is step through their test data and verify each value as I calculate it.  I didn't try using Now and Date, but I ran the calculator at this website:

http://www.usno.navy.mil/USNO/astronomical-applications/data-services/earthview

and they matched.

I'm using your ArcTan and ArcSin functions.
0
 
iscodeAuthor Commented:
All test data is correct except RA = 0  before it goes throug the If statments and there it gets 180 so there must be something wrong in that Arctan fuction because Ysun and Xsun is correct

I am using it like this

Public Function ArcTan(X As Double) As Double
  'Inverse Tangent
    On Error Resume Next
        ArcTan = Atn(X) * (180 / pi)
    On Error GoTo 0
End Function


0
 
mdouganCommented:
Have you defined pi?  I didn't think that VB6 had Pi defined.

        Const pi = 3.14159265358979
Are you calling the DegreesToRadians function for the parameters?

        RA = ArcTan(DegreesToRadians(Ysun) / DegreesToRadians(Xsun))   ' Right Ascension

What exact date/time are you trying to run through, and what "should" the LAT/LONG be?
0
 
iscodeAuthor Commented:
I was missing  the pi ofcourse and it didnt give error on it, thanks so much.

I am using theyr test data from the formula and its correct now :-)
I will try it now with date today and see how it goes
0
 
mdouganCommented:
Cool!  I'm done for the night... I'll see how you're doing tomorrow!

Cheers!
0
 
iscodeAuthor Commented:
I think  this is working perfectly correct enough. Great formula we have there.  I am done for the night too..
sun is already down here so time to look to the moon :-)

we see how tomorrow goes
Cheers
0
 
iscodeAuthor Commented:
Hi
I wonder if we can now use this formula to find moon position lat/lon? I think part of it can help us
because on this website http://stjarnhimlen.se/comp/tutorial.html#8 
at number 8. The Moon's position with higher accuracy. Perturbations
they are using this fundamental arguments:

    Sun's  mean longitude:        Ls     (already computed)
    Moon's mean longitude:        Lm  =  N + w + M (for the Moon)
    Sun's  mean anomaly:          Ms     (already computed)
    Moon's mean anomaly:          Mm     (already computed)
    Moon's mean elongation:       D   =  Lm - Ls
    Moon's argument of latitude:  F   =  Lm - N


But to finish of the sun position I think we need to have at least altitude, azimuth and distance in the formula. We had started it but now I am trying to combine it with this new formula.

As you see on the graph above the worldmap from http://www.jgiesen.de/SME/ they show this for each 24 hour,

So first with the distance, on our previous sample http://www.stargazing.net/kepler/sun.html 

'   Find the Earth - Sun distance
'
'r = 1.00014 - 0.01671 * Cos(g) - 0.00014 * Cos(2 * g)
so... g = mean anomaly of the Sun and that is our MA
r = 1.00014 - 0.01671 * Cos(DegreesToRadians(MA)) - 0.00014 * Cos(DegreesToRadians(2 * MA))

and that seams to be correct with that sample data though I have not found place to compare to real date yet on web.

So to altitude and azimuth,
the altitude (elevation) is the angular distance measured above the horizon,
the azimuth is the angular distance measured along the horizon

We first need to find the local sidereal time LST and then the hour angle of the Sun ha

I am using the test data and it works when I put test_d = -877.04167
but I am not sure what I use now because it is not FRDAY and for alpha we would use RA?

        Dim LST As Double
        LST = CDbl(280.46061837 + (360.98564736629 * test_d) + (-1.9166702079167))
        While LST < 0
            LST = LST + 360
        Wend
        Text2.Text = Format(LST, "###,##0.000000000000")
        'LST = -316320.911064 degs
        '= 119.088936 degs
       
        Dim ha As Double
        ha = LST - RA
        While ha < 0
            ha = ha + 360
        Wend
        Text3.Text = ha
        'ha = 119.088936 - 137.44352degs
        '= -18.354584 degs
        '= 341.64662711177084  (they forgot to 0-360 the value)
0
 
mdouganCommented:
One thing to keep in mind about computing the Altitude and Azimuth of the Sun is that you are computing those values in relation to some other location where you are viewing the sun from.  

So, what we have been calculating so far is where is the sun directly overhead at this precise day/time.  At that place, at that time, the Altitude is 90 degrees and the azimuth is either zero or can't be calculated, because the sun is directly over head.

Now, if you were 1000 miles away from the Sun's GP at that exact moment, then, the sun would appear lower in the sky to you (that would be the altitude), and also , depending on your longitude, it would either appear more to the East or West (or North or South) from where you are.  So, that would be the azimuth.

So, if you are going to calculate the altitude and azimuth, then you have to decide where you are calculating it from.  One of our sample problems gave us a location to use in the example, but how would that work in your application?

I agree with you, that many of the same calculations will probably help us on the moon calculations too.
0
 
iscodeAuthor Commented:
Ofcourse I had that in mind that we are looking at sun from observer view point and he enter his location, select from list name of city or click on map to get lat/lon but in our formula sample we are using this as test Birmingham UK 52.5 latitude, -1.91667 longitude
0
 
iscodeAuthor Commented:
I forgot to say that using gps to have the lat and lon input to the program is another way.
This website is good source for Converting RA and DEC to ALT and AZ
http://www.stargazing.net/kepler/altaz.html 
0
 
iscodeAuthor Commented:
This looks very solid and is working correct on the test data I will try it on real date and let you know how it goes.
0
 
iscodeAuthor Commented:
There are some different but I think its only because of this rounding so it can be tweaked just need more time to format this better or our data is just more acurate :-)
0
 
mdouganCommented:
Yea, I read one comment that you should carry all the decimal positions all the way through the calculations for the greatest accuracy, then, only round before display.  I think that is a great idea.

But, the makers of the samples don't always follow that advice!
0
 
iscodeAuthor Commented:
I was wondering about the calculation for LST we dont take into account in the formula for Time Zone Daylight Time, I have found a good sample how its done here http://63.236.73.107/ShowCode.asp?ID=8371&NoBox=True
and then here is a good sample how to check agains different Time Zone, Daylight Time and Location

0
 
mdouganCommented:
When we've specified a time, so far, we've been specifying Universal Time (UT), this is what all the other websites we've been using to check our results do too.  UT isn't affected by daylight savings time, and because it is basically the time at the zero meridian, it isn't affected by timezones either.  LST, Local Sidereal Time, is like Universal Time, only adjusted for your longitude, so, it also is not affected by daylight savings time or time zones, just longitude.

Time zones are arbitrary national boundries, they're more influenced by politics than geography.  They don't line up neatly on meridians, so, they would be really hard to factor in to any of these calculations.  Same with daylight savings time, the rules are different wherever you live.  Some countries practice it, some do not, and even if they do, they may or may not start it on the same day.

Unless you want to code the rules for every spot on the earth, all you can do is have the user input  time in Universal Time -- so, meaning they are responsible to know how many hours they are ahead of or behind UT and they make the adjustment, or, you can ask them if they are currently under daylight savings time, and what their time zone is and try to adjust for it... that is what that code you posted is doing.  

Most of the sites that do this sort of world mapping don't bother with local time.  But, you can add in the code if you want to!
0
 
iscodeAuthor Commented:
I see but this early comment from you then confuse me :-)
http://www.experts-exchange.com/Programming/Languages/Visual_Basic/Q_24363977.html#24279996

>>
To do it for right now, you would input the date and time that it is right now, however the time has to be UTC - universal time.  For me, I think I'm 4 hours behind UTC, so, the code at the top of your routine would be this...

        y = Year(Now)
        m = Month(Now)
        d = Microsoft.VisualBasic.Day(Now)
        hours = (Hour(Now) + 4) <--//////
        minutes = Minute(Now)
        seconds = Second(Now)
<<
0
 
mdouganCommented:
Well, that's what I was referring to in the last comment when I said this:

"all you can do is have the user input  time in Universal Time -- so, meaning they are responsible to know how many hours they are ahead of or behind UT and they make the adjustment"

When I was using Now, "now" is the time on my computer, which is set to the local time.  Normally, my local time in New York is 5 hours behind Universal Time (usually expressed -5), which would mean that to get to UT, I would add 5 hours to my local time.  But, since we're in daylight savings time, I'm only 4 hours behind UT.  So, by saying Hour(Now) + 4, I was converting my local time to UT.

But, you really shouldn't try to do that as part of your program, because the formula will be different for users all over the world... some people won't have daylight savings time, other people will be -6 or -7 or +3 etc hours behind or ahead of UT.  

So, either, you can ask the user to enter in UT, or, you can ask them to enter in their offset to UT if they know it (in which case at this point, I'd enter -4, but in Winter I'd enter -5).  But, if all you have them enter is their timezone and whether they are in daylight savings time, your program is going to have to be smart enough to calculate what the real UT is.

Keep in mind that LST - Local Sidereal Time is not the same as Local Time... local time is whatever the politicians have agreed the time should be in a particular place, LST is the precise time in relation to where the sun is based on the rotation of the earth and the longitude of the place.  I think there is some place in Canada that is far enough East to be in another timezone from the rest of its province, however, for political reasons, they keep the local time the same as the rest of the province... but, the LST would probably be at least an hour later... meaning, it's so far east, that it is really, say, three in the afternoon according to where the sun is, but their clocks only say 2.

So, if you are having the user enter their local time... it is going to be tricky for you to calculate.  It would be much easier if you just found another way to get the UT.

If you have the
0
 
mdouganCommented:
Forget that last sentence :-)  
0
 
iscodeAuthor Commented:
ok ok :-) but but
Yes this makes things more complicated but very interesting and thanks for this good details.  

Lets say we want to have sunrise and sunset (we are not doing that now...), as we have already lat/lon location from user we would need local time and ofcourse daylight saving and timezone right?

But I really think the ideal way to make a good interface is have the user enter as little information as possible.

So one way to get those information is from Date and Time Properties of the computer system.
The user has already enter those information there and set the timezone and etc.

Another way is to get the lat/lon from IP address but less acurate for pin point location and calculate timezone and etc for that area.

If you were using gps you would have those information from the unit
Just some ideas of posabilities :-)


 
0
 
mdouganCommented:
You should be able to google some formulas for calculating sunrise and sunset.  That would give you the times of sunrise and sunset in UT.  Then, you still have the problem of converting UT to local time.  

Local time on the computer is never really going to be accurate, because it is based on policy, not on geography, so, using the local time is difficult unless you have some source of data that says what each city's policy is.

If you have the users lat/lon, which isn't too hard... even if you have them choose their location from a map, or you get it from the GPS, and you find another way to get the UT (there are time servers where you can make a request to the server and it will return the UT value to you... sort of like the Java applet that we were working with above), then you should be able to calculate all the things you want.

If they have a GPS, the gps reading should also have the UT value, so, it is possible you can get it from there.

But, as far as being able to tell someone, your sunrise is going to be at 6:37am according to the time on your PC... that's not a very easy thing to do :-)
0
 
iscodeAuthor Commented:
Ok lets narrow this down,
we call Now like this

        y = Year(Now)
        m = Month(Now)
        d = Microsoft.VisualBasic.Day(Now)
        hours = (Hour(Now) + 4) <--Here we adjust to UT, this  depents on timezone
        minutes = Minute(Now)
        seconds = Second(Now)

so my question was, we need to have timezone who depends on daylight saving
to do the calculation and get UT.... right?

0
 
mdouganCommented:
Yes, I think so.  I don't know all of the complexities.

If you use "Now", you're going to be taking the time off of the local computer.  The time might be wrong... it is set to whatever the user has it set at.

The time on the computer might be adjusted for daylight savings time... but, how do you know?

The computer might be in Nova Scotia, where they bend the timezone way East so that all of the province shares the same time...

If you're trying to use Now, to get information about the users local Noon, or Sunrise or Sunset, then you have to be able to calculate UT from their local time (say, whatever is on their PC), and I'm suggesting that while you might be able to guess at this and be successful most of the time, it won't be perfect.  That's all!  

When looking at the PC's Now, you need to get their "location" which will usually look something like  "Eastern Daylight Time", and probably has the -4 hour adjustment value for UTC.  As I said, you can probably use that and be successful most of the time, but it is not guaranteed.

:-)

Since this question is getting pretty long, you should probably close it and open up a new question (the amount of points do not matter, just worried about the question becoming so large, it is not useful to anyone else).
0
 
iscodeAuthor Commented:
I just wanted to have it clear that we do need the timezone (and daylight saving if it is avalable) into the calculation to get the UT, it was just your answer at first that confused me :-)

About guarantee for correct time and other information we give that responsability to the user.

I absolut agree about closing this question now and do a continue. I think we made a record in length here but this is so darn fun...., you sure deserv the points and more dont be so humble :-)
I am sure others will find this information useful

Thanks for great work and all the time you put on this with me and see you on the next page

 
0
 
iscodeAuthor Commented:
0
 
mdouganCommented:
Cool, thanks!  See you on the next page!
0
Question has a verified solution.

Are you are experiencing a similar issue? Get a personalized answer when you ask a related question.

Have a better answer? Share it in a comment.

All Courses

From novice to tech pro — start learning today.