def itrace(Y, YOLD, YAPX, YP, BX, BY, BZ, BB, SGN, DS, NSTP):
    """
    Follow a geomagnetic field line until passing its apex
    Cartesian component magnetic field (partial) derivitives steer the trace

    This uses the 4-point Adams formula after initialization.
    First 7 iterations advance point by 3 steps.

    INPUTS:
     YAPX  = Matrix of cartesian coordinates of the three points
             about the apex.  Set in itrace.

    FLDCOMD has geomagnetic field at current trace point:
     BX    = X component (nT)
     BY    = Y component (nT)
     BZ    = Z component (nT)
     BB    = Magnitude   (nT)

    itrace has field line tracing variables determined in findApex:
     NSTP  = Step count.
     Y     = Array containing current tracing point cartesian coordinates.
     YOLD  = Array containing previous tracing point cartesian coordinates.
     SGN   = Determines direction of trace.
     DS    = Step size (arc length in km).

    REFERENCES:
     Stassinopoulos E. G. , Mead Gilbert D., X-841-72-17 (1971) GSFC,
     Greenbelt, Maryland
    """
    YP[0][3] = SGN * BX / BB
    YP[1][3] = SGN * BY / BB
    YP[2][3] = SGN * BZ / BB
    D2 = DS / 2.0
    D6 = DS / 6.0
    D12 = DS / 12.0
    D24 = DS / 24.0
    if NSTP <= 7:
        for I in range(0, 3):
            if NSTP == 1:
                YP[I][0] = YP[I][3]
                YOLD[I] = Y[I]
                YAPX[I][0] = Y[I]
                Y[I] = YOLD[I] + DS * YP[I][0]
            elif NSTP == 2:
                YP[I][1] = YP[I][3]
                Y[I] = YOLD[I] + D2 * (YP[I][1] + YP[I][0])
            elif NSTP == 3:
                Y[I] = YOLD[I] + D6 * (2.0 * YP[I][3] + YP[I][1] + 3.0 * YP
                    [I][0])
            elif NSTP == 4:
                YP[I][1] = YP[I][3]
                YAPX[I][1] = Y[I]
                YOLD[I] = Y[I]
                Y[I] = YOLD[I] + D2 * (3.0 * YP[I][1] - YP[I][0])
            elif NSTP == 5:
                Y[I] = YOLD[I] + D12 * (5.0 * YP[I][3] + 8.0 * YP[I][1] -
                    YP[I][0])
            elif NSTP == 6:
                YP[I][2] = YP[I][3]
                YOLD[I] = Y[I]
                YAPX[I][2] = Y[I]
                Y[I] = YOLD[I] + D12 * (23.0 * YP[I][2] - 16.0 * YP[I][1] +
                    5.0 * YP[I][0])
            elif NSTP == 7:
                YAPX[I][0] = YAPX[I][1]
                YAPX[I][1] = YAPX[I][2]
                Y[I] = YOLD[I] + D24 * (9.0 * YP[I][3] + 19.0 * YP[I][2] - 
                    5.0 * YP[I][1] + YP[I][0])
                YAPX[I][2] = Y[I]
    else:
        for I in range(0, 3):
            YAPX[I][0] = YAPX[I][1]
            YAPX[I][1] = Y[I]
            YOLD[I] = Y[I]
            Y[I] = YOLD[I] + D24 * (55.0 * YP[I][3] - 59.0 * YP[I][2] + 
                37.0 * YP[I][1] - 9.0 * YP[I][0])
            YAPX[I][2] = Y[I]
            for J in range(0, 3):
                YP[I][J] = YP[I][J + 1]
