function y = PSD1D_Function(P,xdata)
% Return the idealised PSD 1D function using the set of parameters P on the
% frequency axis xdata

%---- Supplementary constraint ----%
% (In case solver goes to crazy)
if (P(3)==0)
    A = 1E-99;
else
    A = abs(P(3));
    
    %---- Rescaling units to help solver ----%
    A = A*1E-16;
    
    y = zeros(1,length(xdata))';
    
    %---- Defines the points where the power law is changed ----%
    Node1 = P(1);
    Node2 = P(2);
    
    %---- Continuity relation at the Nodes ----%
    B = A*(power(Node1,P(5)-P(4)));
    C = B*(power(Node2,P(6)-P(5)));
    
    %---- First segment ----%
    index_vec = (xdata<=Node1);
    y(index_vec) = A*power(xdata(index_vec),-P(4));
    
    %---- Second segment ----%
    index_vec = ((xdata>Node1)&(xdata<Node2));
    y(index_vec) = B*(power(xdata(index_vec),-P(5)));
    
    %---- Third segment ----%
    index_vec = (xdata>=Node2);
    y(index_vec) = C*(power(xdata(index_vec),-P(6)));
end