function dy = SatErdeDGL(t,y,m1,Vorf,ae,J2,J3,J4)
% Differentialgleichung

% Vers.2.0  17.12.2013
% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
[mr,mt]=multipol(y,ae,J2,J3,J4);                             % Multipolmomente
%
dy(1) = y(3)/m1;                                             % r
dy(2) = y(4)/(m1*y(1)^2);                                    % theta
dy(3) = y(4)^2/(m1*y(1)^3) - Vorf/(y(1)^2)*(1+mr);           % p_r
dy(4) = Vorf/y(1)*mt;                                        % p_theta
%
dy = dy.';


function [mr,mt]=multipol(y,ae,J2,J3,J4)
% Berechnung der Legendre-Polynome
x = cos(y(2));
%P0 = 1;
%P1 = x;
P2 = 0.5*(3*x^2 - 1);
P3 = 0.5*(5*x^3 - 3*x);
P4 = 0.25*(7*x*P3 - 3*P2);
%
ris = ae/y(1);              % skalierter inverser Abstand
%
mr = 3*J2*ris^2*P2 + 4*J3*ris^3*P3 + 5*J4*ris^4*P4; % radiale Komponente

% Berechnung der Ableitung der Legendre-Polynome
dx = -sin(y(2));
dP2 = dx*3*x;
dP3 = dx*(7.5*x^2 - 1.5);
dP4 = 0.25*(7*(dx*P3 + x*dP3) -3*dP2);
%
mt = J2*ris^2*dP2 + J3*ris^3*dP3 + J4*ris^4*dP4; 
