%Convert ECEF coordinates (e.g. ITRF2014) to local topocentric system (ENU)

clc

%reference ellipsoid definition
ref  = referenceEllipsoid('GRS80');

%Reference coordinates 
fref = fopen('../vstup/CADM_coordinates_2018.txt');                                                                               
tref = textscan(fref, '%f %f %f','MultipleDelimsAsOne',1);                 
fclose(fref);                                                                                
Xref = tref{:,1};                                                                                              
Yref = tref{:,2};                                                                                              
Zref = tref{:,3};

%Output file
out = fopen('../vstup/2018_24_G-enu.txt', 'w')
Xreff = Xref;
Yreff = Yref;
Zreff = Zref;

%Input file  
f01 = fopen('../vstup/2018_24_G-dis.txt');                                                                        
t01 = textscan(f01, '%s %s %f %f %f %s %s %s %s %s %s %s %s %s %s','MultipleDelimsAsOne',2);         
fclose(f01);                                                                                         
D01 = t01{:,1};                                                                                      
T01 = t01{:,2};                                                                                      
X01 = t01{:,3};                                                                                      
Y01 = t01{:,4};                                                                                      
Z01 = t01{:,5};                                                                                      
N01 = t01{:,7};
res = t01{:,6};
i8  = t01{:,8};
i9  = t01{:,9};
i10 = t01{:,10};
i11 = t01{:,11};
i12 = t01{:,12};
i13 = t01{:,13};
i14 = t01{:,14};
i15 = t01{:,15};



hh=1
for i = 1 : 365
    Xref(hh) = Xreff(i); 
    Yref(hh) = Yreff(i);
    Zref(hh) = Zreff(i);
    
    hh=hh+1;
end
      

for i = 1 : length(D01)
    if X01(i) == NaN  %No data
        fprintf(out,'%s %s NaN NaN NaN %s %s %s %s %s %s %s %s %s %s\r\n', ...
            D01{i},...
            T01{i},...            
            res{i},...
            N01{i},...
            i8{i},...
            i9{i},...
            i10{i},...
            i11{i},...
            i12{i},...
            i13{i},...
            i14{i},...
            i15{i});
    else    % ECEF ? ENU transformation
    [refB, refL, refH] = ecef2geodetic(ref, Xref(i), Yref(i), Zref(i));
    [B, L, H] = ecef2geodetic(ref, X01(i), Y01(i), Z01(i));
    [E, N, U] = geodetic2enu(B, L, H, refB, refL, refH ,ref);
    
    fprintf(out,'%s %s %0.4f %0.4f %0.4f %s %s %s %s %s %s %s %s %s %s\r\n', ...
            D01{i},...
            T01{i},...
            E,...
            N,...
            U,...
            res{i},...
            N01{i},...
            i8{i},...
            i9{i},...
            i10{i},...
            i11{i},...
            i12{i},...
            i13{i},...
            i14{i},...
            i15{i});
    end    
end

fclose(out);
disp('Done');
clear

