function [x,y,yprime] = manipulate_geometry(X,Y)

    % Read data
    close all

    % Remove first point
    x = X(2:end);
    y = Y(2:end);

    %%%%% Compute y' and y'' %%%%%%%%%%%%%

        dx = 1; % 1 meter
        Nx = length(x);

        %First derivative 6th order
        x1 = 0.70127127127127;  

        D1x = (1/60*diag(ones(Nx-3,1),3)-9/60*diag(ones(Nx-2,1),2)+45/60*diag(ones(Nx-1,1),1)-45/60*diag(ones(Nx-1,1),-1)+9/60*diag(ones(Nx-2,1),-2)-1/60*diag(ones(Nx-3,1),-3));
        D1x(1:6,1:9) = [-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ...
             715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ...
             x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; ...
             -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, ...
             -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, ...
             -34560/12013*x1+240569/120130, 0, 0, 0; 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, 25920/2711*x1-182261/27110, 0, 0, 0; -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, -17280/5359*x1+710473/321540, 72/5359, 0, 0; 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, 0, 8640/7877*x1, -1296/7877, 144/7877, 0; -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602, -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801];
        D1x(Nx-5:Nx,Nx-8:Nx) = rot90( -D1x(1:6,1:9),2);
        D1x=D1x/dx;


        %Construct 1st derivative operator in x direction
        D1x = sparse(D1x);

        %Compute derivatives
        yprime = D1x*y;
        ybis = D1x*yprime;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                %         % Change points 2 to 9 to get y'(x0) = 0.
                %         % Use 3rd degree polynomial, cubic spline
                %     
                %         x0 = x(1); x1 = x(10);
                %         y0 = y(1); y1 = y(10);
                %         yprime1 = yprime(10);
                %         
                %         A = [x0^3, x0^2, x0, 1; x1^3, x1^2, x1, 1; 3*x0^2, 2*x0, 1, 0; 3*x1^2, 2*x1, 1, 0];
                %         B = [y0;y1;0;yprime1];
                %         v = A\B;
                %         a1 = v(1); a2 = v(2); a3 = v(3); a4 = v(4);
                %     
                %         y(2:9) = a1*x(2:9).^3 + a2*x(2:9).^2 + a3*x(2:9) + a4;
                
        % Change points 2 to 9 to get y'(x0) = 0.
        % Use 4th degree polynomial spline

        x0 = x(1); x1 = x(10);
        y0 = y(1); y1 = y(10);
        yprime1 = yprime(10);
        ybis1 = ybis(10);

        A = [x0^4, x0^3, x0^2, x0, 1; x1^4, x1^3, x1^2, x1, 1; 4*x0^4, 3*x0^2, 2*x0, 1, 0; 4*x1^3, 3*x1^2, 2*x1, 1, 0; 12*x1^2, 6*x1, 2, 0, 0];
        
        B = [y0;y1;0;yprime1;ybis1];
        v = A\B;
        a1 = v(1); a2 = v(2); a3 = v(3); a4 = v(4); a5 = v(5);

        y(2:9) = a1*x(2:9).^4 + a2*x(2:9).^3 + a3*x(2:9).^2 + a4*x(2:9) + a5;
%     
%         % Compute yprime again
%         yprime = D1x*y;
%         
%         % Compute yprime exactly on spline
%         yprime(2:9) = 4*a1*x(2:9).^3 + 3*a2*x(2:9).^2 + 2*a3*x(2:9) + a4;
        
%         % TEST
%         I_C2F = Interp_1to2_4th_order(length(y));
%         y_new = I_C2F*y;
%         x_new = 1:0.5:1999;
%         
%         plot(x,y,'xr');
%         hold on
%         pause
%         plot(X,Y,'ob');
%         pause
% %         plot(x_new,y_new,'*g')
% %         
%         figure
%         plot(x,yprime);
%         hold on
%         pause

%         yprime(1:4) = 3*a1*x(1:4).^2 + 2*a2*x(1:4) + a3;
%         hold on
%         plot(x,yprime,'r');
    


end
