function [lat,lon,hae]=WGS84displace(lat,lon,depth,N,E,D,degrees)
        if nargin < 7
            degrees = true;
        end
        if ~ degrees
            lat = rad2deg(lat);
            lon = rad2deg(lon);
        end

        [x,y,z] = toECEF(lat,lon,depth);

        phi = atan2(z, sqrt(x*x + y*y));

        slon = sin(deg2rad(lon));
        clon = cos(deg2rad(lon));
        sphi = sin(phi);
        cphi = cos(phi);

        x = x -slon * E - clon * sphi * N - clon * cphi * D;
        y = y + clon * E - slon * sphi * N - slon * cphi * D;
        z = z + cphi * N - sphi * D;

        [lat, lon, hae] = toGeodetic(x, y, z);
                                                                                              
        if (D ~= 0)
            hae = depth + D;
        else
            hae = depth;
        end
