function [ dopplertime,rangetime] = getZeroDopplerTime(lineTimeInterval,wavelength,orbit_points ,X,Y,Z)
  numOrbitVec = 20;
  firstVecTime = 0.0;
  secondVecTime = 0.0;
  firstVecFreq = 0.0;
  secondVecFreq = 0.0;
        for i = 0:numOrbitVec-1
            sensorPosition.x = orbit_points(i+1,2);
            sensorPosition.y = orbit_points(i+1,3);
            sensorPosition.z = orbit_points(i+1,4);
            sensorVelocity.x = orbit_points(i+1,5);
            sensorVelocity.y = orbit_points(i+1,6);
            sensorVelocity.z = orbit_points(i+1,7);
            currentFreq = getDopplerFrequency(X,Y,Z, sensorPosition, sensorVelocity, wavelength);
            if (i == 0 || firstVecFreq * currentFreq > 0) 
                firstVecTime = orbit_points(i+1,1);
                firstVecFreq = currentFreq;
             else 
                secondVecTime = orbit_points(i+1,1);
                secondVecFreq = currentFreq;
                break;
            end
        end

        if (firstVecFreq * secondVecFreq >= 0.0) 
            disp( 'NonValidZeroDopplerTime');
        end
        lowerBoundTime = firstVecTime;
        upperBoundTime = secondVecTime;
        lowerBoundFreq = firstVecFreq;
        upperBoundFreq = secondVecFreq;
        diffTime = abs(upperBoundTime - lowerBoundTime);
        absLineTimeInterval = abs(lineTimeInterval);%abs(lineTimeInterval)/86400;
        totalIterations = floor(diffTime/ absLineTimeInterval) + 1;
        numIterations = 0;
        while (diffTime > absLineTimeInterval && numIterations <= totalIterations) 

            midTime = (upperBoundTime + lowerBoundTime) / 2.0;
            [sensorPosition, sensorVelocity] = getPositionVelocity(midTime,orbit_points);
            midFreq = getDopplerFrequency(X,Y,Z, sensorPosition, sensorVelocity, wavelength);

            if (midFreq * lowerBoundFreq > 0.0) 
                lowerBoundTime = midTime;
                lowerBoundFreq = midFreq;
            elseif (midFreq * upperBoundFreq > 0.0) 
                upperBoundTime = midTime;
                upperBoundFreq = midFreq;
            elseif (midFreq == 0) 
                  dopplertime= midTime;
                break;
            end

            diffTime = abs(upperBoundTime - lowerBoundTime);
            numIterations=numIterations+1;
        end


        dopplertime= lowerBoundTime - lowerBoundFreq * (upperBoundTime - lowerBoundTime) / (upperBoundFreq - lowerBoundFreq);
        rangetime = computeSlantRange(X,Y,Z,dopplertime,orbit_points);
end
function Dop=getDopplerFrequency(X,Y,Z, sensorPosition, sensorVelocity, wavelength)
xDiff =  X - sensorPosition.x;
yDiff =  Y - sensorPosition.y;
zDiff =  Z - sensorPosition.z;
distance = sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
Dop=2.0 * (sensorVelocity.x * xDiff + sensorVelocity.y * yDiff + sensorVelocity.z * zDiff) / (distance * wavelength);
end
function rangetime = computeSlantRange(X,Y,Z,dopplerTime,orbitStateVectors)
[sensorPos, ~] = getPositionVelocity(dopplerTime,orbitStateVectors);
xDiff = sensorPos.x -  X;
yDiff = sensorPos.y -  Y;
zDiff = sensorPos.z -  Z;
slantRange = sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
rangetime = slantRange/299792458;
end
function [position,velocity] = getPositionVelocity(time,orbitStateVectors)
nv = 8;
m = 20;
dt = (orbitStateVectors(m,1)-orbitStateVectors(1,1)) / (m-1);

                i0 = max(floor(((time - orbitStateVectors(1,1) )/ dt)) - nv/2 + 1, 0);
                iN = min(i0 + nv - 1, m - 1);
                if iN < m - 1
                i0 = i0;
                else
                i0=iN - nv + 1;
                end
                position.x = 0;
                position.y = 0;
                position.z = 0;
                velocity.x = 0;
                velocity.y = 0;
                velocity.z = 0;
            for i = i0:iN
                orbI.time_mjd = orbitStateVectors(i+1,1);
                orbI.x_pos = orbitStateVectors(i+1,2);
                orbI.y_pos = orbitStateVectors(i+1,3);
                orbI.z_pos = orbitStateVectors(i+1,4);
                orbI.x_vel = orbitStateVectors(i+1,5);
                orbI.y_vel = orbitStateVectors(i+1,6);
                orbI.z_vel = orbitStateVectors(i+1,7);
                weight = 1;
                for j = i0:iN
                    if (j ~= i) 
                        time2 = orbitStateVectors(j+1,1);
                        weight = weight * (time - time2) / (orbI.time_mjd - time2);
                    end
                end
                    position.x = position.x + weight * orbI.x_pos;
                    position.y = position.y + weight * orbI.y_pos;
                    position.z = position.z + weight * orbI.z_pos;
                    velocity.x = velocity.x + weight * orbI.x_vel;
                    velocity.y = velocity.y + weight * orbI.y_vel;
                    velocity.z = velocity.z + weight * orbI.z_vel;
            end
end