Some key reading if you want to understand this stuff is chapters two to five of the IERS Technical Note 36, the IERS Conventions (2010).
It's not just the J2000/FK5 frame (aka the EME2000 frame) that is associated with some epoch date. Every Earth-centered inertial frame has some epoch date. There are two fundamental reasons why this must be the case:
- The Earth's rotation axis is not constant.
- Solar system astronomers are regularly improving their best estimate of what constitutes an inertial frame.
Note that the J2000/FK5 frame is now twice passé. The current best estimate of what constitutes an inertial frame is the International Celestial Reference Frame 2 (ICRF2). It's predecessor, the ICRF, represented a vast improvement over the J2000/FK5 frame. The ICRF2 is even better than the ICRF. The ICRF was supposed to be co-aligned with the J2000/FK5 frame on J2000.0 (12 noon Terrestrial Time on 2000 January 1). It turned out that this was not the case; there's a slight bias between the frames at the epoch. The J2000/FK5 frame also turns out to be rotating a tiny bit, about 3 milliarcseconds/year. Unless you are doing milliarcsecond astronomy, you can ignore that bias and rotation. For most applications, J2000/FK5=ICRF=ICRF2.
The first item, that the Earth's rotation axis is not constant, is important. The Earth's rotation axis precesses with a period of about 26,000 years. Accounting for the change in the precession between the epoch time and the time of interest (e.g., today) yields a transformation from the epoch frame (e.g. J2000) to the mean of date frame. In addition to this long-term precession, the Earth's axis also displays some shorter term variations in where it points. These short-term variations (from ~5.5 days to 18.6 years) are collectively called nutation. Accounting for the Earth's nutation on top of the precession yields the transformation to the true of date frame. Finally, the Earth rotates at about one revolution per sidereal day about this precessed and nutated axis. Applying this rotation of the true of date frame yields the Earth-centered, Earth-fixed frame. A somewhat widely-used name for the result of this process is the Earth RNP (rotation, nutation, and precession) matrix.
Well, almost. Precession and nutation are semi-analytical models, as is the concept of one revolution per sidereal day. There are some things those models just can't capture.
That one revolution per sidereal day is incorrect for two reasons. One is that the Earth's rotation rate is very gradually slowing down. Another is that when you look at the rotation rate very closely, the Earth sometimes rotates faster than nominal, other times slower. There are two key parameters that describe this, dUT1=UT1−UTC and ΔT=TT-UT1. If you care about this detail I suggest you use the latter as it's continuous. dUT1 has discontinuities at the leap seconds. This is a correction that you add when you compute the rotation part of the RNP matrix.
There are some things that the semi-analytical precession and nutation model just don't cover (yet). The Chandler wobble, for example. These are collectively called "polar motion", and can only be observed (and predicted to some extent). Polar motion needs to be applied after computing the RNP matrix. The full result is sometimes called the PRNP (polar motion, rotation, nutation, and precession) matrix. These fine scale variations in the Earth's orientation, along with dUT1 and ΔT, are called the "Earth Orientation Parameters". These are published on a regular basis as "IERS Bulletins A and B". I'll say more about this below.
You do not want to code this on your own. You can get code to do these calculations from a number of places. The best sites are:
International Astronomical Union (IAU) Standards of Fundamental Astronomy (SOFA)
The SOFA code is the "official" version of all of the concepts described above. You can get both Fortran and "Ctran" (Fortran converted to ugly C) versions of the SOFA code from the SOFA website. Also be sure to check out the cookbooks, particular the cookbook on "SOFA Tools for Earth Attitude".
Naval Observatory Vector Astrometry Software (NOVAS)
The US Naval Observatory is responsible for IERS Bulletins A and B. They have their own software, distinct from the SOFA code. The NOVAS software is available in Fortran, C, and Python. The difference in terms of results is negligible, in the microarcseconds. That's the kind of error one would expect from using double precision and performing the same computations a bit differently.
There are a number of others (e.g., JPL Spice, GSFC GMAT, Orekit), but I suggest you go to the source, and that would be either the IAU or the US Naval Observatory.
I mentioned IERS Bulletins A and B a couple of times above. The International Earth Orientation and Reference Systems Service (IERS) is the worldwide organization responsible for defining things such as the ICRF and for determining how the Earth is oriented. (Yes, the acronym doesn't match. It used to before they changed the name but not the acronym.) As far as those technical bulletins are concerned, they just contain numbers (and a tiny bit of text). These numbers are time-tabulated values for the Earth Orientation Parameters. These bulletins are updated monthly.
A couple of final points:
I put the link to IERS Technical Note 36 at the top of this answer. Read it.
Be very careful of time. There are a number of time scales involved in this modeling. A few of them that you will run into are:
TAI - International Atomic Time. Time according to an earthbound physicist who uses an atomic clock at sea level.
TT - Terrestrial Time. Time according to an earthbound astronomer. Physicists and astronomers disagree by 32.184 seconds.
UT1 - Universal Time. Conceptually, what a sundial says the time is, but smoothed to eliminate things such as the equation of time.
UTC - Coordinated Universal Time. That's what the clock on your computer shows if you are using network time protocol. Next time we have a leap second (probably the end of next year), you'll be able to see a minute with 61 seconds in it. UTC ticks at the same rate as TAI and TT but occasionally has leap seconds so as to stay within a second of UT1.
TCB - Barycentric Coordinate Time. A general relativistic timescale that on average ticks faster than clocks on the surface of the Earth.
TDB - Barycentric Dynamical Time. A general relativistic timescale that on average ticks at the same rate as clocks on the surface of the Earth.
GAST - Greenwich Apparent Sidereal Time. If you run into this time scale you are looking at an out-of-date concept for calculating Earth rotation. Use the newer Earth Rotation Angle concept, which relies on UT1.
Update
I didn't answer the title of the question,
How to determine satellite position in J2000 from latitude, longitude and distance from Earth?
This is the easy part. The only tricky aspect is that latitude is almost always geodetic latitude rather than geocentric latitude. I'll assume that latitude $\phi$, longitude $\lambda$, and altitude $h$ are in the WGS84 reference system. See [Department of Defense (2000), "World Geodetic System 1984: Its Definition and Relationships with Local Geodetic Systems" NIMA TR8350.2]
(https://web.archive.org/web/20200101071826/https://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf)
Equations 4-14 and 4-15 in the above reference describe the transformation from latitude $\phi$, longitude $\lambda$, and altitude $h$ are in the WGS84 reference system to Cartesian Earth-centered, Earth-fixed (ECEF) coordinates. The equations below use two key parameters that describe the shape of the Earth (see tables 3.1 and 3.3 of the above reference):
$$
\begin{aligned}
a &= 6378137\ \text{m} && \text{Earth equatorial radius} \\
e^2 &= 6.69437999014\times10^{-3} && \text{Square of Earth eccentricity}
\end{aligned}
$$
First you need to compute the "radius of curvature in the prime vertical" (equation 4-15 in the reference):
$$N = \frac a {\sqrt{1-e^2\sin^2\phi}}$$
Then simply compute the ECEF coordinates via equations 4-14:
$$
\begin{aligned}
x &= (N+h)\cos\phi \cos\lambda \\
y &= (N+h)\cos\phi \sin\lambda \\
z &= ((1-e^2)N+h) \sin\phi
\end{aligned}
$$
