MATLAB: Trouble modelling Earth’s rotation for ‘gravitysp​hericalhar​monic’

aerospacedriftgravitygravitysphericalharmonicMATLABmodelnumerical integrationorbitsimulationverlet

I've been trying to make a simple orbital trajectory simulator using the 'gravitysphericalharmonic' function, and Verlet integration. It worked well at first, but I then tried to take into account the Earth's rotation, and this led to some very significant drifts from the correct trajectory.
To take into account the Earth's rotation, I first converted the satellite's position from cartesian to cylindrical polar coordinates. I then subtracted the angle the Earth had rotated by (from the satellite's argument/theta), and then converted the cylindrical coordinates back into cartesian, before inputting them into the 'gravitysphericalharmonic' function.
I was wondering if anyone could help me figure out what causes the drift?
Edit: I've also tried using a rotation matrix to change reference frames from an absolute one to one relative to the Earth, but I got the same results.
Orbital plot ignoring the Earth's rotation:
Orbital plot trying to take into account the Earth's rotation:

Best Answer

Found the problem - the accelerations I was using were in the Earth reference frame, not the absolute one. I just had to convert them to the absolute reference frame and everything fixed itself.