def planet_approx(JD, planet_flag):
    """Approximate orbital elements of each of the planets

    This function provides the classical orbital elements and the rates of
    change for any of the planets in the range of 3000 BC to 3000 AD. A higher
    accuracy model if provided over the range of 1800 - 2050 AD.

    Parameters
    ----------
    jd : float
        Julian date for the epoch
    planet_flag : int
        Number 0 to 8 for each of the planets (Pluto forever)

    Returns
    -------
    a0 : float
        Semimajor axis in AU at epoch
    adot : float
        Rate of change of semimajor axis in AU/century
    e0 : float
        Eccentricty of orbit at epoch
    edot :  float
        Rate of change of eccentricy per century
    inc0 : float
        inclination in degrees at epoch
    incdot : float
        Rate of change of inclination in degrees/century
    meanL0 : float
        Mean longitude in degrees at epoch
    meanLdot : float
        Rate of change of mean longitude in degrees/century
    lonperi0 : float
        Longitude of perihelion in degrees at epoch
    lonperidot : float
        Rate of change of Longitude of perihelion in degrees/century
    raan0 : float
        Longitude of the ascending node in degrees at epoch
    raandot : float
        Rate of RAAN in degrees/century
    b, c, f, s : floats
        Additional parameters for the outer planets when using the coarse model

    See Also
    --------
    planet_coe : Use this function to compute the orbital elements

    Notes
    -----
    The orbital elements are defined with respect to the mean ecliptic and
    equinox of J2000. The fundamental plane is the ecliptic, or the orbit of
    the Earth about the sun. The fundamental direction is aligned with the
    first point of Ares or the vernal equinox.

    Rather than using this function directly, simply use planet_coe to get the
    classical orbital elements wrt J2000 ecliptic.

    Author
    ------
    Shankar Kulumani		GWU		skulumani@gwu.edu   23 June 2017

    References
    ----------

    .. [1] https://ssd.jpl.nasa.gov/?planet_pos
    .. [2] MEEUS, Jean H.. Astronomical Algorithms. Willmann-Bell,
    Incorporated, 1991.

    Examples
    --------
    An example of how to use the function


    """
    JD_1800AD = 2378497.0
    JD_2050AD = 2469808.0
    JD_3000BC = 625674.0
    JD_3000AD = 2816788.0
    if JD > JD_1800AD and JD < JD_2050AD:
        b = 0.0
        c = 0.0
        f = 0.0
        s = 0.0
        if planet_flag == 0:
            a0 = 0.38709927
            adot = 3.7e-07
            e0 = 0.20563593
            edot = 1.906e-05
            inc0 = 7.00497902
            incdot = -0.00594749
            meanL0 = 252.2503235
            meanLdot = 149472.67411175
            lonperi0 = 77.45779628
            lonperidot = 0.16047689
            raan0 = 48.33076593
            raandot = -0.12534081
        elif planet_flag == 1:
            a0 = 0.72333566
            adot = 3.9e-06
            e0 = 0.00677672
            edot = -4.107e-05
            inc0 = 3.39467605
            incdot = -0.0007889
            meanL0 = 181.9790995
            meanLdot = 58517.81538729
            lonperi0 = 131.60246718
            lonperidot = 0.00268329
            raan0 = 76.67984255
            raandot = -0.27769418
        elif planet_flag == 2:
            a0 = 1.00000261
            adot = 5.62e-06
            e0 = 0.01671123
            edot = -4.392e-05
            inc0 = -1.531e-05
            incdot = -0.01294668
            meanL0 = 100.46457166
            meanLdot = 35999.37244981
            lonperi0 = 102.93768193
            lonperidot = 0.32327364
            raan0 = 0.0
            raandot = 0.0
        elif planet_flag == 3:
            a0 = 1.52371034
            adot = 1.847e-05
            e0 = 0.0933941
            edot = 7.882e-05
            inc0 = 1.84969142
            incdot = -0.00813131
            meanL0 = -4.55343205
            meanLdot = 19140.30268499
            lonperi0 = -23.94362959
            lonperidot = 0.44441088
            raan0 = 49.55953891
            raandot = -0.29257343
        elif planet_flag == 4:
            a0 = 5.202887
            adot = -0.00011607
            e0 = 0.04838624
            edot = -0.00013253
            inc0 = 1.30439695
            incdot = -0.00183714
            meanL0 = 34.39644051
            meanLdot = 3034.74612775
            lonperi0 = 14.72847983
            lonperidot = 0.21252668
            raan0 = 100.47390909
            raandot = 0.20469106
        elif planet_flag == 5:
            a0 = 9.53667594
            adot = -0.0012506
            e0 = 0.05386179
            edot = -0.00050991
            inc0 = 2.48599187
            incdot = 0.00193609
            meanL0 = 49.95424423
            meanLdot = 1222.49362201
            lonperi0 = 92.59887831
            lonperidot = -0.41897216
            raan0 = 113.66242448
            raandot = -0.28867794
        elif planet_flag == 6:
            a0 = 19.18916464
            adot = -0.00196176
            e0 = 0.04725744
            edot = -4.397e-05
            inc0 = 0.77263783
            incdot = -0.00242939
            meanL0 = 313.23810451
            meanLdot = 428.48202785
            lonperi0 = 170.9542763
            lonperidot = 0.40805281
            raan0 = 74.01692503
            raandot = 0.04240589
        elif planet_flag == 7:
            a0 = 30.06992276
            adot = 0.00026291
            e0 = 0.00859048
            edot = 5.105e-05
            inc0 = 1.77004347
            incdot = 0.00035372
            meanL0 = -55.12002969
            meanLdot = 218.45945325
            lonperi0 = 44.96476227
            lonperidot = -0.32241464
            raan0 = 131.78422574
            raandot = -0.00508664
        elif planet_flag == 8:
            a0 = 39.48211675
            adot = -0.00031596
            e0 = 0.2488273
            edot = 5.17e-05
            inc0 = 17.14001206
            incdot = 4.818e-05
            meanL0 = 238.92903833
            meanLdot = 145.20780515
            lonperi0 = 224.06891629
            lonperidot = -0.04062942
            raan0 = 110.30393684
            raandot = -0.01183482
        else:
            print('Incorrect planet flag should be between 0 and 8')
    elif JD > JD_3000BC and JD < JD_3000AD:
        b = 0.0
        c = 0.0
        f = 0.0
        s = 0.0
        if planet_flag == 0:
            a0 = 0.38709843
            adot = 0.0
            e0 = 0.20563661
            edot = 2.123e-05
            inc0 = 7.00559432
            incdot = -0.00590158
            meanL0 = 252.25166724
            meanLdot = 149472.67486623
            lonperi0 = 77.45771895
            lonperidot = 0.15940013
            raan0 = 48.33961819
            raandot = -0.12214182
        elif planet_flag == 1:
            a0 = 0.72332102
            adot = -2.6e-07
            e0 = 0.00676399
            edot = -5.107e-05
            inc0 = 3.39777545
            incdot = 0.00043494
            meanL0 = 181.9797085
            meanLdot = 58517.8156026
            lonperi0 = 131.76755713
            lonperidot = 0.05679648
            raan0 = 76.67261496
            raandot = -0.27274174
        elif planet_flag == 2:
            a0 = 1.00000018
            adot = -3e-08
            e0 = 0.01673163
            edot = -3.661e-05
            inc0 = -0.00054346
            incdot = -0.01337178
            meanL0 = 100.46691572
            meanLdot = 35999.37306329
            lonperi0 = 102.93005885
            lonperidot = 0.3179526
            raan0 = -5.11260389
            raandot = -0.24123856
        elif planet_flag == 3:
            a0 = 1.52371243
            adot = 9.7e-07
            e0 = 0.09336511
            edot = 9.149e-05
            inc0 = 1.85181869
            incdot = -0.00724757
            meanL0 = -4.56813164
            meanLdot = 19140.29934243
            lonperi0 = -23.91744784
            lonperidot = 0.45223625
            raan0 = 49.71320984
            raandot = -0.26852431
        elif planet_flag == 4:
            a0 = 5.20248019
            adot = -2.864e-05
            e0 = 0.0485359
            edot = 0.00018026
            inc0 = 1.29861416
            incdot = -0.00322699
            meanL0 = 34.33479152
            meanLdot = 3034.90371757
            lonperi0 = 14.27495244
            lonperidot = 0.18199196
            raan0 = 100.29282654
            raandot = 0.13024619
        elif planet_flag == 5:
            a0 = 9.54149883
            adot = -3.065e-05
            e0 = 0.05550825
            edot = -0.00032044
            inc0 = 2.49424102
            incdot = 0.00451969
            meanL0 = 50.07571329
            meanLdot = 1222.11494724
            lonperi0 = 92.86136063
            lonperidot = 0.54179478
            raan0 = 113.63998702
            raandot = -0.25015002
        elif planet_flag == 6:
            a0 = 19.18797948
            adot = -0.00020455
            e0 = 0.0468574
            edot = -1.55e-05
            inc0 = 0.77298127
            incdot = -0.00180155
            meanL0 = 314.20276625
            meanLdot = 428.49512595
            lonperi0 = 172.43404441
            lonperidot = 0.09266985
            raan0 = 73.96250215
            raandot = 0.05739699
        elif planet_flag == 7:
            a0 = 30.06952752
            adot = 6.447e-05
            e0 = 0.00895439
            edot = 8.18e-06
            inc0 = 1.7700552
            incdot = 0.000224
            meanL0 = 304.22289287
            meanLdot = 218.46515314
            lonperi0 = 46.68158724
            lonperidot = 0.01009938
            raan0 = 131.78635853
            raandot = -0.00606302
        elif planet_flag == 8:
            a0 = 39.48686035
            adot = 0.00449751
            e0 = 0.24885238
            edot = 6.016e-05
            inc0 = 17.1410426
            incdot = 5.01e-06
            meanL0 = 238.96535011
            meanLdot = 145.18042903
            lonperi0 = 224.09702598
            lonperidot = -0.00968827
            raan0 = 110.30167986
            raandot = -0.00809981
        else:
            print('Incorrect planet flag should be between 0 and 8')
    else:
        raise ValueError('Julian Date is outside of the range 3000BC-3000AD')
    return (a0, adot, e0, edot, inc0, incdot, meanL0, meanLdot, lonperi0,
        lonperidot, raan0, raandot, b, c, f, s)
