%% ZFG / ul 4 / 3.5.2009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;
format long g;
close all;

rho = pi/180;

%% ZADANE PARAMETRY %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = 2400.0000; %zadana vyska v [m]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

GM  = 3986005e+8;
C2  = -1082.63e-6;
C4  = -2.37091222e-6;
a   = 6378137;
omg = 7292115e-11;
Rk  = 6371000.79;
E   = 521854.0097;
gmA = 9.7803267715;
gmB = 9.8321863685;
f4  = 0.0000232955287;
FId = (0:2.5:90);
FI = (0:2.5:90)*rho;

%%%%%%%%%%%%%%%%%%%%%%

q = ((omg^2)*(a^3))/GM;
ah = a + H;
qh = ((omg^2)*(ah^3))/GM;

%% PRIKLAD 1 / CLAIRAUTUV SFEROID %%
Uco = (GM/a)*(1 - C2/2 + q/2);

for i = 1:length(FI)
    Ro = a;
    delta = 1;
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    while delta > 0.001
        R = (GM/Uco)*(1 + ((a/Ro)^2)*C2*P2 + ((1/3)*q*((a/Ro)^(-3)))*(1 - P2));
        delta = abs(R-Ro);
        Ro = R;
    end 
    Rc1(i) = R;
end

Uch = (GM/ah)*(1 - C2/2 + qh/2);

for i = 1:length(FI)
    Ro = ah;
    delta = 1;
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    while delta > 0.001
        Rh = (GM/Uch)*(1 + ((ah/Ro)^2)*C2*P2 + ((1/3)*qh*((ah/Ro)^(-3)))*(1 - P2));
        delta = abs(Rh-Ro);
        Ro = Rh;
    end 
    Rc2(i) = Rh;
end

rozdilC = Rc2 - Rc1;

for i=1:(length(FI)-1)
   difC(i) =rozdilC(i+1)-rozdilC(i);
end
zmenaC = [0, difC];
gradC = [difC, 0];
disp('PRUVODIC CLAIRAUTOVA SFEROIDU')
    disp('¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨')
    fprintf('SIRKA          H = 0 m       H =% 5.0d m     rozdil [m]', H)
    disp(':')
       for i=1:length(FI)
        fprintf('%4.0f°%13.4f%13.4f%10.4f%8.4f\n',FId(i), Rc1(i), Rc2(i), rozdilC(i), gradC(i))
       end
disp('*******************************************************')

figure(1)
plot(FId,Rc1,'r','LineWidth',2)
hold on
plot(FId,Rc2,'b','LineWidth',2)
grid on
title('Prubeh hladinove plochy Clairautova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('delka pruvodice [m]')
legend('H = 0 m','H = 2400 m',1)
hold off

figure(2)
plot( FId,rozdilC, 'g', 'LineWidth', 2)
grid on
title('Sbihavost hladinovych ploch Clairautova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('rozdil pruvodicu hladinovych ploch Clairautova sferoidu [m]')

figure(3)
plot(FId, zmenaC, 'y', 'LineWidth', 2)
grid on
title('Zmena sbihavosti hladinovych ploch Clairautova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('zmena sbihavosti [m]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% PRIKLAD 1 / HELMERTUV SFEROID %%
Uho = (GM/a)*(1 - C2/2 + (3/8)*C4 + q/2);

for i = 1:length(FI)
    Ro = a;
    delta = 1;
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    P4 = (35/8)*((sin(FI(i)))^4) - (15/4)*((sin(FI(i)))^2) + 3/8;
    while delta > 0.001
        RH = (GM/Uho)*(1 + ((q/3)*((a/Ro)^(-3)))*(1 - P2) + ((a/Ro)^2)*C2*P2 + ((a/Ro)^4)*C4*P4);
        delta = abs(RH-Ro);
        Ro = RH;
    end 
    RH1(i) = RH;
end

UhH = (GM/ah)*(1 - C2/2 + (3/8)*C4 + qh/2);

for i = 1:length(FI)
    Ro = ah;
    delta = 1;
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    P4 = (35/8)*((sin(FI(i)))^4) - (15/4)*((sin(FI(i)))^2) + 3/8;
    while delta > 0.001
        RHh = (GM/UhH)*(1 + ((qh/3)*((ah/Ro)^(-3)))*(1 - P2) + ((ah/Ro)^2)*C2*P2 + ((ah/Ro)^4)*C4*P4);
        delta = abs(RHh-Ro);
        Ro = RHh;
    end 
    RH2(i) = RHh;
end

rozdilH = RH2 - RH1;

for i=1:(length(FI)-1)
   difH(i) =rozdilH(i+1)-rozdilH(i);
end
zmenaH = [0, difH];
gradH = [difH, 0];

disp('PRUVODIC HELMERTOVA SFEROIDU')
    disp('¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨')
    fprintf('SIRKA          H = 0 m       H =% 5.0d m     rozdil [m]', H)
    disp(':')
       for i=1:length(FI)
        fprintf('%4.0f°%13.4f%13.4f%10.4f%8.4f\n',FId(i), RH1(i), RH2(i), rozdilH(i), gradH(i))
    end
disp('*******************************************************')

figure(4)
plot(FId,RH1,'r','LineWidth',2)
hold on
plot(FId,RH2,'b','LineWidth',2)
grid on
title('Prubeh hladinove plochy Helmertova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('delka pruvodice [m]')
legend('H = 0 m','H = 2400 m',1)
hold off

figure(5)
plot( FId,rozdilH, 'g', 'LineWidth', 2)
grid on
title('Sbihavost hladinovych ploch Helmertova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('rozdil pruvodicu hladinovych ploch Helmertova sferoidu [m]')

figure(6)
plot(FId, zmenaH, 'y', 'LineWidth', 2)
grid on
title('Zmena sbihavosti hladinovych ploch Helmertova sferoidu')
xlabel('sfericka sirka [^o]')
ylabel('zmena sbihavosti [m]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% PRIKLAD 2 %%

%%hladinova koule
Uk = GM/Rk + (1/3)*(omg^2)*(Rk^2); %normalni tihovy potencial

for i = 1:length(FI)               %tihove zrychleni
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    gk = GM/(Rk^2) + (omg^2)*Rk*P2 - (omg^2)*Rk*((cos(FI(i)))^2);
    gK(i) = gk;
end

%%hladinovy elipsoid
b = sqrt(a^2 - E^2);
f2 = ((gmB - gmA)/gmA) - f4;

Uel = (GM/E)*atan(E/b) + (1/3)*(omg^2)*(a^2);

for i = 1:length(FI)
    gel = gmA*(1 + f2*((sin(FI(i)))^2) + f4*((sin(FI(i)))^4));
    gEL(i) = gel;
end

%%Clairautuv sferoid
Ucl = (GM/a)*(1 - C2/2 + q/2);
Rcl = GM/Ucl;

for i = 1:length(FI)
    P2 = (3/2)*((sin(FI(i)))^2) - 1/2;
    gcl = (GM/(Rcl^2))*(1 - (4/3)*q + (C2 + (4/3)*q)*P2);
    gCL(i) = gcl;
end

%%Helmertuv sferoid
Uhe = (GM/a)*(1 - C2/2 + (3/8)*C4 + q/2);

for i = 1:length(FI)
   ghe = 9.78030*(1 + 0.005302*((sin(FI(i)))^2) - 0.000007*((sin(2*FI(i)))^2) );
   gHE(i) = ghe;
end

disp('NORMALNI TIHOVY POTENCIAL')
    disp('¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨')
    fprintf('Clairaut          Helmert           koule             elipsoid')
    disp(':')
    fprintf('%16.7f  %16.7f  %16.7f  %16.7f\n',Ucl, Uhe, Uk, Uel)
    disp('*******************************************************')

disp('NORMALNI TIHOVE ZRYCHLENI')
    disp('¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨')
    fprintf('SIRKA      Clairaut          Helmert           koule             elipsoid')
    disp(':')
       for i=1:10
        fprintf('%4.0f°%16.4f  %16.4f  %16.4f  %16.4f\n',FId(i), gCL(i)*1e5, gHE(i)*1e5, gK(i)*1e5, gEL(i)*1e5)
    end
disp('*******************************************************')

figure(7)
plot(FId, gCL, '--b', 'LineWidth', 5)
hold on
plot(FId, gHE, 'r', 'LineWidth', 3)
plot(FId, gK, 'g', 'LineWidth', 2)
plot(FId, gEL, ':y', 'LineWidth', 2)
grid on
title('Prubeh normalniho tihoveho zrychleni')
xlabel('zemepisna sirka [^o]')
ylabel('normalni tihove zrychleni [m*s^-^2]')
legend('Clairautuv sferoid','Helmertuv sferoid', 'Hladinova koule', 'Hladinovy elipsoid',2)
hold off

%%KONEC%%