clear all
p_elipsoidy_GRS80;

%dane
%punkt C
a_fi_d = 53;
a_fi_m = 38;
a_fi_s = 0;

a_la_d = 21;
a_la_m = 22;
a_la_s = 0;

%punkt D
b_fi_d = 53;
b_fi_m = 21;
b_fi_s = 51.9190;

b_la_d = 21;
b_la_m = 9;
b_la_s = 25.8088;

%obliczenia dla 1992
%zamiana na radiany
a_fi = stopnie_d_to_rad(stopnie_to_stopnie_d( a_fi_d, a_fi_m, a_fi_s ));
a_la = stopnie_d_to_rad(stopnie_to_stopnie_d( a_la_d, a_la_m, a_la_s ));

b_fi = stopnie_d_to_rad(stopnie_to_stopnie_d( b_fi_d, b_fi_m, b_fi_s ));
b_la = stopnie_d_to_rad(stopnie_to_stopnie_d( b_la_d, b_la_m, b_la_s ));

%wyznaczenie poludnika strefowego dla ukladu 1992
L_0 = 19;
L_0 = stopnie_d_to_rad(L_0);

%wyznaczenie wsp na płaszczyznie gaussa 
[a_x_gk,a_y_gk] = bl_gk(a_fi,a_la,L_0,a,e2);
[b_x_gk,b_y_gk] = bl_gk(b_fi,b_la,L_0,a,e2);

%zamiana na wspołrzędne w układzie 1992
[a_x_1992, a_y_1992] = cechowanie_1992(a_x_gk, a_y_gk);
[b_x_1992, b_y_1992] = cechowanie_1992(b_x_gk, b_y_gk);

%obliczenie odwrotne
[a_x_gk_test, a_y_gk_test] = odcechowanie_1992(a_x_1992, a_y_1992);
[b_x_gk_test, b_y_gk_test] = odcechowanie_1992(b_x_1992, b_y_1992);

[a_fi_test, a_la_test] = gk_bl(a_x_gk_test, a_y_gk_test, L_0, a, e2);
[b_fi_test, b_la_test] = gk_bl(b_x_gk_test, b_y_gk_test, L_0, a, e2);

%zamiana na stopnie
[a_fi_d_test, a_fi_m_test, a_fi_s_test] = rad_to_degres(a_fi_test);
[a_la_d_test, a_la_m_test, a_la_s_test] = rad_to_degres(a_la_test);

[b_fi_d_test, b_fi_m_test, b_fi_s_test] = rad_to_degres(b_fi_test);
[b_la_d_test, b_la_m_test, b_la_s_test] = rad_to_degres(b_la_test);

%raport po 1992
disp([' ']);
disp(['________________________________________________________________________________']);
disp(['Raport obliczeniowy programu GeoKał produkcji Kałowego Softwaru. /Michał Kałuski']);
disp(['Układy GK - przeliczenie współrzędnych geodezyjnych na prostokątne w układzie PUWG1992 oraz PUWG2000']);

disp([' ']);
disp(['Zadanie 1']);

disp([' ']);
disp(['Dane wejściowe:']);

disp(['c_fi = ',num2str(a_fi_d), ' d ' , num2str(a_fi_m), ' m ', num2str(a_fi_s, 8), ' s']);
disp(['c_la = ',num2str(a_la_d), ' d ' , num2str(a_la_m), ' m ', num2str(a_la_s, 8), ' s']);
disp([' ']);
disp(['d_fi = ',num2str(b_fi_d), ' d ' , num2str(b_fi_m), ' m ', num2str(b_fi_s, 8), ' s']);
disp(['d_la = ',num2str(b_la_d), ' d ' , num2str(b_la_m), ' m ', num2str(b_la_s, 8), ' s']);

disp([' ']);
disp(['Wyniki 1992:']);

disp(['c_x_gk = ',num2str(a_x_gk),]);
disp(['c_y_gk = ',num2str(a_y_gk),]);
disp([' ']);
disp(['d_x_gk = ',num2str(b_x_gk),]);
disp(['d_y_gk = ',num2str(b_y_gk),]);
disp([' ']);

disp(['c_x_1992 = ',num2str(a_x_1992),]);
disp(['c_y_1992 = ',num2str(a_y_1992),]);
disp([' ']);
disp(['d_x_1992 = ',num2str(b_x_1992),]);
disp(['d_y_1992 = ',num2str(b_y_1992),]);

disp([' ']);
disp(['Kontrola:']);
disp(['c_fi_odwrotne = ',num2str(a_fi_d_test), ' d ' , num2str(a_fi_m_test), ' m ', num2str(a_fi_s_test, 8), ' s']);
disp(['c_la_odwrotne = ',num2str(a_la_d_test), ' d ' , num2str(a_la_m_test), ' m ', num2str(a_la_s_test, 8), ' s']);
disp([' ']);
disp(['d_fi_odwrotne = ',num2str(b_fi_d_test), ' d ' , num2str(b_fi_m_test), ' m ', num2str(b_fi_s_test, 8), ' s']);
disp(['d_la_odwrotne = ',num2str(b_la_d_test), ' d ' , num2str(b_la_m_test), ' m ', num2str(b_la_s_test, 8), ' s']);


disp([' ']);

%obliczenia dla 2000
%poludnik zerowy musi byc w stopniach dziesietnych zeby zadzialal IF
a_la_dz = rad_to_stopnie_d(a_la);
b_la_dz = rad_to_stopnie_d(b_la);

%wyznaczenie południka osiowego i numeru strefy
[ a_nr_strefy, a_L_0 ] = poszukiwanie_strefy_2000( a_la_dz );
[ b_nr_strefy, b_L_0 ] = poszukiwanie_strefy_2000( b_la_dz );

% do dalszych funkcji potrzebujemy juz poludnika zerowego ale w radianach
a_L_0_d = a_L_0;
b_L_0_d = b_L_0;
a_L_0 = stopnie_d_to_rad(a_L_0);
b_L_0 = stopnie_d_to_rad(b_L_0);

%wyznaczenie wsp na plaszczyznie gaussa
[a_x_gk,a_y_gk] = bl_gk(a_fi,a_la,a_L_0,a,e2);
[b_x_gk,b_y_gk] = bl_gk(b_fi,b_la,b_L_0,a,e2);

%zamiana na wspołrzędne w układzie 2000
[a_x_2000, a_y_2000] = cechowanie_2000(a_x_gk, a_y_gk, a_nr_strefy);
[b_x_2000, b_y_2000] = cechowanie_2000(b_x_gk, b_y_gk, b_nr_strefy);

%obliczenie odwrotne
[ a_x_gk_test, a_y_gk_test ] = odcechowanie_2000( a_x_2000, a_y_2000, a_nr_strefy );
[ b_x_gk_test, b_y_gk_test ] = odcechowanie_2000( b_x_2000, b_y_2000, b_nr_strefy );

[a_fi_test, a_la_test] = gk_bl(a_x_gk_test, a_y_gk_test, a_L_0, a, e2);
[b_fi_test, b_la_test] = gk_bl(b_x_gk_test, b_y_gk_test, b_L_0, a, e2);

%zamiana na stopnie
[a_fi_d_test, a_fi_m_test, a_fi_s_test] = rad_to_degres(a_fi_test);
[a_la_d_test, a_la_m_test, a_la_s_test] = rad_to_degres(a_la_test);

[b_fi_d_test, b_fi_m_test, b_fi_s_test] = rad_to_degres(b_fi_test);
[b_la_d_test, b_la_m_test, b_la_s_test] = rad_to_degres(b_la_test);

%raport po 2000
disp(['Wyniki 2000:']);

disp(['poludnik osiowy punktu C = ',num2str(a_L_0_d), ' d']);
disp(['numer strefy punktu C = ',num2str(b_nr_strefy)]);
disp([' ']);
disp(['poludnik osiowy punktu D = ',num2str(b_L_0_d), ' d']);
disp(['numer strefy punktu D = ',num2str(a_nr_strefy)]);
disp([' ']);

disp(['c_x_gk = ',num2str(a_x_gk),]);
disp(['c_y_gk = ',num2str(a_y_gk),]);
disp([' ']);
disp(['d_x_gk = ',num2str(b_x_gk),]);
disp(['d_y_gk = ',num2str(b_y_gk),]);
disp([' ']);

disp(['c_x_2000 = ',num2str(a_x_2000),]);
disp(['c_y_2000 = ',num2str(a_y_2000),]);
disp([' ']);
disp(['d_x_2000 = ',num2str(b_x_2000),]);
disp(['d_y_2000 = ',num2str(b_y_2000),]);

disp([' ']);
disp(['Kontrola:']);
disp(['c_fi_odwrotne = ',num2str(a_fi_d_test), ' d ' , num2str(a_fi_m_test), ' m ', num2str(a_fi_s_test, 8), ' s']);
disp(['c_la_odwrotne = ',num2str(a_la_d_test), ' d ' , num2str(a_la_m_test), ' m ', num2str(a_la_s_test, 8), ' s']);
disp([' ']);
disp(['d_fi_odwrotne = ',num2str(b_fi_d_test), ' d ' , num2str(b_fi_m_test), ' m ', num2str(b_fi_s_test, 8), ' s']);
disp(['d_la_odwrotne = ',num2str(b_la_d_test), ' d ' , num2str(b_la_m_test), ' m ', num2str(b_la_s_test, 8), ' s']);

%zadanie 2
%wyznaczenie sasiedniego pasa odwzorowawczego
if b_nr_strefy == 5
    b_nr_strefy_sasiedni = 6;
elseif a_nr_strefy == 8
    b_nr_strefy_sasiedni = 7;
else
    b_nr_strefy_sasiedni = b_nr_strefy + 1;
end; 

if b_nr_strefy_sasiedni == 5 
    b_L_0 = 15;
elseif b_nr_strefy_sasiedni == 6
    b_L_0 = 18;
elseif b_nr_strefy_sasiedni == 7
    b_L_0 = 21;
elseif b_nr_strefy_sasiedni == 8
    b_L_0 = 24;
end;

b_L_0_d = b_L_0;
b_L_0 = stopnie_d_to_rad(b_L_0);

%wyznaczenie wsp na plaszczyznie gaussa
[b_x_gk,b_y_gk] = bl_gk(b_fi,b_la,b_L_0,a,e2);

%zamiana na wspołrzędne w układzie 2000
[b_x_2000, b_y_2000] = cechowanie_2000(b_x_gk, b_y_gk, b_nr_strefy_sasiedni);

%raport końcowy
disp([' ']);
disp(['Zadanie 2']);
disp([' ']);

disp(['Wyniki:']);
disp(['nowy południk osiowy = ', num2str(b_L_0_d)]);
disp(['nowa strefa = ', num2str(b_nr_strefy_sasiedni)]);
disp(['d_x_2000 w pasie sasiednim = ',num2str(b_x_2000),]);
disp(['d_y_2000 w pasie sasiednim = ',num2str(b_y_2000),]);
disp([' ']);
disp(['Współrzędne w układzie 1992 z oczywistych względów są takie same jak w Zadaniu 1']);

disp([' ']);
disp(['________________________________________________________________________________']);