function Map = Do_Virtual_Map(G,Power_law,varargin)
%Do_Virtual_map() Create a synthetic map according to a parametrised PSD
%   Grid is the Grid object for the calculations
%   Power_law is a vector of length 8 with the power law, see the function
%   2023 version from M. Lejean and J. Degallaix


p = inputParser;
p.FunctionName = 'Create a virtual map';

% Check if the first argument is a grid
addRequired(p,'G',@(x)isa(x, 'Grid'));

% Check if the PSD law is given
addRequired(p,'Power_law',@(x)isnumeric(x))

% Check if the phase is given
p.addParameter('phase',[],@(x)isnumeric(x));

% Check if the frequency interval for the PSD_1D generation is given
p.addParameter('bandwidth',[],@(x)isnumeric(x));

p.parse(G,Power_law,varargin{:})

debug = true;

if ~isempty(p.Results.phase)
    Phase = p.Results.phase;
else
    M = zeros(G.Num_point);
    for i = 1 : G.Num_point
        M(i,:) = pi - 2*pi*rand(1,G.Num_point); %generate random number for the row;
        M(:,i) = -M(i,:); % take the same number and put it in the antisymmetric slot
    end
    Phase = exp(1i*M);
end

FFT_Rad_2D = sqrt(G.D2_FFT_X.^2 + G.D2_FFT_Y.^2);

for pp = G.Half_num_point+1:G.Num_point
    FFT_Radius_Axis(pp-G.Half_num_point) = FFT_Rad_2D(pp,pp);
end

FFT_Radius_Axis_h = G.Axis_FFT(G.Half_num_point+1:end); % will take the radius along the diagonal
df = mean(diff(FFT_Radius_Axis_h));

% Draw the PSD 1D on the new frequency vector
PSD_1D  = PSD1D_Function(Power_law,FFT_Radius_Axis);
PSD_1D(1) = PSD_1D(2);
PSD_1D(isnan(PSD_1D)) = 0;

if debug
    figure(201)
    loglog(FFT_Radius_Axis,PSD_1D)
end

PSD_2D =  interp1(FFT_Radius_Axis,PSD_1D,FFT_Rad_2D,'linear',0);


if debug
    figure(202)
    imagesc(FFT_Rad_2D)
    imagesc(log(PSD_2D))
end

ASD = sqrt(PSD_2D.*G.Length^2)/G.Step;
Map = ifft2(ifftshift(ASD.*Phase),'symmetric');

end
