In this example we will implement the J2 Perturbation (first-order) propagator which accounts for secular variations in the orbit elements due to Earth oblateness. We do not consider atmospheric drag or solar or lunar gravitational forces.It uses two-line element (TLE) sets as an input. Additionally as an example we wil use the implemented propagator to plot orbits of RADUGA satellites.
% TLE Data Source. http://www.celestrak.com/ January 20, 2012 clc; clear all; Earth3DPlot(1); mu = 398600; % Earth’s gravitational parameter [km^3/s^2] R_earth = 6378; % Earth radius [km] J2 = 0.0010836; we = 360*(1 + 1/365.25)/(3600*24); % Earth's rotation [deg/s] fname = 'Raduga.txt'; % TLE file name % Open the TLE file and read TLE elements fid = fopen(fname, 'rb'); while ~feof(fid) % Reading TLE File L1 = fscanf(fid,'%23c%*s',1); L2 = fscanf(fid,'%d%6d%*c%5d%*3c%*2f%f%f%5d%*c%*d%5d%*c%*d%d%5d',[1,9]); L3 = fscanf(fid,'%d%6d%f%f%f%f%f%f%f%f',[1,9]); epoch = L2(1,4)*24*3600; % Epoch Date and Julian Date Fraction Db = L2(1,5); % Ballistic Coefficient i = L3(1,3); % Inclination [deg] RAAN = L3(1,4); % Right Ascension of the Ascending Node [deg] e = L3(1,5)/1e7; % Eccentricity omega = L3(1,6); % Argument of periapsis [deg] M = L3(1,7); % Mean anomaly [deg] n = L3(1,8); % Mean motion [Revs/day] % Orbital parametres a = (mu/(n*2*pi/(24*3600))^2)^(1/3); % Semi-major axis [km] T = 2*pi*sqrt(a^3/mu)/60; % Period in [min] rp = a*(1-e); h = (mu*rp*(1 + e))^0.5; % Angular momentum E = keplerEq(M,e,eps); theta = acos((cos(E) -e)/(1 - e*cos(E)))*180/pi; % [deg] True anomaly % J2 Propagator dRAAN = -(1.5*mu^0.5*J2*R_earth^2/((1-e^2)*a^3.5))*cosd(i)*180/pi; domega = dRAAN*(2.5*sind(i)^2 - 2)/cosd(i); eps = 1E-9; dt = 5; % time step [sec] ti = 0; j = 1; while(ti <= T); E = 2*atan(tand(theta/2)*((1-e)/(1+e))^0.5); M = E - e*sin(E); t0 = M/(2*pi)*T; t = t0 + dt; M = 2*pi*t/T; E = keplerEq(M,e,eps); theta = 2*atan(tan(E/2)*((1+e)/(1-e))^0.5)*180/pi; RAAN = RAAN + dRAAN*dt ; omega = omega + domega*dt; [R V] = Orbital2State( h, i, RAAN, e,omega,theta); % Considering Earth's rotation fi_earth = we*ti; Rot = [cosd(fi_earth), sind(fi_earth),0;... -sind(fi_earth),cosd(fi_earth),0;0,0,1]; R = Rot*R; ti = ti+dt; Rv(j,:) = R'; j = j+1; end scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.'); end fclose(fid); hold off; % Ground Track % Open the TLE file and read TLE elements EarthTopographicMap( 2,820,420); fid = fopen(fname, 'rb'); while ~feof(fid) [ h, i, RAAN, e,omega,theta ] = TLE2OE( fid ); [Rv, alfa,delta ] = J2PropagR( h, i, RAAN, e,omega,theta) ; scatter(alfa,delta,'.'); end fclose(fid); figure(2); legend('RADUGA 32', 'RADUGA-1 4', 'RADUGA-1 5', 'RADUGA-1 7', 'RADUGA-1M 1'); text(270,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',... 'HorizontalAlignment','left','FontSize',14 ); title('RADUGA satellites ground track');
Input TLE File: Source www.celestrak.com
Raduga.txt RADUGA 32 1 23448U 94087A 13017.58278277 .00000000 00000-0 10000-3 0 5863 2 23448 12.7840 37.1920 0005679 124.1377 235.9487 1.00272544 83823 RADUGA-1 4 1 25642U 99010A 13017.71286665 .00000130 00000-0 10000-3 0 6727 2 25642 11.7705 50.8546 0001080 106.9153 253.1213 1.00296963 50866 RADUGA-1 5 1 26477U 00049A 13019.59753669 -.00000103 00000-0 10000-3 0 5404 2 26477 9.3910 51.0039 0002612 354.7657 5.2534 1.00203902 45393 RADUGA-1 7 1 28194U 04010A 13018.70907761 .00000086 00000-0 10000-3 0 6657 2 28194 6.8757 64.2408 0001484 24.4802 335.5290 1.00365778 32290 RADUGA-1M 1 1 32373U 07058A 13018.90229375 -.00000043 00000-0 00000+0 0 4653 2 32373 0.0263 70.5780 0002023 210.8419 231.8461 1.00272457 19621
Small Satellites 
Where did you get the modeling equations for the J2 Propagator? Is there a description or explanation anywhere?
Orbital Mechanics for Engineering Students Second Edition Aerospace Engineering. This book is a good resource !
You can find all equations you need for J2 there .
Hope it helps !