开发者

How can I convert these VB calculations to Objective-c?

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];
0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜