I have an older program that I wrote some time ago in VB then late updated for VB.Net. The program calculates the ephemerides of comets based on orbital elements downloaded from the Minor Planet Center website. I am now trying to make an iPhone/iPad application using the same idea and the same calculations. Unfortunately I am too new at objective-c that trying to convert the extensive calculations is proving quite difficult for me. I am hoping that someone can help me out and at least give me some pointers.
In particular, I have no idea how to accomplish this in objective-c:
Dim d1 As Date = DateSerial(yr_peri%, mo_peri%, dy_peri%)
dateser# = d1.ToOADate()
dateser# = dateser# + frac_day_peri#
jdperi# = 2415018.5 + dateser#
and this:
' Find Julian Day Number of start date & time
Dim d2 As Date = DateSerial(yr%, mo%, dd%)
Dim inputdate# = d2.ToOADate()
jdp# = inputdate + 2415018.5 ' use this day for phenom (rise/set) calcs
jd# = jdp# + fd# ' add fractional days
Here is the full VB version of the calculations in case it helps to understand the context:
Public Sub calc(ByVal cname$, ByVal qq#, ByVal e#, ByVal w#, ByVal om#, ByVal I#, ByVal yr_peri%, ByVal mo_peri%, ByVal day_peri#, ByVal abs_mag#, ByVal mag_n#)
Dim ep#
Dim dy_peri%, frac_day_peri#, dateser#, jdperi#
Dim yr%, mo%, dd%, fd#
Dim jdp#, jd#, interval#
Dim twi%, p_date#, ut#, er%
Dim t#, tt#, f#, p#, g#, q#, h#, r#
Dim a1#, a2#, b1#, b2#, c1#, c2#
Dim sdt#, a#, n#, period#, m#, e0#
Dim mrs#, frs%, ers0#, drs#, jrs%, mrs1#, e2#
Dim tv#, v#, r1#, r2#, w0#, s#, s0#
Dim x#, y#, z#, ls#, ms#, es#, th#, vs#, c#, th2000#
Dim rs#, xs#, ys#, zs#, xx#, yy#, radeg#, rad#, ra#, ras$
Dim dl#, dlm#, sd#, dec#, ds$
Dim cospsi#, psi#, elong#, cosbeta#, beta#, phase#
Dim ha#, ang#, az#, sinalt#, alt#, m1#, vl#
Dim Lati$, Longi$
Dim outpt$
'----------------------------------------------
'** AA=ASTRONOMICAL ALGORITHMS BY JEAN MEEUS **
'----------------------------------------------
counter% = 0
ep# = 23.439291 ' year 2000, AA page 214
cname$ = Trim(cname$)
'Get orbital elements
MainForm.grpEphemerides.Text = "Comet: " & cname$ & " - " + when$ ' Name
'qq# = Val(comet$(ind%, 1)) ' perihelion dist
'e# = Val(comet$(ind%, 2)) ' eccentricity (1=parabola)
'w# = Val(comet$(ind%, 3)) ' arg of perihelion - little omega
'om# = Val(comet$(ind%, 4)) ' long of ascending node - big omega
'I# = Val(comet$(ind%, 5)) ' inclination
'yr_peri% = Val(comet$(ind%, 6)) ' year
'mo_peri% = Val(comet$(ind%, 7)) ' month
'day_peri# = Val(comet$(ind%, 8)) ' day
'abs_mag# = Val(comet$(ind%, 9)) ' magnitude coef #1
'mag_n# = Val(comet$(ind%, 10)) ' magnitude coef #2
' Find JD of perihelion
dy_peri% = Int(day_peri#) ' day
frac_day_peri# = day_peri# - dy_peri% ' fractional day
Dim d1 As Date = DateSerial(yr_peri%, mo_peri%, dy_peri%)
dateser# = d1.ToOADate()
dateser# = dateser# + frac_day_peri#
jdperi# = 2415018.5 + dateser#
' Start date/time
yr% = Val(MainForm.txtYear.Text)
mo% = Val(MainForm.txtMonth.Text)
dd% = Val(MainForm.txtDay.Text)
fd# = Val(MainForm.txtHour.Text) / 24 + Val(MainForm.txtMin.Text) / 1440 + Val(MainForm.txtSec.Text) / 86400
lat# = Val(MainForm.txtLat.Text)
lon# = Val(MainForm.txtLon.Text)
' Find Julian Day Number of start date & time
Dim d2 As Date = DateSerial(yr%, mo%, dd%)
Dim inputdate# = d2.ToOADate()
jdp# = inputdate + 2415018.5 ' use this day for phenom (rise/set) calcs
jd# = jdp# + fd# ' add fractional days
interval# = Val(MainForm.txtInterval.Text)
' Determine if fixed times or rise/set phenomena
If MainForm.rdbFixedTimes.Checked = True Then
twi% = 0 ' use fixed times
Else
twi% = 1 ' use twilight times
interval# = Int(interval#) ' Set integer days (fraction doesn't make sense)
If interval# < 1 Then interval# = 1 ' make sure interval is at least one day
End If
'-----------
' Loop here
'-----------
While counter% <= Val(MainForm.txtNumSteps.Text)
If twi% = 1 Then
'Serial date of this calculation
p_date# = jdp# - 2415018.5
Call phenom(p_date#, ut#, er%) ' calcs for rise/set, returns er% (among other stuff)
jd# = Int(p_date#) + 2415018.5 + ut# / 24 + 1 / 2880
If er% = 1 Then ' er%=1 means no phenom (ex: sunrise) that day
ser_date# = jd# - 2415018.5
outpt$ = Format$(ser_date#, "mm/dd/yy") + " " + "No Rise/Set Phenomena on this Date"
GoTo No_phenom
End If
End If
' Days from perihelion of comet
t# = jd# - jdperi#
' Time (in Julian Centuries (36525 days)) since 1900 Jan 0.5 ET
' AA page 151
tt# = (jd# - 2451545.0#) / 36525.0!
'-------------- CALC F,G,H,P,Q,R (PP214 IN AA) ---------------
f# = Cos(dr# * om#)
p# = -Sin(dr# * om#) * Cos(dr# * I#)
g# = Sin(dr# * om#) * Cos(dr# * ep#)
q# = Cos(dr# * om#) * Cos(dr# * I#) * Cos(dr# * ep#) - Sin(dr# * I#) * Sin(dr# * ep#)
h# = Sin(dr# * om#) * Sin(dr# * ep#)
r# = Cos(dr# * om#) * Cos(dr# * I#) * Sin(dr# * ep#) + Sin(dr# * I#) * Cos(dr# * ep#)
'---------------- CALC A1,A2,B1,B2,C1,C2 (PP214 AA) -----------------
Call rec2pol(p#, f#, a1#, a2#)
Call rec2pol(q#, g#, b1#, b2#)
Call rec2pol(r#, h#, c1#, c2#)
'--------------- SIDEREAL TIME AT GREENWICH (PP84 AA) ----------------
' In degrees
sdt# = 280.46061837 + 360.98564736629 * (jd# - 2451545) + 0.000387933 * tt# * tt# - tt# * tt# * tt# / 38710000
sdt# = sdt# - 360 * Int(sdt# / 360)
'------------------- CALC ECCENTRIC & TRUE ANOMALY --------------------
If e# < 1 Then ' elliptical motion (pp214 AA)
a# = qq# / (1 - e#) ' semi major axis (AU)
n# = 0.9856076686 / a# / Sqrt(a#) ' mean motion (deg/day)
period# = 360 / n# / 365.25
MainForm.lblOrbitInf.Text = "Semimajor axis: " + Format$(a#, "#.000") + " AU Period: " + Format$(period#, "#.00") + " Years"
m# = n# * t# ' mean anomaly
e0# = rd# * e# ' "modified" eccentricity, page 187)
'ecc anomaly (page 195)
mrs# = m# * dr#
frs% = Sign(mrs#)
mrs# = Abs(mrs#) / (2 * pi#)
mrs# = (mrs# - Int(mrs#)) * 2 * pi# * frs%
If mrs# < 0 Then mrs# = mrs# + 2 * pi#
frs% = 1
If mrs# > pi# Then frs% = -1
If mrs# > pi# Then mrs# = 2 * pi# - mrs#
ers0# = pi# / 2
drs# = pi# / 4
For jrs% = 1 To 53
mrs1# = ers0# - e# * Sin(ers0#)
ers0# = ers0# + drs# * Sign(mrs# - mrs1#)
drs# = drs# / 2
Next jrs%
ers0# = ers0# * frs%
e2# = ers0# * rd#
'alpha# = (1 - e#) / (4 * e# + .5)
'beta# = m# / (8 * e# + 1)
'signbeta% = Sgn(beta#)
'cube# = beta# + signbeta% * Sqr(beta# ^ 2 + alpha# ^ 3)
'signz% = Sgn(cube#)
'zz# = signz% * (Abs(cube#)) ^ (1 / 3)
'ss0# = zz# - alpha# / 2
'ss# = ss0# - .078 * ss0# ^ 5 / (1 + e#)
'e1# = m# + e# * (3 * ss# - 4 * ss# ^ 3)
'loop1% = 0 ' counter
1040:
'e2# = e1# + (m# + e0# * Sin(dr# * e1#) - e1#) / (1 - e# * Cos(dr# * e1#))
'If Abs(e2# - e1#) < .0000000001 Or loop1% > 1000 Then ' e2#:ecc anom found
tv# = Sqrt((1 + e#) / (1 - e#)) * Tan(dr# * e2# / 2)
v# = 2 * rd# * Atan(tv#) ' true anomaly
r1# = a# * (1 - e# * Cos(dr# * e2#)) ' dist from sun AU
GoTo 2120
'Else
' e1# = e2# ' new ecc anomaly
' loop1% = loop1% + 1
' GoTo 1040
'End If
Else ' parabolic motion (pp225 AA)
MainForm.lblOrbitInf.Text = "This comet has a parabolic orbit-semimajor axis and period not determined"
w0# = 0.0364911624 * t# / qq# / Sqrt(qq#)
s# = 0
2110:
s0# = (2 * s# * s# * s# + w0#) / 3 / (s# * s# + 1)
If Abs(s# - s0#) > 0.0000000001 Then
s# = s0#
GoTo 2110
End If
s# = s0#
v# = 2 * rd# * Atan(s#) ' true anomaly
r1# = qq# * (1 + s# * s#) ' distance from sun (AU)
End If
2120:
v# = v# - 360 * Int(v# / 360) ' true anomaly
'r1#= dist from sun in au (above)
r2# = r1# * 92955807.0# ' Dist from sun (miles)
'-------- CALC HELIOCENTRIC RECT COORD OF COMET (PP215 AA) ----------
x# = r1# * a2# * Sin(dr# * (a1# + w# + v#))
y# = r1# * b2# * Sin(dr# * (b1# + w# + v#))
z# = r1# * c2# * Sin(dr# * (c1# + w# + v#))
'-------- CALC GEOCENTRIC RECT COORD OF SUN (PP151,152 AA) ---------
ls# = 280.46645 + 36000.76983 * tt# + 0.0003032 * tt# * tt#
ls# = ls# - 360 * Int(ls# / 360)
ms# = 357.5291 + 35999.0503 * tt# - 0.0001559 * tt# * tt# - 0.00000048 * tt# * tt# * tt#
ms# = ms# - 360 * Int(ms# / 360)
es# = 0.016708617 - 0.000042037 * tt# - 0.0000001236 * tt# * tt#
c# = (1.9146 - 0.004817 * tt# - 0.000014 * tt# * tt#) * Sin(dr# * ms#) + (0.019993 - 0.000101 * tt#) * Sin(2 * dr# * ms#) + 0.00029 * Sin(3 * dr# * ms#)
th# = ls# + c# ' TH is Sun's true longitude
vs# = ms# + c# ' VS is Sun's true anomaly
' pp152 AA
th2000# = th# - 0.01397 * (yr% + mo% / 12 + dd% / 365 - 2000) ' Sun true long referred to 2000
' sun's radius vector
rs# = (1.000001018 * (1 - es# * es#)) / (1 + es# * Cos(dr# * vs#))
' rectangular coord of sun (pp159 AA)
xs# = rs# * Cos(dr# * th2000#) ' Xsun
ys# = rs# * Sin(dr# * th2000#) * Cos(dr# * ep#) ' Ysun
zs# = rs# * Sin(dr# * th2000#) * Sin(dr# * ep#) ' Zsun
'------------- CALC GEOCEN开发者_JAVA技巧TRIC RA & DEC (PP 119 IN AFFC) --------------
xx# = xs# + x#
yy# = ys# + y#
Call rec2pol(xx#, yy#, radeg#, rad#)
ra# = radeg# / 15 ' RA of comet (2000) in hours
' 0 means time
ras$ = dec2ddmm(ra#, 0) ' ra in string form
'dist from comet to earth (au)
dl# = Sqrt((xs# + x#) ^ 2 + (ys# + y#) ^ 2 + (zs# + z#) ^ 2)
dlm# = dl# * 92955807.0# ' distance from earth in miles
sd# = (zs# + z#) / dl#
'DEC of comet (2000)
dec# = rd# * Atan(sd# / Sqrt(-sd# * sd# + 1))
' 1 means angle
ds$ = dec2ddmm(dec#, 1) ' dec in string form
' Calc elongation
cospsi# = (rs# * rs# + dl# * dl# - r1# * r1#) / (2 * rs# * dl#)
psi# = -Atan(cospsi# / Sqrt(-cospsi# * cospsi# + 1)) + pi# / 2
elong# = psi# * rd#
elong# = elong# - 360 * Int(elong# / 360)
' Calc phase
cosbeta# = (r1# * r1# + dl# * dl# - rs# * rs#) / (2 * r1# * dl#)
beta# = -Atan(cosbeta# / Sqrt(-cosbeta# * cosbeta# + 1)) + pi# / 2
phase# = beta# * rd#
phase# = phase# - 360 * Int(phase# / 360)
'------------ CALCULATE ALTITUDE & AZIMUTH (AA PP89) --------------
ha# = sdt# - lon# - radeg#
yy# = Sin(dr# * ha#)
xx# = (Cos(dr# * ha#) * Sin(dr# * lat#) - Tan(dr# * dec#) * Cos(dr# * lat#))
' rectangular to polar
Call rec2pol(xx#, yy#, ang#, rad#)
az# = ang# + 180
az# = az# - 360 * Int(az# / 360)
sinalt# = Sin(dr# * lat#) * Sin(dr# * dec#) + Cos(dr# * lat#) * Cos(dr# * dec#) * Cos(dr# * ha#)
alt# = rd# * Atan(sinalt# / Sqrt(1 - sinalt# * sinalt#))
'-------------- CALC EST MAGNITUDES & ORBITAL VELOCITY (AA PP 216) ---------------
m1# = abs_mag# + 5 * Log(dl#) / Log(10) + 2.5 * mag_n# * Log(r1#) / Log(10)
If e# < 1 Then
vl# = 42.1219 * Sqrt((1 / r1#) - (1 / (2 * a#))) ' elliptical orbital velocity, page 223
Else
vl# = 42.1219 * Sqrt(1 / r1#) ' parabolic orital velocity (my approximation)
End If
'------------- Stuff for output form -------------
Dim dt As DateTime = DateTime.FromOADate(dateser#)
MainForm.lblJDperi.Text = "Date of perihelion: " + dt.ToString("dd MMM yyyy hh:mm:ss") + " UTC JD: " + Format$(jdperi#, "#.00000")
Lati$ = " "
Longi$ = " "
If lat# > 0 Then Lati$ = "N " Else Lati$ = "S "
If lon# > 0 Then Longi$ = "W" Else Longi$ = "E"
MainForm.lblLoc.Text = "Observing Location: " + Format$(lat#, "+00.0000;-00.0000") + Chr(176) + Lati$ + Format$(lon#, "+000.0000;-000.0000") + Chr(176) + Longi$
'Serial date of this calculation
ser_date# = jd# - 2415018.5
' put the main info into a string
dt = DateTime.FromOADate(ser_date#)
outpt = dt.ToString("MM/dd/yy HH:mm:ss") & " " & ras$ & " " & ds$ & " " & Format$(alt#, "+00.0;-00.0") & Chr(176) & " " & Format$(az#, "000.0") & Chr(176)
outpt = outpt & " " & Format$(r1#, "00.0000") & " " & Format$(dl#, "00.0000") & " " & Format$(m1#, "+00.0;-00.0") + " " & Format$(elong#, "000.0") & Chr(176)
No_phenom:
' put the string into the list box
MainForm.lstPosition.Items.Add(outpt)
' put the "extra" data into the array
extra#(counter%, 0) = t# ' time to perih
extra#(counter%, 1) = vl# ' speed
extra#(counter%, 2) = phase# ' phase
extra#(counter%, 3) = v# ' true anomaly
' Increment counter
counter% = counter% + 1
' New JD
jd# = jd# + interval#
jdp# = jdp# + interval#
End While ' do again until done
End Sub
I appreciate any help or direction anyone can give on converting this function to objective-c.
Regards, Keith
When working with dates in Objective-C you will often need to work with more than 1 class unlike other languages such as C#, Java and VB.NET. For creating a date from a string or getting a string from a date you will want to use the NSDateFormatter. Now when attempting to modify or create a date with "components" (year,month,day etc.) you will want to use NSDateComponents along with the correct NSCalendar. The component solution looks like it will assist you with the DateSerial part. The line jdp# = inputdate + 2415018.5 should be convertible (if it is seconds) to
jdp = [NSDate dateWithTimeInterval:2415018.5 sinceDate:inputdate];
![Interactive visualization of a graph in python [closed]](https://www.devze.com/res/2023/04-10/09/92d32fe8c0d22fb96bd6f6e8b7d1f457.gif)
加载中,请稍侯......
精彩评论