function [llh, xyz] = compute_ATS_trajectory(ats, zda, option)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  WGS-84 Constants
WGS84.semi_major_axis   = 6378137; 
WGS84.f_inv             = 298.257223563;
WGS84.c_light           = 299792458;
WGS84.pi                = 3.1415926535898;
WGS84.gravitation       = 3.986005e-14;
WGS84.earth_rate        = 7.292115147e-5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% adjust the zda time (used for the ATS time tag) by adding 86400 seconds
% (as required) until the zda time 'agrees' with the GGA time from the RTK
% receiver (which accounts for the day of the week).  This assumes that the
% ATS/ZDA data collection starts within 2 hours of the RTK data collection!
% while( abs(temp_data.time(1)-data.data_1.zda.zda_GPS_timeSS(1)) > 7200)
%     data.data_1.zda.zda_GPS_timeSS = data.data_1.zda.zda_GPS_timeSS + 86400;
% end
zda.data.zda_GPS_timeSS = zda.data.zda_GPS_timeSS + option.data.GPSday_start*86400;

% interpolate to establish the zda time at the time that the ATS
% measurements were made, based on the common pc clock time between the two
% files (and already accounting for the 83ms latency for the ATS
% measurements).
% "interp_GPS_time" now contains the best estimate for the GPS time that
% would have been recorded for each ATS record
interp_GPS_time = interp1(zda.data.zda_pctime, zda.data.zda_GPS_timeSS, ats.data.ats_pctime);


% compute the lat/lon/height of the ATS phase center
lph_Orig_ATS = [ats.data.lon_rad, ...
        ats.data.lat_rad, ...
        ats.data.IH_ht];

start = floor(zda.data.zda_GPS_timeSS(1));
stop  = ceil(zda.data.zda_GPS_timeSS(end));

% time_array contains a continuous vector from the start time to the stop
% time with a data spacing equal to the GPS data rate (typically 1Hz)
time_array = [start : option.data.data_rate : stop]';

% limit the ATS data based on time duration and ATS mode (sync only)
ats_index = find( (interp_GPS_time>=start)&(interp_GPS_time<=stop)&(ats.data.ats_state==69));
ats_gps_time = interp_GPS_time(ats_index);

% interpolate the data for every epoch based on the data rate
interp_vang   = interp1(interp_GPS_time(ats_index), ats.data.ats_vang(ats_index), time_array);
interp_hang   = interp1(interp_GPS_time(ats_index), ats.data.ats_hang(ats_index), time_array);
interp_range  = interp1(interp_GPS_time(ats_index), ats.data.ats_range(ats_index), time_array);

% find the intersection between the time_array data and the RTK data
% [c,file_ats_index,file_rtk_index] = intersect(time_array, temp_data.time);


% create a vector for elevation, angle, slope distance from the ATS
% trajectory data file corresponding to the GPS times from the GGA data
eas = [(90-interp_vang)*pi/180.0, ...
        interp_hang*pi/180.0, ...
        interp_range];

% convert the angles and ranges to lambda/phi/h based on the original
% coordinates for the ATS phase center
[ats.lph] = eas2lph(eas, lph_Orig_ATS, WGS84.semi_major_axis, WGS84.f_inv);
[templph.lph] = eas2lph(eas, lph_Orig_ATS, WGS84.semi_major_axis, WGS84.f_inv);

% adjust the longitude to be between -pi and +pi, not 0 and 2*pi
temp = find(ats.lph(:,1)>pi);
ats.lph(temp,1) = -1*(2*pi-ats.lph(temp,1));

[llh(:,1)] = time_array;
[llh(:,2)] = ats.lph(:,1);
[llh(:,3)] = ats.lph(:,2);
[llh(:,4)] = ats.lph(:,3);


% find indices of all elements that are not flagged as NaN
%not_NaN = find((isnan(data.data_1.ats.lph(:,1)))==0);
not_NaN = find(time_array>=0);

% print the ATS trajectory data to file
% print_ATS_lph(data, temp_data.time(not_NaN), ...
%     data.data_1.ats.lph(not_NaN,2), ...
%     data.data_1.ats.lph(not_NaN,1), ...
%     data.data_1.ats.lph(not_NaN,3), ...
%     'lph');

%temp.lph = data.data_1.ats.lph;
templph.lph(:,3) = templph.lph(:,3) + ats.data.v_offset;

[xyz_temp] = lph2xyz(templph.lph, WGS84.semi_major_axis, WGS84.f_inv);
[xyz(:,1)] = time_array;
[xyz(:,2)] = xyz_temp(:,1);
[xyz(:,3)] = xyz_temp(:,2);
[xyz(:,4)] = xyz_temp(:,3);





return;