Attribute VB_Name = "Module1" Option Explicit Public Sub OverviewNotes() ' Notes ' These functions allow for the translation between three coordinate ' systems - the equatorial, the galactic and the orthogonal (xyz). ' Decimal degrees is the common positional format. Other utilities ' translate longitudinal decimal degrees back to right ascension format, ' latitudinal decimal degrees back to declination format, or ' distance measurements to light years and/or parsecs. ' Distance is required for orthogonal (xyz) 2-D or 3-D plotting ' ' 2-D plotting, in either the x-y or y-z planes can be ' done within a spreadsheet. 3-D plotting can be done ' with Euler and/or VRML rendering. ' ' The system of functions and related external programs ' are shown schematically as follows: ' ' Extended-object-boundaries=>AngleBetw2Objs=>Size ' | ' | Parallax2lightyrs=>| | ' Distance=>| Parallax2parsecs =>| | ' measure | Parsecs2light_yrs=>|------->--| ' | RVel2Mparsecs =>| | ' v ' Galactic l,b,dist,size ' ^ | ' RA=>Convert_RA2DecDeg---| | | ' |=>Equ2gal_coord_convert | ' Dec=>Convert_DMS2DecDeg-| | v ' | Gal2ortho_coord_convert ' RA <=\ v | ' Convert_DecDeg2RADMS<=Gal2equ_coord_convert | ' Dec<=/ v ' |-->-Star charting software(equ)---<--2-D-----Plotting ' | options ' | | ' v v ' 2-D Spreadsheet 3-D ' (x-y ortho, size) | ' / ' VRML_proto_writer--------<------------/ ' (text file) ' | ' ------------------------------- ' | | | | ' v v v v ' OpenGL DirectXL Netscape MS-Explorer ' ' If plotting in 3D with VRML and MS-Explorer, installation ' of Parallel Graphics enhanced VRML rendered plug-in for ' MS-Explorer is recommended. The native VRML plotting ' of MS-Explorer is too primitive for most renderings. ' ' A variety of catalogue input and output formats are supported. ' Handling of various text catalogue formats can be ' easily extended modified by the end user. ' ' For some functions, multiple coordinates are returned ' in a string separated by spaces. A generic string ' parsing function - GetCSWord - is included for decomposing ' these multi-variable strings into separate variables. ' ' These functions are written in Visual Basic for Applications ' for Windows users of MS-Excel for MS-Access. Excel users can ' install these functions in the Visual Basic Scripting Editor window ' and then access each function using the Excel "functions" button. ' The functions will appear under the "User defined functions" list ' and operate as built-in Excel functions or MS-Access ' modules. MS-Excel does not support polar or 3-D plotting. ' ' Robust run-time error trapping caused by improper ' pairing of source text with input formats is not provided. ' ' Author: Kurt A. Fisher 2003 fisherkaNO@SPAMcsolutions.net ' End Sub Public Function Gal2equ_coord_convert(ByVal l_dec_tmp As Double, ByVal b_dec_tmp As Double, ByVal RA_NGP_dec_tmp As Double, ByVal Dec_NGP_dec_tmp As Double, ByVal galLong_NCP_dec_tmp As Double) As Variant ' Receives galactic latitude ("b") and longitude ("l") in decimal form ' Returns equatorial right ascension and declination in ' decimal form as a string of "RA" Space "Dec". ' Three constants relating to the position of the galactic north pole and the ' celestral north pole in the galactic system are set according to ' the epoch applicable to the galactic latitude and longitude. ' ' RA_ngp_dec Dec_ngp_dec galLong_CP_dec ' 2000 192.85948 27.12825 122.932 ' 1950 192.25 27.4 123 ' ' Usage and test data: ' 1st source quadrant objects ' M57 ' Gal2equ_coord_convert(63.17,13.98,192.85948,27.12825,122.932) ' 283.394079491403 33.029750249075 ' 2nd source quadrant objects ' Association Cepheus OB 6 Epoch 2000 (Simbad query Ass Cep OB 6-) ' Gal2equ_coord_convert(105.7,0.12,192.85948,27.12825,122.932) ' 337.49304230001 57.999348445433 ' Association Cepheus OB 5 Epoch 2000 (Simbad query Ass Cep OB 5-) ' Gal2equ_coord_convert(118.44,+4.73,192.85948,27.12825,122.932) ' 1.24167085734632 67.1955701001172 ' 3rd source quadrant objects ' Duffet-Smith test problem in epoch 1950 at p. 44 ' Gal2equ_coord_convert(232.247778,+51.122222,192.25,27.4,123) ' 155.249925119455 10.0530878846083 ' Association Carina OB 1 Epoch 2000 (Simbad query Ass Car OB 1-) ' Gal2equ_coord_convert(286.53,-0.55,192.85948,27.12825,122.932) ' 159.498759821857 -59.1046928396689 ' 4th quadrant objects ' IC2602 tet Car Cluster ' Gal2equ_coord_convert(289.62,-4.89,192.85948,27.12825,122.932) ' 160.796946998307 -64.3945991521751 ' ' Algorithm sources: ' Declination: Benny, J. and Merrifeld, M. 1998. Galactic Astronomy. ' Princeton Univ. Press. p. 31. ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 94. ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. p. 44. ' Right ascension: Duffett-Smith ' ' Author: Kurt A. Fisher 2003 Dim l_dec As Double Dim b_dec As Double Dim RA_ngp_dec As Double ' RA of the north galactic pole Dim Dec_ngp_dec As Double ' Dec of the north galactic pole Dim l_CP_dec As Double ' Galactic longitude of the celestial pole Dim A, B, C, D, E, f, g, H, x, y As Double ' Temporary working variables Dim sinDecGP As Double Dim cosDecGP As Double Dim Dec_result_dec As Double Dim RA_result_dec As Double ' get galactic coordinates and constants l_dec = l_dec_tmp b_dec = b_dec_tmp RA_ngp_dec = RA_NGP_dec_tmp Dec_ngp_dec = Dec_NGP_dec_tmp l_CP_dec = galLong_NCP_dec_tmp ' create reused constants sinDecGP = Sin(Deg2Rad(Dec_ngp_dec)) cosDecGP = Cos(Deg2Rad(Dec_ngp_dec)) ' convert to equatorial declination Duffett steps 2-4 ' compute equation components ' sinDecGP A = sinDecGP ' sin_b B = Sin(Deg2Rad(b_dec)) ' cosDecGP C = cosDecGP ' cos_b D = Cos(Deg2Rad(b_dec)) ' cos(l_CP-l) E = Cos(Deg2Rad(l_CP_dec - l_dec)) ' set equatorial declination ' as sine f = (A * B) + (C * D * E) ' convert to radians g = ArcSin((f)) ' convert to degrees Dec_result_dec = Rad2Deg((g)) ' reset working variables A = 0 B = 0 C = 0 D = 0 E = 0 f = 0 g = 0 H = 0 x = 0 y = 0 ' convert right ascension ' compute the y intermediary coordinate Duffett step 4 ' cos_b A = Cos(Deg2Rad(b_dec)) ' cos(l+90-l_CP) or l-33 B = Cos(Deg2Rad(l_dec + 90 - l_CP_dec)) ' evalute y intermediary coordinate y = A * B ' reset working variables A = 0 B = 0 ' compute the x intermediary coordinate Duffett step 5 ' sin_b A = Sin(Deg2Rad(b_dec)) ' cos_Dec_GP B = cosDecGP ' cos_b C = Cos(Deg2Rad(b_dec)) ' sin_GP D = sinDecGP ' Sin (l + 90 - l_CP) Or l - 33 E = Sin(Deg2Rad(l_dec + 90 - l_CP_dec)) ' evalute x intermediary coordinate x = (A * B) - (C * D * E) ' determine the inverse of the tangent Duffett step 6-7a g = Atn((y / x)) H = Rad2Deg((g)) ' Convert radians to degrees ' adjust quadrant Duffet step 7b 'using Meeus ambiguity resolution (Meeus at p. 9) If x < 0 Then H = H + 180 ' Add ascending node of equatorial plane or + 192 Duffet step 8 H = H + RA_ngp_dec ' Adjust quadrant if addition exceeds 360 deg If H > 360 Then H = H - 360 RA_result_dec = H ' return result Gal2equ_coord_convert = Trim(Str(RA_result_dec)) & " " & Trim(Str(Dec_result_dec)) End Function Public Function Equ2gal_coord_convert(ByVal RA_dec_tmp As Double, ByVal Dec_dec_tmp As Double, ByVal RA_NGP_dec_tmp As Double, ByVal Dec_NGP_dec_tmp As Double, ByVal galLong_NCP_dec_tmp As Double) As Variant ' Receives in equatorial right ascension and declination in decimal form ' Returns galactic latitude ("b") and longitude ("l") in ' decimal form as a string of "l" Space "b". ' Three constants relating to the position of the galactic north pole and the ' celestral north pole in the galactic system are set according to ' the epoch applicable to the galactic latitude and longitude. ' ' RA_ngp_dec Dec_ngp_dec galLong_CP_dec ' 2000 192.85948 27.12825 122.932 ' 1950 192.25 27.4 123 ' ' Usage and test data: ' 1st source quadrant objects ' Association Cepheus OB 5 Epoch 2000 (Simbad query Ass Cep OB 5-) ' Equ2gal_coord_convert(1.242, +67.20,192.85948,27.12825,122.932) ' 118.44,4.7343 ' 2nd source quadrant objects ' Duffet-Smith test problem in epoch 1950 at p. 44 ' Equ2gal_coord_convert(155.250, +10.053,192.25,27.4,123) ' 232.248,51.1222 ' Association Carina OB 1 Epoch 2000 (Simbad query Ass Car OB 1-) ' Equ2gal_coord_convert(159.498, -59.105,192.85948,27.12825,122.932) ' 286.53,-.550458 ' 3rd source quadrant objects ' M101 ' Equ2gal_coord_convert(210.80,54.35,192.85948,27.12825,122.932) ' ' 4th quadrant objects ' Association Cepheus OB 6 Epoch 2000 (Simbad query Ass Cep OB 6-) ' Equ2gal_coord_convert(337.493, +57.99,192.85948,27.12825,122.932) ' 105.7,.11199 ' ' Algorithm sources: ' Declination: Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 94. ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. p. 43. ' Right ascension: Duffett-Smith ' ' Author: Kurt A. Fisher 2003 Dim l_dec As Double Dim b_dec As Double Dim RA_ngp_dec As Double ' RA of the north galactic pole Dim Dec_ngp_dec As Double ' Dec of the north galactic pole Dim l_CP_dec As Double ' Galactic longitude of the celestial pole Dim A, B, C, D, E, f, g, H, x, y As Double ' Temporary working variables Dim sinDecGP As Double Dim cosDecGP As Double Dim b_result_dec As Double Dim l_result_dec As Double Dim RA_dec As Double Dim Dec_dec As Double ' get equalatorial coordinates and constants RA_dec = RA_dec_tmp Dec_dec = Dec_dec_tmp RA_ngp_dec = RA_NGP_dec_tmp Dec_ngp_dec = Dec_NGP_dec_tmp l_CP_dec = galLong_NCP_dec_tmp ' create reused constants sinDecGP = Sin(Deg2Rad(Dec_ngp_dec)) cosDecGP = Cos(Deg2Rad(Dec_ngp_dec)) ' convert to galactic latitude Duffett steps 2-4 ' compute equation components ' cos_Dec A = Cos(Deg2Rad((Dec_dec))) ' cosDecGP B = cosDecGP ' cos(l_CP-l) C = Cos(Deg2Rad(RA_dec - RA_ngp_dec)) 'less 192 deg Epoch 2000 ' sin_b D = Sin(Deg2Rad((Dec_dec))) ' sinDecGP E = sinDecGP ' set galactic declination Duffett Step 3 ' as sine f = (A * B * C) + (D * E) ' convert to radians g = ArcSin((f)) ' convert to degrees Duffett Step 4 b_result_dec = Rad2Deg((g)) ' reset working variables A = 0 B = 0 C = 0 D = 0 E = 0 f = 0 g = 0 H = 0 x = 0 y = 0 ' convert right ascension ' compute the y intermediary coordinate Duffett step 5 ' sinDec_dec A = Sin(Deg2Rad((Dec_dec))) ' sin_l_result_dec B = Sin(Deg2Rad(b_result_dec)) ' sin_GP C = sinDecGP ' evalute y intermediary coordinate y = A - (B * C) ' reset working variables A = 0 B = 0 C = 0 ' compute the x intermediary coordinate Duffett step 6 ' cos_Dec_dec A = Cos(Deg2Rad((Dec_dec))) ' Sin (RA - RA_ngp) ' or less 192 B = Sin(Deg2Rad(RA_dec - RA_ngp_dec)) ' cos_GP C = cosDecGP ' evalute x intermediary coordinate x = A * B * C ' determine the inverse of the tangent Duffett step 8a g = Atn((y / x)) H = Rad2Deg((g)) ' Convert radians to degrees ' adjust quadrant Duffet step 8b 'using Meeus ambiguity resolution (Meeus at p. 9) If x < 0 Then H = H + 180 ' Add ascending node of galactic plane or l_CP_dec-90 or + 33 Duffet step 9 H = H + (l_CP_dec - 90) ' Adjust quadrant if addition exceeds 360 deg If H > 360 Then H = H - 360 ' Correction 11/2005 If H < 0 Then H = H + 360 l_result_dec = H ' return result Equ2gal_coord_convert = Trim(Str(l_result_dec)) & " " & Trim(Str(b_result_dec)) End Function Public Function Gal2ortho_coord_convert(ByVal l_dec_tmp As Double, ByVal b_dec_tmp As Double, ByVal Distance_tmp As Double, ByVal l_size_arcmin_tmp As Double, ByVal b_size_arcmin_tmp As Double) As Variant ' Receives galactic latitude ("b"), longitude ("l"), and distance ' in decimal form; "l" object size and "b" object size in arcminutes; and a ' plotting scalar. Returns x,y,z coordinates and object diameters in xyz units, ' in decimal form, as a string. ' ' Used for plotting objects in 2-D spreadsheet and 3-D VRML. ' ' Size "0" designates a point object for plotting purposes. ' Note for purposes of plotting objects in the galactic plane ' that the position of the first and fourth and second and ' third quadrants are reversed. To generate a plot that ' corresponds to the visual appearance of the celestial sky, ' the sign of the x coordinate would be reversed. ' This distinction is not addressed in this function, ' but is left to the plotting utility. ' ' Usage and test data: ' ' Eight quadrant check ' 1st source quadrant stars ' +z ' Name l b Epoch RA Dec parallax ' ALF HER/ 35.53/+27.82/2000/258.66/+14.39/8.53 ' ? Gal2ortho_coord_convert(35.53,+27.82,100,10,10) ' x y z x_size y_size ' 71.97495 51.39611 46.66953 .29088 .29088 ' -z ' BET CAP/ 29.15/-26.37/2000/305.25/-14.78/9.48 ' ? Gal2ortho_coord_convert(29.15,-26.37,100,10,10) ' 78.24708 43.64124 -44.41661 .29088 .29088 ' ' 2nd source quadrant stars ' +z ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/26.38 ' ? Gal2ortho_coord_convert(142.85,+51.01,100,10,10) ' -50.14962 37.9967 77.72557 .29088 .29088 ' -z ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/16.36 ' ? Gal2ortho_coord_convert(127.11,-27.10,100,10,10) ' -53.71073 70.99256 -45.55449 .29088 .29088 ' ' 3rd source quadrant stars ' +z ' ALF GEM/187.44/+22.48/2000/113.65/+31.89/63.27 ' ? Gal2ortho_coord_convert(187.44,+22.48,100,10,10) ' -91.62337 -11.96484 38.23609 .29088 .29088 ' -z ' ALF ORI/199.79/ -8.96/2000/ 88.79/ +7.41/7.63 ' ? Gal2ortho_coord_convert(199.79,-8.96,100,10,10) ' -92.94578 -33.44421 -15.57448 .29088 .29088 ' ' 4th source quadrant stars ' +z ' ALF VIR/316.11/+50.84/2000/201.30/-11.16/12.44 ' ? Gal2ortho_coord_convert(316.11,+50.84,100,10,10) ' 45.50958 -43.77956 77.53855 .29088 .29088 ' -z ' ALF ARA/340.76/ -8.83/2000/262.96/-49.88/13.46 ' ? Gal2ortho_coord_convert(340.76,-8.83,100,10,10) ' 93.29566 -32.56203 -15.35032 .29088 .29088 ' ' Algorithm sources: ' X,Y,Z conversion ' Meeus, J. 1988. Astronomical Formulae for Calculators. 4th ed. ' Willmann-Bell. p. 85 (rectangular coordinates of the sun) ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 172. ' Jan Gullberg. 1997. Formulas Transforming Spherical Coordinates ' into Orthogonal Coordinates. In Mathematics: From the Birth ' of Numbers. W.W. Norton. p. 584-585. ' MS Neatcode 1997 Lat/Lon_XYZ example function ' ' Size conversion ' Author's ' 1/2 S = Abs( R * tan(rads(1/2*angle)) ) ' S = 2 * Abs( R * tan(rads(1/2*angle)) ) ' ' Author: Kurt A. Fisher 2003 Dim l_dec As Double Dim b_dec As Double Dim R As Double ' Distance in arbitrary units Dim l_size_arcsec As Double ' Set to 0 for stellar point objects Dim b_size_arcsec As Double ' Set to 0 for stellar point objects Dim l_size_arcdeg As Double ' Set to 0 for stellar point objects Dim b_size_arcdeg As Double ' Set to 0 for stellar point objects Dim x As Double ' X result Dim y As Double ' Y result Dim z As Double ' Z result Dim s_x As Double ' Size in xyz in x dimension Dim s_y As Double ' Size in xyz in y dimension ' get galactic coordinates l_dec = l_dec_tmp b_dec = b_dec_tmp R = Distance_tmp l_size_arcdeg = l_size_arcmin_tmp / 60 ' convert arcsecs to degrees b_size_arcdeg = b_size_arcmin_tmp / 60 ' convert arcsecs to degrees ' convert X,Y,Z coordinates x = R * Sin(Deg2Rad((90 - b_dec))) * Cos(Deg2Rad((l_dec))) y = R * Sin(Deg2Rad((90 - b_dec))) * Sin(Deg2Rad((l_dec))) z = R * Cos(Deg2Rad((90 - b_dec))) ' convert distance s_x = Abs(2 * R * Tan(Deg2Rad(l_size_arcdeg / 2))) s_y = Abs(2 * R * Tan(Deg2Rad(b_size_arcdeg / 2))) ' return result Gal2ortho_coord_convert = Trim(Str(x)) & " " & Trim(Str(y)) & " " & Trim(Str(z)) & " " & Trim(Str(s_x)) & " " & Trim(Str(s_y)) End Function Function Equ2Horizon_coord_convert(ByVal objRA_dec_tmp As Double, ByVal objDec_dec_tmp As Double, ByVal strTime_tmp As String, ByVal dayDate_tmp As Date, ByVal decLongitude_tmp As Double, ByVal decLatitude_tmp As Double) As Variant ' Receives in equatorial right ascension and declination in decimal form ' Returns local horizon system coordinates azimuth ("az") and altitude ("a") in ' decimal form as a string of "az" Space "a". ' Azimuth is reckoned from North meridian instead of Meuss's use of South meridian. ' Note in navigation applications - azimuth - is the "to bearing." ' Used to build narrative report on current status based on relationship between ' rising time, current time and setting time, based on object's current altitude. ' Uses: TimeLocalSidereal ' Source: Built around ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. pp. 35-37, Secs 25 and 24. ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. pp. 94-95. ' Author: Kurt A. Fisher 10/2005 ' ' Usage and test data: ' Duffet-Smith example at page 35 ' Where 18.539167 RA equals 278.087505 degrees ' ? Equ2Horizon_coord_convert(278.087505, 23.219444, "14:36:51.67", #4/22/1980#, -64#, 52#) ' 283.271728215061 19.3337758286988 ' ? Equ2Horizon_coord_convert(347.31801,-6.64924,"19:21:00",#4/10/1987#,-77.065556,38.921389) ' 248.088414619519 15.1718465489576 ' 1st quad, high declination, Andromeda Galaxy M31 from SLC, UT ' ? Equ2Horizon_coord_convert(10.675014, 41.266667, "21:02:00", #10/30/2005#, -111.890833, 40.756389) ' 36.5827323942514 4.68595609131246 ' Checks against http://www.roman-britain.org/astronomy/astro.htm# ' 1st quad, ultra high declination, from SLC, UT ' Equ2Horizon_coord_convert(10.675014, 81.266667, "21:02:00", #10/30/2005#, -111.890833, 40.756389) ' 8.4317970029173 35.0853510131691 ' Checks against http://www.roman-britain.org/astronomy/astro.htm# ' 2nd quadrant source objects, above horizon ' Mars from SLC ' ? Equ2Horizon_coord_convert(45.30962, 16.224, "06:00:00", #10/31/2005#, -111.890833, 40.756389) ' 127.761658515905 55.813214328665 ' Checks against MPC epheremis ' 2nd quadrant source objects, extreme near south celestial pole horizon ' ? Equ2Horizon_coord_convert(150#, -80#, "06:00:00", #10/31/2005#, -111.890833, 40.756389) ' 169.135817674973 -47.0008156136534 ' Checks against http://www.roman-britain.org/astronomy/astro.htm# ' 2nd quadrant source objects, below horizon ' Saturn ' ? Equ2Horizon_coord_convert(133.44501, 17.96344, "05:00:00", #10/31/2005#, -111.890833, 40.756389) ' 48.4491184842782 -15.5963969322819 ' Checks against MPC epheremis ' 3rd source quadrant objects, above horizon ' Hypothetical object from SLC ' ? Equ2Horizon_coord_convert(345# , -45#, "06:00:00", #10/31/2005#, -111.890833, 40.756389 ) ' 202.493216228253 -0.63932541425704 ' Checks against http://www.roman-britain.org/astronomy/astro.htm# ' 3rd source quadrant objects, extreme below horizon ' Hypothetical object from SLC ' ? Equ2Horizon_coord_convert(345# , -85#, "06:00:00", #10/31/2005#, -111.890833, 40.756389 ) ' 183.362777115236 -36.5009520227761 ' Checks against http://www.roman-britain.org/astronomy/astro.htm# ' 4th quadrant objects, below horizon ' Sun test at night from SLC, Utah - ' ? Equ2Horizon_coord_convert(215.50685, -14.13643, "05:00:00", #10/31/2005#, -111.890833, 40.756389) ' 303.463075847706 -50.9790772544255 ' Checks against MPC epheremis ' Define variables Dim A, B, C, D, E, f, g, H, x, y As Double ' Temporary working variables Dim strOut As String Dim objRA_dec As Double Dim objDecDeg As Double Dim strTime As String Dim dayDate As Date Dim decLongitude As Double Dim decLatitude As Double Dim dblHourAngle As Double Dim sinAltitude As Double Dim cosAzimuth As Double Dim dblAzimuth As Double Dim dblAltitude As Double ' Get input variables objRA_dec = objRA_dec_tmp objDecDeg = objDec_dec_tmp strTime = strTime_tmp dayDate = dayDate_tmp decLongitude = decLongitude_tmp decLatitude = decLatitude_tmp ' Find value of hour angle H per Duffett-Smith p. 35 and Meuss p. 92 ' Find the decimal hour time at Greenwich and convert to decimal hours ' Find the sidereal hour angle and divide by 15 degrees per hour B = (CDbl(GetCSWord(TimeLocalSidereal(strTime, dayDate, decLongitude), 2)) / 15#) f = (objRA_dec / 15#) dblHourAngle = B - f ' Test if the angle is less than zero, adjust boundary to ' a range of 0-24 hours ' If greater or equal than 24, subtract 24 and add one to the day If dblHourAngle < 0 Then dblHourAngle = dblHourAngle + 24 End If If dblHourAngle > 24 Then dblHourAngle = dblHourAngle - 24 End If If dblHourAngle = 24 Then dblHourAngle = 0 End If ' Multiply the hour angle to convert it to degrees (Duffett-Smith step 3, Sec. 25) dblHourAngle = dblHourAngle * 15# ' Find the sine of the altitude - Duffett-Smith step 4 sinAltitude = CDbl((Sin(Deg2Rad(objDecDeg)) * Sin(Deg2Rad(decLatitude))) + (Cos(Deg2Rad(objDecDeg)) * Cos(Deg2Rad(decLatitude)) * Cos(Deg2Rad(dblHourAngle)))) ' Find a from the sine's inverse - Duffett-Smith's step 5 dblAltitude = Rad2Deg(ArcSin(sinAltitude)) ' Stop ' Find the cosine of the azimuth using the computed altitude - Duffett-Smith's step 6 ' Compute the numerator C = CDbl(Sin(Deg2Rad(objDecDeg)) - (Sin(Deg2Rad(decLatitude)) * Sin(Deg2Rad(dblAltitude)))) ' Compute the denominator D = CDbl((Cos(Deg2Rad(decLatitude)) * Cos(Deg2Rad(dblAltitude)))) ' Find the cosine of the azimuth cosAzimuth = CDbl(C / D) ' Find the inverses of the sine and cosine to reduce to angles - Duffett-Smith's step 7 dblAzimuth = Rad2Deg(ArcCos(cosAzimuth)) ' Adjust the azimuth value for 360 deg rotation per Duffett-Smith's step 8 ' If sign of the sine of HourAngle is postive or zero, adjust the value by subtracting ' it from 360 degrees ' Test the sin of the Hourangle E = CDbl(Sin(Deg2Rad(dblHourAngle))) If Sgn(E) >= 0 Then dblAzimuth = 360# - dblAzimuth End If ' Return the result. Equ2Horizon_coord_convert = CStr(Trim(dblAzimuth)) & " " & CStr(Trim(dblAltitude)) End Function Public Function Convert_RA2DecDeg(ByVal RA_tmp As String, ByVal s_inAngFormat_tmp As String, ByVal s_outAngFormat_tmp As String) As Variant ' Converts RA to a decimal degree ' 360 deg = 24 hours ' 15 deg = 1 hour ' 15 deg = 60 minutes ' 1/4 deg = 1 minute ' 1/4 deg = 60 seconds ' 1/4 * 1/60 deg = 1 second ' 0.004167 deg = 1 second ' ' s_intAngform is the input format of the right ascension ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' "HHMMSS.SS" ' "HH:MM:SS" ' "HHh MMm SSs" ' ' s_outAngform is the desired output format, as either string or double, ' consisting of valid types: ' "string" or "double" ' ' Usage and test data: ' FK5 data ' Name l b Epoch RA_d Dec_d RA_hms Dec_dms ' ALF HER/ 35.53/+27.82/2000/258.66/+14.39/17 14 38.86 +14 23 25.2 ' Convert_RA2DecDeg("171439","HHMMSS","string") yields 258.662513 ' BET CAP/ 29.15/-26.37/2000/305.25/-14.78/20 21 00.68 -14 46 52.9 ' Convert_RA2DecDeg("20 21 01","HH MM SS","string") yields 305.254167 ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/11 03 43.67 +61 45 03.7 ' Convert_RA2DecDeg("11:03:44","HH:MM:SS","string") yields 165.933348 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_RA2DecDeg("01h 09m 44s","HHh MMm SSs","string") yields 017.433348 ' ALF GEM/187.44/+22.48/2000/113.65/+31.89/07 34 35.86 +31 53 17.8 ' Convert_RA2DecDeg("073435.86","HH MM SS.SS","string") yields 113.650012 ' Convert_RA2DecDeg("073436","HHMMSS","string") yields 113.645845 ' BET ORI/209.24/-25.25/2000/ 78.63/ -8.20/05 14 32.27 -08 12 05.9 ' Convert_RA2DecDeg("051432","HHMMSS","string") yields 078.633344 ' ' Source: Author ' Author: Kurt A. Fisher 2003 ' initialize variables Dim s_RA As String Dim s_DecDeg As String Dim s_inAngFormat As String Dim s_outAngFormat As String Dim iDeg As Integer 'Degrees Dim iMin As Integer 'Minutes Dim dSec As Integer 'Seconds Dim A As Double 'Temporary working variables ' get function parameters s_RA = RA_tmp s_inAngFormat = s_inAngFormat_tmp s_outAngFormat = s_outAngFormat_tmp ' translate the input format ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' "HHMMSS.SS" ' "HH:MM:SS" ' "HHh MMm SSs" ' This format includes HH° MM' SS" Select Case s_inAngFormat Case "HHMMSS" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 3, 2)) dSec = CDbl(Mid(s_RA, 5, 2)) Case "HH MM SS" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 4, 2)) dSec = CDbl(Mid(s_RA, 7, 2)) Case "HH MM SS.SS" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 4, 2)) dSec = CDbl(Mid(s_RA, 7, 5)) Case "HHMMSS.SS" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 3, 2)) dSec = CDbl(Mid(s_RA, 5, 4)) Case "HH:MM:SS" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 4, 2)) dSec = CDbl(Mid(s_RA, 7, 2)) Case "HHh MMm SSs" iDeg = CInt(Mid(s_RA, 1, 2)) iMin = CInt(Mid(s_RA, 5, 2)) dSec = CDbl(Mid(s_RA, 9, 2)) Case Else ' graceful return on unsupported type without error message iDeg = 0 iMin = 0 dSec = 0 End Select ' Convert to decimal degrees A = (iDeg * 15) + (iMin * 0.25) + (dSec * 0.004167) ' Return the result in requested format Select Case s_outAngFormat Case "string" ' Make string formatted output uniform size by adding leading zeros Convert_RA2DecDeg = CStr(Format(A, "00#.######")) Case "double" Convert_RA2DecDeg = CDbl(A) Case Else ' Gracefully handle error condition by returning double Convert_RA2DecDeg = CDbl(A) End Select End Function Public Function Convert_DMS2DecDeg(ByVal DecDeg_tmp As String, ByVal inAngFormat_tmp As String, ByVal outAngFormat_tmp As String) As Variant ' Converts sDDMMSS formated angles to ' decimal angle format (and not to decimal hour format) ' Used to convert text data from astronomical catalogues ' into numeric data types. ' inAngFormat is the input format, consisting of valid types: ' "sDDMMSS" ' Sign is optional ' "sDD MM SS" ' "sDD MMm SS.S" ' High precision format Dec ' "sDDd MMm SSs" ' This format includes sDD° MM' SS" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format ' "sDDDd MMm SSs" ' outAngform is the desired output format, as either string or double, ' consisting of valid types: ' "string" ' Two position padded angle +61 45 36 ' "string3" ' Three position padded angle +061 45 36 ' "double" ' ' Usage and test data: ' FK5 data ' Name l b Epoch RA_d Dec_d RA_hms Dec_dms ' ALF HER/ 35.53/+27.82/2000/258.66/+14.39/17 14 38.86 +14 23 25.2 ' Convert_DMS2DecDeg("+142325","sDDMMSS","string") yields +014.390278 ' BET CAP/ 29.15/-26.37/2000/305.25/-14.78/20 21 00.68 -14 46 52.9 ' Convert_DMS2DecDeg("-14 45 53","sDD MM SS","double") yields -14.764722 ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/11 03 43.67 +61 45 03.7 ' Convert_DMS2DecDeg("+61 45 03.7","sDD MM SS.S","string") yields +61.751111 ' Convert_DMS2DecDeg("+61 45 03.7","sDD MM SS.S","string3") yields +061.751111 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_DMS2DecDeg("+0010944s","sDDDMMSS","string3") yields +001.162222 ' ALF GEM/187.44/+22.48/2000/113.65/+31.89/07 34 35.86 +31 53 17.8 ' Convert_DMS2DecDeg("031 53 17.8","sDDD MM SS.S","string3") yields +031.052222 ' BET ORI/209.24/-25.25/2000/ 78.63/ -8.20/05 14 32.27 -08 12 05.9 ' Convert_DMS2DecDeg("+005d 14m 32s","sDDDd MMm SSs","double") yields 5.242222 ' ' Sources: ' Getz, Ken. 1997. VB5 Developer's Guide. ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. Press. p. 44. ' ' Author: Kurt Fisher 2003 ' initialize variables Dim s_DecDeg As String Dim s_inAngFormat As String Dim s_outAngFormat As String Dim iSign As Integer 'Angle sign Dim iDeg As Integer 'Degrees Dim iMin As Integer 'Minutes Dim dSec As Integer 'Seconds Dim A As Double 'Temporary working variables ' get function parameters s_DecDeg = DecDeg_tmp s_inAngFormat = inAngFormat_tmp s_outAngFormat = outAngFormat_tmp ' translate the input format ' "sDDMMSS" ' Sign is optional ' "sDD MM SS" ' "sDD MMm SS.S" ' High precision format ' "sDDd MMm SSs" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format ' "sDDDd MMm SSs" Select Case s_inAngFormat Case "sDDMMSS" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 6, 2)) Case "sDD MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 8, 2)) Case "sDD MM SS.S" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 8, 4)) Case "sDD MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 8, 2)) Case "sDDd MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 10, 2)) Case "sDDDMMSS" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 2)) Case "sDDD MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) Case "sDDD MM SS.S" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 4)) Case "sDDDd MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 7, 2)) dSec = CDbl(Mid(s_DecDeg, 11, 2)) Case Else ' graceful return on unsupported type without error message iDeg = 0 iMin = 0 dSec = 0 End Select If iDeg < 0 Then iSign = -1 iDeg = Abs(iDeg) Else iSign = 1 End If 'Determine the decimal value A = Round(((iDeg + (iMin / 60) + (dSec / 3600)) * iSign), 6) ' Return the result in requested format Select Case s_outAngFormat Case "string" ' Make string formatted output uniform size by adding plus sign If iSign = 1 Then Convert_DMS2DecDeg = "+" & CStr(Format(A, "0#.######")) Else ' if sign is minus Convert_DMS2DecDeg = CStr(Format(A, "0#.######")) End If Case "string3" ' Make string formatted output uniform size by adding plus sign If iSign = 1 Then Convert_DMS2DecDeg = "+" & CStr(Format(A, "00#.######")) Else ' if sign is minus Convert_DMS2DecDeg = CStr(Format(A, "00#.######")) End If Case "double" Convert_DMS2DecDeg = CDbl(A) End Select End Function Public Function Convert_DecDeg2RADMS(ByVal DecDeg_tmp As Double, ByVal s_outAngFormat_tmp As String) As String ' Converts decimal number angles to HHMMSS or sDDMMSS format ' Used to convert decimal angle data into formatted text data ' for use in astronomical catalogues ' outAngFormat is the output string format, consisting of valid types: ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' High precision format RA ' "HHMMSS.SS" ' High precision format RA w/o RA ' "HH° MMm SS.SS" ' e.g. 12° 51' 23.36" ' "HH:MM:SS" ' "HHh MMm SSs" ' "DDD.DD" ' Leave as decimal for galactic longitude ' "sDD.DD" ' Leave as decimal for galactic latitude ' "sDDMMSS" ' "sDD MM SS" ' "sDD MM SS.S" ' High precision format Dec ' "sDDMMSS.S" ' High precision format Dec w/o spaces ' "sDD° MMm SS.S" ' -12° 51' 23.5" Use Alt+248 to generate ° symbol ' "sDDd MMm SSs" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format Dec ' "sDDDd MMm SSs" ' ' Exceptions: ' There is no error trapping to prevent a negative Dec ' being feed to the RA hours conversion function. ' A errorneous negative dec will be translated as ' a postive RA and reported without error ' Usage and testing ' FK5 data ' Name l b Epoch RA_d Dec_d RA_hms Dec_dms ' ALF HER/ 35.53/+27.82/2000/258.66/+14.39/17 14 38.86 +14 23 25.2 ' Convert_DecDeg2RADMS(258.66,"HHh MMm SS.SS") yields 17h 14' 38.40" ' BET CAP/ 29.15/-26.37/2000/305.25/-14.78/20 21 00.68 -14 46 52.9 ' Convert_DecDeg2RADMS(-14.78,"sDD° MMm SS.S") yields -14º 46' 48.0" ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/11 03 43.67 +61 45 03.7 ' Convert_DecDeg2RADMS(165.93,"HHMMSS") yields 110343 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Used galactic coordinate catalogue formatting ' Convert_DecDeg2RADMS(27.1139827,"DDD.DD") yields 027.11 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Used galactic coordinate catalogue formatting ' Convert_DecDeg2RADMS(-27.101981,"sDD.DD") yields -27.10 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_DecDeg2RADMS(35.62,"sDDMMSS") yields +353712 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_DecDeg2RADMS(35.62,"sDDMMSS") yields +353712 ' ALF GEM/187.44/+22.48/2000/113.65/+31.89/07 34 35.86 +31 53 17.8 ' Convert_DecDeg2RADMS(113.65,"HH MM SS.SS") yields 07 34 36.00 ' BET ORI/209.24/-25.25/2000/ 78.63/ -8.20/05 14 32.27 -08 12 05.9 ' Convert_DecDeg2RADMS(-8.20,"sDD MM SS.S") yields -08 12 00.0 ' ALF VIR/316.11/+50.84/2000/201.30/-11.16/13 25 11.58 -11 09 40.8 ' Convert_DecDeg2RADMS(201.30,"HHh MMm SSs") yields 13h 25m 12s ' ALF ARA/340.76/ -8.83/2000/262.96/-49.88/17 31 50.49 -49 52 34.1 ' Convert_DecDeg2RADMS(-49.88,"sDDDd MMm SSs") yields -049d 52m 48s ' ALF CMA/227.23/ -8.89/2000/101.29/-16.72/06 45 08.92 -16 42 58.0 ' Convert_DecDeg2RADMS(101.29,"HH:MM:SS") yields 06:45:10 ' initialize variables Const RndError As Double = 0.00000001 'Adjusts for rounding error in binary Dim chrDegree As String chrDegree = Chr(186) Dim chrSingleQuote As String chrSingleQuote = Chr(39) Dim chrDoubleQuote As String chrDoubleQuote = Chr(34) Dim s_outAngFormat As String Dim iSign As Integer 'Sign Dim dDeg As Double 'Degrees Dim dDeg_sav As Double ' Degrees Dim dMin As Double 'Minutes Dim dSec As Double 'Seconds Dim sSign As String 'Sign Dim sDeg As String 'Degrees Dim sMin As String 'Minutes Dim sSec As String 'Seconds Dim sSec1p As String 'Seconds precision 1 digit Dim sSec2p As String 'Seconds precision 2 digit Dim A, B, C As Double 'Temporary working variables ' get function parameters dDeg = DecDeg_tmp dDeg_sav = DecDeg_tmp s_outAngFormat = s_outAngFormat_tmp ' If dDeg < 0# Then If Sgn(dDeg) < 0# Then iSign = -1 dDeg = Abs(dDeg) Else iSign = 1 End If ' If output is in RA hours, convert decimal degrees to ' decimal hours If InStr(1, s_outAngFormat, "HH", vbTextCompare) > 0 Then ' dDeg = Convert_DecDeg2HDec(dDeg) A = dDeg / 15 ' Note there is no error trapping to prevent a negative RA ' being feed to this hours function. ' Otherwise, no adjustment is necessary Else A = dDeg ' Temporary degrees storage End If ' Determine Hours-Deg, Min, Sec dDeg = Int(A + RndError) 'Get degrees or hours dMin = (A - dDeg) * 60# 'Get minutes & seconds dSec = (dMin - Int(dMin + RndError)) * 60# 'Get seconds dMin = Int(dMin + RndError) 'Get minutes ' Format the output string ' preformat with leading zeros If iSign < 0# Then sSign = "-" Else sSign = "+" ' if output format is three position decimal angle, ' change leading zero to three positions If InStr(1, s_outAngFormat, "DDD", vbTextCompare) > 0 Then sDeg = CStr(Format(Int(dDeg), "000")) Else sDeg = CStr(Format(Int(dDeg), "00")) End If sMin = CStr(Format(Int(dMin), "00")) sSec = CStr(Format(Round(dSec), "00")) sSec1p = CStr(Format(dSec, "00.0")) sSec2p = CStr(Format(dSec, "00.00")) ' translate the output format ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' High precision format RA ' "HHMMSS.SS" ' High precision format RA w/o spaces ' "HH° MMm SS.SS" ' e.g. 12° 51' 23.36" ' "HH:MM:SS" ' "HHh MMm SSs" ' ' "HH° MMm SSs" ' e.g. 12° 51' 23" ' "DDD.DD" ' Leave as decimal for galactic longitude ' "sDD.DD" ' Leave as decimal for galactic latitude ' "sDDMMSS" ' "sDD MM SS" ' "sDD MM SS.S" ' High precision format Dec ' "sDDMMSS.S" ' High precision format Dec matching ' "sDD° MMm SS.S" ' e.g. -12° 51' 23.5" ' "sDDd MMm SSs" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format Dec ' "sDDD° MMm SS.S" ' e.g. -121° 51' 23.5" ' "sDDDd MMm SSs" Select Case s_outAngFormat Case "HHMMSS" Convert_DecDeg2RADMS = sDeg & sMin & sSec Case "HH MM SS" Convert_DecDeg2RADMS = sDeg & " " & sMin & " " & sSec Case "HH MM SS.SS" Convert_DecDeg2RADMS = sDeg & " " & sMin & " " & sSec2p Case "HHMMSS.SS" Convert_DecDeg2RADMS = sDeg & sMin & sSec2p Case "HH° MMm SS.SS" Convert_DecDeg2RADMS = sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec2p & chrDoubleQuote Case "HH:MM:SS" Convert_DecDeg2RADMS = sDeg & ":" & sMin & ":" & sSec Case "HHh MMm SSs" Convert_DecDeg2RADMS = sDeg & "h " & sMin & "m " & sSec & "s" Case "HH° MMm SSs" Convert_DecDeg2RADMS = sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec & chrDoubleQuote Case "DDD.DD" ' No error trapping for passing negative galactic longitude Convert_DecDeg2RADMS = CStr(Format(Abs(dDeg_sav), "000.00")) Case "sDD.DD" Convert_DecDeg2RADMS = sSign & CStr(Format(Abs(dDeg_sav), "00.00")) Case "sDDMMSS" Convert_DecDeg2RADMS = sSign & sDeg & sMin & sSec Case "sDD MM SS" Convert_DecDeg2RADMS = sSign & sDeg & " " & sMin & " " & sSec Case "sDD MM SS.S" Convert_DecDeg2RADMS = sSign & sDeg & " " & sMin & " " & sSec1p Case "sDDMMSS.S" Convert_DecDeg2RADMS = sSign & sDeg & sMin & sSec1p Case "sDD° MMm SS.S" Convert_DecDeg2RADMS = sSign & sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec1p & chrDoubleQuote Case "sDDd MMm SSs" Convert_DecDeg2RADMS = sSign & sDeg & "d " & sMin & "m " & sSec & "s" Case "sDDDMMSS" Convert_DecDeg2RADMS = sSign & sDeg & sMin & sSec Case "sDDD MM SS" Convert_DecDeg2RADMS = sSign & sDeg & " " & sMin & " " & sSec Case "sDDD MM SS.S" Convert_DecDeg2RADMS = sSign & sDeg & " " & sMin & " " & sSec1p Case "sDDD° MMm SS.S" Convert_DecDeg2RADMS = sSign & sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec1p & chrDoubleQuote Case "sDDDd MMm SSs" Convert_DecDeg2RADMS = sSign & sDeg & "d " & sMin & "m " & sSec & "s" Case Else ' graceful return on unsupported type without error message Convert_DecDeg2RADMS = "00 00 00" End Select End Function Public Function Convert_ReformatTxtCoor(ByVal sCoord_tmp As String, ByVal s_inAngFormat_tmp As String, ByVal s_outAngFormat_tmp As String, Optional s_IAU_acryon As String, Optional s_IAU_specifier As String) As String ' Reformats text coordinates from one text format to another ' E.g. - "+67° 33 22.2" to "+67 33 22" ' Used to reformat text coordinate data in astronomical catalogues ' inAngFormat and outAngFormat are the in and output string formats, ' consisting of valid types: ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' High precision format RA ' "HHMMSS.SS" ' "HH° MMm SS.SS" ' e.g. 12° 51' 23.36" ' "HH:MM:SS" ' "HHh MMm SSs" ' This format includes HH° MM' SS" ' "arc cHHMMSS.SS" ' For IAU RA formats ' "HH° MMm SSs" ' ' "HH MM.M" ' "DDD.DD" ' Leave as decimal for galactic longitude ' "sDD.DD" ' Leave as decimal for galactic latitude ' "sDDMMSS" ' "sDDMMSS.SS" ' For IAU declination formats ' "sDD MM SS" ' "sDD MM SS.S" ' High precision format Dec ' "sDD° MMm SS.S" ' -12° 51' 23.5" Use Alt+248 to generate ° symbol ' "sDDd MMm SSs" ' This format includes Dd° MM' SS" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format Dec ' "sDDD° MMm SSs" ' "sDDDd MMm SSs" ' ' Exceptions: ' It is assumed that the user will only use RA-type input formats ' with similar RA-output formats. There is no error trapping or ' logic control to prevent a user from sending an RA-type input to ' an Declination-like output format. ' ' There is no error trapping to prevent a negative Dec ' being feed to the RA hours conversion function. ' Known errors: ' When converting from "DDD.DDDD" to "DDD.DD" or "sDD.DDD" to "sDD.DD" ' there may be rounding errors of 1 second. Recommend using the ' Convert_DecDeg2RADMS function instead. E.g. - ' Convert_DecDeg2RADMS(Val("27.1139827"),"DDD.DD") ' Runtime bug for declinations between -0.5 and -0.00. ' Bug fixed 1/12/2006 ' ' Usage and testing ' FK5 data ' Name l b Epoch RA_d Dec_d RA_hms Dec_dms ' ALF HER/ 35.53/+27.82/2000/258.66/+14.39/17 14 38.86 +14 23 25.2 ' Convert_ReformatTxtCoor("258.66","DDD.DD","HH° MMm SS.SSs") yields 17º 14' 38.40" ' BET CAP/ 29.15/-26.37/2000/305.25/-14.78/20 21 00.68 -14 46 52.9 ' Convert_ReformatTxtCoor("-14.78","sDD.DD","sDD° MMm SS.S") yields -14º 46' 48.0" ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/11 03 43.67 +61 45 03.7 ' Convert_ReformatTxtCoor("165.93","DDD.DD","HHMMSS") yields 110343 ' ALF UMA/142.85/+51.01/2000/165.93/+61.75/11 03 43.67 +61 45 03.7 ' IAU RA format ' Convert_ReformatTxtCoor("165.93", "DDD.DD", "acr cHHMMSS.SS", "J") yields J110343 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_ReformatTxtCoor("27.1139827", "Dec DDD.DD", "DDD.DD") yields 027.11 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_ReformatTxtCoor("-27.1139827","sDD.DD","sDD.DD") yields -27.10 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_ReformatTxtCoor("35.62","sDD.DD","sDDMMSS") yields +353712 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' Convert_ReformatTxtCoor("35.62","sDD.DD","sDDMMSS.S") yields +353712.0 ' BET AND/127.11/-27.10/2000/ 17.43/+35.62/01 09 43.92 +35 37 14.0 ' IAU Dec format ' Convert_ReformatTxtCoor("35.62","sDD.DD","sDDMMSS.S spec","","(QSO=195)") yields +353712.0 (QSO=195) ' ALF GEM/187.44/+22.48/2000/113.65/+31.89/07 34 35.86 +31 53 17.8 ' Convert_ReformatTxtCoor("113.65","DDD.DD","HH MM SS.SS")yields 07 43 06.00 ' BET ORI/209.24/-25.25/2000/ 78.63/ -8.20/05 14 32.27 -08 12 05.9 ' Convert_ReformatTxtCoor("-8.20","sDD.DD","sDD MM SS.S") yields -08 12 00.0 ' ALF VIR/316.11/+50.84/2000/201.30/-11.16/13 25 11.58 -11 09 40.8 ' Convert_ReformatTxtCoor("201.30","DDD.DD","HHh MMm SSs") yields 13h 25m 02s ' ALF ARA/340.76/ -8.83/2000/262.96/-49.88/17 31 50.49 -49 52 34.1 ' Convert_ReformatTxtCoor("-49.88","sDD.DD","sDDDd MMm SSs") yields -049d 52m 48s ' ALF CMA/227.23/ -8.89/2000/101.29/-16.72/06 45 08.92 -16 42 58.0 ' Convert_ReformatTxtCoor(101.29,"HH:MM:SS") yields 06:45:10 ' initialize variables Const RndError As Double = 0.00000001 'Adjusts for rounding error in binary Dim chrDegree As String chrDegree = Chr(186) Dim chrSingleQuote As String chrSingleQuote = Chr(39) Dim chrDoubleQuote As String chrDoubleQuote = Chr(34) Dim s_inAngFormat As String Dim s_outAngFormat As String Dim s_DecDeg As String ' holding for string coordinate input Dim iSign As Integer 'Sign Dim iDeg As Integer 'Degrees Dim iMin As Double 'Minutes Dim dSec As Double 'Seconds Dim sSign As String 'Sign Dim sAcronym As String 'Acronym portion of IAU format Dim sDeg As String 'Degrees Dim sMin As String 'Minutes Dim sSec As String 'Seconds Dim sSec1p As String 'Seconds precision 1 digit Dim sSec2p As String 'Seconds precision 2 digit Dim A, B, C As Double 'Temporary working variables ' get function parameters s_DecDeg = sCoord_tmp s_inAngFormat = s_inAngFormat_tmp s_outAngFormat = s_outAngFormat_tmp sAcronym = s_IAU_acryon ' e.g. "QSO J" (equatorial 2000), "B" or "G" (galactic) ' translate the RA-like input format ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' "HHMMSS.SS" ' "HH° MMm SS.SSs" ' "HH:MM:SS" ' "HHh MMm SSs" ' This format includes HH° MM' SS" ' "HH° MMm SSs" ' Duplicate input for symmetry with output ' "acr HHMMSS.SS" ' IAU RA ' translate the Declination-like input format ' "DDD.DD" ' Leave as decimal for galactic longitude ' "sDD.DD" ' Leave as decimal for galactic latitude ' "sDDMMSS" ' "sDD MM SS" ' "sDD MM SS.S" ' High precision format Dec ' "sDD° MMm SS.S" ' -12° 51' 23.5" Use Alt+248 to generate ° symbol ' "sDDd MMm SSs" ' This format includes Dd° MM' SS" ' "sDDDMMSS" ' "DDMMSS.SS" ' For IAU format ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format Dec ' "sDDD° MMm SS.S" ' "sDDDd MMm SSs" ' "sDD MM" ' Select Case s_inAngFormat ' translate the RA-like input format Case "HHMMSS" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 3, 2)) dSec = CDbl(Mid(s_DecDeg, 5, 2)) Case "HH MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 2)) Case "HH MM SS.SS" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 5)) Case "HHMMSS.SS" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 3, 2)) dSec = CDbl(Mid(s_DecDeg, 5, 4)) Case "HH° MMm SS.SSs" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 5)) Case "HH:MM:SS" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 2)) Case "HHh MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) Case "HH° MMm SSs" ' duplicate format kept for symmetry without output iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) Case "arc cHHMMSS.SS" ' IAU RA format In tested iDeg = CInt(Mid(s_DecDeg, 6, 2)) iMin = CInt(Mid(s_DecDeg, 8, 2)) dSec = CDbl(Mid(s_DecDeg, 10, 2)) Case "HH MM.M" ' iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 4, 3)) dSec = 0# ' Translate the Declination-like input format Case "sDDMMSS" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 6, 2)) Case "sDD MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 8, 2)) Case "sDDMMSS.S" ' IAU dec format iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 6, 4)) Case "sDD MM SS.S" ' Tested iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 8, 4)) Case "sDD° MMm SS.S" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 4)) Case "sDDd MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) Case "sDDDMMSS" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 2)) Case "sDDD MM SS" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 2)) Case "sDDD MM SS.S" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 6, 2)) dSec = CDbl(Mid(s_DecDeg, 9, 4)) Case "sDDD° MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 7, 2)) dSec = CDbl(Mid(s_DecDeg, 11, 2)) Case "sDDDd MMm SSs" iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 7, 2)) dSec = CDbl(Mid(s_DecDeg, 11, 2)) Case "sDD MM" iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = 0# ' Translate a decimal degree input format Case "DDD.DD" ' Tested s_DecDeg = Convert_DecDeg2RADMS(Val(s_DecDeg), "HHMMSS") iDeg = CInt(Mid(s_DecDeg, 1, 2)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 6, 2)) Case "Dec DDD.DD" ' Tested s_DecDeg = Convert_DecDeg2RADMS(Val(s_DecDeg), "sDDDMMSS") iDeg = CInt(Mid(s_DecDeg, 1, 4)) iMin = CInt(Mid(s_DecDeg, 5, 2)) dSec = CDbl(Mid(s_DecDeg, 7, 2)) Case "sDD.DD" ' Tested s_DecDeg = Convert_DecDeg2RADMS(Val(s_DecDeg), "sDDMMSS") iDeg = CInt(Mid(s_DecDeg, 1, 3)) iMin = CInt(Mid(s_DecDeg, 4, 2)) dSec = CDbl(Mid(s_DecDeg, 6, 2)) Case Else ' graceful return on unsupported type without error message iDeg = 0 iMin = 0 dSec = 0 End Select ' 1/12/2006 To Correct runtime error in negative declinationss ' between -0.5 and -0.00 declination ' If Sgn(iDeg) < 0 Then ' iSign = -1 ' iDeg = Abs(iDeg) 'Else ' iSign = 1 ' End If If Mid(s_DecDeg, 1, 1) = "-" Then iSign = -1 iDeg = Abs(iDeg) Else iSign = 1 End If ' Format the output string ' preformat with leading zeros If iSign < 0# Then sSign = "-" Else sSign = "+" ' if output format is three position decimal angle, ' change leading zero to three positions If InStr(1, s_outAngFormat, "DDD", vbTextCompare) > 0 Then sDeg = CStr(Format(Int(iDeg), "000")) Else sDeg = CStr(Format(Int(iDeg), "00")) End If sMin = CStr(Format(Int(iMin), "00")) sSec = CStr(Format(Round(dSec), "00")) sSec1p = CStr(Format(dSec, "00.0")) sSec2p = CStr(Format(dSec, "00.00")) ' translate the output format ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' High precision format RA ' "HHMMSS.SS" ' "HH° MMm SS.SSs" ' e.g. 12° 51' 23.36" ' "HH:MM:SS" ' "HHh MMm SSs" ' ' "HH° MMm SSs" ' e.g. 12° 51' 23" ' "acr HHMMSS.SS" ' IAU RA ' translate the Declination-like output format ' "DDD.DD" ' Leave as decimal for galactic longitude ' "sDD.DD" ' Leave as decimal for galactic latitude ' "sDDMMSS" ' "sDDMMSS.SS" ' IAU declination ' "sDD MM SS" ' "sDD MM SS.S" ' High precision format Dec ' "sDD° MMm SS.S" ' e.g. -12° 51' 23.5" ' "sDDd MMm SSs" ' "sDDDMMSS" ' "sDDD MM SS" ' "sDDD MM SS.S" ' High precision format Dec ' "sDDD° MMm SS.S" ' e.g. -121° 51' 23.5" ' "sDDDd MMm SSs" Select Case s_outAngFormat ' translate the RA-like output format Case "HHMMSS" Convert_ReformatTxtCoor = sDeg & sMin & sSec Case "HH MM SS" Convert_ReformatTxtCoor = sDeg & " " & sMin & " " & sSec Case "HH MM SS.SS" Convert_ReformatTxtCoor = sDeg & " " & sMin & " " & sSec2p Case "HHMMSS.SS" Convert_ReformatTxtCoor = sDeg & sMin & sSec2p Case "HH° MMm SS.SSs" ' Tested Convert_ReformatTxtCoor = sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec2p & chrDoubleQuote Case "HH:MM:SS" Convert_ReformatTxtCoor = sDeg & ":" & sMin & ":" & sSec Case "HHh MMm SSs" Convert_ReformatTxtCoor = sDeg & "h " & sMin & "m " & sSec & "s" Case "HH° MMm SSs" Convert_ReformatTxtCoor = sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec & chrDoubleQuote Case "acr cHHMMSS.SS" ' IAU RA format tested Convert_ReformatTxtCoor = sAcronym & sDeg & sMin & sSec2p ' translate the Declination-like output format Case "DDD.DD" ' tested Convert_ReformatTxtCoor = CStr(Format(Val(Convert_DMS2DecDeg(sDeg & sMin & sSec, "sDDMMSS", "string")), "000.00")) Case "sDD.DD" ' ' tested Convert_ReformatTxtCoor = sSign & CStr(Format(Val(Convert_DMS2DecDeg("+" & sDeg & sMin & sSec, "sDDMMSS", "string")), "00.00")) Case "sDDMMSS" Convert_ReformatTxtCoor = sSign & sDeg & sMin & sSec Case "sDD MM SS" ' Tested Convert_ReformatTxtCoor = sSign & sDeg & " " & sMin & " " & sSec Case "sDD MM SS.S" Convert_ReformatTxtCoor = sSign & sDeg & " " & sMin & " " & sSec1p Case "sDDMMSS.S" ' IAU dec format without specifier Convert_ReformatTxtCoor = sSign & sDeg & sMin & sSec1p Case "sDDMMSS.S spec" ' IAU dec format with specifier Convert_ReformatTxtCoor = sSign & sDeg & sMin & sSec1p & " " & s_IAU_specifier Case "sDD° MMm SS.S" Convert_ReformatTxtCoor = sSign & sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec1p & chrDoubleQuote Case "sDDd MMm SSs" Convert_ReformatTxtCoor = sSign & sDeg & "d " & sMin & "m " & sSec & "s" Case "sDDDMMSS" Convert_ReformatTxtCoor = sSign & sDeg & sMin & sSec Case "sDDD MM SS" Convert_ReformatTxtCoor = sSign & sDeg & " " & sMin & " " & sSec Case "sDDD MM SS.S" Convert_ReformatTxtCoor = sSign & sDeg & " " & sMin & " " & sSec1p Case "sDDD° MMm SS.S" Convert_ReformatTxtCoor = sSign & sDeg & chrDegree & " " & sMin & chrSingleQuote & " " & sSec1p & chrDoubleQuote Case "sDDDd MMm SSs" Convert_ReformatTxtCoor = sSign & sDeg & "d " & sMin & "m " & sSec & "s" Case Else ' graceful return on unsupported type without error message Convert_ReformatTxtCoor = "00 00 00" End Select End Function Function DecDeg2arcmins(ByVal Deg_dec_tmp As Double) As Double ' Utility to quickly convert the units of the diameter of an extended object ' from degrees to arcmins ' Exceptions: No error trapping for negative values. Dim dDeg As Double 'Degrees Dim dMin As Double 'Minutes Dim iSign As Integer Dim A As Double 'Temporary working variables Const RndError As Double = 0.00000001 'Adjusts for rounding error in binary ' get function parameters dDeg = Deg_dec_tmp If dDeg < 0 Then iSign = -1 dDeg = Abs(dDeg) Else iSign = 1 End If A = dDeg ' Temporary degrees storage ' Determine Hours-Deg, Min, Sec dDeg = Int(A + RndError) 'Get degrees or hours dMin = (A - dDeg) * 60# 'Get minutes & seconds DecDeg2arcmins = Round(((dDeg * 60) + dMin) * iSign, 0) End Function Function DecDeg2arcsecs(ByVal Deg_dec_tmp As Double) As Double ' Utility to quickly convert the units of the diameter of an extended object ' from degrees to arcseconds ' Exceptions: No error trapping for negative values. Dim dDeg As Double 'Degrees Dim dMin As Double 'Minutes Dim dSec As Double 'Seconds Dim iSign As Integer Dim A As Double 'Temporary working variables Const RndError As Double = 0.00000001 'Adjusts for rounding error in binary ' get function parameters dDeg = Deg_dec_tmp If dDeg < 0 Then iSign = -1 dDeg = Abs(dDeg) Else iSign = 1 End If A = dDeg ' Temporary degrees storage ' Determine Hours-Deg, Min, Sec dDeg = Int(A + RndError) 'Get degrees or hours dMin = (A - dDeg) * 60# 'Get minutes & seconds dSec = (dMin - Int(dMin + RndError)) * 60# 'Get seconds dMin = Int(dMin + RndError) 'Get minutes DecDeg2arcsecs = Round(((dDeg * 3600) + (dMin * 60) + dSec) * iSign, 0) End Function Function DecArcsecs2DecDeg(ByVal Arcsecs_tmp As Double) As Double ' Utility to quickly convert declination arcsecs to decimal degrees ' There are 3600 arcsecs in a declination degree DecArcsecs2DecDeg = Arcsecs_tmp / 3600# End Function Function RAArcsecs2DecDeg(ByVal Arcsecs_tmp As Double) As Double ' Utility to quickly convert right ascension arcsecs to decimal degrees ' There are 3600 arcsecs in a degree. ' There are 24 hours of right ascension in 360 degrees or 1 hour = 15 degrees ' Circumference = 360 degrees * 60 minutes * 60 seconds declination ' Circumference = 24 hours * 60 minutes * 60 seconds right ascension ' Ratio is 15:1 ' therefore, 1 minute right ascension = 15 minutes declination or decimal arcmins ' and 1 arcsec right ascension = 15 arcsecs declination or decimal arcsecs RAArcsecs2DecDeg = Arcsecs_tmp * 15# / (3600#) End Function Public Function Rad2arcsec(ByVal Rads As Double) As Double ' Receives radians; returns arcsecs ' A_arcsec= A_radians * 360/2pi() degs per rad * ( 60 * 60 ) arcsecs per degree ' The constants reduce to 206265 Rad2arcsec = 206265 * Rads End Function Function Parallax2parsecs(ByVal Plex_millarcsec As Double) As Variant ' Converts parallax measurements in millarcsecs to parsecs ' Used in measuring distances to intra-galactic objects within 800 parsecs ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 396 Parallax2parsecs = 1# / (Plex_millarcsec / 1000) End Function Function Parallax2light_yrs(ByVal Plex_millarcsec As Double) As Variant ' Converts parallax measurements in millarcsecs to parsecs ' Used in measuring distances to intra-galactic objects within 800 parsecs ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 396 Parallax2light_yrs = 3.2616 / (Plex_millarcsec / 1000) End Function Function Parsecs2light_yrs(ByVal Parsecs As Double) As Variant ' Converts parallax measurements in millarcsecs to parsecs ' Used in measuring distances to intra-galactic objects within 800 parsecs ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 396 Parsecs2light_yrs = Parsecs / 3.2616 End Function Function Light_yrs2Parsecs(ByVal Light_yrs As Double) As Variant ' Converts light years to parsecs ' Used in measuring distances to intra-galactic objects within 800 parsecs ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 396 Light_yrs2Parsecs = Light_yrs * 3.2616 End Function Function RVel2Mparsecs(ByVal RVel As Double) As Variant ' Converts redshifts to distance in kiloparsecs ' using Hubble's constant to determine distance of ' extra-galactic objects ' v = Hd ' D = v / H ' Hubble Constant per deVaucouleurs 3rd Ref. Cat. of Bright Galaxies ' 100 km/s/ MPC ' Hubble Constant per Allan Sandage of the Carnegie Institutions ' 50 km/s / MPC ' Hubble Constant per the HST Project Team 2001ApJ 553, 47F ' 72 km/s -+ 8 / MPC RVel2Mparsecs = Round((RVel / 72), 4) End Function Public Function DistHubbleNowComoving(ByVal z_tmp As Double, ByVal H_0_tmp As Double) As String ' Returns distance in Mpc and billion light-years from z ' Working variables Const C As Double = 299792 ' kilometers / second Dim z, v, H_0, D_mpc, D_blyrs As Double ' Get input z = z_tmp H_0 = H_0_tmp ' Compute result v = C * z D_mpc = v / H_0 D_blyrs = (D_mpc / 1000#) * 3.26 ' Return result DistHubbleNowComoving = CStr(Trim(D_mpc)) & " " & CStr(Trim(D_blyrs)) End Function Public Function DistKpc2PowerofTen(ByVal Dist_kpc_tmp As Double, ByVal Round_int_tmp As Integer) As String ' Converts kiloparsec to a Power of Ten Analogy ' Returns multiplier and log10 exponent ' Source: Ned Wright's Cosmology Calculator gives the following constants ' 1 Gly = 1,000,000,000 light years or 9.461*10^26 cm. ' 1 Gyr = 1,000,000,000 years. ' 1 Mpc = 1,000,000 parsecs = 3.08568*10^24 cm, or 3,261,566 light years. ' Working Variables Dim Dist_kpc, DistanceMultipler_cm, Distance_Log10Exp_cm, Round_int As Double Dim A, B, C, D As Double Dim S As String Dim k As Double ' Constants k = 3.08568 * 10 ^ 21 ' parsecs to cm ' Get input Dist_kpc = Dist_kpc_tmp Round_int = Round_int_tmp ' Compute ' Find distance in cm A = Dist_kpc * k ' Find the exponent B = Log10(A) C = CDbl(Int(A)) Distance_Log10Exp_cm = CDbl(Right(Trim(CStr(C)), 2)) ' Find the multiplier DistanceMultipler_cm = A / (10 ^ Distance_Log10Exp_cm) ' Round mulitplier DistanceMultipler_cm = NumTrunc(DistanceMultipler_cm, Round_int) ' Return result DistKpc2PowerofTen = Trim(CStr(A)) & " " & Trim(CStr(DistanceMultipler_cm)) & " " & Trim(CStr(Distance_Log10Exp_cm)) End Function Public Function phyDistModulus2Parsecs(ByVal Apparent_mag_tmp As Double, ByVal Absolute_mag_tmp As Double) As Double ' Returns a distance estimate in parsecs from the distance modulus. ' Source: Inglis, M. 2003. Observer's Guide to Stellar Evolution: The Birth, Life, and Death of Stars. Springer. ISBN 1-85233-465-7 ' At p. 16 ' log(d) = (appMag -absMag + 5 ) / 5 ' Test ' Inglis example at p.16; Abs_mag = 6; App_mag=16 answer = 1000 parsecs ' ? DistModulus2Parsecs(16,6) ' 1000 ' Working variables Dim App_mag, Abs_mag, A, B, C As Double ' Get input App_mag = Apparent_mag_tmp Abs_mag = Absolute_mag_tmp ' Find the log of the distance A = (App_mag - Abs_mag + 5) / 5 phyDistModulus2Parsecs = 10 ^ A End Function Public Function phtJohnson2mv(ByVal dblB As Double, ByVal dbl_uncB As Double, ByVal dblV As Double, ByVal dbl_uncV As Double, ByVal intPrecision As Integer) As String ' Converts Johnson UBV photometry B-V color magnitudes to the ' visual magnitude system with uncertainty handling ' Then returns the Johnson Color Index (B-V) ' Based on Stanton, R. 1999. Visual Magnitudes and the Average Observer: ' the SS Cygnii Field Experiment. JAAVSO 27:97. ' ' m_v = V + 0.210(B-V) ' ' where ' m=v is the visual magnitude system ' B is the Johnson UBV B filter magnitude ' V is the Johnson UBV B filter magnitude ' and CI (color index) = B-V ' ' m_v observations have an inherent uncertainty of 0.2 mags; not 0.1 mags, for ' experienced observers and 0.3 mags for inexperienced observers. ' ' This implementation defaults to 0.2 mag uncertainty of observation. ' ' Inputs: ' B, uncertainty of B, V, uncertainty of V from Johnson UBV ' data obtained, usually, by photoelectric means. ' Precision formatting for result. ' Outputs: ' Returns a string with the m_v (magnitude in visual system) and ' the +- uncertainty as positive double ' Useage ' ? phtJohnson2mv(1.3,0.1,2.3,0.1,1) yields "2 .5" Dim A As Double 'working variables Dim C As Double 'working variables Dim D As Double 'working variables Dim E As Double 'working variables Dim f As Double 'working variables Dim g As Double 'working variables Dim B As Double Dim v As Double Dim uncB As Double Dim uncV As Double ' get variables B = CDbl(dblB) ' syntax - coerce to double to prevent runtime error v = CDbl(dblV) uncB = Abs(dbl_uncB) uncV = Abs(dbl_uncV) ' compute the magnitude results C = v + 0.21 * (B - v) ' compute the uncertainty ' uncertainty of first part 0.21*(B - V) D = UncertDif(B, uncB, v, uncV, "uncertainty", intPrecision) E = 0.21 * (B - v) ' uncertainty of second part V + (0.21 * (B - V)) f = UncertSum(v, uncV, E, D, "uncertainty", intPrecision) ' add effect of 0.2 mag uncertainty of average visual observer ' use )+m_v with uncertainty. g = UncertSum(C, f, 1, 0.2, "uncertainty", intPrecision) ' Return the result phtJohnson2mv = Str(NumTrunc(C, intPrecision)) & " " & Str(NumTrunc(g, intPrecision)) End Function Public Function phtPhotoElect2mv(ByVal dblB As Double, ByVal dbl_uncB As Double, ByVal dblV As Double, ByVal dbl_uncV As Double, ByVal intPrecision As Integer) As String ' Converts generic photoelectric UBV photometry B-V color magnitudes to the ' visual magnitude system with uncertainty handling ' This is a generic equation, although each study, photoelectric detector or catalogue ' may provide a specific conversion equation. ' Photoelectric magnitude measurements have the lowest standard errors (<0.10). Photographic ' magnitude estimates have higher error rates (>0.3) ' Based on against tests on the US Naval Obs. 1974 catalogue: ' Stanton, R. 1981. Photoelectric Measuraes of AAVSO Comparison Star Sequences - II. ' JAAVSO 10:1S. ' Stanton, R. 1999. Visual Magnitudes and the Average Observer: ' the SS Cygnii Field Experiment. JAAVSO 27:97. ' ' m_v = V_e + 0.182*(B_e-V_e)-0.15 ' ' where ' m=v is the visual magnitude system ' B is the photoelectric B filter magnitude ' V is the photoelectric V filter magnitude ' and CI (color index) = B-V ' ' m_v observations have an inherent uncertainty of 0.2 mags; not 0.1 mags, for ' experienced observers and 0.3 mags for inexperienced observers. ' ' This implementation defaults to 0.2 mag uncertainty of observation. ' ' Inputs: ' B, uncertainty of B, V, uncertainty of V from photoelectric UBV ' data obtained, usually, by photoelectric means. ' Precision formatting for result. ' Outputs: ' Returns a string with the m_v (magnitude in visual system) and ' the +- uncertainty as positive double ' Useage ' ? phtPhotoElect2mv(1.3,0.01,2.3,0.01,1) yields "1.9 .2" Dim A As Double 'working variables Dim C As Double 'working variables Dim D As Double 'working variables Dim E As Double 'working variables Dim f As Double 'working variables Dim g As Double 'working variables Dim B As Double Dim v As Double Dim uncB As Double Dim uncV As Double ' get variables B = CDbl(dblB) ' syntax - coerce to double to prevent runtime error v = CDbl(dblV) uncB = Abs(dbl_uncB) uncV = Abs(dbl_uncV) ' compute the magnitude results C = v + (0.182 * (B - v)) - 0.15 ' compute the uncertainty ' uncertainty of first part 0.182*(B - V) D = UncertDif(B, uncB, v, uncV, "uncertainty", intPrecision) E = (0.182 * (B - v)) - 0.15 ' uncertainty of second part V + (0.182 * (B - V) - 0.15) f = UncertSum(v, uncV, E, D, "uncertainty", intPrecision) ' add effect of 0.2 mag uncertainty of average visual observer ' use )+m_v with uncertainty. g = UncertSum(C, f, 1, 0.2, "uncertainty", intPrecision) ' Return the result phtPhotoElect2mv = Str(NumTrunc(C, intPrecision)) & " " & Str(NumTrunc(g, intPrecision)) End Function Public Function phtTycho2mv(ByVal dblB As Double, ByVal dbl_uncB As Double, ByVal dblV As Double, ByVal dbl_uncV As Double, ByVal intPrecision As Integer) As String ' Converts Tycho UBV photometry B-V color magnitudes to the ' visual magnitude system with uncertainty handling ' Based on: ' Stanton, R. 1999. Visual Magnitudes and the Average Observer: ' the SS Cygnii Field Experiment. JAAVSO 27:97. ' Bessell, M. July, 2000. The Hipparcos and Tycho Photometric System Passbands. ' PASP 112:961. ' Hog E. et al. 2000. The Tycho-2 Catalogue of the 2.5 Million Brightest Stars ' Astron. Astrophys. 355, L27 =2000A&A...355L..27H ' The Tycho-2 catalogue notes that magnitude standard error is: ' Photometric std. errors (3) on VT ' VT < 9 mag 0.013 mag ' all stars 0.10 mag ' Note 7 to the Tycho-2 readme file states: ' Blank when no magnitude is available. Either BTmag or VTmag is ' always given. Approximate Johnson photometry may be obtained as: ' V = VT - 0.09 * (BT - VT) ' B - V = 0.85 * (BT - VT) ' Consult Sect 1.3 of Vol 1 of "The Hipparcos and Tycho Catalogues", ' ESA SP-1200, 1997, for details. ' Conversion equation from Stanton 1999: ' m_v = V_t + 0.089(B_t-V_t) ' ' where ' m=v is the visual magnitude system ' B is the Tycho-2 UBV B filter magnitude ' V is the Tycho-2 UBV B filter magnitude ' and CI (color index) = B-V ' ' m_v observations have an inherent uncertainty of 0.2 mags; not 0.1 mags, for ' experienced observers and 0.3 mags for inexperienced observers. ' ' This implementation defaults to 0.2 mag uncertainty of observation. ' ' Inputs: ' B, uncertainty of B, V, uncertainty of V from Johnson UBV ' data obtained, usually, by photoelectric means. ' Precision formatting for result. ' Outputs: ' Returns a string with the m_v (magnitude in visual system) and ' the +- uncertainty as positive double ' Useage ' ? phtTycho2mv(1.3,0.1,2.3,0.1,1) yields "2.2 .5" Dim A As Double 'working variables Dim C As Double 'working variables Dim D As Double 'working variables Dim E As Double 'working variables Dim f As Double 'working variables Dim g As Double 'working variables Dim B As Double Dim v As Double Dim uncB As Double Dim uncV As Double ' get variables B = CDbl(dblB) ' syntax - coerce to double to prevent runtime error v = CDbl(dblV) uncB = Abs(dbl_uncB) uncV = Abs(dbl_uncV) ' compute the magnitude results C = v + 0.089 * (B - v) ' compute the uncertainty ' uncertainty of first part 0.21*(B - V) D = UncertDif(B, uncB, v, uncV, "uncertainty", intPrecision) E = 0.089 * (B - v) ' uncertainty of second part V + (0.21 * (B - V)) f = UncertSum(v, uncV, E, D, "uncertainty", intPrecision) ' add effect of 0.2 mag uncertainty of average visual observer ' use )+m_v with uncertainty. g = UncertSum(C, f, 1, 0.2, "uncertainty", intPrecision) ' Return the result phtTycho2mv = Str(NumTrunc(C, intPrecision)) & " " & Str(NumTrunc(g, intPrecision)) End Function Public Function phtTycho2Johnson(ByVal dblB As Double, ByVal dbl_uncB As Double, ByVal dblV As Double, ByVal dbl_uncV As Double, ByVal intPrecision As Integer) As String ' Converts Tycho UBV photometry B-V color magnitudes to the ' Johnson UBV with uncertainty handling ' Based on: ' Bessell, M. July, 2000. The Hipparcos and Tycho Photometric System Passbands. ' PASP 112:961. ' Hog E. et al. 2000. The Tycho-2 Catalogue of the 2.5 Million Brightest Stars ' Astron. Astrophys. 355, L27 =2000A&A...355L..27H ' The Tycho-2 catalogue notes that magnitude standard error is: ' Photometric std. errors (3) on VT ' VT < 9 mag 0.013 mag ' all stars 0.10 mag ' Note 7 to the Tycho-2 readme file states: ' Blank when no magnitude is available. Either BTmag or VTmag is ' always given. Approximate Johnson photometry may be obtained as: ' V = VT - 0.09 * (BT - VT) ' B - V = 0.85 * (BT - VT) ' Consult Sect 1.3 of Vol 1 of "The Hipparcos and Tycho Catalogues", ' ESA SP-1200, 1997, for details. ' Conversion equation from Tycho-2 catalogue: ' m_v = V_t + 0.09(B_t-V_t) ' ' where ' m=v is the visual magnitude system ' B is the Tycho-2 UBV B filter magnitude ' V is the Tycho-2 UBV B filter magnitude ' and CI (color index) = B-V ' ' m_v observations have an inherent uncertainty of 0.2 mags; not 0.1 mags, for ' experienced observers and 0.3 mags for inexperienced observers. ' ' This implementation defaults to 0.2 mag uncertainty of observation. ' ' Inputs: ' B_Tycho, uncertainty of B_Tycho, V_Tycho, uncertainty of V_Tycho ' data obtained, usually, by photoelectric means. ' Precision formatting for result. ' Outputs: ' Returns a string with the Johnson V (magnitude in Johnson system) and ' the +- uncertainty as positive double ' Useage ' ? phtTycho2Johnson(1.3,0.1,2.3,0.1,1) yields "2.2 .5" Dim A As Double 'working variables Dim C As Double 'working variables Dim D As Double 'working variables Dim E As Double 'working variables Dim f As Double 'working variables Dim g As Double 'working variables Dim B As Double Dim v As Double Dim uncB As Double Dim uncV As Double ' get variables B = CDbl(dblB) ' syntax - coerce to double to prevent runtime error v = CDbl(dblV) uncB = Abs(dbl_uncB) uncV = Abs(dbl_uncV) ' compute the magnitude results C = v + 0.09 * (B - v) ' compute the uncertainty ' uncertainty of first part 0.09*(B - V) D = UncertDif(B, uncB, v, uncV, "uncertainty", intPrecision) E = 0.09 * (B - v) ' uncertainty of second part V - (0.09 * (B - V)) ' Note - no uncertainty change for constant multiplication f = UncertDif(v, uncV, E, D, "uncertainty", intPrecision) ' add effect of 0.2 mag uncertainty of average visual observer ' use )+m_v with uncertainty. ' G = UncertSum(C, F, 1, 0.2, "uncertainty", 2) ' Return the result phtTycho2Johnson = Str(NumTrunc(C, intPrecision)) & " " & Str(NumTrunc(f, intPrecision)) End Function Public Function phtTycho2JohnsonCI(ByVal dblB As Double, ByVal dbl_uncB As Double, ByVal dblV As Double, ByVal dbl_uncV As Double, ByVal intPrecision As Integer) As String ' Converts Tycho UBV photometry B-V color magnitudes to the ' Johnson UBV with uncertainty handling, but returns the ' Johnson CI (B-V) instead of Johnson V ' Based on: ' Bessell, M. July, 2000. The Hipparcos and Tycho Photometric System Passbands. ' PASP 112:961. ' Hog E. et al. 2000. The Tycho-2 Catalogue of the 2.5 Million Brightest Stars ' Astron. Astrophys. 355, L27 =2000A&A...355L..27H ' The Tycho-2 catalogue notes that magnitude standard error is: ' Photometric std. errors (3) on VT ' VT < 9 mag 0.013 mag ' all stars 0.10 mag ' Note 7 to the Tycho-2 readme file states: ' Blank when no magnitude is available. Either BTmag or VTmag is ' always given. Approximate Johnson photometry may be obtained as: ' V = VT - 0.09 * (BT - VT) ' B - V = 0.85 * (BT - VT) ' Consult Sect 1.3 of Vol 1 of "The Hipparcos and Tycho Catalogues", ' ESA SP-1200, 1997, for details. ' Conversion equation from Tycho-2 catalogue: ' m_v = V_t + 0.09(B_t-V_t) ' m_ci = 0.85 * (B_t - V_t) ' ' where ' m=v is the visual magnitude system ' B is the Tycho-2 UBV B filter magnitude ' V is the Tycho-2 UBV B filter magnitude ' and CI (color index) = B-V ' ' m_v observations have an inherent uncertainty of 0.2 mags; not 0.1 mags, for ' experienced observers and 0.3 mags for inexperienced observers. ' ' This implementation defaults to 0.2 mag uncertainty of observation. ' ' Inputs: ' B_Tycho, uncertainty of B_Tycho, V_Tycho, uncertainty of V_Tycho ' data obtained, usually, by photoelectric means. ' Precision formatting for result. ' Outputs: ' Returns a string with the Johnson CI (magnitude in Johnson system) and ' the +- uncertainty as positive double ' Useage ' ? phtTycho2JohnsonCI(1.3,0.1,2.3,0.1,1) yields "-0.9 .2" Dim A As Double 'working variables Dim C As Double 'working variables Dim D As Double 'working variables Dim E As Double 'working variables Dim f As Double 'working variables Dim g As Double 'working variables Dim B As Double Dim v As Double Dim uncB As Double Dim uncV As Double ' get variables B = CDbl(dblB) ' syntax - coerce to double to prevent runtime error v = CDbl(dblV) uncB = Abs(dbl_uncB) uncV = Abs(dbl_uncV) ' compute the magnitude results ' m_ci = 0.85 * (B_t - V_t) C = 0.85 * (B - v) ' compute the uncertainty ' uncertainty of first part 0.85 * (B - V) ' Note - No reduction uncertainty for multiplication by constant D = UncertDif(B, uncB, v, uncV, "uncertainty", intPrecision) ' Return the result phtTycho2JohnsonCI = Str(NumTrunc(C, intPrecision)) & " " & Str(NumTrunc(D, intPrecision)) End Function Public Function TimeJulianDay(ByVal Year_tmp As Integer, ByVal Month_tmp As Integer, ByVal Day_tmp As Integer, ByVal UTC_Hour_tmp As Integer, _ ByVal UTC_Min_tmp As Integer, ByVal UTC_Sec_tmp As Double, ByVal IsGregorian_tmp As Integer) As String ' Returns julian day number for a gregorian calender date and UTC time, if IsGregorian=1. ' Returns julian day number for a julian day and UTC time, if IsGregorian = 0. ' Gregorian calendar input date is the default. ' Returns string containing the decimal Julian day and an IAU formatted Julian decimal day. ' Future enchancements ' Inverse function JD decimal day to Gregorian Days HH:MM:SS.SS ' Need HH:MM:SS.SS to decimal day ' Usage and test data: ' ' Example 1: ' Meeus example at p. 61: ' Date Sputnik launched 1957 Oct. 4.81 ' ? TimeJulianDay(1957,10,4,19,30,0,1) yields Julian Day 2436116.3125 ' Example 2: ' Meeus example at p. 61: ' January 27, 333 at noon. ' ? TimeJulianDay(333,1,27,12,0,0,0) yields Julian Day 1842713 ' Algorithm sources: ' Julian Day: Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 61. ' See ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. ' ' Author: Kurt A. Fisher 2004 ' Initialize variables Dim A As Double ' Meeus intermediate variable Dim B As Double ' Meeus intermediate variable Dim intY, intM, intD, intH, intMin As Integer Dim bolGregorian As Boolean Dim dblSec, dblDay, dblJD As Double ' Working variables Dim S As String ' Output string Dim varMonthNames(12) As String varMonthNames(1) = "Jan" varMonthNames(2) = "Feb" varMonthNames(3) = "Mar" varMonthNames(4) = "Apr" varMonthNames(5) = "May" varMonthNames(6) = "Jun" varMonthNames(7) = "Jul" varMonthNames(8) = "Aug" varMonthNames(9) = "Sep" varMonthNames(10) = "Oct" varMonthNames(11) = "Nov" varMonthNames(12) = "Dec" ' Transfer interface input ' Get M/D/Y UTC-Hour/UTC-Min/UTC-Sec intY = Year_tmp intM = Month_tmp intD = Day_tmp intH = UTC_Hour_tmp intMin = UTC_Min_tmp dblSec = UTC_Sec_tmp Select Case IsGregorian_tmp Case 0 bolGregorian = False Case 1 bolGregorian = True Case Else bolGregorian = True End Select ' Compute decimal_day of the month from time, including hours and secs. dblDay = CDbl(intD) + ((CDbl((intH * 60 * 60)) + CDbl((intMin * 60)) + dblSec) / (24# * 60# * 60#)) ' Build the IAU formatted date S = CStr(intY) & varMonthNames(intM) & Format(dblDay, "00.00000") ' Evaluate M - Meeus p. 61 Select Case intM Case 1 intY = intY - 1 intM = intM + 12 Case 2 intY = intY - 1 intM = intM + 12 Case Is > 2 intY = intY intM = intM Case Else ' Assume > 2 intY = intY intM = intM End Select ' Determine A - Meeus p. 61 A = Int(intY / 100#) ' Determine B - Meeus p. 61 Select Case bolGregorian Case True B = 2# - A + Int(A / 4#) Case False B = 0# Case Else ' Assume is Gregorian calendar date B = 2# - A + Int(A / 4#) End Select ' Determine Julian Day Meeus p. 61 dblJD = Int(365.25 * (CDbl(intY) + 4716#)) + Int(30.6001 * (CDbl(intM) + 1)) + dblDay + B - 1524.5 ' Format output as Julian decimal day plus IAU formatted string S = CStr(dblJD) & " " & S ' Return result TimeJulianDay = S End Function Public Function TimeLocalSiderealST(ByVal strUTCtmp As String, ByVal dayDate As Date, ByVal decLongitude As Double) As String ' Computes mean, not apparent, sidereal time using a low precision formula from ' ______. Oct. 1993. Finding Sidereal Time. Sky & Telescope. p. 72. ' Inputs: 24 hour clock time as Universal Time ' Date as date ' West longitude as negative number ' Returns as string the RA as string and decimal degrees. ' ST article claims accuracy within 0.04 hours or two minutes ' Uses Convert_DecDeg2RADMS ' For the Meeus example problem, this alogrithm yields RA 08 35 16.98. ' The Meeus higher precision formula yields RA 08 34 57.09. ' ' This ST algorithm uses the west longitudSe only. ' Example: Find LST at Greenwich on 4/10/1987 at UT 19:21:00. ' ? TimeLocalSiderealST("19:21:00",#04/10/1987#,0) ' 083516.98 ' Find LST at SLC, UT 111 West longitude on 9/24/2005 at UT 05:59:59. ' TimeLocalSiderealST("05:59:59", #9/24/2005#,-111) ' 224610.64 22.7696213888889 Dim D, t, L, LST As Double Dim strLST As String ' Convert the time value to decimal hours t = CDbl(DatePart("h", TimeValue(strUTCtmp))) t = t + (CDbl((DatePart("n", TimeValue(strUTCtmp)))) / 60#) t = t + (CDbl((DatePart("s", TimeValue(strUTCtmp)))) / 3600#) ' Get longitude ' Convert to positive for West longitude L = -1# * CDbl(decLongitude) ' Get the Day of Year D = CDbl(DatePart("y", dayDate)) ' Compute local sidereal time LST = 6.61 + (1.003 * t) + (0.0657 * D) - (L / 15#) If LST >= 24# Then LST = LST - 24# End If ' Need to convert LST to a string in format hh mm ss.s LST = LST * 15 strLST = Convert_DecDeg2RADMS(LST, "HHMMSS.SS") & " " & CStr(LST) ' Return result TimeLocalSiderealST = strLST End Function Public Function TimeLocalSidereal(ByVal strTimetmp As String, ByVal dayDate As Date, decLongitude As Double) As String ' Computes mean, not apparent, sidereal time using a higher precision formula from ' Modified from Meeus. 1998. Sidereal Time at Greenwich. Chap. 12. In ' Astronomical Algorithms. ' Inputs: 24 hour clock time as Universal Time e.g. "18:21:45" or "18:21:45.67" ' Date as date ' Longitude as number - West is negative ' Returns as string the RA and Sidereal Hour Angle (SHA). ' Uses Meuss T factor at page 88. ' Uses Convert_DecDeg2RADMS and TimeJulianDay ' ' Examples ' Meeus - this implementation gives the same result as Meeus example problems: ' Find LST at Greenwich on 4/10/1987 at UT 00:00:00. ' ? TimeLocalSidereal("00:00:00",#04/10/1987#,0) ' 131046.37 197.693195 ' Find LST at Greenwich on 4/10/1987 at UT 19:21:00. ' ? TimeLocalSidereal("19:21:00",#04/10/1987#,0) ' 083457.09 128.737873 ' Peter Duffet-Smith's Practical Astronomy with Your Calculator, Sec. 12, p. 17: ' TimeLocalSidereal("14:36:51.67", #4/22/1980#, 0) ' 044005.23 70.021788 ' Find LST at SLC, UT 111 West longitude on 9/24/2005 at UT 05:59:59. ' TimeLocalSidereal("05:59:59", #9/24/2005#, -111) ' 224840.79 342.169968 Dim intDay, intMonth, intYear, intHour, intMinute As Integer Dim dblSec, L As Double Dim intIsGregorian As Integer Dim LST, LST_dec, LST_dec_adj As Double ' Local sidereal time in degrees Dim strLST As String ' Output build string Dim JD_c As Double ' Julian Day - computed Dim t As Double Dim JD_0h As Double Dim strTimetmp_adj As String ' Get input; set variables ' strTimetmp_adj special handles time formats with decimal seconds, e.g. "12:29:19.87" strTimetmp_adj = Left(strTimetmp, 8) intYear = CInt(DatePart("yyyy", dayDate)) intMonth = CInt(DatePart("m", dayDate)) intDay = CInt(DatePart("d", dayDate)) intHour = CInt(DatePart("h", TimeValue(strTimetmp_adj))) intMinute = CInt(DatePart("n", TimeValue(strTimetmp_adj))) dblSec = CDbl(Right(strTimetmp, Len(strTimetmp) - 6)) ' Note seconds are a double, not an integer ' CDbl(DatePart("s", TimeValue(strTimetmp))) L = CDbl(decLongitude) intIsGregorian = 1 ' Date is gregorian, not Julian ' Get the Julian Day ' at the moment of UTC JD_c = CDbl(GetCSWord((TimeJulianDay(intYear, intMonth, intDay, intHour, intMinute, dblSec, intIsGregorian)), 1)) ' at the zero hour JD_0h = CDbl(GetCSWord((TimeJulianDay(intYear, intMonth, intDay, 0, 0, 0, intIsGregorian)), 1)) ' Compute T factor t = (JD_0h - 2451545#) / 36525# ' Compute the instanteous local mean sidereal time ' Meeus formula 12.4 - yields instanteous LST in degrees - not just LST at 0h UTC LST_dec = 280.46061837 + (360.98564736629 * (JD_c - 2451545#)) + (0.000387933 * (t ^ 2)) - ((t ^ 3) / 38710000#) ' Typically yields a large number of degrees, e.g. 1677831.19946226 ' Parse out the fractional part of the neartest 360 degree circle for the current day sidereal time ' Use iteration method If LST_dec >= 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop End If If LST_dec < 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj <= 360# LST_dec_adj = LST_dec_adj + 360# Loop LST_dec_adj = LST_dec_adj - 360# End If ' Add time for the geographic longitude LST_dec_adj = LST_dec_adj + (L) ' LST_dec_adj is in degrees already ' Adjust for 24 hours rounding If CDbl(LST_dec_adj) > 360# Then LST_dec_adj = LST_dec_adj - 360# End If If CDbl(LST_dec_adj) = 360# Then LST_dec_adj = 0# End If ' Note: This error check may be redunant or deprecated If CDbl(LST_dec_adj) < 0# Then LST_dec_adj = LST_dec_adj + 360# End If ' Trap runtime error where LST is exactly 360 If CStr(LST_dec_adj) = "360" Then LST_dec_adj = 0# End If ' Trap runtime error in Convert_DecDeg2RADMS where LST is 359.9999999 or 179.9999999, 89.999999 or higher precision LST_dec_adj = NumTrunc(LST_dec_adj, 6) ' Build the output string strLST = Convert_DecDeg2RADMS(LST_dec_adj, "HHMMSS.SS") & " " & CStr(LST_dec_adj) ' Return result TimeLocalSidereal = strLST End Function Public Function TimeRA2DecHours(ByVal RA_tmp As String, ByVal s_inAngFormat_tmp As String) As Double ' Receives RA time in input format of "HHMMSS.SS" ' and outputs decimal hours ' Uses RADMS2Decdeg to convert string to decimal degrees and then divides by 15 ' After Duffett-Smith at Sec. 7, p. 10. ' Input: Alternatively uses all input formats of RADMS2Decdeg ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' "HHMMSS" ' "HH:MM:SS" ' "HHh MMm SSs" ' Example Duffett-Smith Sec. 7, p. 10: ' TimeRA2DecHours("183127", "HHMMSS") ' 18.5241672666667 Dim dblTime As Double dblTime = Convert_RA2DecDeg(RA_tmp, s_inAngFormat_tmp, "double") dblTime = dblTime / 15# TimeRA2DecHours = dblTime End Function Public Function TimeLStandT2UTC(ByVal dtDate_tmp As Date, ByVal intOffset_tmp As Integer) As String ' Receives MS-Windows system date-time (US) and converts it to ' a UTC date time ' Example: ' ? TimeLStandT2UTC(Now(),-6) ' ? TimeLStandT2UTC(#09/27/2005 16:02:18#,-6) ' 09/27/2005 22:04:21 ' initialize variables Dim dtDate As Date Dim intDay, intMonth, intYear, intHour, intMinute, intSec As Integer Dim strTime As String Dim dblTime As Double Dim intOffset As Integer Dim dtDay As Date Dim strUTC As String ' Output build string ' Get input; set variables dtDate = dtDate_tmp intOffset = intOffset_tmp intYear = CInt(DatePart("yyyy", dtDate)) intMonth = CInt(DatePart("m", dtDate)) intDay = CInt(DatePart("d", dtDate)) intHour = CInt(DatePart("h", dtDate)) intMinute = CInt(DatePart("n", dtDate)) intSec = CInt(DatePart("s", dtDate)) ' Compute the zero-UTC hour day dtDay = DateSerial(intYear, intMonth, intDay) ' Build 24 hour local time string strTime = Format(intHour, "00") & ":" & Format(intMinute, "00") & ":" & Format(intSec, "00") ' Convert to decimal time dblTime = TimeRA2DecHours(strTime, "HH:MM:SS") ' Add the offset dblTime = dblTime - CDbl(intOffset) ' If greater or equal than 24, subtract 24 and add one to the day If dblTime >= 24 Then dblTime = dblTime - 24 dtDay = dtDay + 1 End If ' Build the output string ' strUTC = CStr(dtDay) & " " & TimeDecHours2RA(dblTime, "HH:MM:SS") strUTC = CStr(Format(dtDay, "mm/dd/yyyy")) & " " & TimeDecHours2RA(dblTime, "HH:MM:SS") ' Return TimeLStandT2UTC = strUTC End Function Public Function TimeDecHours2RA(ByVal DecHours_tmp As Double, ByVal s_outAngFormat_tmp As String) ' Receives time in input format of decimal hours ' and outputs RA formatted hours ' Multiples decimal hours by 15 and uses Convert_DecDeg2RADMS to ' convert to RA formatted hours ' After Duffett-Smith at Sec. 8, p. 11. ' Input: Alternatively uses all output formats of Convert_DecDeg2RADMS ' "HHMMSS" ' "HH MM SS" ' "HH MM SS.SS" ' "HHMMSS.SS" ' "HH:MM:SS" ' "HHh MMm SSs" ' Example Duffett-Smith Sec. 8, p. 11: ' TimeDecHours2RA(18.5241672666667, "HH MM SS.SS") ' 18 31 27.00 Dim dblTime As Double Dim outString As String dblTime = CDbl(DecHours_tmp * 15#) outString = Convert_DecDeg2RADMS(dblTime, s_outAngFormat_tmp) TimeDecHours2RA = outString End Function Public Function TimeGST(ByVal strTimetmp As String, ByVal dayDate As Date) As String ' Yields current GST (Greenwich Sidereal Time) at a specific UTC time and date ' Inputs: 24 hour clock time as Universal Time ' Date as date ' Uses TimeLocalSidereal ' Examples: ' TimeGST("16:02:18.23",#09/27/2005#) ' 162828.63 247.11930765782 ' TimeGST("16:02:18.23", #01/27/2005#) ' 003025.68 7.60699764930177 ' Initialize variables and get input Dim strTime As String Dim dtDay As Date Dim strGreenwichLST As String strTime = strTimetmp dtDay = dayDate strGreenwichLST = TimeLocalSidereal(strTime, dtDay, 0) ' Return the result TimeGST = strGreenwichLST End Function Public Function TimeGHA_Aries(ByVal strTimetmp As String, ByVal dayDate As Date) ' Yields current GHA (Greenwich Hour Angle to the First Point in Aries as time and angle. ' Includes a LST check sum for a specific UTC time and date ' Outputs: GHA as decimal time, GHA, backcheck LST for First Point in Aries, backcheck SHA to First Point in Aries (s/b = 0 ) ' Inputs: 24 hour clock time as Universal Time ' Date as date ' Uses TimeGST, TimeLocalSidereal ' Examples: ' TimeGHA_Aries("16:02:18.23", #9/27/2005#) ' GHA_RAhrs GHA LST LST_angle_decdeg ' 073131.37 112.880693 000000.00 0 ' TimeGHA_Aries("16:02:18.23", #1/27/2005#) ' 232934.32 352.393003 000000.00 0 ' Initialize variables and get input Dim strTime, strGST, strGHA As String Dim dtDay As Date Dim dblGHA_Aries As Double Dim dblGHA_Aries_dechours As Double Dim dblLongitude_Aries As Double ' For checking using TimeLocalSidereal Dim strLST_check As String strTime = strTimetmp dtDay = dayDate ' Get the current Greenwich Sidereal Time and Greenwich Hour Angle ' Returns string in form RA, decdeg (hour angle) strGST = TimeGST(strTime, dtDay) ' Compute the greenwich hour angle and decimal hours to the First Point in Aries dblGHA_Aries = 360# - CDbl(GetCSWord(strGST, 2)) dblGHA_Aries_dechours = dblGHA_Aries / 15# ' Compute the longitude of the GHA_Aries for back checking purposes If dblGHA_Aries <= 180# Then dblLongitude_Aries = dblGHA_Aries Else dblLongitude_Aries = dblGHA_Aries - 360# End If ' Back check to LST strLST_check = TimeLocalSidereal(strTime, dtDay, dblLongitude_Aries) ' Build the output string strGHA = TimeDecHours2RA(dblGHA_Aries_dechours, "HHMMSS.SS") strGHA = strGHA & " " & CStr(dblGHA_Aries) & " " & strLST_check TimeGHA_Aries = strGHA End Function Public Function ConvertGHA2Longitude(ByVal dblAngle_tmp As Double) As Double ' Receives a Greenwich Hour Angle as a decimal angle between 0 and 360 ' Return an east west longitude. Positive is east; negative is west Dim dblAngle, dblAngle_adj As Double ' Get input dblAngle = dblAngle_tmp If CDbl(dblAngle) = 0# Or CDbl(dblAngle) = 360# Then dblAngle_adj = 0# End If If CDbl(dblAngle) > 0# And CDbl(dblAngle) < 180# Then dblAngle_adj = CDbl(dblAngle) End If If CDbl(dblAngle) > 180# And CDbl(dblAngle) < 360# Then dblAngle_adj = CDbl(dblAngle) - 360# End If If CDbl(dblAngle) = 180# Then dblAngle_adj = 180# End If ' Error trap - if out of range gracefully return 0 If CDbl(dblAngle) < 0# Or CDbl(dblAngle) > 360# Then dblAngle_adj = 0# End If ConvertGHA2Longitude = dblAngle_adj End Function Function PositionPrecessionEpochtoDate(ByVal RA_dec_tmp As Double, ByVal Dec_dec_tmp As Double, ByVal PM_RA_milliarcsec_tmp As Double, ByVal PM_Dec_milliarcsec_tmp As Double, ByVal Epoch_Start_Jdec_tmp As Double, ByVal Epoch_Final_Jdec_tmp As Double) As Variant ' Receives RA and declination position in decimal form, proper motions in milliarcsec ' Receives initial, reference and final equinox of date epochs in Julian decimal format ' Returns reference and final positions in decimal form as a string of "RA" Space "Dec" and "Epoch" ' and supplement information, e.g. - ' ' # Contents ' 1 RA final in decimal degrees ' 2 Declination final in decimal degrees ' 3 Epoch final in decimal epochs ' 4 Annual variation RA in RA arcsecs due to precession only ' 5 Annual variation Dec in arcsecs due to precession only ' 6 Annual variation RA in RA arcsecs due to precession with proper motions ' 7 Annual variation Dec in arcsecs due to precession with proper motions ' 8 Total variation RA in RA arcsecs over T years ' 9 Total variation Dec in arcsecs over T years ' 10 Total variation RA in decimal degrees over T years ' 11 Total variation Dec in decimal degrees over T years ' 12 T years, the interval between the starting and final epoch ' Function addresses the fact that the vernal equinox regresses about 50" per year along ' ecliptic, and the plane of the ecliptic rotates around the North Pole at about 47" per ' century. ' ' Usage and test data: ' ' Example 1: ' Meeus example at pp. 132-133: ' Reduce alp Leo at Epoch 2000.0 to Epoch 1978.0 where ' alp Leo J10h08m22.3s +11°58'02" or 152.093° +11.967°, and ' alp Leo proper motions are RA -0.0169 RA arcsecs, Dec 0.006 arcsecs ' ' ? PositionPrecessionEpochtoDate(152.091674,11.967,-1.69,6.0,2000.0,1978.0) yields ' 151.79920038843 12.0752035767259 1978 3.20752121712505 -17.7120398278677 3.19062121712505 -17.7060398278677 -70.1936667767512 389.532876213089 -0.292473611569797 0.108203576725858 -22 ' ' ? Convert_DecDeg2RADMS(151.79920038843,"HH MM SS.SS") ' 10 07 11.81 ' Meeus gives 10h 07m 12.1s ' ' ? Convert_DecDeg2RADMS(12.0752035767259,"sDD MM SS") ' +12 04 31 ' Meeus gives +12° 04' 31" ' ' Example 2: ' Duffet Smith at p. 57: ' Reduce star at Epoch 1950.0 to Epoch 1979.5 where ' star J09h10m43s +14°23'25" or 137.679167° +14.390278° and ' there are no proper motions ' ' ? PositionPrecessionEpochtoDate(137.679167,14.390278,0,0,1950.0,1979.5) yields ' 138.0963328752 14.2714146737041 1979.5 3.39389186603259 -14.5053550734007 3.39389186603259 -14.5053550734007 100.119810047961 -427.907974665321 0.417165875199839 -0.118863326295923 29.5 ' ' ? Convert_DecDeg2RADMS(138.0963328752,"HH MM SS.SS") ' 09 12 23.12 ' Duffet-Smith gives 09h 12m 20s ' ' ? Convert_DecDeg2RADMS(14.2714146737041,"sDD MM SS") ' +14 16 17 ' Duffet-Smith gives 14° 16' 08" ' ' The reason for the deviation in the Duffet-Smith result is ' Meeus (1998) uses slightly different values for m,n, and n' ' than Duffet-Smith (1988) ' ' Algorithm sources: ' Precession: Low precision precession model per: ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. p. 132-133. ' See ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. p. 56. ' ' Author: Kurt A. Fisher 2004 Dim RA_start_dec As Double ' Right ascension position at starting epoch in decimal degrees Dim Dec_start_dec As Double ' Declination position at starting epoch in decimal degrees Dim PM_RA_arcsec As Double ' Proper motion in right ascension arcsecs Dim PM_Dec_arcsec As Double ' Proper motion in arcsecs Dim Epoch_Start_Jdec As Double ' Julian starting epoch in decimal form e.g. "1910.9283" or "1950.0" Dim Epoch_Final_Jdec As Double ' Julian final epoch in decimal form e.g. "2000.0" or "2004.5921" ' Usual order of these epochs on a timeline are: ' Starting Epoch -> Reference Epoch -> Final Epoch Dim t As Double ' Number of years between ' Starting Epoch -> Reference Epoch ' or ' Starting Epoch -> Final Epoch Dim RA_final_dec As Double ' Right ascension position at final epoch in decimal degrees Dim Dec_final_dec As Double ' Declination position at final epoch in decimal degrees Dim m_dec As Double ' Working intermediate variable Dim m_arcsec As Double ' Working intermediate variable Dim n_arcsec_dec As Double ' Working intermediate variable Dim n_dec As Double ' Working intermediate variable Dim n_arcsec_RA As Double ' Working intermediate variable Dim delta_RA_arcsec As Double ' Working accrual of right ascension precession in RA arcsecs Dim delta_Dec_arcsec As Double ' Working accrual of declination precession in arcsecs Dim T_reference As Double ' Working variable Dim A, B, C, D, E, f, g, H, x, y As Double ' Temporary working variables Dim S As String ' get input variables RA_start_dec = RA_dec_tmp Dec_start_dec = Dec_dec_tmp PM_RA_arcsec = PM_RA_milliarcsec_tmp / 1000# ' Convert millarcsecs to arcsecs PM_Dec_arcsec = PM_Dec_milliarcsec_tmp / 1000# ' Convert millarcsecs to arcsecs Epoch_Start_Jdec = Epoch_Start_Jdec_tmp Epoch_Final_Jdec = Epoch_Final_Jdec_tmp ' create reused constants ' initialize variables ' Graceful handle no proper motion data If IsNull(PM_RA_arcsec) Then PM_RA_arcsec = 0 If IsNull(PM_Dec_arcsec) Then PM_RA_arcsec = 0 ' Epoch 2000 is the assumed reference year and final epoch if ' final epoch is not set. If IsNull(Epoch_Final_Jdec) Then Epoch_Final_Jdec = 2000 ' Find the interval between the start and final dates t = Epoch_Final_Jdec - Epoch_Start_Jdec T_reference = 2000# - Epoch_Start_Jdec ' Convert input units to units applicable to this algorithm ' Step 1: Find annual precession ' Find the m and n prime components, converting equation to handle decimal arcsecs ' See equations Meeus at p. 132; Duffet at 56 ' Epoch_JDec_tmp is the final or ending epoch. ' These values are negative because the vernal equinox regresses ' about 50" per year along ecliptic. m_arcsec = 3.07496 + (0.00186 * T_reference) ' Is in RA arcseconds n_arcsec_dec = 20.0431 - (0.0085 * T_reference) n_arcsec_RA = n_arcsec_dec / 15 ' Debug.Print m_arcsec, n_arcsec, T, "m,n,T" ' Accrue the total _annual_ precession delta_RA_arcsec = m_arcsec + n_arcsec_RA * Sin(Deg2Rad(RA_start_dec)) * Tan(Deg2Rad(Dec_start_dec)) delta_Dec_arcsec = n_arcsec_dec * Cos(Deg2Rad(RA_start_dec)) A = delta_RA_arcsec ' store annual precession for output B = delta_Dec_arcsec ' store annual precession for output ' Debug.Print delta_RA_arcsec, delta_Dec_arcsec, "Annual Variation, no PM" ' Step 2: Add proper motion to the _annual_ precession delta_RA_arcsec = delta_RA_arcsec + PM_RA_arcsec delta_Dec_arcsec = delta_Dec_arcsec + PM_Dec_arcsec C = delta_RA_arcsec ' store annual precession with PM for output D = delta_Dec_arcsec ' store annual precession with PM for output ' Debug.Print PM_RA_arcsec, PM_Dec_arcsec, "Proper motions" ' Debug.Print delta_RA_arcsec, delta_Dec_arcsec, "Annual Variation with PM" ' Step 3: Multiply the number of years times the annual change to ' obtain total change delta_RA_arcsec = delta_RA_arcsec * t ' in right ascension arcsecs delta_Dec_arcsec = delta_Dec_arcsec * t E = delta_RA_arcsec ' store total variation in arcsecs for output f = delta_Dec_arcsec ' store total variation in arcsecs for output g = RAArcsecs2DecDeg(delta_RA_arcsec) ' store total variation in decimal degrees for output H = DecArcsecs2DecDeg(delta_Dec_arcsec) ' store total variation in decimal degrees for output ' Debug.Print delta_RA_arcsec, delta_Dec_arcsec, T, "Total variation for T years" ' Step 4: Subtract total change from initial positions to obtain ' position at final epoch and reference epoch ' while converting arcsec accrual of change to decimal degrees RA_final_dec = RA_start_dec + RAArcsecs2DecDeg(delta_RA_arcsec) Dec_final_dec = Dec_start_dec + DecArcsecs2DecDeg(delta_Dec_arcsec) ' Debug.Print RA_start_dec, RA_final_dec, RAArcsecs2DecDeg(delta_RA_arcsec) ' Debug.Print Dec_start_dec, Dec_final_dec, DecArcsecs2DecDeg(delta_Dec_arcsec) ' build the output string S = CStr(RA_final_dec) & " " & CStr(Dec_final_dec) & " " & CStr(Epoch_Final_Jdec) S = S & " " & CStr(A) & " " & CStr(B) & " " & CStr(C) S = S & " " & CStr(D) & " " & CStr(E) & " " & CStr(f) S = S & " " & CStr(g) & " " & CStr(H) & " " & CStr(t) ' return result PositionPrecessionEpochtoDate = S End Function Public Function PosEastWestHorizon(ByVal dblLongitude_tmp As Double) As String ' Returns the geographical longitude (unadjusted for atmospheric refraction) ' of the east and west horizons of an observer at geographic longitude in ' decimal form. ' Input: Longitude in decimal degrees. Positive is east; negative is west. ' Range -180 to plus 180. Invalid range > 180 or < - 180 ' Output: ' Returns string rising (eastern) longitude, rising GHA, zenith-OP longitude, ' zenith GHA, setting (western) longitude, setting GHA ' Examples - ' ? PosEastWestHorizon(0) ' 90 90 0 0 -90 270 ' Rising Long 90 ' Rising GHA 90 ' OP Long 0 ' Zenith GHA 0 ' Setting Long -90 ' Setting GHA 270 ' ? PosEastWestHorizon(-118.2) ' -28.2 331.8 -118.2 241.8 151.8 151.8 Dim dblLongitude, L, dblGHA As Double Dim dblRisingGHA, dblSettingGHA, dblRisingLongitude, dblSettingLongitude As Double Dim strOut As String ' Get input dblLongitude = dblLongitude_tmp L = dblLongitude ' Convert longitude back to Greenwich Hour Angle If L = 0 Then dblGHA = 0 End If If L > 0 And L < 180 Then dblGHA = L End If If L = 180 Or L = -180 Then dblGHA = 180 End If If L > -180 And L < 0 Then dblGHA = 360 - Abs(L) End If ' Error trap case If Abs(L) > 180 Then ' Handle error entry gracefully by returning zero dblGHA = 0 End If ' Compute eastern horizon - add 90 deg dblRisingGHA = dblGHA + 90# ' Check and adjust for rotations greater than 360 dblRisingGHA = Adj360Rotation(dblRisingGHA) ' Compute western horizon - subtract 90 deg dblSettingGHA = dblGHA - 90# ' Check and adjust for rotations greater than 360 dblSettingGHA = Adj360Rotation(dblSettingGHA) ' Convert GHA for the eastern horizon back to longitude dblRisingLongitude = ConvertGHA2Longitude(dblRisingGHA) ' Convert GHA for the western horizon back to longitude dblSettingLongitude = ConvertGHA2Longitude(dblSettingGHA) ' Build the output string strOut = Trim(Str(dblRisingLongitude)) & " " & Trim(Str(dblRisingGHA)) & " " strOut = strOut & Trim(Str(dblLongitude)) & " " & Trim(Str(dblGHA)) & " " strOut = strOut & Trim(Str(dblSettingLongitude)) & " " & Trim(Str(dblSettingGHA)) ' Return result PosEastWestHorizon = strOut ' Debug.Print "Rising Long", dblRisingLongitude ' Debug.Print "Rising GHA", dblRisingGHA ' Debug.Print "OP Long", L ' Debug.Print "Zenith GHA", dblGHA ' Debug.Print "Setting Long", dblSettingLongitude ' Debug.Print "Setting GHA", dblSettingGHA End Function Public Function PositionRiseSetZenithLST(ByVal strTime_tmp As String, ByVal dayDate_tmp As Date, decLongitude_tmp As Double) As String ' Computes from longitude and UTC, ' the local sidereal time and GHA of the zenith ' the local sidereal time of the rising eastern horizon without refraction and GHA ' the local sidereal time of the setting western horizon without refraction and GHA ' See TimeLocalSidereal for other restrictions on input ' Inputs: 24 hour clock time as Universal Time e.g. "18:21:45" or "18:21:45.67" ' Date as date ' Longitude as number - West is negative ' Returns as string the RA and Sidereal Hour Angle (SHA) and longitude of the ' rising horizon, observing point and setting horizon ' Examples ' PositionRiseSetZenithLST("04:30:45", #10/14/2005#, -111.87) ' 043434.44 68.6435 -21.87 223434.44 338.6435 -111.87 163434.44 248.6435 158.13 Dim strOut As String Dim strTime As String Dim dayDate As Date Dim decLongitude As Double Dim strZenith As String Dim strRising As String Dim strSetting As String Dim strPosEastWestHorizon As String strTime = strTime_tmp dayDate = dayDate_tmp decLongitude = decLongitude_tmp ' Find the rising and setting longitude strPosEastWestHorizon = PosEastWestHorizon(decLongitude) ' Get the zenith local sideral time and longitude strZenith = Trim(TimeLocalSidereal(strTime, dayDate, decLongitude)) & " " & Trim(CStr(decLongitude)) ' Get the rising local sideral time and rising longitude strRising = Trim(TimeLocalSidereal(strTime, dayDate, GetCSWord(strPosEastWestHorizon, 1))) & " " & Trim(GetCSWord(strPosEastWestHorizon, 1)) ' Get the setting local sideral time and setting longitude strSetting = Trim(TimeLocalSidereal(strTime, dayDate, GetCSWord(strPosEastWestHorizon, 5))) & " " & Trim(GetCSWord(strPosEastWestHorizon, 5)) strOut = Trim(strRising) & " " & Trim(strZenith) & " " & Trim(strSetting) PositionRiseSetZenithLST = strOut ' Debug.Print strPosEastWestHorizon ' Debug.Print strRising ' Debug.Print strZenith ' Debug.Print strSetting End Function Public Function TimeDurationSidereal2Hour24(ByVal DecHours_tmp As Double) As Double ' Converts a duration in decimal hours and sidereal time units to ' local 24 hour time units ' After Duffett-Smith and http://scienceworld.wolfram.com/astronomy/LocalSideralTime.html ' accessed 9/27/2005 ' Example: ' ? TimeDurationSidereal2Hour24(23.9344695931105) ' 24 TimeDurationSidereal2Hour24 = DecHours_tmp * 1.0027379093 End Function Public Function TimeDuration24Hour2Sidereal(ByVal DecHours_tmp As Double) As Double ' Converts a duration in decimal hours and 24 hour time units to ' sidereal time units ' After Duffett-Smith and http://scienceworld.wolfram.com/astronomy/LocalSideralTime.html ' accessed 9/27/2005 ' ? TimeDuration24Hour2Sidereal(24#) ' 23.9344695931105 TimeDuration24Hour2Sidereal = DecHours_tmp / 1.0027379093 End Function Public Function PositionObjectRiseZenithSetMeuss(ByVal objRA_dec_tmp As Double, ByVal objDecDeg_tmp As Double, ByVal strTime_tmp As String, ByVal dayDate_tmp As Date, ByVal decLongitude_tmp As Double, ByVal decLatitude_tmp As Double, altCenterBody_arcmin_tmp) As String ' Given an observing point and the right ascension and declination ' returns the rising, transit and setting times for the object ' and the local hour angle _during a day of Universal Time Coordinated_ ' and not during a day of local time. 'Inputs - ' dayDate_tmp is the current UTC day, not the current local day. (See Meeus at p. 102) ' altCenterBody_deg ' Adjustment for top of object rising before center ' is -34' or 0.5667 degrees for stars and planets ' is -50' or 0.8333 for Sun 'Outputs - A space separated string containing ' Type of computations - (Word #1) ' "Circumpolar-always-visible" ' "Rises-sets" ' "Never-rises-sets" ' For type "Rises-sets" provides triplet stings of rising, transit and setting ' UTC Times Day and Time pairs (Words #2-3,4-5,6-7) ' Hour Angles (Words #8,9,10) ' Azimuth-Altitude pairs (Words #11-12, 13-14, 15-16) ' Use GetCSWord to decode by position ' After Meeus, Astronomical Algorithms, at Chap. 15 (1998) ' Implementation notes - ' This implementation does not follow the Meeus convention of using postive ' longitudes for west of Greenwich. Input using the IAU convention that ' longitudes west of Greenwish are negative. This implementation handles adjusting ' the sign of the longitude accordingly. ' ' This implementaton does not follow the Meuss convention of expressing ' the Hour Angle westwards from the south meridian. Before output, the ' Meuss Hour Angles conventions are converted to a U.S. convention of east is positive ' and west is negative. ' ' Duffet-Smith at sec. 33 is used to compute the rising and setting azimuth. ' ' This implementation gives approximate results per Note 3 in Meeus at page. 104. ' This implementation uses only the first approximately of m values. ' A more accurate implemented m value is not used. ' ' Results are not adjusted for apparent positions, e.g. atmospheric refraction. ' Example Test Data: ' Duffett-Smith at sec 33, p 52 ' PositionObjectRiseZenithSetMeuss(23.655556, 21.7, "00:00:00", #8/24/1980#, 64, 30, -34) ' This Meeus implementation does not yield the Duffett-Smith example results. ' Therefore, this implementation was tested against the MPC Ephemris generator. ' This Meeus implementation was confirmed against the MPC application. ' Therefore, the Duffett-Smith example was rejected. ' Meuss at 103 ' PositionObjectRiseZenithSetMeuss(41.73129, 18.44092, "00:00:00", #3/20/1988#, -71.0833, 42.3333, -34) ' Rises-Sets 03/20/1988 12:26:09 03/20/1988 19:40:18 03/20/1988 02:54:26 ' Hour Angles: 108.023594386663 -0.807882976226566 -108.653713339117 ' Rise Az-H: 64.6651303020785 -0.226569783894222 ' Transit Az-H: 0 66.0977602069409 ' Set Az-H: 295.334869697921 -0.646028580737908 ' Checks out. Rise-set times and Hour angles agrees with Meuss to first m value approximation ' Test against hypothetical current zenithal object ' Note: Yields 90 degree hour angles for only object on celestial equator ' For objects at higher declinations, can yield hour angles greater than 180. ' In northern hemisphere, can yield hour angles lower than 90 degrees ' PositionObjectRiseZenithSetMeuss(30.0157, 0, "00:00:00", #10/19/2005#, -111.87, 41, -34) ' Rises-Sets 10/19/2005 01:34:16 10/19/2005 07:37:16 10/19/2005 13:40:16 ' Hour Angles: 90.6863285822195 -0.312988872815549 -91.3123063278506 ' Rise Az-H: 90 -0.51797342432719 ' Transit Az-H: 180 48.9990165851848 ' Set Az-H: 270 -0.990372882455052 ' Checks out ' Test for Mars at SLC on 10/20/2005 against MPC Ephemris generator ' MPC - Rises between 1:45-2:00 UTC at az 67.1572-69.6041; ' Sets 15:45-16:00 UTC at 292.0648-294.5466 ' PositionObjectRiseZenithSetMeuss(48.87458, 16.55497, "00:00:00", #10/20/2005#, -111.87, 41, -34) ' Rises-Sets 10/20/2005 01:45:37 10/20/2005 08:48:46 10/20/2005 15:51:55 ' Hour Angles: 105.715374866755 -0.361924143651954 -106.439223154059 ' Rise Az-H: 67.8184436931826 -0.516333831064078 ' Transit Az-H: 0 65.5529717673036 ' Set Az-H: 292.181556306817 -1.01954637334685 ' Checks out ' Test for Moon at SLC on 10/20/2005 against MPC Ephemris generator ' MPC sets between 17:30 and 17:35 ' ? PositionObjectRiseZenithSetMeuss(64.59965, 25.18278, "00:00:00", #10/20/2005#, -111.87, 41, -34) ' Rises-Sets 10/20/2005 02:11:31 10/20/2005 09:51:40 10/20/2005 17:31:49 ' Hour Angles: 114.947840965639 -0.4049779438472 -115.757796853333 ' Rise Az-H: 55.6808072050928 -0.510942817247003 ' Transit Az-H: 0 74.1791941672002 ' Set Az-H: 304.319192794907 -1.01088856107392 ' Checks out ' For input variable aliasing Dim strOut As String Dim objRA_dec As Double Dim objDecDeg As Double Dim strTime As String Dim dayDate As Date Dim decLongitude As Double Dim decLatitude As Double Dim altCenterBody_deg As Double ' Adjustment for top of object rising before center ' is -34' or 0.5667 degrees for stars and planets ' is -50' or 0.8333 for Sun Dim B, C, D, E, f, g, H, I, J, L As Double ' Working variables Dim strTemp As String ' Working variables ' For Meuss rising setting times and hour angles Dim decdegGHAGreenwich As Double Dim bolVisibleTest As Boolean Dim decLongitude_adj As Double Dim intDay, intMonth, intYear, intHour, intMinute As Integer Dim strComputType As String Dim m_rising_dayfrac As Double Dim m_transit_dayfrac As Double Dim m_setting_dayfrac As Double Dim dblH0 As Double Dim strDateDay As String Dim strTransit_UTC As String Dim strRising_UTC As String Dim strSetting_UTC As String Dim LST_dec_adj As Double ' Working variable Dim decdegTransitGST As Double Dim decdegRisingGST As Double Dim decdegSettingGST As Double Dim decdegHourAngle_Transit As Double Dim decdegHourAngle_Rising As Double Dim decdegHourAngle_Setting As Double ' For Duffet-Smith rising, transit, setting azimuths Dim A_rising, A_transit, A_setting As Double Dim strAzRising As String Dim strAzTransit As String Dim strAzSetting As String Dim strh_Rising As String Dim strh_Transit As String Dim strh_Setting As String Dim h_rising As Double Dim h_transit As Double Dim h_setting As Double Dim strOutStore As String ' Get input variables objRA_dec = objRA_dec_tmp objDecDeg = objDecDeg_tmp strTime = strTime_tmp dayDate = dayDate_tmp decLongitude = decLongitude_tmp decLatitude = decLatitude_tmp altCenterBody_deg = altCenterBody_arcmin_tmp / 60# ' Convert arcmins to degrees ' Build the target day as a string for later use in output intYear = CInt(DatePart("yyyy", dayDate)) intMonth = CInt(DatePart("m", dayDate)) intDay = CInt(DatePart("d", dayDate)) strDateDay = Format(intMonth, "00") & "/" & Format(intDay, "00") & "/" & Format(intYear, "0000") ' Compute the LST at Greenwhich at 0h UTC on the target date ' in degrees strTemp = TimeLocalSidereal("00:00:00", dayDate, 0) decdegGHAGreenwich = CDbl(GetCSWord(strTemp, 2)) ' Express right ascension of object in degrees. ' objRA_dec already is in decimal degrees ' Express longitude in degrees and invert the sign per Meuss convention. decLongitude_adj = -1# * decLongitude ' Compute quantity H - a parameter used in computing ' the m values of the rising and setting times ' Test quantity H for range ' Meuss test on "second member" at page 102, Eq. 15.1 C = CDbl(-1# * Tan(Deg2Rad(objDecDeg)) * Tan(Deg2Rad(decLatitude))) ' Test and filter for circumpolar and below horizon objects If C > 1 Then strComputType = "Never-rises-or-sets" bolVisibleTest = False ' Store values for output string strRising_UTC = "00/00/0000 00:00:00" strTransit_UTC = "00/00/0000 00:00:00" strSetting_UTC = "00/00/0000 00:00:00" decdegHourAngle_Rising = 0 decdegHourAngle_Transit = 0 decdegHourAngle_Setting = 0 strAzRising = "0" strAzTransit = "0" strAzSetting = "0" strh_Rising = "0" strh_Transit = "0" strh_Setting = "0" End If If C < -1 Then strComputType = "Circumpolar-always-visible" bolVisibleTest = False strRising_UTC = "00/00/0000 00:00:00" strTransit_UTC = "00/00/0000 00:00:00" strSetting_UTC = "00/00/0000 00:00:00" decdegHourAngle_Rising = 0 decdegHourAngle_Transit = 0 decdegHourAngle_Setting = 0 strAzRising = "0" strAzTransit = "0" strAzSetting = "0" strh_Rising = "0" strh_Transit = "0" strh_Setting = "0" End If If C >= -1 And C <= 1 Then bolVisibleTest = True strComputType = "Rises-Sets" End If ' If object is visible ' Find the quantity H0 in degrees If bolVisibleTest = True Then ' Meuss Eq. 15.1 D = CDbl(Sin(Deg2Rad(altCenterBody_deg)) - (Sin(Deg2Rad(objDecDeg)) * Sin(Deg2Rad(decLatitude)))) E = CDbl((Cos(Deg2Rad(objDecDeg)) * Cos(Deg2Rad(decLatitude)))) dblH0 = Rad2Deg(ArcCos(D / E)) End If ' Compute m values for rising transit and setting ' m is a fractional part of a day If bolVisibleTest = True Then ' Compute the transit m value ' Meuss Eq. 15.2 m_transit_dayfrac = CDbl((objRA_dec + decLongitude_adj - decdegGHAGreenwich) / 360#) m_rising_dayfrac = CDbl(m_transit_dayfrac - (dblH0 / 360#)) m_setting_dayfrac = CDbl(m_transit_dayfrac + (dblH0 / 360#)) ' and adjust range of m so it is between 0 and 1 ' Meuss p. 102 If m_transit_dayfrac > 1 Then m_transit_dayfrac = m_transit_dayfrac - 1 End If If m_transit_dayfrac < 0 Then m_transit_dayfrac = m_transit_dayfrac + 1 End If If m_rising_dayfrac > 1 Then m_rising_dayfrac = m_rising_dayfrac - 1 End If If m_rising_dayfrac < 0 Then m_rising_dayfrac = m_rising_dayfrac + 1 End If If m_setting_dayfrac > 1 Then m_setting_dayfrac = m_setting_dayfrac - 1 End If If m_setting_dayfrac < 0 Then m_setting_dayfrac = m_setting_dayfrac + 1 End If ' Store the m values in output formats strTransit_UTC = strDateDay & " " & TimeDecHours2RA((m_transit_dayfrac * 24#), "HH:MM:SS") strRising_UTC = strDateDay & " " & TimeDecHours2RA((m_rising_dayfrac * 24#), "HH:MM:SS") strSetting_UTC = strDateDay & " " & TimeDecHours2RA((m_setting_dayfrac * 24#), "HH:MM:SS") End If 'If bolVisibleTest = True for compute m values ' Compute the Hour Angle for the object at rising, setting transit ' Meuss at p. 103 If bolVisibleTest = True Then ' Find the rising Greenwich Sidereal Time expressed in degrees decdegRisingGST = CDbl(decdegGHAGreenwich + (360.985647 * m_rising_dayfrac)) ' Round to less than 360 degs if greater than 360 (omitted in Meuss) LST_dec_adj = decdegRisingGST Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop decdegRisingGST = LST_dec_adj decdegHourAngle_Rising = decdegRisingGST - decLongitude_adj - objRA_dec ' Find the transit Greenwich Sidereal Time expressed in degrees decdegTransitGST = CDbl(decdegGHAGreenwich + (360.985647 * m_transit_dayfrac)) ' Round to less than 360 degs if greater than 360 (omitted in Meuss) LST_dec_adj = decdegTransitGST Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop decdegTransitGST = LST_dec_adj decdegHourAngle_Transit = decdegTransitGST - decLongitude_adj - objRA_dec ' Find the setting Greenwich Sidereal Time expressed in degrees decdegSettingGST = CDbl(decdegGHAGreenwich + (360.985647 * m_setting_dayfrac)) ' Round to less than 360 degs if greater than 360 (omitted in Meuss) LST_dec_adj = decdegSettingGST Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop decdegSettingGST = LST_dec_adj decdegHourAngle_Setting = decdegSettingGST - decLongitude_adj - objRA_dec End If ' If bolVisibleTest for Compute the Hour Angle ' Compute the rising, transit and setting azimuth using Duffett-Smith If bolVisibleTest = True Then D = CDbl(Sin(Deg2Rad(objDecDeg)) / Cos(Deg2Rad(decLatitude))) A_rising = Rad2Deg(ArcCos(CDbl(D))) A_setting = 360# - A_rising strAzRising = CStr(A_rising) strAzSetting = CStr(A_setting) ' Handle special case of transit value in northern vs southern portions of local horizon coordinate system If A_rising < 90# Then strAzTransit = CStr(0#) Else strAzTransit = CStr(180#) End If End If ' If bolVisibleTest for Compute rising and setting azimuth ' Compute the rising, transit and setting altiude using Meuss Eq. 13.6 at 93 If bolVisibleTest = True Then f = CDbl(Sin(Deg2Rad(objDecDeg)) * Sin(Deg2Rad(decLatitude))) g = CDbl(Cos(Deg2Rad(objDecDeg)) * Cos(Deg2Rad(decLatitude))) ' Rising case I = Cos(Deg2Rad(decdegHourAngle_Rising)) h_rising = Rad2Deg(ArcSin(CDbl(f + (g * I)))) ' Transit case I = Cos(Deg2Rad(decdegHourAngle_Transit)) h_transit = Rad2Deg(ArcSin(CDbl(f + (g * I)))) ' Setting case I = Cos(Deg2Rad(decdegHourAngle_Setting)) h_setting = Rad2Deg(ArcSin(CDbl(f + (g * I)))) strh_Rising = CStr(h_rising) strh_Transit = CStr(h_transit) strh_Setting = CStr(h_setting) End If ' If bolVisibleTest for Compute rising and setting altitude ' Store rising, transit and setting times and hour angle values for later output. strOutStore = strComputType & " " & Trim(strRising_UTC) & " " & Trim(strTransit_UTC) & " " & Trim(strSetting_UTC) & " " ' Adjust to reverse Meuss sign convention strOutStore = strOutStore & Trim(CStr(-decdegHourAngle_Rising)) & " " & Trim(CStr(-decdegHourAngle_Transit)) & " " & Trim(CStr(-decdegHourAngle_Setting)) & " " ' Add the rising, transit, setting azimuth and altitude strOutStore = strOutStore & Trim(strAzRising) & " " & Trim(strh_Rising) & " " & Trim(strAzTransit) & " " & Trim(strh_Transit) & " " strOutStore = strOutStore & Trim(strAzSetting) & " " & Trim(strh_Setting) & " " PositionObjectRiseZenithSetMeuss = strOutStore End Function Public Function EclipticObliquity(ByVal JulianDay_dbl_tmp As Double) As Double ' Returns a low precision obliquity of the ecliptic on a given Julian Day, ' using the algorithm of Duffett-Smith. ' Source: Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. p. 41, Sec. 27. ' Example - Duffett-Smith at p. 41, obliquity of ecliptic on 1/1/1980.0 ' Find obliquity of the ecliptic on Julian Date 1980 Jan. 0.0 ' which is Gregorian dated 1/1/1980 12:00:00 UTC ' or Julian Day 2444238.50 ' ? EclipticObliquity(2444238.5) ' 23.4418933780795 ' Define variables Dim JulianDay_dbl As Double Dim t As Double Dim deltaObliquity_arcsec As Double Dim Obliquity_dec_deg As Double Dim Obliquity_arcsec As Double ' Get input variables JulianDay_dbl = JulianDay_dbl_tmp ' Convert Julian Day to Julian centuries since Julian 1/1/2000 at 0:00 (or noon UTC) ' This is value T t = CDbl((JulianDay_dbl - 2451545#) / 36525#) ' Apply Duffett-Smith alogrithm using T to find the change in obliquity in arcsec Obliquity_arcsec = CDbl((46.815 * t) - (0.0006 * t ^ 2) + (0.00181 * t ^ 3)) ' Apply the change in the obliquity to obliquity at ' Julian 1/1/2000 at 0:00 (or noon UTC) Obliquity_dec_deg = 23.439292 - CDbl(Obliquity_arcsec / 3600#) ' Return the result. EclipticObliquity = Obliquity_dec_deg End Function Public Function TimeGHA2Object(ByVal objRA_dec_tmp As Double, ByVal strTime_tmp As String, ByVal dayDate_tmp As Date) As String ' Gives angle between object after last passage of the meridian. ' Returns GHA as ' a positive westward angle - 0 to 360 degrees ' a time in hours ' as a terresterial longitude of object with east as positive and ' west as negative ' The relationship between an object's hour angle, the local sidereal time ' and the object's position is ' LHA = LST - R.A. ' For Greenwich sidereal time, the relationship is: ' GHA = GST - R.A. ' An object's Greenwich hour angle, when converted to a -+ range between ' 0 and 180 degrees gives the object's current longitude. When combined with ' the object's declination, this identifies the point on the Earth over ' which the celestial object is located. This is used in celestial ' navigation with a sextant to determine the observer's current position. ' Source: Built around ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. pp. 35-37, Secs 25 and 24. ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. pp. 92-95. ' Author: Kurt A. Fisher 10/2005 ' Uses TimeLocalSidereal and ConvertLHA2Longitude ' Inputs: ' Adapted from PositionObjectRiseZenithSetMeuss except that inputs ' decLongitude_tmp = 0 and decLatitude_tmp = 0 and altCenterBody_deg = 0 ' are omitted and set to zero and - ' objDecDeg_tmp - ' declination is omitted. ' Output: ' 1) Input string as day ' 2) Input UTC time ' 3) Greenwich Sidereal Time as degrees e.g 261.183465 ' 4) Greenwich Sidereal Time as decimal hours e.g 17.412231 ' 5) Greenwich Hour angle Object to Greenwich as degree angle e.g 343.09596 ' 6) Greenwich Hour angle Object (time since last transit) to Greenwich as decimal hours e.g 22.873064 ' 7) Greenwich Hour angle Object (time since last transit) to Greenwich as RA hours e.g 225223.03 ' 8) Greenwich Local Hour angle Object to Greenwich as longitude angle degrees (angle to next transit) e.g 16.90404 ' 9) Greenwich Local Hour angle Object to Greenwich as as longitude decimal hours (time to next transit) e.g 1.126936 ' 10) Greenwich Local Hour angle Object to Greenwich as as longitude RA hours (time to next transit) e.g 010736.97 ' Testing: ' Duffet-Smith example at page 35 ' Where 18.539167 RA equals 278.087505 degrees ' ? TimeGHA2Object(278.087505, "14:36:51.67", #11/2/2005#) ' 11/02/2005 14:36:51.67 261.183465 17.412231 343.09596 22.873064 225223.03 16.90404 1.126936 010736.97 ' Input time and date 11/02/2005 14:36:51.67 ' Greenwich Sidereal Time as degrees 261.183465 ' Greenwich Sidereal Time as decimal hours 17.412231 ' Greenwich Hour Angle Object to Greenwich as degree angle 343.09596 ' Greenwich Hour Angle Object (time since last transit) to Greenwich as decimal hours 22.873064 ' Greenwich Hour Angle Object (time since last transit) to Greenwich as RA hours 225223.03 ' Local Hour Angle Object to Greenwich as longitude angle degrees 16.90404 ' Local Hour Angle Object to Greenwich as longitude decimal hours 1.126936 ' Local Hour Angle Object to Greenwich as longitude RA hours 010736.97 ' which Checks against http://www.roman-britain.org/astronomy/astro.htm# ' Typical usage extracts the Greenwich hour angle as either a degree, decimal hour or RA hour, e.g. ' ? GetCSWord(TimeGHA2Object(278.087505, "14:36:51.67", #11/2/2005#), 5) ' ? GetCSWord(TimeGHA2Object(278.087505, "14:36:51.67", #11/2/2005#), 6) ' ? GetCSWord(TimeGHA2Object(278.087505, "14:36:51.67", #11/2/2005#), 7) ' For input variable aliasing Dim strOut As String Dim strOutStore As String Dim objRA_dec As Double Dim objDecDeg As Double Dim strTime As String Dim dayDate As Date Dim decLongitude As Double Dim decLatitude As Double Dim altCenterBody_deg As Double ' Adjustment for top of object rising before center ' is -34' or 0.5667 degrees for stars and planets ' is -50' or 0.8333 for Sun Dim B, C, D, E, f, g, H, I, J, L As Double ' Working variables Dim strTemp As String ' Working variables ' Greenwich sidereal time as degrees Dim decdegGSTGreenwich As Double Dim decLongitude_adj As Double Dim intDay, intMonth, intYear, intHour, intMinute As Integer Dim strDateDay As String Dim LST_dec_adj As Double ' Working variable Dim LST_dec As Double Dim decdegHourAngle As Double Dim sgnHourAngle As String ' Get input variables objRA_dec = objRA_dec_tmp ' objDecDeg = objDecDeg_tmp strTime = strTime_tmp dayDate = dayDate_tmp decLongitude = 0# ' Forced to Greenwich longitude and latitude decLatitude = 0# ' Forced to Greenwich longitude and latitude altCenterBody_deg = 0# / 60# ' Forced to zero - not setting or rising, is transit ' Convert arcmins to degrees ' Build the target day as a string for later use in output intYear = CInt(DatePart("yyyy", dayDate)) intMonth = CInt(DatePart("m", dayDate)) intDay = CInt(DatePart("d", dayDate)) strDateDay = Format(intMonth, "00") & "/" & Format(intDay, "00") & "/" & Format(intYear, "0000") ' Compute the current LST at Greenwich at the current UTC on the target date ' in degrees strTemp = TimeLocalSidereal(strTime, dayDate, 0) decdegGSTGreenwich = CDbl(GetCSWord(strTemp, 2)) ' Express right ascension of object in degrees. ' objRA_dec already is in decimal degrees ' Express longitude in degrees and invert the sign per Meuss convention. ' decLongitude_adj = -1# * decLongitude ' Find the Greenwich Hour Angle of the object decdegHourAngle = decdegGSTGreenwich - objRA_dec - decLongitude_adj ' Round to less than 360 degs if greater than 360 or less than 0 LST_dec = decdegHourAngle ' Use iteration method If LST_dec >= 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop End If If LST_dec < 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj <= 360# LST_dec_adj = LST_dec_adj + 360# Loop LST_dec_adj = LST_dec_adj - 360# End If decdegHourAngle = LST_dec_adj ' Build the output string ' Store the hour angle sign as string for use in output formatting If Sgn(ConvertLHA2Longitude(decdegHourAngle)) = -1 Then sgnHourAngle = "-" Else sgnHourAngle = "+" End If ' Debug.Print "Input time and date", strDateDay, strTime ' Debug.Print "Greenwich Sidereal Time as degrees", decdegGSTGreenwich ' Debug.Print "Greenwich Sidereal Time as decimal hours", decdegGSTGreenwich / 15# ' Debug.Print "Greenwich Hour angle Object to Greenwich as degree angle", decdegHourAngle ' Debug.Print "Greenwich Hour angle Object (time since last transit) to Greenwich as decimal hours", decdegHourAngle / 15# ' Debug.Print "Greenwich Hour angle Object (time since last transit) to Greenwich as RA hours", TimeDecHours2RA(decdegHourAngle / 15#, "HHMMSS.SS") ' Debug.Print "Local Hour angle Object to Greenwich as longitude angle degrees (angle to next transit)", ConvertLHA2Longitude(decdegHourAngle) ' Debug.Print "Local Hour angle Object to Greenwich as longitude decimal hours (time to next transit)", ConvertLHA2Longitude(decdegHourAngle) / 15# ' Debug.Print "Local Hour angle Object to Greenwich as longitude RA hours (time to next transit)", sgnHourAngle & TimeDecHours2RA(ConvertLHA2Longitude(decdegHourAngle) / 15#, "HHMMSS.SS") ' Store input date/time and Greenwich sidereal times for output. strOutStore = Trim(strDateDay) & " " & Trim(strTime) & " " & Trim(CStr(decdegGSTGreenwich)) & " " & Trim(CStr(decdegGSTGreenwich / 15#)) & " " ' Store Greenwich Hour Angle (time since last transit) as degree angle, decimal hours, and RA hours strOutStore = strOutStore & Trim(CStr(decdegHourAngle)) & " " & Trim(CStr(decdegHourAngle / 15#)) & " " & Trim(TimeDecHours2RA(decdegHourAngle / 15#, "HHMMSS.SS")) & " " ' Store Local Hour Angle (time to next transit) as degree angle, decimal hours, and RA hours strOut = strOutStore & Trim(CStr(ConvertLHA2Longitude(decdegHourAngle))) & " " & Trim(CStr(ConvertLHA2Longitude(decdegHourAngle) / 15#)) & " " & sgnHourAngle & Trim(TimeDecHours2RA(ConvertLHA2Longitude(decdegHourAngle) / 15#, "HHMMSS.SS")) TimeGHA2Object = strOut End Function Public Function ConvertLHA2Longitude(ByVal dblAngle_tmp As Double) As Double ' Receives a Local Hour Angle as a decimal angle between 0 and 360 ' Return an east west longitude. Positive is east; negative is west Dim dblAngle, dblAngle_adj As Double ' Get input dblAngle = dblAngle_tmp If CDbl(dblAngle) = 0# Or CDbl(dblAngle) = 360# Then dblAngle_adj = 0# End If If CDbl(dblAngle) > 0# And CDbl(dblAngle) < 180# Then dblAngle_adj = -1 * CDbl(dblAngle) End If If CDbl(dblAngle) > 180# And CDbl(dblAngle) < 360# Then dblAngle_adj = 360# - CDbl(dblAngle) End If If CDbl(dblAngle) = 180# Then dblAngle_adj = 180# End If ' Error trap - if out of range gracefully return 0 If CDbl(dblAngle) < 0# Or CDbl(dblAngle) > 360# Then dblAngle_adj = 0# End If ConvertLHA2Longitude = dblAngle_adj End Function Public Function ConvertLongLatDMS2Dec(ByVal strDirection_tmp As String, ByVal intDeg_tmp As Integer, ByVal intMin_tmp As Integer, ByVal dblSec_tmp As Double) As Double ' Converts a longitude or latitude entered with a "N", "S", "E", "W" and DMS values ' into a signed decimal degree, e.g - ' W 43 12 8.0 = -43.385556 ' Uses Convert_DMS2DecDeg ' Define variables Dim intDeg, intMin As Integer Dim dblSec As Double Dim strDirection As String ' Working variables Dim strSign As String Dim A As Double Dim strTemp As String ' Get input strDirection = strDirection_tmp intDeg = intDeg_tmp intMin = intMin_tmp dblSec = dblSec_tmp ' Convert the sign Select Case strDirection Case "South", "S", "s", "West", "W", "w" strSign = "-" Case "North", "N", "n", "East", "E", "e" strSign = "+" Case Else ' Gracefully handle error by returning a + strSign = "+" End Select strTemp = strSign & Trim(Format(intDeg, "000")) & " " & Trim(Format(intMin, "00")) & " " & Trim(Format(NumTrunc(dblSec, 1), "00.0")) ConvertLongLatDMS2Dec = Convert_DMS2DecDeg(strTemp, "sDDD MM SS.S", "double") End Function Public Function TimeLHA2Object(ByVal objRA_dec_tmp As Double, ByVal strTime_tmp As String, ByVal dayDate_tmp As Date, ByVal decLongitude_tmp As Double) As String ' Gives angle between object after last passage of the meridian. ' Returns LHA as a positive westward angle - 0 to 360 degrees and ' as a meridian longitude with east as positive and west as negative ' Gives angle between object after last passage of the meridian. ' Returns GHA as ' a positive westward angle - 0 to 360 degrees ' a time in hours ' as a terresterial longitude of object with east as positive and ' west as negative ' The relationship between an object's hour angle, the local sidereal time ' and the object's position is ' LHA = LST_greenwich - R.A. - Longitude ' or: ' LHA = LST - R.A. ' Source: Built around ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. pp. 35-37, Secs 25 and 24. ' Meeus, J. 1998. Astronomical Algorithms. 2nd ed. ' Willmann-Bell. pp. 92-95. ' Author: Kurt A. Fisher 10/2005 ' Uses TimeLocalSidereal and ConvertLHA2Longitude ' Inputs: ' Adapted from PositionObjectRiseZenithSetMeuss except that inputs ' decLatitude_tmp = 0 and altCenterBody_deg = 0 ' are omitted and set to zero and - ' objDecDeg_tmp ' declination is omitted. ' Output: ' 1) Input string day ' 2) Input string UTC time ' 3) Local Sidereal Time as degrees ' 4) Local Sidereal Time as decimal hours ' 5) Local Hour Angle Meridian to Object degree angle ' 6) Local Hour Angle Meridian to Object (time since last transit) as decimal hours ' 7) Local Hour Angle Meridian to Object (time since last transit) as RA hours ' 8) Local Hour Angle Meridian to Object as longitude angle degrees ' 9) Local Hour Angle Meridian to Object as longitude decimal hours ' 10) Local Hour Angle Meridian to Object as longitude RA hours ' Testing: ' Object off SLC local meridian ' ? TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87) ' Input time and date 11/02/2005 14:36:51.67 ' Local Sidereal Time as degrees 149.313465 ' Local Sidereal Time as decimal hours 9.954231 ' Local Hour Angle Meridian to Object degree angle 231.22596 ' Local Hour Angle Meridian to Object (time since last transit) as decimal hours 15.415064 ' Local Hour Angle Meridian to Object (time since last transit) as RA hours 152454.23 ' Local Hour Angle Meridian to Object as longitude angle degrees 128.77404 ' Local Hour Angle Meridian to Object as longitude decimal hours 8.584936 ' Local Hour Angle Meridian to Object as longitude RA hours +083505.77' which Checks against http://www.roman-britain.org/astronomy/astro.htm# ' Typical usage extracts the local hour angle as either a degree, decimal hour or RA hour, e.g. ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 5) ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 6) ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 7) ' or to extract LHA as longitude degree, decimal hour or RAMS formatted time ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 8) ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 9) ' ? GetCSWord(TimeLHA2Object(278.087505, "14:36:51.67", #11/2/2005#, -111.87), 10) ' For input variable aliasing Dim strOut As String Dim strOutStore As String Dim objRA_dec As Double Dim objDecDeg As Double Dim strTime As String Dim dayDate As Date Dim decLongitude As Double Dim decLatitude As Double Dim altCenterBody_deg As Double ' Adjustment for top of object rising before center ' is -34' or 0.5667 degrees for stars and planets ' is -50' or 0.8333 for Sun Dim B, C, D, E, f, g, H, I, J, L As Double ' Working variables Dim strTemp As String ' Working variables ' Local sidereal time as degrees Dim decdegLST As Double Dim decLongitude_adj As Double Dim intDay, intMonth, intYear, intHour, intMinute As Integer Dim strDateDay As String Dim LST_dec_adj As Double ' Working variable Dim LST_dec As Double Dim decdegHourAngle As Double Dim sgnHourAngle As String ' Get input variables objRA_dec = objRA_dec_tmp ' objDecDeg = objDecDeg_tmp strTime = strTime_tmp dayDate = dayDate_tmp decLongitude = decLongitude_tmp decLatitude = 0# ' Not used altCenterBody_deg = 0# / 60# ' Forced to zero - not setting or rising, is transit ' Convert arcmins to degrees ' Build the target day as a string for later use in output intYear = CInt(DatePart("yyyy", dayDate)) intMonth = CInt(DatePart("m", dayDate)) intDay = CInt(DatePart("d", dayDate)) strDateDay = Format(intMonth, "00") & "/" & Format(intDay, "00") & "/" & Format(intYear, "0000") ' Compute the current LST at the target UTC and date ' in degrees strTemp = TimeLocalSidereal(strTime, dayDate, decLongitude) decdegLST = CDbl(GetCSWord(strTemp, 2)) ' Express right ascension of object in degrees. ' objRA_dec already is in decimal degrees ' Express longitude in degrees and _do not_ invert the sign per Meuss convention. ' Compare to PositionObjectRiseZenithSetMeuss decLongitude_adj = 1# * decLongitude ' Find the Local Hour Angle of the object decdegHourAngle = decdegLST - objRA_dec ' Debug.Print decdegHourAngle ' Round to less than 360 degs if greater than 360 or less than 0 LST_dec = decdegHourAngle ' Use iteration method If LST_dec >= 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj >= 360# LST_dec_adj = LST_dec_adj - 360 Loop End If If LST_dec < 0 Then LST_dec_adj = LST_dec Do While LST_dec_adj <= 360# LST_dec_adj = LST_dec_adj + 360# Loop LST_dec_adj = LST_dec_adj - 360# End If decdegHourAngle = LST_dec_adj ' Build the output string If Sgn(ConvertLHA2Longitude(decdegHourAngle)) = -1 Then sgnHourAngle = "-" Else sgnHourAngle = "+" End If ' Debug.Print "Input time and date", strDateDay, strTime ' Debug.Print "Local Sidereal Time as degrees", decdegLST ' Debug.Print "Local Sidereal Time as decimal hours", decdegLST / 15# ' Debug.Print "Local Hour Angle Meridian to Object degree angle", decdegHourAngle ' Debug.Print "Local Hour Angle Meridian to Object (time since last transit) as decimal hours", decdegHourAngle / 15# ' Debug.Print "Local Hour Angle Meridian to Object (time since last transit) as RA hours", TimeDecHours2RA(decdegHourAngle / 15#, "HHMMSS.SS") ' Debug.Print "Local Hour Angle Meridian to Object as longitude angle degrees", ConvertLHA2Longitude(decdegHourAngle) ' Debug.Print "Local Hour Angle Meridian to Object as longitude decimal hours", ConvertLHA2Longitude(decdegHourAngle) / 15# ' Debug.Print "Local Hour Angle Meridian to Object as longitude RA hours", sgnHourAngle & TimeDecHours2RA(ConvertLHA2Longitude(decdegHourAngle) / 15#, "HHMMSS.SS") ' Store input date/time and Greenwich sidereal times for output. strOutStore = Trim(strDateDay) & " " & Trim(strTime) & " " & Trim(CStr(decdegLST)) & " " & Trim(CStr(decdegLST / 15#)) & " " ' Store Greenwich Hour Angle (time since last transit) as degree angle, decimal hours, and RA hours strOutStore = strOutStore & Trim(CStr(decdegHourAngle)) & " " & Trim(CStr(decdegHourAngle / 15#)) & " " & Trim(TimeDecHours2RA(decdegHourAngle / 15#, "HHMMSS.SS")) & " " ' Store Local Hour Angle (time to next transit) as degree angle, decimal hours, and RA hours strOut = strOutStore & Trim(CStr(ConvertLHA2Longitude(decdegHourAngle))) & " " & Trim(CStr(-ConvertGHA2Longitude(decdegHourAngle) / 15#)) & " " & Trim(TimeDecHours2RA(ConvertLHA2Longitude(decdegHourAngle) / 15#, "HHMMSS.SS")) TimeLHA2Object = strOut End Function Public Function TimeJulianCenturiesPeriod(ByVal Year_tmp As Integer, ByVal Month_tmp As Integer, ByVal Day_tmp As Integer, ByVal UTC_Hour_tmp As Integer, _ ByVal UTC_Min_tmp As Integer, ByVal UTC_Sec_tmp As Double, ByVal Year2_tmp As Integer, ByVal Month2_tmp As Integer, ByVal Day2_tmp As Integer, _ ByVal UTC_Hour2_tmp As Integer, ByVal UTC_Min2_tmp As Integer, ByVal UTC_Sec2_tmp As Double, ByVal IsGregorian_tmp As Integer) As String ' Accepts two times and dates and returns the duration between the two in Julian centuries. ' Uses TimeJulianDay ' After Meeus (1998) at p. 134 ' Test and usage ' IsGregorian_tmp is normally set to 1; 0 is for Julian dates ' Meeus (198) at p. 135 ' Find the Julian centuries between J2000.0 and Nov. 13.19 on 2028 ' ? TimeDecHours2RA(0.19 * 24,"HHMMSS") = 043336 ' ? TimeJulianCenturiesPeriod(2028, 11, 13, 4, 33, 36, 2000, 1, 1, 12, 0, 0, 1) ' 0.288670499657767 ' Define variables Dim intY, intM, intD, intH, intMin As Integer Dim dblSec As Double Dim intY2, intM2, intD2, intH2, intMin2 As Integer Dim dblSec2 As Double Dim dblJulianDay1 As Double Dim dblJulianDay2 As Double Dim IsGregorian_tmp2 As Integer ' Get inputs ' First date intY = Year_tmp intM = Month_tmp intD = Day_tmp intH = UTC_Hour_tmp intMin = UTC_Min_tmp dblSec = UTC_Sec_tmp ' Second date intY2 = Year2_tmp intM2 = Month2_tmp intD2 = Day2_tmp intH2 = UTC_Hour2_tmp intMin2 = UTC_Min2_tmp dblSec2 = UTC_Sec2_tmp ' Gregorian flag IsGregorian_tmp2 = IsGregorian_tmp ' Get the Julian Day for the two dates dblJulianDay1 = CDbl(GetCSWord(TimeJulianDay(intY, intM, intD, intH, intMin, dblSec, IsGregorian_tmp2), 1)) dblJulianDay2 = CDbl(GetCSWord(TimeJulianDay(intY2, intM2, intD2, intH2, intMin2, dblSec2, IsGregorian_tmp2), 1)) ' Compute and return the result TimeJulianCenturiesPeriod = (dblJulianDay1 - dblJulianDay2) / 36525 End Function Public Function TimeDecimalYear(ByVal Year_tmp As Integer, ByVal Month_tmp As Integer, ByVal Day_tmp As Integer, ByVal UTC_Hour_tmp As Integer, _ ByVal UTC_Min_tmp As Integer, ByVal UTC_Sec_tmp As Double) As String ' Accepts a Gregorian calendar time string and date and returns the calendar ' date in decimal year format - ' 1993.39382 is Year 1993 plus 0.39382 fractional part of a year ' Computes based on sideral year or 365.2564 days in a year ' Define variables Dim intY, intM, intD, intH, intMin As Integer Dim dblSec As Double Dim dblYearLength As Double dblYearLength = 365.2564 ' Sideral year or 365.2564 days in a year Dim dblJulianDayBegYear As Double ' Julian Day of the start of the year Dim dblJulianDayDate As Double ' Julian Day of the target date Dim C As Double ' Result - decimal year ' Get inputs intY = Year_tmp intM = Month_tmp intD = Day_tmp intH = UTC_Hour_tmp intMin = UTC_Min_tmp dblSec = UTC_Sec_tmp ' Find the time interval in days between the start of the year and the target day and time dblJulianDayBegYear = CDbl(GetCSWord(TimeJulianDay(intY, 0, 0, 0, 0, 0, 1), 1)) dblJulianDayDate = CDbl(GetCSWord(TimeJulianDay(intY, intM, intD, intH, intM, dblSec, 1), 1)) ' Get the date part express in centuries C = (dblJulianDayDate - dblJulianDayBegYear) / (dblYearLength) C = C + CDbl(intY) TimeDecimalYear = C End Function Public Function SchaeferLimitingMag(ByVal dblAperature_inch As Double, ByVal dblMagnification As Double, ByVal dblNELM As Double, ByVal dblAge As Double, ByVal strTelescopeType As String, ByVal dblCleanliness As Double, ByVal dblColorIndex As Double, ByVal dblZenithDist_deg As Double, ByVal dblExtinctionCoeff As Double, ByVal dblSeeing_arsec As Double, ByVal dblExperience As Double) As String ' Returns limiting magnitude for stellar point sources based on characteristics of ' observer, sky brightness, telescope and magnification, etc. ' ' Per ' Schaefer, Bradley E. Oct. 2002. Astronomy and the limits of vision. Vistas in Astronomy. 36(4):311-361 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1993VA.....36..311S ' ' Schaefer, B.E. Feb. 1990. Telescopic Limiting Magnitude. PASP 102:212-229 ' http://adsbit.harvard.edu/cgi-bin/nph-iarticle_query?bibcode=1990PASP..102..212S ' ' Schaefer, B. March 1989, "How faint can you see?", Sky & Telescope, 77(3):332 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1989S%26T....77..332S ' ' Schaefer, B. Nov. 1989, "Your telescope's limiting magnitude", Sky & Telescope 78(5):552 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1989S%26T....78..522S ' ' Schaefer, B. E., Nov. 1989. The Visual Recovery of Halley's Comet. 78(5):525 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1989S%26T....78..525S ' ' Clark, R. April 1994, "How faint can you see?", Sky & Telescope, 87(4):106 ' ' Garstang, R. H. 2000. Limiting visual magnitude and night sky brightness. Memorie della Societa Astronomia Italiana 71:83 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2000MmSAI..71...83G ' ' Garstang, R.H May. 1999. Vision thresholds revisited (abstract). American Astronomical Society, 194th AAS Meeting, #09.16 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1999AAS...194.0916G ' ' Garstang, R.H. April 1999. New Formula for Optimum Magnification and Telescopic Limiting Magnitude. Journal of the Royal Astronomical Society of Canada. 93:80 ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1999JRASC..93...80G ' ' Knoll, H. A. et al. Aug. 1946. Visual thresholds of steady point sources of light in fields of brightness from dark to daylight. 1946JOSA...36..480K ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1946JOSA...36..480K ' ' Blackwell, H. R. Nov. 1946. Constant thresholds of the human eye. 1946JOSA...36..624B ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1946JOSA...36..624B ' ' Hecht, S. et al. 1947. Visibility of Lines and Squares. 1947JOSA...37..500H ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1947JOSA...37..500H ' ' Kumnick, L.S. 1954. Pupillary psychosensory restitution and aging. 1954JOSA...44..735K ' http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1954JOSA...44..735K ' ' Schaefer's algorithm models star brightness (I*) based on the product of a series of ' corrective factors: ' ' I* = L * F_b * F_c * F_t * F_p * F_a * F_r * F_sc * F_c * F_s ' ' where ' ' L = is the response of the human eye to light ' F_b = transforms the binocular eye perception equation of Knoll (1946) and Hecht (1947) to monocular viewing ' F_e = extinction factor ' F_t = transmission factor of the telescope ' F_p = correction factor for the observer's exit pupil size based on the observer's age ' F_m = correction factor for the dimming of the image resulting from magnification ' F_r = extended object size correction factor ' F_sc = correction factor for the Stiles-Crawford effect in the human eye ' F_c = correction factor for the color of the star observed ' F_s = correction factor for the experience of observer ' Test: ' Schaefer, S&T. Nov. 1989 at 525. The O'Meara 1985 recovery of Halley's Comet ' Aperature 24 Magnification 549 NELM 8.2 Age 28 ' TelescopeType Reflector Cleanliness 9 StarColorIndex 0.7 ' Zenith_dist_degrees 20 ExtinctionCoeff 0.1 SeeingDiskDiameter 0.7 ' Experience 9 yields TLM of 19.2 ' ? SchaeferLimitingMag(24,549,8.2,28,"Reflector",9,0.7,20,0.1,0.7,9) ' TLM Dia. Exit Pupil Dia. Eye Pupil ' 19.2683027088087 1.1103825136612 6.73090864766789 ' Initialize variables Dim z, CC, KK, k, x, FC, FD, DS, D, FO, DE, AG, DP, MG, TH, SE, FB As Double Dim FE, FT, FP, FA, FM, FR, CI, CL, KV, FS, MZ, BS, XX, B, I, dblIS, M, EXPER As Double Dim Fl As Double Dim TT, strOut As String ' Translate interface variables to Schaefer variables D = dblAperature_inch ' Aperature inches MG = dblMagnification ' Magnification (Power X) MZ = dblNELM ' Visual magnitude limit (NELM) AG = dblAge ' Age TT = strTelescopeType ' Telescope Type Reflector | Refractor | Schmidt Cassegrain CL = dblCleanliness ' Cleanliness index 0 - 9 Dirty - Clean | Default 4 CI = dblColorIndex ' Color Index (B-V) ' A0 type:CI = 0.0 | B7 type: CI=-0.11 | G2: CI = 0.63 | M2 type star: CI = 1.85 z = dblZenithDist_deg ' Zenith Distance degrees KV = dblExtinctionCoeff ' Extinction Coefficient SE = dblSeeing_arsec ' Seeing - arcseconds EXPER = dblExperience ' Experience 0 - 9 Novice - Expert | Default 4 ' Do Schaefer limiting magnitude computation ' Initial variable set D = D * 25.4 ' Inches input converted to mm ' Set F_t - the transmission factor of the telescope Select Case TT ' Evaluate telescope type Case "Reflector" DS = 0 Fl = (0.99 ^ 4) Case "Refractor" DS = 0.15 * D Fl = (0.88 ^ 2) Case "Schmidt Cassegrain" DS = 0.35 * D Fl = ((0.99 * 0.88) ^ 2) Case Else ' Other values - gracefully handle by resetting to type reflector. ' Debug.Print "Error bad telescope type" DS = 0 Fl = (0.99 ^ 4) End Select ' Run algorithm ' Set constants z = z / 57.296 ' CONVERT TO RADIANS CC = 1.58 * (10 ^ -10) ' NIGHT-VISION CONSTANT KK = 0.0126 ' NIGHT-VISION CONSTANT k = 1.2 * KV ' CORRECT KV TO 5100 A x = 1# / Cos(z) ' AIR MASS ' Determine correction factors FD = 1# - ((DS / D) * (DS / D)) ' ALLOW FOR OBSTRUCTION FO = 0.99 ^ 4 ' EYEPIECE (4 COATED AIR-GLASS SURFACES) DE = 7# * Exp(-AG * AG / 20000#) ' DIAM EYE PUPIL IN MM DP = D / MG ' DIAM OF EXIT PUPIL TH = 2# * SE * MG ' APPARENT DIAMETER OF SEEING DISK (ARC-SEC) CORRECTIONS FB = Sqr(2#) ' BINOCULAR VISION FE = 10 ^ (0.4 * k * x) ' ATMOSPHERIC EXTINCTION ' Runtime error was here due to parsing of denominator FT = 1# / (Fl * FD * FO - (0.01 * (9# - CL))) ' TOTAL TRANSMISSION OF TELESCOPE FP = 1# ' LIGHT OUTSIDE PUPIL If (DE < DP) Then FP = (DP * DP) / DE / DE End If FA = (DE / D) * (DE / D) ' LIGHT-COLLECTING AREA FM = MG * MG ' SPREAD OF SKY PHOTONS FR = 1# ' POINT SOURCE CAN If (TH > 900) Then FR = Sqr(TH / 900) ' APPEAR EXTENDED End If FC = 10 ^ (0.4 * ((CI / 2) - 1#)) ' COLOR OF STAR FS = 1# ' OBSERVER'S SENSITIVITY ' Determine sky brightness ' CALCULATE SKY BRIGHTNESS If (MZ >= (7# - k)) Then BS = 54 ' BEST POSS. SKY BRIGHTNESS FS = 10 ^ (0.4 * (7 - k - MZ)) ' & GOOD EYESIGHT Else XX = 0.2 * (8.68 - k - MZ) ' FS ASSUMED = 1 BS = 39.7 * (((10 ^ XX) - 1#) ^ 2#) ' SKY BRIGHTNESS FOR MZ End If BS = BS * ((z * z * 0.5) + 1#) ' ZENITH HAS DARKEST SKY B = BS / (FB * FT * FP * FA * FM * FC) ' BACKGROUND BRIGHTNESS IN TELESCOPE ' Determine limiting magnitude ' CALCULATE LIMITING MAGNITUDE ' Apply Hecht's equation. See Schaefer (Feb 1990). I = CC * ((1# + Sqr(KK * B)) ^ 2#) ' Hecht (JOSA,v37,p59,1947) ' Apply Schaefer's corrective factor model. See Schaefer (Feb 1990). dblIS = I * FB * FE * FT * FP * FA * FR * FC * FS ' FOR NO SCOPE, NO AIR M = -16.57 - (2.5 * Log(dblIS) / Log(10#)) ' INTENSITY TO V MAG M = M + ((EXPER - 6#) * 0.16) ' EMPIRICAL EXPERIENCE CORRECTION ' Assemble output string strOut = M & " " & DP & " " & DE ' Return result SchaeferLimitingMag = strOut End Function Public Function DiffractionDiskAiryAngularDia_arcsec(ByVal D_mm As Double, ByVal lambda As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the size-diameter of the Airy disk. It is the diameter of the diffraction disk to the ' first diffraction dark ring and contains 84% of a point source's light energy. ' Default wavelength is yellow green 550 nm light. ' Other standard wavelengths commonly used are ' 460nm (Dawes blue-white star assumption and ' 656 nm - Rutten. ' 500 nm - Suitter. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' See Kitchin, C. 2003. 2ed. Telescopes and Techniques. Springer. ISBN 1-85233-725-7 ' at page 33 ' Argyle, R. (ed). 2004. Observing and Measuring Visual Double Stars. Springer. ISBN 1-85233-558-0 ' at page 85, Eq. 10.1 ' Test data ' Berry at p. 7-8. A 200mm (8 inch) scope at 656 nm wavelength is 0.16 arcsecs ' ? DiffractionDiskAiryAngularDia_arcsec (200,656) ' 1.650780048 ' Standard definition of Dawes limit is 4.56 arcsec at 1", which is ' one-half the diameter - the radius of the Airy disk. ' Dawes limit is empirical and based on equal magnitude blue/white stars ' so a wavelength of 460 nm is used. ' ? DiffractionDiskAiryAngularDia_arcsec (25.4,460) ' 9.11463921259842 ' ? 9.11463921259842 / 2# ' 4.55731960629921 ' Notes ' The traditional formula is: ' A_radians = 2.44 * lamba / D ' A_arcsecs = A_radians * 360/2pi() * 60 * 60 ' convert radians to arcsecs ' ' All of the above reduces to: ' A_arcsec = 244/ D_mm ' _conditioned on blue-white stars emitting white-light at 460 nm. Dim A_rad, D_meter, lambda_meter As Double 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (2.44 * lambda_meter) / D_meter DiffractionDiskAiryAngularDia_arcsec = Rad2arcsec(A_rad) End Function Public Function DiffractionDiskAiryLinearDia_microns(ByVal F_ratio As Double, ByVal lambda As Double) As Double ' Receives focal ratio of a scope, a scalar, and light wavelength in nanometers and returns ' the linear size(diameter) of the Airy Disk. It is the diameter of the diffraction disk to the ' first diffraction dark ring and contains 84% of a point source's light energy. ' Default wavelength is yellow green 550 nm light. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' Suiter, H.R. 1994. Star Testing Astronomical Telescopes: A Manual for Optical Evaluation and Adjustment. Willmann-Bell. ISBN 0943396441 ' at page 11 ' Test data ' Suiter 1994 at 11 - a 6" f/8 = 5.5 microns ; 12" f/4 = 5.5 microns ' Suiter generally uses 550nm light in _Star Testing_ ' ? DiffractionDiskAiryLinearDia_microns(4,550) ' 5.368 ' ? DiffractionDiskAiryLinearDia_microns(8,550) ' 10.736 Dim A_meters As Double Dim lambda_meter As Double 'Convert input units to same basis - meters lambda_meter = (lambda * (10# ^ -9#)) A_meters = (2.44 * lambda_meter) * F_ratio DiffractionDiskAiryLinearDia_microns = A_meters * (10 ^ 6) ' microns per meter End Function Public Function DiffractionDiskFWHMAngularDia_arcsec(ByVal D_mm As Double, ByVal lambda As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the size-diameter of the FWHM disk. It is the diameter of the diffraction disk to the ' 50% energy diameter. ' Default wavelength is yellow green 550 nm light. ' Other standard wavelengths commonly used are ' 460nm (Dawes blue-white star assumption and ' 656 nm - Rutten. ' 500 nm - Suitter. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 9 ' Test data ' Dim A_rad, D_meter, lambda_meter As Double 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (1.02 * lambda_meter) / D_meter DiffractionDiskFWHMAngularDia_arcsec = Rad2arcsec(A_rad) End Function Public Function DiffractionDiskFWHMLinearDia_microns(ByVal F_ratio As Double, ByVal lambda As Double) As Double ' Receives focal ratio of a scope, a scalar, and light wavelength in nanometers and returns ' the linear size(diameter) of the Full-Width-Half-Maximum disk. It is the diameter of the diffraction disk to the ' diameter that contains 50% of a point source's light energy. ' Default wavelength is yellow green 550 nm light. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 9 ' Test data ' Berry at 9; f/10 optical system yields a FWHM of 6.7 microns ' ? DiffractionDiskFWHMLinearDia_microns(10,656) ' 6.6912 Dim A_meters As Double Dim lambda_meter As Double 'Convert input units to same basis - meters lambda_meter = (lambda * (10# ^ -9#)) A_meters = (1.02 * lambda_meter) * F_ratio DiffractionDiskFWHMLinearDia_microns = A_meters * (10 ^ 6) ' microns per meter End Function Public Function ResolutionSidgwickTheorAngularRadius_arcsec(ByVal D_mm As Double, ByVal lambda As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the separation resolution test per Sidgwick's theoretical criteria. ' Default wavelength is yellow green 550 nm light, unlike Dawes criteria that is ' based on 460nm. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' See Kitchin, C. 2003. 2ed. Telescopes and Techniques. Springer. ISBN 1-85233-725-7 ' at page 33 ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 38 Eq. d and 45 Eq. f. ' Test data ' Sidgwick's theoretical resolution is typically stated as 5.45 arcsecs for 1" of aperature ' yellow green 550 nm light ' ? ResolutionSidgwickTheorAngularRadius_arcsec(25.4, 550) ' 5.44896909448819 ' Notes ' The traditional formula is: ' A_radians = 1.22 * lamba / D ' A_arcsecs = A_radians * 360/2pi() * 60 * 60 ' convert radians to arcsecs ' ' All of the above reduces to: ' A_arcsec = 138/ D_mm ' _conditioned on yellow green 550nm light. Dim A_rad, D_meter, lambda_meter As Double 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (1.22 * lambda_meter) / D_meter ResolutionSidgwickTheorAngularRadius_arcsec = Rad2arcsec(A_rad) End Function Public Function ResolutionAiryAngularRadius_arcsec(ByVal D_mm As Double, ByVal lambda As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the separation resolution test of the radius of the Airy disk. It is the radius of the diffraction disk to the ' first diffraction dark ring and contains 42% of a point source's light energy. ' The Airy disk radius at 460nm light is the same as the Dawes Limit resolution test. ' Default wavelength is blue/white 460 nm light. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' See Kitchin, C. 2003. 2ed. Telescopes and Techniques. Springer. ISBN 1-85233-725-7 ' at page 33 ' Argyle, R. (ed). 2004. Observing and Measuring Visual Double Stars. Springer. ISBN 1-85233-558-0 ' at page 86 ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 45, Eq. g. ' Test data ' Airy's ratio is typically stated as 4.56 arcsecs for 1" of aperature ' Using blue/white light at 460 nm ' ? ResolutionAiryAngularRadius_arcsec(25.4, 460) ' 4.55731960629921 ' Notes ' The traditional formula is: ' A_radians = 1.22 * lamba / D ' A_arcsecs = A_radians * 360/2pi() * 60 * 60 ' convert radians to arcsecs ' ' All of the above reduces to, per Kitchin, at 500nm: ' A_arcsec = 122/ D_mm ' _conditioned_ on using 500nm as the light wavelength, e.g. - ' ? ResolutionAiryAngularRadius_arcsec(150, 500) ' 0.838811 ' 122 / 150# ' 0.813333333333333 ' Other web variations, that are based on slightly different wavelenghts, include: ' A_arcsec = 127/ D_mm; A_arcsec = 126/ D_mm ' Other authors, using 460nm light, summarize this rule as: ' A_arcsec = 116/ D_mm ' _conditioned_ on using 460nm as the light wavelength, e.g. - ' ? ResolutionAiryAngularRadius_arcsec(25.4, 460) ' 4.55731960629921 ' 116 / 25.4 ' 4.56692913385827 ' Other web published variants include: ' A_arcsec = 115.8/ D_mm Dim A_rad, D_meter, lambda_meter As Double 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (1.22 * lambda_meter) / D_meter ResolutionAiryAngularRadius_arcsec = Rad2arcsec(A_rad) End Function Public Function ResolutionAiryLinearRadius_microns(ByVal F_ratio As Double, ByVal lambda As Double) As Double ' Receives focal ratio of a scope, a scalar, and light wavelength in nanometers and returns ' a resolution test in the form of the linear size(radius) of the Airy Disk. It is the radius of the diffraction disk to the ' first diffraction dark ring and contains 42% of a point source's light energy. ' Default wavelength is blue/white 460 nm light. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' The Airy disk test is strictly only valid for two matched blue/white stars ' e.g. - blue/white light near 460nm. ' Test data - Dim A_meters As Double Dim lambda_meter As Double 'Convert input units to same basis - meters lambda_meter = (lambda * (10# ^ -9#)) A_meters = (1.02 * lambda_meter) * F_ratio ResolutionAiryLinearRadius_microns = A_meters * (10 ^ 6) ' microns per meter End Function Public Function ResolutionDawesLimit_arcsec(ByVal D_mm As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the separation resolution test of the Dawes Limit where 460nm blue/white light is used. ' It is the radius of the diffraction disk to the ' first diffraction dark ring and contains 42% of a point source's light energy. ' The Airy disk radius at 460nm light is the same as the Dawes Limit resolution test. ' Only valid if wavelength is blue/white 460 nm light. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' See Kitchin, C. 2003. 2ed. Telescopes and Techniques. Springer. ISBN 1-85233-725-7 ' at page 33 ' Argyle, R. (ed). 2004. Observing and Measuring Visual Double Stars. Springer. ISBN 1-85233-558-0 ' at page 86 ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 45, Eq. g. ' Dawes limit is usually stated as 4.5" per inch at 460nm ' Other authors, using 460nm light, summarize this rule as: ' A_arcsec = 116/ D_mm ' _conditioned_ on using 460nm as the light wavelength, e.g. - ' ? ResolutionAiryAngularRadius_arcsec(25.4, 460) ' 4.55731960629921 ' 116 / 25.4 ' 4.56692913385827 ' Other web published variants include: ' A_arcsec = 115.8/ D_mm; Covington's _Astrophotography_: 115 / D_mm using 455nm Dim A_rad, D_meter, lambda_meter, lambda As Double ' Force light wavelength setting lambda = 460# ' 460 nm 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (1.22 * lambda_meter) / D_meter ResolutionDawesLimit_arcsec = Rad2arcsec(A_rad) End Function Public Function ResolutionRayleighLimit_arcsec(ByVal D_mm As Double, ByVal lambda As Double) As Double ' Receives aperature diameter in millimeters and light wavelength in nanometers and returns ' the separation resolution test of the radius of the Airy disk as seen with ' yellow-green light at 550 nm. It is the radius of the diffraction disk to the ' first diffraction dark ring and contains 42% of a point source's light energy. ' Default wavelength is green/yellow 550 nm light. ' The human eye is most sensitve at 550nm. ' Photographic film is most sensitive at 400nm. (per Covington _Astrophotography per Amateurs_). ' Rayleigh's limit is the same as the Dawes limit - it's just at a different wavelength - ' closer to the usual Johnson V visual band. ' Source: Berry, R & Burnell, J. 2005. 2d. Handbook of Astronomical Processing. (HAIP). Willman-Bell. ' at page 8 ' See Kitchin, C. 2003. 2ed. Telescopes and Techniques. Springer. ISBN 1-85233-725-7 ' at page 33 ' Argyle, R. (ed). 2004. Observing and Measuring Visual Double Stars. Springer. ISBN 1-85233-558-0 ' at page 86 ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 45, Eq. g. ' Test data ' The traditional Rayleigh test is stated as 5.45" per for 1" of aperature ' Using green-yellow light at 550 nm ' ? ResolutionRayleighLimit_arcsec(25.4,550) ' 5.44896909448819 ' Covington's _Astrophotography for Amateurs_ at p. 83 ' A 100mm scope at 550nm should have a Rayleigh limit of 1.4 arcsecs ' ResolutionRayleighLimit_arcsec(100, 550) ' 1.38403815 ' Notes ' The traditional formula is: ' A_radians = 1.22 * lamba / D ' A_arcsecs = A_radians * 360/2pi() * 60 * 60 ' convert radians to arcsecs ' All of the above reduces to at 550nm: ' A_arcsec = 138/ D_mm ' ? ResolutionRayleighLimit_arcsec(25.4,550) ' 5.44896909448819 ' 138 / 25.4# ' 5.43307086614173 Dim A_rad, D_meter, lambda_meter As Double 'Convert input units to same basis - meters D_meter = D_mm / 1000# lambda_meter = (lambda * (10# ^ -9#)) A_rad = (1.22 * lambda_meter) / D_meter ResolutionRayleighLimit_arcsec = Rad2arcsec(A_rad) End Function Public Function ResolutionAngular2Linear_lpmm(ByVal ResolutionAngular_arcsec As Double, ByVal Fl_mm As Double) As Double ' Converts angular resolution into linear resolutions in form of line pairs per millimeter ' Use Rayleigh's limit as input. For film, film is more sensitive at 400nm per Covington. ' Source: Rutten (2002) Telescope Optics. at 211-212. ' Test data ' Rutten at 212 A 1000 fl telescope with 1.14" resolution yields 180 lp/mm ' ? ResolutionAngular2Linear_lpmm(1.14,1000) ' 180.934210526316 ResolutionAngular2Linear_lpmm = 206265 / (ResolutionAngular_arcsec * Fl_mm) End Function Public Function ResolutionLinear2Angular_arcsec(ByVal ResolutionLinear_lpmm As Double, ByVal Fl_mm As Double) As Double ' Converts linear resolution into angular resolutions in form of line pairs per millimeter to arcsecs ' Use Rayleigh's limit as input. For film, film is more sensitive at 400nm per Covington. ' Source: Rutten (2002) Telescope Optics. at 211-212. ' Test data ' Rutten at 212 A 1000 fl telescope with 1.14" resolution yields 180 lp/mm ' ? ResolutionLinear2Angular_arcsec(180, 1000) ' 1.14591666666667 ResolutionLinear2Angular_arcsec = 206265 / (ResolutionLinear_lpmm * Fl_mm) End Function Public Function ResolutionLinearTheor_lpmm(ByVal F_ratio As Double, ByVal lambda As Double) As Double ' Receives focal ratio of a scope, a scalar, and light wavelength in nanometers and returns ' the linear resolution in terms of a linear height millimeters. ' Default wavelength is yellow green 550 nm light. ' Photographic linear resolution is based on 400nm light per Covington. ' Rutten's method does not agree with Covington (1999) at 84-85 and Table 6.3. ' Source: Rutten (2002) Telescope Optics. at 211-212. ' Test data ' Rutten at 212 A f/10 telescope at 550nm yields 180 lp/mm ' ? ResolutionLinear_lpmm(10, 550) ' 181.818181818182 Dim A_meters As Double Dim lambda_meter As Double 'Convert input units to same basis - meters lambda_meter = (lambda * (10# ^ -9#)) A_meters = 1# / (F_ratio * lambda_meter) ResolutionLinearTheor_lpmm = A_meters / (10 ^ 3) ' convert to millmeters End Function Public Function ResolutionLinearFilmPractical_lpmm(ByVal F_ratio As Double) As Double ' Receives focal ratio of a scope, a scalar, and returns ' the linear resolution in terms of a line pairs / mm. ' Default wavelength is yellow green 550 nm light. ' Photographic linear resolution is based on 400nm light per Covington. ' Is Covington's empirical formula. ' Source: Covington (1999). Astrophotography for the Amateur. At. 84. ' See Table 6.3 at 85. ' Test data ' Covington - an f/4 scope has a practical resolution power of 250 lp/mm ' ? ResolutionLinearFilmPractical_lpmm(4) ' 250 ResolutionLinearFilmPractical_lpmm = 1000# / F_ratio ' is empircally estimate End Function Public Function ScaleDetector2LinePairs_per_mm(ByVal EFl_mm_tmp As Double, ByVal Detector_maj_axis_microns_tmp As Double) As String ' Receives detector size and returns a practical lines pairs per mm ' Also receives the focal length and returns image scale in terms of arcsecs covered by a line pair ' Used to estimate photographic resolution ' Source: Covington (1999) at 219 for lines per millimeter ' Covington uses line pairs per mmm and lines per mm interchangeably ' There are two pixels for each line pair, so the linear height of two pixels ' gives the linear height of a line. This is converted into a a millimeter basis. ' After ScalePixel2Arcsecs ' Sample - a 2000 fl scope with 6 micron detectors ' ? ScaleLinePairs_per_mm(2000, 6) ' 83.3333333333333 1.23768 ' Working variables Dim Detector_maj_axis_microns, Detector_maj_axis_meter, Detector_maj_axis_mm As Double Dim Fl_mm, Fl_meter As Double Dim A, B As Double Dim LinePairs_mm As Double Dim LinePairs_arcsec As Double 'Get and convert input units to basis Fl_mm = EFl_mm_tmp Detector_maj_axis_microns = Detector_maj_axis_microns_tmp Detector_maj_axis_mm = Detector_maj_axis_microns * (10 ^ -3) ' Double the detector size to create line pair Detector_maj_axis_mm = 2 * Detector_maj_axis_mm ' Find lp/mm - Take the inverse of the detector size to obtain the number of line pairs per mm LinePairs_mm = 1# / Detector_maj_axis_mm ' Find the arcsecs covered by a line pair B = ScaleHeight2Image_decdeg(Fl_mm, Detector_maj_axis_mm) LinePairs_arcsec = B * 3600# ' Convert to arcsecs ' Return result ' Line pairs per mm Arcsecs covered by a line pair ScaleDetector2LinePairs_per_mm = CStr(LinePairs_mm) & " " & CStr(LinePairs_arcsec) End Function Public Function ScaleLinePairs_per_mm2Detector(ByVal EFl_mm_tmp As Double, ByVal LinePairs_mm_tmp As Double) As String ' Receives a film type detector size and returns the comparable for detector size ' in square microns ' Used to equate film resolution characteristics into a common based detector size ' as CCD pixels. With common detector basis, exposure computations can be done ' with a single equation set ' Used to estimate photographic resolution, but note this function is termed ' "scale" not "resolution". Resolution involves characteristics of the ' diffraction disk size of incoming light and not detector size alone. ' Usually, expert books refer to this estimator as "resolution." ' Source: After ScaleDetector2LinePairs_per_mm; ScalePixel2Arcsecs ' Dependencies: ScalePixel2Arcsecs ' Sample - a 2000 fl scope at 83 line pairs per mm ' ? ScaleLinePairs_per_mm2Detector(2000, 83.333333) ' 6.000000024 0.61884000247536 ' Working variables Dim Detector_maj_axis_microns, Detector_maj_axis_meter, Detector_maj_axis_mm As Double Dim Fl_mm, Fl_meter As Double Dim A, B As Double Dim LinePairs_mm As Double Dim LinePairs_arcsec As Double 'Get and convert input units to basis Fl_mm = EFl_mm_tmp LinePairs_mm = LinePairs_mm_tmp ' Find line size in mm - Take the inverse of the line pairs to obtain size of the line in mm Detector_maj_axis_mm = 1# / LinePairs_mm ' Halve the line size find the detector size Detector_maj_axis_mm = Detector_maj_axis_mm / 2 ' Convert the detector size in mm or microns Detector_maj_axis_microns = Detector_maj_axis_mm / (10 ^ -3) ' Find the arcsecs covered by the detector B = ScalePixel2Arcsecs(Fl_mm, Detector_maj_axis_microns) ' Return result ' Detector (pixel) size microns Arcsecs per pixel ScaleLinePairs_per_mm2Detector = CStr(Detector_maj_axis_microns) & " " & CStr(B) End Function Public Function ScaleImage2HeightLinear(ByVal Fl_mm As Double, ByVal lambda_decdegs As Double) As Double ' Receives focal length of a scope and object true size in decimal degrees and returns ' the linear size (diameter) of an image, usually at prime focus, in millimeters. ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 52, Eq. a, 346-347 ' Test data ' Sidgwick at 52 - a 4ft fl solar scope produces an image of a 0.53 deg Sun of 0.4" ' ? ScaleImage2HeightLinear (4 * 12 * 25.4, 0.53) ' ? 11.2778987473669 / 25.4 ' 0.444011761707358 ' Covington (1999) at 75, Table 6.2 1250 fl Moon 1800" is 11mm on film ' ? ScaleImage2HeightLinear(1250, 1800 / 3600) ' 10.9083078249646 ' Covington (1999) at 75, Table 6.2 3000 fl Jupiter 40" is 0.6mm on film ' ? ScaleImage2HeightLinear(3000, 40 / 3600) ' 0.581776417331443 ' Dim A_meters As Double Dim lambda_radians, Fl_meter As Double 'Get and convert input units to basis Fl_meter = Fl_mm / 1000# lambda_radians = Deg2Rad(lambda_decdegs) A_meters = lambda_radians * Fl_meter ScaleImage2HeightLinear = A_meters * 1000 ' meter convert to mm End Function Public Function ScaleHeight2Image_decdeg(ByVal Fl_mm As Double, ByVal Height_mm_tmp As Double) As Double ' Receives focal length of a scope and detector or image height at prime focus and returns ' the angular size (diameter) of the true object in decimal degrees. ' Primarily used to determine detector TFOV sizes - including CCD pixel and film sizes. ' Sidgwick1971a: Sidgwick, J.B. 1971 (3ed). Amateur Astronomer's Handbook. Dover. 1971QB44.S558...… ' at 52, Eq. a, 346-347 ' Covington (1999). Astrophotography for the Amateur at 75 ' Berry, HAIP (2005, 2d) at 29. ' Test data ' Sidgwick at 52 - a 4ft fl solar scope produces an image of a 0.53 deg Sun of 0.4" ' ? ScaleHeight2Image_decdeg(4 * 12 * 25.4, 11.2778987473669) ' 0.530039040538159 ' Berry, HAIP (2005, 2d) at 29. The true field of a 17 micron detector is 3.5" ' Note 17 microns is 17^10-3 millimeters - the height input units ' ? ScaleHeight2Image_decdeg(1000, 17 * 10 ^ -3) ' 0.0009741 ' ? 0.0009741 * 3600 ' 3.50676 Dim Height_meter As Double Dim Fl_meter As Double 'Get and convert input units to basis Fl_meter = Fl_mm / 1000# Height_meter = Height_mm_tmp / 1000# ScaleHeight2Image_decdeg = 57.3 * (Height_meter / Fl_meter) ' meter convert to mm End Function Public Function ScalePhysicalDistArcsec(ByVal ObjDist_tmp As Double, ByVal ObjDiameter_arcsec_tmp As Double) As Double ' Receives the distance to an object in astronomical units (pc, kpc, Mpc, light-years, etc) and ' its major axis diameter and returns its physical size in those same units. Uses the ' sine rule. ' Source: ' Test data ' Working variables Dim A As Double Dim ObjDist As Double Dim ObjDiameter_arcsec As Double Dim ObjDiameter_rads As Double 'Get and convert input units to basis ObjDist = ObjDist_tmp ObjDiameter_arcsec = ObjDiameter_arcsec_tmp ' Convert to degrees and radians A = Deg2Rad((ObjDiameter_arcsec / 3600#)) ' Return result ScalePhysicalDistArcsec = ObjDist * A End Function Public Function ScalePixel2Arcsecs(ByVal EFl_mm_tmp As Double, ByVal Detector_maj_axis_microns_tmp As Double) As Double ' Receives focal length and detector size and returns the pixel scale of the image. ' Used for computing pixel scale of image ' Source: Berry, HAIP (2005) at 29 ' Dependencies: ScaleHeight2Image_decdeg ' Test data ' Berry, HAIP (2005, 2d) at 29. The true field of a 17 micron detector is 3.5" ' Note 17 microns is 17^10-3 millimeters - the height input units ' ? ScalePixel2Arcsecs(1000, 17) ' 3.50676 ' Working variables Dim Detector_maj_axis_microns, Detector_maj_axis_meter, Detector_maj_axis_mm As Double Dim Fl_mm, Fl_meter As Double Dim A, B As Double 'Get and convert input units to basis Fl_mm = EFl_mm_tmp Detector_maj_axis_microns = Detector_maj_axis_microns_tmp ' Detector_maj_axis_meter = Detector_maj_axis_microns * (10 ^ -6#) ' Detector_maj_axis_mm = Detector_maj_axis_meter / 1000# Detector_maj_axis_mm = Detector_maj_axis_microns * (10 ^ -3) B = ScaleHeight2Image_decdeg(Fl_mm, Detector_maj_axis_mm) ScalePixel2Arcsecs = B * 3600# ' Convert to arcsecs End Function Public Function ScaleEFL2Plate(ByVal EFl_mm_tmp As Double) As Double ' Receives effective focal length of a telescopic system and returns ' the plate scale of the image made at prime focus ' Source: Edmund Scientific. 1998. Photography with Your Telescope at 32. ' Test data ' Edmund at 32: A scope with 7.4 inches EFL has a plate scale of 1100" per mm ' ? ScaleEFL2Plate(+7.4*25.4) ' 1097.38774207278 ' Create working variables and get input Dim EFL_mm As Double EFL_mm = EFl_mm_tmp ScaleEFL2Plate = 206265# / EFL_mm ' return plate scale in arcsecs End Function Public Function ScalePlate2EFL_mm(ByVal PlateScale_arcsec_mm_tmp As Double) As Double ' Receives effective focal length of a telescopic system and returns ' the plate scale of the image made at prime focus ' Source: Edmund Scientific. 1998. Photography with Your Telescope at 32. ' Test data ' Edmund at 32: A scope with 7.4 inches EFL has a plate scale of 1100" per mm ' ? ScalePlate2EFL_mm(1100) ' ? 187.513636363636 / 25.4 ' 7.38242662848961 ' Create working variables and get input Dim PlateScale_arcsec_mm As Double PlateScale_arcsec_mm = PlateScale_arcsec_mm_tmp ScalePlate2EFL_mm = 206265# / PlateScale_arcsec_mm ' return plate scale in arcsecs End Function Public Function ScaleMoonObjSizeArcsec2Km(ByVal Obj_size_arcsec_tmp As Double, ByVal Obj_Sellong_tmp As Double, ByVal Obj_Sellat_tmp As Double, ByVal MoonDist_km_tmp As Double) As String ' Receives lunar feature arcsec size, object selenographic long and lat, lunar distance, and ' returns object longitude and latitude size in kilometers ' Source: ' Astronomy Dept. of the U. of Michigan. 2000. Lunar Features and ' Mountain Heights. (Webpage) ' http://www.astro.lsa.umich.edu/Course/Labs/lunar_features/lunar_intro.html ' Test data ' ? MoonSizeArcsec2Km(60,0,0,384500) ' 111.846411170097 111.846411170097 ' Working variables Dim Obj_size_arcsec, Obj_Sellong, Obj_Sellat, MoonDist_km As Double Dim Obj_Sellong_rad, Obj_Sellat_rad As Double Dim Objsize_on_Moonkm As Double Dim Objsize_on_MoonLatkm As Double Dim Objsize_on_MoonLongkm As Double ' Get input Obj_size_arcsec = Obj_size_arcsec_tmp Obj_Sellong = Obj_Sellong_tmp Obj_Sellong_rad = Deg2Rad(Abs(Obj_Sellong)) Obj_Sellat = Obj_Sellat_tmp Obj_Sellat_rad = Deg2Rad(Abs(Obj_Sellat)) MoonDist_km = MoonDist_km_tmp ' Evaluate ' What is the physical size of a one arcsec TFOV on the Moon ' if it where located at selenographic 0 long 0 lat Objsize_on_Moonkm = (1 / 206265#) * MoonDist_km ' What is the physical size of a one arcsec TFOV on the Moon ' if it where located at selenographic N long N lat Objsize_on_MoonLatkm = Objsize_on_Moonkm / Cos(Obj_Sellat_rad) Objsize_on_MoonLongkm = Objsize_on_Moonkm / Cos(Obj_Sellong_rad) ' Scale sizes to observed object arcsec size Objsize_on_MoonLatkm = Objsize_on_MoonLatkm * Obj_size_arcsec Objsize_on_MoonLongkm = Objsize_on_MoonLongkm * Obj_size_arcsec ' Return result ScaleMoonObjSizeArcsec2Km = CStr(Objsize_on_MoonLongkm) & " " & CStr(Objsize_on_MoonLatkm) End Function Public Function TelesDepthofFocus(ByVal F_ratio_tmp As Double, ByVal lambda As Double) As Double ' Receives focal ratio and returns in millimeters, the current scope's depth of focus ' using Sinnott's formula Dim F_ratio, FocuserTravel As Double Dim lambda_meter As Double ' Get and convert input units to same basis - meters lambda_meter = (lambda * (10# ^ -9#)) F_ratio = F_ratio_tmp FocuserTravel = 4 * lambda_meter * (F_ratio ^ 2) TelesDepthofFocus = FocuserTravel * (10 ^ 3) ' Convert to millimeters End Function Public Function TelesDiffractionLimitedFieldSize(ByVal F_ratio_tmp As Double) As Double ' Receives focal ratio and returns ' the size of the diffraction free limited field in millmeters ' per Sinnott's equation ' Source: Sinnott, Roger W. May 1991. Focus and Collimination: How Critical? Sky&Telescope pp528-531. Dim F_ratio As Double Dim A As Double ' Get input F_ratio = F_ratio_tmp A = 0.0007 * (F_ratio ^ 3) ' field-width inches - Sinnott's empirical equation TelesDiffractionLimitedFieldSize = A * 25.4 ' convert to millimeters End Function Public Function TelesComaFreeFieldSize_mm(ByVal F_ratio_tmp As Double) As Double ' Receives focal ratio and returns ' the size of the coma free limited field in arcsecs ' per Covington's equation ' Source: Covington 1999. Astrophotography for Amateurs. At 72. ' Test data ' Covington Table 6.1 at p. 73 ' f/5 =12.5; f/8 =32; f/10 = 50 ' ? TelesComaFreeFieldSize_mm(5) ' 12.5 ' ? Print TelesComaFreeFieldSize_mm(10) ' 50 Dim F_ratio, A As Double ' Get and convert input units to same basis - meters F_ratio = F_ratio_tmp A = (F_ratio ^ 2) / 2 ' field-width inches - Sinnott's empirical equation TelesComaFreeFieldSize_mm = A ' convert to millimeters End Function Public Function ScaleHalfAngleFormula(ByVal Objwidth_tmp As Double, ByVal ObjDist_tmp As Double) As Double ' Receives distance to object and object size in same units and ' returns the angular size of the object in arcsecs ' Object distance is usually the focal length of an eyepiece or objective ' Imaging scaling functions are simplified version of the half-angle ' equation used for long distances, i.e. - 1% accurate above 1000mm ' Input: Object distance and width must be in same units ' Source: Covington 1999. Astrophotography for Amateurs. At 74. ' Test data ' Covington at 75, Table 6.2 ' The Moon has a linear height of 3.6mm in a 400mm fl scope. ' What is it's angular size in arcsecs? ' ? ScaleHalfAngleFormula(3.6, 400) ' 1856.37072578913 ' Usage example: ' What is the field of view of a telescope with a 1200mm fl ' and a tube-size of 1 1/4"? ' ? ScaleHalfAngleFormula(1.25*25.4,1200) ' ? 5457.10466263391 / 3600 ' in degrees ' 1.5158624062872 ' degrees ' Working variables Dim Objwidth, ObjDist, A, B, C As Double ' Get input and convert Objwidth = Objwidth_tmp ObjDist = ObjDist_tmp ' Compute A = Objwidth / (2 * ObjDist) B = 2 * Atn(A) C = Rad2Deg(B) ScaleHalfAngleFormula = C * 3600# ' Convert output to arcsecs End Function Public Function AtmosphericRefractionSimple_arcsecs(ByVal ObjAltitude_decdeg_tmp) As Double ' Receives the altitude of an object (not the zenithal distance) ' and returns a simple refraction estimate in arcsecs ' Limitations: Good to zenithal distance 70 deg; altitude 30 deg ' Is based on the simple estimate formula and not the more complex Comstock formula ' including temperature and pressure. ' Assumes a constant k of 58.2, which is for sea level ' Source: Kitchin (2003). Telescopes and Techniques. At 105. ' Sidgwick. (1971 3ed). Amateur Astronomer's Handbook. At 453 ' Test data ' Sidgwick. At z = 40, Pulkova gives 47.9"; Bessel 48.0" ' ? AtmosphericRefractionSimple_arcsecs(50) ' 48.8355985345177 ' Difference is attributable to higher k value used by Kitchin 58.2 vs. Sidgwick 57.0 ' Working variables Dim ObjAltitude_decdeg, z As Double ' Get input ObjAltitude_decdeg = ObjAltitude_decdeg_tmp ' Compute and return result z = 90# - ObjAltitude_decdeg z = Deg2Rad(z) AtmosphericRefractionSimple_arcsecs = 58.2 * Tan(z) End Function Public Function AtmosphericRefractionComplex_arcsecs(ByVal ObjAltitude_decdeg_tmp, ByVal TemperatureF_tmp As Double, ByVal BarometricPressure_in_tmp As Double) As Double ' Receives the altitude of an object (not the zenithal distance), temperature and pressure, ' and returns a simple refraction estimate in arcsecs ' Limitations: Good to zenithal distance 70 deg; altitude 30 deg ' Is based on the a more complex Comstock formula including temperature and pressure. ' Source: Sidgwick. (1971 3ed). Amateur Astronomer's Handbook. At 453 ' Test data ' Sidgwick. At z = 40, Pulkova gives 47.9"; Bessel 48.0" ' ? AtmosphericRefractionComplex_arcsecs(50,32,29) ' 48.6183194836803 ' Working variables Dim ObjAltitude_decdeg, TemperatureF, BarometricPressure_in Dim z, A As Double ' Get input ObjAltitude_decdeg = ObjAltitude_decdeg_tmp TemperatureF = TemperatureF_tmp BarometricPressure_in = BarometricPressure_in_tmp ' Compute z = 90# - ObjAltitude_decdeg z = Deg2Rad(z) A = (983 * BarometricPressure_in) / (460 + TemperatureF) ' Return result AtmosphericRefractionComplex_arcsecs = A * Tan(z) End Function Public Function TelesPeriodicErrorCumulative_arcsecs(ByVal Timesecs_tmp As Double, ByVal Err_Velocity_arcsecsec_tmp As Double, ByVal Err_Acceleration_arcsecsec_tmp As Double) As Double ' Gets periodic errors in velocity and accelerations and total time ' and returns an estimate of the total periodic error at time t ' Source: Berry (2005 2ed) HAIP at 96-97 ' Test data ' An mount with an error of velocity of 0.03 arcsecs/sec and acceleration of 0.0008 ' ? TelesPeriodicErrorCumulative_arcsecs(33,0.03,0.0008/15#) ' 1.01904 ' Working Variables Dim A As Double Dim Err_Velocity_arcsecsec, Err_Acceleration_arcsecsec, Timesecs As Double ' Get input Timesecs = Timesecs_tmp Err_Velocity_arcsecsec = Err_Velocity_arcsecsec_tmp Err_Acceleration_arcsecsec = Err_Acceleration_arcsecsec_tmp ' Compute A = Err_Velocity_arcsecsec * Timesecs + (0.5 * Err_Acceleration_arcsecsec * (Timesecs ^ 2)) ' Return result TelesPeriodicErrorCumulative_arcsecs = A End Function Public Function astrometryGalacticArmEquiAngularSpiral(ByVal theta_deg_tmp As Double, ByVal pitch_angle_deg_tmp As Double, ByVal R_arm0_kparsecs_tmp As Double, ByVal z_warp_angle_dec_tmp As Double, ByVal R_0_kparsecs_tmp As Double, ByVal z_0_kparsecs_tmp As Double) As String ' Generates orthogonal coordinates of the Milky Way Spiral Arms. ' Receives: ' theta - the galactic-centric rotation angle in degrees of this Arm ' beta - pitch_angle - the pitch angle of this Arm ' R_0arm_kparsecs - the galactic-center distance to the Arm at theta = 0 ' R_0kparsecs - the Sol-galactic distance (used to adjust to Sol-centered ortho coordinates ' gamma - z-warp-angle is the rotation angle between the mean galactic plane and the Arm. ' z_0 is the offset between Sol and mean galactic plane = 30 parsecs ' Returns ' r = distance from galactic center to spiral point ' x,y,z ortho in Sol-based galactic coordinate system ' x,y,z ortho measured from the galactic center ' Limitations: ' z-warp-angle is valid only between theta 0 and 270 degrees. ' After 270 degrees, z_warp result is forced to z=0 ' Test data: Perseus Arm at -10.3 kpc from galactic center and -1.7 kpc from Orion Arm ' ? phyGalacticArmEquiAngularSpiral(180, 12, 5.318285035, 8.5, 3, -0.03) ' ? 10.268947300293 -10.268947300293 1.25754181110059E-15 -0.03 -1.76894730029299 1.25754181110059E-15 -0.03 ' Source: ' General form of the equiangular spiral is: ' r = rho_0 * e^(theta*beta) where beta is the pitch-angle for the spiral. ' Once the current r for the spiral is known, the ortho coordinates are found with: ' x(theta) = r * cos(theta) ' y(theta) = r * sin(theta) ' For the z-warp angle, it is similar to reducing ecliptic coordinates using the 23.5 Earth ecliptic angle ' z(theta) = sin(theta) * sin(z_warp) ' Then the coordinate is moved relative to the Earth's position above the mean galactic plane. ' Create working variables Dim A, B, C, R As Double Dim X_ortho, Y_ortho, Z_ortho, X_ortho_galcen, Y_ortho_galcen, Z_ortho_galcen Dim strOut As String Dim theta_deg, beta_deg, pitch_angle_deg, R_arm0_kparsecs, R_0_kparsecs, z_warp_angle_dec, z_0_kparsecs As Double Dim theta_rads, beta_rads, pitch_angle_rads, z_warp_angle_rads As Double ' Get the input theta_deg = theta_deg_tmp pitch_angle_deg = pitch_angle_deg_tmp R_arm0_kparsecs = R_arm0_kparsecs_tmp R_0_kparsecs = R_0_kparsecs_tmp z_warp_angle_dec = z_warp_angle_dec_tmp z_0_kparsecs = z_0_kparsecs_tmp ' Convert to radians theta_rads = Deg2Rad(theta_deg) pitch_angle_rads = Deg2Rad(pitch_angle_deg) z_warp_angle_rads = Deg2Rad(z_warp_angle_dec) ' Find r (distance spiral origin to point on spiral at angle theta R = R_arm0_kparsecs * Exp(theta_rads * pitch_angle_rads) ' Find the galactic center x-y ortho coordinates in the galactic centered coordinate system X_ortho_galcen = R * Cos(theta_rads) Y_ortho_galcen = R * Sin(theta_rads) ' Programming note: Handle runtime error limit in z-warp sin equation (270 degrees) by forcing ' result to zero If theta_deg <= 270 Then ' Y_ortho_galcen scaled to tilt of arm plane Z_ortho_galcen = Y_ortho_galcen * Sin(z_warp_angle_rads) Else Z_ortho_galcen = 0 End If ' Add the z-offset for Sol position above galactic plane Z_ortho = Z_ortho_galcen + z_0_kparsecs ' Find the x-y ortho coordinates in the Sol-centered galactic coordinate sytem X_ortho = X_ortho_galcen + R_0_kparsecs Y_ortho = Y_ortho_galcen Z_ortho = Z_ortho_galcen ' Return the result ' Distance gc->obj Sol-based galactic coordinates Galactic center based ortho coordinates astrometryGalacticArmEquiAngularSpiral = CStr(R) & " " & CStr(X_ortho) & " " & CStr(Y_ortho) & " " & CStr(Z_ortho) & " " & CStr(X_ortho_galcen) & " " & CStr(Y_ortho_galcen) & " " & CStr(Z_ortho_galcen) End Function Public Function PositionAngle(ByVal Outlier_RA1_dec_tmp As Double, ByVal Outlier_Dec1_dec_tmp As Double, ByVal CentralPoint_RA2_dec_tmp As Double, ByVal CentralPoint_Dec2_dec_tmp As Double) As String ' Meeus Chapter 17 function at p. 116 ' Equatorial positions of two objects, what is their position angle ' Similar to the terresterial navigation "bearing" ' Position Angle is measured counterclockwise from North ' (Terresterial bearings are measured clockwise from North) ' Receives to two RA and Dec positions and returns a string containing - ' PositionAngle ' Terresterial navigation To Bearing ' Terresterial navigation From Bearing ' Precision - Is untested for very small angles ' Test and Usage: ' Assume the celestial sphere viewed from the inside, ' with RA1h, Dec0d as the origin surrounded by stars ' at RA1h and +-Dec15d intervals (RA1h and RA0h). ' Converted to degrees, the points on the celestial equator are ' from left to right: 30,0;15,0 and 0,0. ' ? PositionAngle(15, 15, 15, 0) ' Direct above ' 0 0 180 ' ? PositionAngle(30, 15, 15, 0) ' 1st Quad ' 44.0070271956363 315.992972804364 135.992972804364 ' ? PositionAngle(30, 0, 15, 0) ' 90 degrees - on the left ' 90 270 90 ' ? PositionAngle(30, -15, 15, 0) ' 2nd Quad ' 134.007027195636 225.992972804364 45.9929728043637 ' ? PositionAngle(15, -15, 15, 0) ' 180 degrees Directly below ' 180 180 0 ' ? PositionAngle(0, -15, 15, 0) ' 3rd Quad ' 224.007027195636 135.992972804364 315.992972804364 ' ? PositionAngle(0, 0, 15, 0) ' 270 degrees - on the right ' 270 90 270 ' ? PositionAngle(0, 15, 15, 0) ' 4th Quad ' 314.007027195636 45.9929728043637 225.992972804364 ' Define variables Dim RA1_dec As Double Dim Dec1_dec As Double Dim RA2_dec As Double Dim Dec2_dec As Double Dim A, B, C, D, E, f, g, H, x, y As Double ' Temporary working variables Dim tanP As Double ' Tangent of the position angle Dim Angle_result_dec As Double Dim dblPositionAngle_result As Double Dim dblBearingTo_result As Double Dim dblBearingFrom_result As Double Dim Angle_dec As Double ' Temporary working variables Dim Angle_dec_adj ' Temporary working variables ' get equilatorial or galactic coordinates of two objects RA1_dec = Outlier_RA1_dec_tmp Dec1_dec = Outlier_Dec1_dec_tmp RA2_dec = CentralPoint_RA2_dec_tmp Dec2_dec = CentralPoint_Dec2_dec_tmp ' Evaluate Meeus Eq. at p. 116 ' Evaluate the numerator ' Find the difference between the two right ascensions ' convert to sine A = RA1_dec - RA2_dec y = Sin(Deg2Rad(CDbl(A))) ' Evaluate the first term of the demoninator C = Cos(Deg2Rad(Dec2_dec)) * Tan(Deg2Rad(Dec1_dec)) ' Evaluate the second term of the demoninator D = Sin(Deg2Rad(Dec2_dec)) * Cos(Deg2Rad(CDbl(A))) x = C - D ' Find the tan of P ' Except error case at 90 and 270 degrees If (x <> 0) Then tanP = y / x ' Find the arctangent ' Compute the raw degrees E = Rad2Deg(Atn(tanP)) ' Find and adjust to the position angle coordinate system ' Evaluate normal cases ' Note: System adjusted to rotated coordinate system ' where North is 0 degrees and East is 270 degrees ' Adjust for quadrant If (y > 0) And (x > 0) Then f = Abs(E) ' 1st Quad If (y > 0) And (x < 0) Then f = Abs(E) + 90# ' 2nd Quad If (y < 0) And (x > 0) Then f = Abs(E) + 270# ' 4th Quad If (y < 0) And (x < 0) Then f = Abs(E) + 180# ' 3rd Quad ' Adjust for special cases If (y = 0) And (x > 0) Then f = 0# ' Directly above or origin If (y = 0) And (x < 0) Then f = 180# ' Directly below End If ' Evaluate special cases of 90 and 270 degrees If x = 0 Then ' Adjust for special cases If (y >= 0) And (x = 0) Then f = 90# ' Directly above or origin If (y < 0) And (x = 0) Then f = 270# ' Directly below End If dblPositionAngle_result = f ' Find the equivalent terresterial "To" bearing ' 360 less position angle is the "To bearing" g = 360 - dblPositionAngle_result ' Is always a positive number ' Round to less than 360 degs if greater than 360 Angle_dec = g ' Use iteration method If Angle_dec >= 0 Then Angle_dec_adj = Angle_dec Do While Angle_dec_adj >= 360# Angle_dec_adj = Angle_dec_adj - 360 Loop End If dblBearingTo_result = Angle_dec_adj ' Find the equivalent terresterial "From" bearing Angle_dec_adj = 0 g = 0 ' Subtract 180 degrees to the "To Bearing" g = dblBearingTo_result - 180 ' Round to less than 360 degs if greater than 360 or less than zero Angle_dec = g ' Use iteration method If Angle_dec >= 0 Then Angle_dec_adj = Angle_dec Do While Angle_dec_adj >= 360# Angle_dec_adj = Angle_dec_adj - 360 Loop End If If Angle_dec < 0 Then Angle_dec_adj = Angle_dec Do While Angle_dec_adj <= 360# Angle_dec_adj = Angle_dec_adj + 360# Loop Angle_dec_adj = Angle_dec_adj - 360# End If dblBearingFrom_result = Angle_dec_adj ' Return result PositionAngle = Trim(CStr(dblPositionAngle_result)) & " " & Trim(CStr(dblBearingTo_result)) & " " & Trim(CStr(dblBearingFrom_result)) End Function Public Function astrometryDistBet2Points3D(ByVal x1_tmp As Double, ByVal y1_tmp As Double, ByVal z1_tmp As Double, ByVal x2_tmp As Double, ByVal y2_tmp As Double, ByVal z2_tmp As Double) As Double ' Receives the orthogonal coords of two points in three dimensions and ' returns the distance between those two points ' =SQRT((x2-x1)^2 + (y2-y1)^2 + (z1-z2)^2 ) ' Test data ' ? DistBet2Points3D(0, 0, 0, 1, 1, 0) ' 1.4142135623731 ' ? DistBet2Points3D(0, 0, 0, 1, 1, 1) ' 1.73205080756888 ' Create working variables Dim x1, y1, z1, x2, y2, z2 As Double ' Get input x1 = x1_tmp y1 = y1_tmp z1 = z1_tmp x2 = x2_tmp y2 = y2_tmp z2 = z2_tmp ' Return result astrometryDistBet2Points3D = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2 + (z2 - z1) ^ 2) End Function Public Function astrometrySizeDriftRate(ByVal DecDeg_tmp As Double, ByVal DriftTime_secs_tmp As Double) As String ' Receives the declination of an object and a number of seconds of drift. ' Returns the size of the object in arcsecs and the celestial sphere drift rate at that declination ' Drift rate of a star along right ascension ' with your telescope's star drive turned off. ' v = drift rate of star across field of view on right ascension ' reference Line ' d = true field of view (or object size) in arcseconds ' D = 1 or infinity for celestial sphere ' v = ( D x cosine(radians(abs(declination)) ' * ( 360 deg x 60 min x 60 secs) / (24 hrs x 60 min x 60 secs) ' v = 15 x cosine(radians(abs(declination)) (Eq. 2.0) ' T = observed seconds to transit field of view ' d = v x T (distance = rate x time) (Eq. 3.0) ' Test data ' At the celestial equator the drift rate is approx. 15 arcsecs per second ' ? astrometrySizeDriftRate(1,1) ' 15 15 ' Working variables Dim DecDeg, DriftTime_secs, v, t, D, Size_arcsecs, Size_arcmins ' Get input DecDeg = DecDeg_tmp DriftTime_secs = DriftTime_secs_tmp v = 15 * Cos(Deg2Rad(Abs(DecDeg))) D = v * DriftTime_secs ' Return result astrometrySizeDriftRate = CStr(D) & " " & CStr(v) End Function Public Function RotFigureXY(ByVal R_tmp As Double, ByVal angNewdeg_tmp As Double, ByVal angOlddeg_tmp As Double) As String ' Receives the distance between the origin and a point, ' the old rotation angle of the point, usually, zero, in degrees ' and the new rotation angle, in degrees ' Returns the ortho x,y coordinates of the rotated point. ' Working variables Dim R, angNew, angOld, angNew_rads, angOld_rads, x, y As Double ' Get input R = R_tmp angNew = angNewdeg_tmp angOld = angOlddeg_tmp angNew_rads = Deg2Rad(angNew) angOld_rads = Deg2Rad(angOld) x = (Cos(angNew_rads) * Cos(angOld_rads)) - (Sin(angNew_rads) * Sin(angOld_rads)) y = (Sin(angNew_rads) * Cos(angOld_rads)) + (Cos(angNew_rads) * Sin(angOld_rads)) x = R * x y = R * y 'Return result RotFigureXY = CStr(Trim(x)) & " " & CStr(Trim(y)) End Function Public Function Adj360Rotation(ByVal dblDegdec_tmp As Double) As Double ' Utility function that returns a degree angle adjusted to ' limit of 360 where degrees of return is greater than 360. Dim dblDegdec, dblDegdec_adj As Double ' Get input dblDegdec = CDbl(dblDegdec_tmp) If CDbl(dblDegdec) = 360# Or CDbl(dblDegdec) = -360# Or CDbl(dblDegdec) = 0# Then dblDegdec_adj = 0# End If If CDbl(dblDegdec) > 0# And CDbl(dblDegdec) < 360# Then dblDegdec_adj = CDbl(dblDegdec) End If If CDbl(dblDegdec) < 0# And CDbl(dblDegdec) > -360# Then dblDegdec_adj = dblDegdec + 360# End If ' For > 360, parse out the fractional part of the neartest 360 degree circle If CDbl(dblDegdec) > 360# Then dblDegdec_adj = CDbl(dblDegdec) - CDbl(CInt(dblDegdec / 360#) * 360#) End If If CDbl(dblDegdec) < -360# Then dblDegdec_adj = CDbl(dblDegdec) - CDbl(CInt(dblDegdec / 360#) * 360#) dblDegdec_adj = dblDegdec_adj + 360 End If Adj360Rotation = dblDegdec_adj End Function Public Function Atan2Angle360(ByVal X_ortho_tmp As Double, ByVal Y_ortho_tmp As Double) As Double ' Receives an ortho coordinates and returns and angle in degrees between 0 and 360 ' Test data ' ? Atan2Angle360(1, 1) ' 45 ' ? Atan2Angle360(-1, 1) ' 135 ' ? Atan2Angle360(-1, -1)' 225 ' ? Atan2Angle360(1, -1) ' 275 ' Create working variables Dim X_ortho, Y_ortho, A As Double ' Get input X_ortho = X_ortho_tmp Y_ortho = Y_ortho_tmp A = Atan2(Y_ortho, X_ortho) If A < 0 Then A = A + (2 * PI()) End If Atan2Angle360 = Rad2Deg(A) End Function Function AngleBetw2Objs(ByVal RA1_dec_tmp As Double, ByVal Dec1_dec_tmp As Double, ByVal RA2_dec_tmp As Double, ByVal Dec2_dec_tmp As Double, ByVal s_outAngFormat_tmp As String) As Variant ' Returns angular distance between two objects in degrees ' Sources may be any angular coordinate system - equatorial or galactic ' ' Dependencies: Deg2arcmins, Deg2arcsecs ' ' Exceptions: ' Per Duffett, does not work within 10 arcmins of 0 deg or 180 deg ' Negative sizes are not error trapped. ' ' s_outAngleFormat formats the return value. Valid options are: ' "degdec" ' "arcmins" ' "arcsecs" ' ' Coordinates can be either equalatorial or galactic ' Duffett-Smith, P. 1988. Practical Astronomy with ' your calculator. 3rd ed. Cambridge Univ. ' Press. p. 51 ' ' Usage: ' Duffett-Smith example - ' Object 1 beta Ori ' BET ORI/209.24/-25.25/2000/ 78.63/ -8.20/4.22 ' Object 2 alpha CMa ' ALF CMA/227.23/ -8.89/2000/101.29/-16.72/379.21 ' ' AngleBetw2Objs(209.24,-25.25,227.23,-8.89,"degdec") yields 23.6770335744856 ' AngleBetw2Objs(209.24,-25.25,227.23,-8.89,"arcmins") yields 1421 ' AngleBetw2Objs(209.24,-25.25,227.23,-8.89,"arcsecs") yields 85237 ' ' Author: Kurt A. Fisher 2003 Dim RA_dec As Double Dim Dec_dec As Double Dim RA1_dec As Double Dim Dec1_dec As Double Dim RA2_dec As Double Dim Dec2_dec As Double Dim A, B, C, D, E, f, g, H As Double ' Temporary working variables Dim Angle_result_dec As Double Dim s_outAngFormat As String ' get equalatorial or galactic coordinates of two objects RA1_dec = RA1_dec_tmp Dec1_dec = Dec1_dec_tmp RA2_dec = RA2_dec_tmp Dec2_dec = Dec2_dec_tmp s_outAngFormat = s_outAngFormat_tmp ' evaluate equation components A = Sin(Deg2Rad((Dec1_dec))) B = Sin(Deg2Rad((Dec2_dec))) C = Cos(Deg2Rad((Dec1_dec))) D = Cos(Deg2Rad((Dec2_dec))) E = Cos(Deg2Rad(RA1_dec - RA2_dec)) ' evaluate equation yielding cos angular distance in degrees f = (A * B) + (C * D * E) ' convert to degrees and report result g = ArcCos((f)) H = Rad2Deg((g)) Select Case s_outAngFormat Case "degdec" AngleBetw2Objs = H Case "arcmins" AngleBetw2Objs = DecDeg2arcmins((H)) Case "arcsecs" AngleBetw2Objs = DecDeg2arcsecs((H)) Case Else ' gracefully handle by invalid option by returning decimal degrees AngleBetw2Objs = H End Select End Function Public Function AngleBetw2ObjsSmall(ByVal deltaA_arcsecs_tmp As Double, ByVal deltaD_arcsecs_tmp As Double, Dec1_decdeg_tmp As Double, Dec2_decdeg_tmp As Double) As Double ' Compute small angular separation distances ( < 1 degree or 10' ) where differences are expressed in arcminutes or arcseconds ' Note on input units: ' RA1_arcsecs - RA2_arcsecs must be expressed in decimal degrees, ' e.g. 15 * RA_arcsecs = A_arcsecs ' Sources: ' Duffett-Smith (3d 1998) at Sec. 32, p. 51 ' Meeus (1998) Eq. 17.2 at p. 109 function for small distances ' This formula uses are variation of the Pythagorean formula directly on arcsec ' angles ' Inputs: ' deltaA_arcsecs_tmp Difference in arcsecs between objects the right ascensions ' delta_D_arcsecs_tmp Difference in arcsecs between the objects in declination ' Dec1_decdeg_tmp Declination of the first object in decimal degrees ' Dec2_decdeg_tmp Declination of the second object in decimal degrees ' Define variables Dim deltaA_arcsecs, deltaD_arcsecs, DecAverage_rads, Dec1_decdeg, Dec2_decdeg As Double Dim deltaA_rads, deltaD_rads, DecDiff_rads Dim A, B, C, D, E As Double ' Working variables ' Get input deltaA_arcsecs = deltaA_arcsecs_tmp deltaD_arcsecs = deltaD_arcsecs_tmp Dec1_decdeg = Dec1_decdeg_tmp Dec2_decdeg = Dec2_decdeg_tmp ' Find the average declination in arcsecs A = CDbl((Dec1_decdeg + Dec2_decdeg) / 2#) ' Convert to radians DecAverage_rads = Deg2Rad(CDbl(A)) ' Find the cosine of the average declination difference B = Cos(CDbl(DecAverage_rads)) ' Find the true degree-arcsec size of the deltaA's at this declination C = B * deltaA_arcsecs ' Apply the Pythagorean formula to the arcsecs dimensions directly D = Sqr(C ^ 2 + deltaD_arcsecs ^ 2) AngleBetw2ObjsSmall = CDbl(D) End Function Function PI() As Double ' Generic dervied math functions from numerous published sources ' and Microsoft 1997 Neatcode library PI = Atn(1) * 4 End Function Function Deg2Rad(ByVal x As Double) As Double ' Degrees to radians Deg2Rad = x / 180 * PI() End Function Function Rad2Deg(ByVal x As Double) As Double ' Radians to Degrees Rad2Deg = x / PI() * 180 End Function Function ArcSin(ByVal x As Double) As Double ' Inverse Sine If x = 1 Then ArcSin = PI() / 2 ElseIf x = -1 Then ArcSin = -PI() / 2 Else ArcSin = Atn(x / Sqr(-x * x + 1)) End If End Function Function ArcCos(ByVal x As Double) As Double ' Inverse Cosine If x = 1 Then ArcCos = 0 ElseIf x = -1 Then ArcCos = -PI() Else ' Corrects formula error in MS Neatcode library ArcCos = Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1) End If End Function Public Function Atan2(ByVal y As Double, ByVal x As Double) As Double ' Four quadrant arctangent ' Source: http://www.vb-helper.com/howto_atan2.html Dim theta As Double If (Abs(x) < 0.0000001) Then If (Abs(y) < 0.0000001) Then theta = 0# ElseIf (y > 0#) Then theta = 1.5707963267949 Else theta = -1.5707963267949 End If Else theta = Atn(y / x) If (x < 0) Then If (y >= 0#) Then theta = 3.14159265358979 + theta Else theta = theta - 3.14159265358979 End If End If End If Atan2 = theta End Function Static Function Log10(ByVal x) ' From Visual Basic Help If x >= 0 Then Log10 = Log(x) / Log(10#) End If If x < 0 Then Log10 = -1# * (Log(Abs(x)) / Log(10#)) End If End Function Public Function NumTrunc(ByVal dblN, ByVal intPrecision) As Double ' Returns a truncated precision number ' Limited to integer size Dim A, B, C As Double ' Temporary working variables C = 10 ^ intPrecision A = dblN * C B = Int(A) NumTrunc = B / C End Function Public Function FracPartNumber(ByVal N_tmp As Double) As String ' Receives number; returns the parts to the left and right of the decimal point Dim strLeft, strRight, strSign As String Dim n, R As Double Dim intSign As Integer ' Get the input n = N_tmp ' Get the sign of the number: Negative -1; Postive or Zero 0 or 1 intSign = Sgn(n) ' Make the number positive n = Abs(n) ' Add small amount to right side to assure a decimal point exists n = n + 0.000000000001 ' Build the sign string If intSign >= 0 Then strSign = "" Else strSign = "-" End If ' Get the left side strLeft = strSign & Left(CStr(n), InStr(1, CStr(n), ".", 1)) & "0" ' Get the right side R = CDbl("0." & Right(CStr(n), Len(CStr(n)) - InStr(1, CStr(n), ".", 1))) - 0.000000000001 strRight = CStr(R) ' Return the result FracPartNumber = strLeft & " " & strRight End Function Public Function UncertSum(ByVal dblA As Double, ByVal dbluncA As Double, ByVal dblB As Double, ByVal dbluncB As Double, ByVal strReturnType As String, ByVal intPrecision As Integer) As Double ' This function computes scientific uncertainty with an upper and lower bound. ' This utility provides the sum operation. ' Source: Gullberg, J. 1997. Mathematics: From the Birth of Numbers. W.W. Norton. pp. 166-167. ' Inputs: ' dblA ' value A ' dbluncA ' +- uncertainty of value A ' dblB ' value B ' dbluncB ' +- uncertainty of value B ' Outputs: ' Returns a collection containing valuea and +- uncertainty as positive double ' Usage ' ? UncertSum (4.3760,0.0005,2.3560,0.0005,"value",4) yields 6.732 ' ? UncertSum (4.3760,0.0005,2.3560,0.0005," ",4) yields 0.001 Dim C, D, E, f As Double 'working variables ' Gracefully handle negative input runtime error dbluncA = Abs(dbluncA) dbluncB = Abs(dbluncB) 'Value - sum or difference C = dblA + dblB 'Lower bound ' D = (dblA - dbluncA) + (dblB - dbluncB) ' D = (dblA + dblB) - (dbluncA + dbluncB) ' Therefore ' D = C - (dbluncA + dbluncB) ' Therefore, lower bound is ' E = -(dbluncA + dbluncB) 'Upper bound ' D = (dblA + dbluncA) + (dblB + dbluncB) ' D = (dblA + dblB) + (dbluncA + dbluncB) ' Therefore ' D = C + (dbluncA + dbluncB) ' Therefore, upper bound is E = (dbluncA + dbluncB) If strReturnType = "value" Then C = NumTrunc(C, intPrecision) UncertSum = C Else E = NumTrunc(E, intPrecision) UncertSum = E End If End Function Public Function UncertDif(ByVal dblA As Double, ByVal dbluncA As Double, ByVal dblB As Double, ByVal dbluncB As Double, ByVal strReturnType As String, ByVal intPrecision As Integer) As Double ' This function computes scientific uncertainty with an upper and lower bound. ' This utility provides the difference operation. ' Source: Gullberg, J. 1997. Mathematics: From the Birth of Numbers. W.W. Norton. pp. 166-167. ' Inputs: ' dblA ' value A ' dbluncA ' +- uncertainty of value A ' dblB ' value B ' dbluncB ' +- uncertainty of value B ' Outputs: ' Returns a collection containing valuea and +- uncertainty as positive double ' Usage ' ? UncertDif (4.3760,0.0005,2.3560,0.0005,"value",4) yields 2.02 ' ? UncertDif (4.3760,0.0005,2.3560,0.0005," ",4) yields 0.001 Dim C, D, E, f As Double 'working variables ' Gracefully handle negative input runtime error dbluncA = Abs(dbluncA) dbluncB = Abs(dbluncB) 'Value - sum or difference C = dblA - dblB 'Lower bound ' D = (dblA - dbluncA) - (dblB + dbluncB) ' D = (dblA - dblB) - (dbluncA + dbluncB) ' Therefore ' D = C - (dbluncA + dbluncB) ' Therefore, lower bound is ' E = -(dbluncA + dbluncB) 'Upper bound ' D = (dblA + dbluncA) - (dblB - dbluncB) ' D = (dblA - dblB) + (dbluncA + dbluncB) ' Therefore ' D = C + (dbluncA + dbluncB) ' Therefore, upper bound is E = (dbluncA + dbluncB) If strReturnType = "value" Then C = NumTrunc(C, intPrecision) UncertDif = C Else E = NumTrunc(E, intPrecision) UncertDif = E End If End Function Public Function UncertProd(ByVal dblA As Double, ByVal dbluncA As Double, ByVal dblB As Double, ByVal dbluncB As Double, ByVal strReturnType As String, ByVal intPrecision As Integer) As Double ' This function computes scientific uncertainty with an upper and lower bound. ' This utility provides the product operations. ' Source: Gullberg, J. 1997. Mathematics: From the Birth of Numbers. W.W. Norton. pp. 166-167. ' Inputs: ' dblA ' value A ' dbluncA ' +- uncertainty of value A ' dblB ' value B ' dbluncB ' +- uncertainty of value B ' Outputs: ' Returns a collection containing value and +- uncertainty as positive double ' ? UncertProd(4.3760,0.0005,2.3560,0.0005,"value",4) yields 10.0398 ' ? UncertProd(4.3760,0.0005,2.3560,0.0005," ",4) yields 0.0033 Dim C, D, E, f As Double 'working variables Dim colResult As Collection 'working variables ' Gracefully handle negative input runtime error dbluncA = Abs(dbluncA) dbluncB = Abs(dbluncB) 'Value - product C = dblA * dblB 'Lower bound 'D = C - ((dblA - dbluncA) * (dblB - dbluncB)) 'Upper bound E = ((dblA + dbluncA) * (dblB + dbluncB)) - C If strReturnType = "value" Then C = NumTrunc(C, intPrecision) UncertProd = C Else E = NumTrunc(E, intPrecision) UncertProd = E End If End Function Public Function UncertQuot(ByVal dblA As Double, ByVal dbluncA As Double, ByVal dblB As Double, ByVal dbluncB As Double, ByVal strReturnType As String, ByVal intPrecision As Integer) As Double ' This function computes scientific uncertainty with an upper and lower bound. ' This utility provides the quotient operations. ' Source: Gullberg, J. 1997. Mathematics: From the Birth of Numbers. W.W. Norton. pp. 166-167. ' Inputs: ' dblA ' value A ' dbluncA ' +- uncertainty of value A ' dblB ' value B ' dbluncB ' +- uncertainty of value B ' Constraints ' dbluncB may not be zero ' abs(dblB) - abs(dbluncB) may not be zero ' Outputs: ' Returns a collection containing value and +- uncertainty as positive double ' Useage ' ? UncertQuot(4.3760,0.0005,2.3560,0.0005,"value",4) yields 1.8573 ' ? UncertQuot(4.3760,0.0005,2.3560,0.0005," ",4) yields 0.0006 Dim C, D, E, f As Double 'working variables Dim colResult As Collection 'working variables ' Gracefully handle negative input runtime error dbluncA = Abs(dbluncA) dbluncB = Abs(dbluncB) ' Error trapping for zero values ' Gracefully handle error conditions If dblB = 0 Then dblB = 0.0000000001 f = dblB - dbluncB If f = 0 Then f = 0.0000000001 'Value - product C = dblA / dblB 'Lower bound 'D = C - ((dblA - dbluncA) / (dblB + dbluncB)) 'Upper bound E = ((dblA + dbluncA) / f) - C If strReturnType = "value" Then C = NumTrunc(C, intPrecision) UncertQuot = C Else E = NumTrunc(E, intPrecision) UncertQuot = E End If End Function Function GetCSWord(ByVal S, ByVal Indx As Integer) 'Returns the nth word in a specific field separated by spaces. ' Source: Dev Ashish's Access Web 2003 ' http://www.mvps.org/access/strings/str0003.htm Dim WC As Integer, Count As Integer Dim SPos As Integer, EPos As Integer WC = CountCSWords(S) If Indx < 1 Or Indx > WC Then GetCSWord = Null Exit Function End If Count = 1 SPos = 1 For Count = 2 To Indx SPos = InStr(SPos, S, " ") + 1 Next Count EPos = InStr(SPos, S, " ") - 1 If EPos <= 0 Then EPos = Len(S) GetCSWord = Trim(Mid(S, SPos, EPos - SPos + 1)) End Function Function CountCSWords(ByVal S) As Integer 'Counts the words in a string that are separated by spaces. ' Source: Dev Ashish's Access Web 2003 ' http://www.mvps.org/access/strings/str0003.htm Dim WC As Integer, Pos As Integer If VarType(S) <> 8 Or Len(S) = 0 Then CountCSWords = 0 Exit Function End If WC = 1 Pos = InStr(S, " ") Do While Pos > 0 WC = WC + 1 Pos = InStr(Pos + 1, S, " ") Loop CountCSWords = WC End Function Public Function PadFixedLengthString(ByVal S_tmp As String, ByVal intLength_tmp As Integer) As String ' Returns a space-padded, left justified string for use in making ' astronomical catlogues. ' intLength is maximum length of padded string ' Intialize working variables Dim S, t As String ' Temporary working variables Dim intLength As Integer ' Get input S = S_tmp intLength = intLength_tmp ' Build string t = Space(50) ' Trim source string first - left and right t = Trim(S) & t PadFixedLengthString = Mid(t, 1, intLength_tmp) End Function Public Function MakeCielExportCat(ByVal PlotString_tmp As String, ByVal RA_s_tmp As String, ByVal Dec_s_tmp As String, ByVal Specifier_type_tmp As String, ByVal PlotSize_tmp As String, ByVal Descrip1_tmp As String, ByVal Descrip2_tmp As String, ByVal Descrip3_tmp As String) As String ' Creates a string that can be saved as a Ciel nebulae size text catalogue ' file for use with Cartes de Ciel ' ' Author: K. Fisher fisherka@csolutions.net 2/2004 ' Ciel text catalogue input fields are: ' Field StartPos Length VBA utils var Example ' Id 1 8 PlotString_tmp m_v - as needed ' 9 1 n/a ' 10 1 Catalogue formatter "J" ' RAh 11 2 RA and Dec decimal transformed ' RAm 13 2 to form "JHHMMSS.SS+DDMMSS.S " ' RAs 15 5 ' DE+ 20 1 ' DEd 21 2 ' DEm 23 2 ' DEs 25 4 ' 29 1 n/a ' Specifier 30 5 Specifier_type "*i*" - Simbad type specifier ' 35 1 n/a ' Size 36 5 PlotSize_tmp "5____" ' 41 1 n/a ' Descrip1 42 15 Descrip1 Id or Name ' 57 1 n/a ' Descrip2 58 15 Descript2 Alternate name or other info ' 73 1 n/a ' Descrip3 74 10 Descript3 Spectral type ' Create working variables Dim S, t, U As String ' Temporary working variables Dim PlotString As String Dim RA_s As String Dim Dec_s As String Dim Specifier_type As String Dim PlotSize As String Dim Descrip1 As String Dim Descrip2 As String Dim Descrip3 As String ' Get input PlotString = PlotString_tmp RA_s = RA_s_tmp Dec_s = Dec_s_tmp Specifier_type = Specifier_type_tmp PlotSize = PlotSize_tmp Descrip1 = Descrip1_tmp Descrip2 = Descrip2_tmp Descrip3 = Descrip3_tmp ' Build padded string S = PadFixedLengthString(PlotString, 8) ' Make the aui format coordinates substring t = Convert_ReformatTxtCoor((RA_s), "HH MM SS.SS", "acr cHHMMSS.SS", "J", "") U = Convert_ReformatTxtCoor((Dec_s), "sDD MM SS.S", "sDDMMSS.S", "", "") ' Complete building string S = (S) & " " & (t) & (U) S = (S) & " " & PadFixedLengthString(Specifier_type, 5) S = (S) & " " & PadFixedLengthString(PlotSize, 5) S = (S) & " " & PadFixedLengthString(Descrip1, 15) S = (S) & " " & PadFixedLengthString(Descrip2, 15) S = (S) & " " & PadFixedLengthString(Descrip3, 10) ' Return result MakeCielExportCat = (S) End Function Public Function XMLstrEscapeSpecialChars(ByVal strIn As String) As String ' Replaces special characters in a string with escaped XSL-XML compatible characters ' Source: Andrew H. Watt. 2003. Teach Yourself XML in 10 Minutes. Sams Pub. p. 27 ' Note: MS-IE XML-XSL processor errors on unique character combinations can be particularly vexsome. ' Typically, a xml-xsl file has to be processed, and any individual line errors have to examined individually ' and the character string that presents a problem for the processor has to be adjusted manually. ' Working variables Dim S, charSpec As String ' Get input S = strIn ' Replace special characters < < > > ' ' " " & & ' Also replace other character references - colon : semi-colon ; ' backslash \ forwardslash / degree sign ° ' Programming note - can't replace semicolon ' S = Replace(S, ";", ";", 1, -1, 1) ' Programming note - replace ; and & characters first to avoid runtime errors S = Replace(S, "&", "&", 1, -1, 1) S = Replace(S, "<", "<", 1, -1, 1) S = Replace(S, ">", ">", 1, -1, 1) S = Replace(S, "'", "'", 1, -1, 1) ' Double quote S = Replace(S, Chr(34), """, 1, -1, 1) S = Replace(S, ":", ":", 1, -1, 1) ' Backslash \ S = Replace(S, "\", "\", 1, -1, 1) ' Forwardslash / S = Replace(S, "/", "/", 1, -1, 1) ' Programming note - MS-IE xml will not pass the standard xml-html character entity ' Degree sign S = Replace(S, "°", "°", 1, -1, 1) ' Programming note - MS-IE xml aborts on this unique character string - special spectral type ' "C..." - Replacement is not always successful and may have to be done manually S = Replace(S, "...", "C . . .", 1, -1, 1) ' Return result XMLstrEscapeSpecialChars = S End Function