Есть код на матлабе, нужно конвертировать в другой язык (С, С++, Паскаль)

clc
clear
 
mu=398.6*10^3;
delta=66.07*10^3;
 
Rz=6371;
i(1)=deg2rad(58); 
omega(1)=deg2rad(40); 
U=70; 
 
hp = 374 + Rz; 
ha = 374 + Rz; 
 
Sb = (2*23)/(2*2600);
ro = 1.0*10^-10; 
r_m = 6800;
 
del_r = -4*pi*Sb*ro*(r_m^2); 
del_e = 12*(pi^2)*Sb*ro*(r_m^2);
del_vy = 2*pi*Sb*ro*sqrt(mu*r_m);
del_T = -12*(pi^2)*Sb*sqrt((r_m^2)/mu);
 
t(1)=0; 
dU=0.1;
 
lambda1(1)=0.1;  
lambda2(1)=0.1; 
e(1)=sqrt(lambda1(1)^2+lambda2(1)^2);
 
 
q=1;   
 for U=0:0.1:6*pi+0.1
    if e(q)<=0.9
         
        if lambda1(q)>0
            w(q)=atan(lambda2(q)/lambda1(q));
        end
        if lambda1(q)<0
            w(q)=atan((lambda2(q)/lambda1(q))+pi);
        end
        if (lambda1(q)==0)&&(lambda2(q)>0)
            w(q)=pi/2;
        end
        if (lambda1(q)==0)&&(lambda2(q)<0)
            w(q)=-pi/2;
        end
 
        th(q)=U-w(q);
        a = (hp + ha)/2;
        p(q) = a/(1 - e(q))^2;        
        r(q)= (p(q)/(1 + e(q)*cos(U))) + del_r; 
 
x(q)=r(q)*(cos(U)*cos((omega(q)))-sin(U)*sin((omega(q)))*cos((i(q))));
        y(q)=r(q)*(cos(U)*sin((omega(q)))+sin(U)*cos((omega(q)))*cos((i(q)))) + del_vy;
        
z(q)=r(q)*sin(U)*sin((i(q)));
               
T(q)=-((mu*delta)/r(q)^4)*sin(2*U)*sin((i(q)))^2;        S(q)=((mu*delta)/r(q)^4)*(3*sin(U)^2*sin((i(q)))^2-1);
W(q)= 0
 
 
         
        j(q)=(1+(r(q)^3/(mu*p(q))*W(q)*cot((i(q)))*sin(U)))^(-1);
        
 
        p(q+1)=p(q)+dU*((2*j(q)/mu)*r(q)^3*T(q));
        lambda1(q+1)=lambda1(q)+dU*((j(q)*r(q)^2/mu)*(T(q)*(1+r(q)/p(q))*cos(U)+S(q)*sin(U)+r(q)/p(q)*(T(q)*lambda2(q)-W(q)*lambda1(q)*cot((i(q)))*sin(U))));
        lambda2(q+1)=lambda2(q)+dU*((j(q)*r(q)^2/mu)*(T(q)*(1+r(q)/p(q))*sin(U)-S(q)*cos(U)+r(q)/p(q)*(T(q)*lambda1(q)+W(q)*lambda2(q)*cot((i(q)))*sin(U))));
        omega(q+1)=omega(q)-dU*(((j(q)*r(q)^3)/(mu*p(q)))*(sin(U)/sin((i(q))))*W(q));
        i(q+1)=i(q)-dU*(((j(q)*r(q)^3)/(mu*p(q)))*cos(U)*W(q));
        t(q+1)=t(q)+dU*((j(q)*r(q)^2)/sqrt(mu*p(q)));
        e(q+1)=sqrt(lambda1(q+1)^2+lambda2(q+1)^2) + del_e; 
    else
break
end    
q=q+1;
 end
figure(1)
[X,Y,Z]=sphere(36);
X=X*6371;
Y=Y*6371;
Z=Z*6371;
mesh(X,Y,Z);
hold on;
plot3(x,z,y);
xlim([-10000 10000]);
ylim([-10000 10000]);
zlim([-10000 10000]);
xlabel('x');
ylabel('z');
zlabel('y');
 
 figure(2)
 U=0:0.1:6*pi+0.1;
 polar(U,r)
  
 figure(3)
 subplot(3,2,1)
 plot(U,p(1:length(U))),grid;
 xlabel('U');
 ylabel('p');
 subplot(3,2,2)
 plot(U,w(1:length(U))),grid;
 xlabel('U');
 ylabel('w');
 subplot(3,2,3)
 plot(U,omega(1:length(U))),grid;
 xlabel('U');
 ylabel('omega');
 subplot(3,2,4)
 plot(U,i(1:length(U))),grid;
 xlabel('U');
 ylabel('i');
 subplot(3,2,5)
 plot(U,t(1:length(U))),grid;
 xlabel('U');
 ylabel('t');
 subplot(3,2,6)
 plot(U,e(1:length(U))),grid;
 xlabel('U');
 ylabel('e');
 
figure(4)
plot(U,W(1:length(U))),grid
xlabel('U');
ylabel('W');
 
 figure(5)
plot(U,T(1:length(U))),grid
xlabel('U');
ylabel('T');
figure (6)
plot(U,S(1:length(U))),grid
xlabel('U');
ylabel('S');