%
%--------------------------------------------------------------------------
% FTE_FFT2Dsurface_map.m
%
% Test example that uses reads a Finesse mirror map and 
% plots the information is various ways.
%
% This script requires `testmap.txt'.
%
% Part of the SimTools package
% Andreas Freise  12.02.2011
%--------------------------------------------------------------------------
%

% Description: Example for plotting a Finesse mirror map
% Keywords: finesse, example, simulation, mirror, map, PSD, FFT

% make sure map is centered around center of mirror

mapfile='testmap.txt';

[map]=FT_read_surface_map(mapfile);

[ypoints, xpoints]=size(map.data);
points=min(xpoints,ypoints);
map.data=map.data(1:points,1:points);
map.nan=map.nan(1:points,1:points);
map.notnan=map.notnan(1:points,1:points);

xsize=points*map.xstep;  
ysize=points*map.ystep;

%%% check whether we can use offsets here, problems with fftshift I think
[mygrid] = FT_init_grid(points, points, xsize, ysize, 0, 0);

% window only works for square grids at the moment
% (distance from border of map to mirror has been determined by hand)
[window]= FT_hanning_window2D(points,14);

% FFT with hanning window
myfft1=fft2(map.data.*window.data);

% plot windowed mirror map
map.data=map.data.*window.data;
FT_plot_virgo_map(map,0);

fftdata=abs(fftshift(myfft1)).^2;
% normalise for PSD
psd2d=fftdata*map.xstep/window.s2;

%{
%% FFT without window
myfft1=fft2(map.data);
fftdata=abs(fftshift(myfft1)).^2;
psd2d=fftdata*map.xstep/points.^2;
%}

%% compute linear psd
psd1d=zeros(ceil(points/2),1);
psdcount=zeros(ceil(points/2),1);

% average 2d into 1d
centerp=ceil(points/2);
for i=1:centerp
  for j=1:points
     r=window.R(i,j);
     if (r<centerp)
       idx=round(r);
       psd1d(idx)=psd1d(idx)+psd2d(i,j);
       psdcount(idx)=psdcount(idx)+1;
     end
  end
end
idx=find(psdcount~=0);
psd1d(idx)=2*psd1d(idx)./psdcount(idx);


% plot 2D fft ------------------------------------------------------------------
figure(111)
plotlsd=log10(sqrt(psd2d));
surf(mygrid.fft_X,mygrid.fft_Y,plotlsd,'Edgecolor','none');
view(2);
axis square;
axis tight;

% Set the y and x scales to linear and attach axis labels
set(gca, 'YScale', 'lin');
set(gca, 'XScale', 'lin');
xlabel('kx ');
ylabel('ky ');
zlabel('log10(Height) [nm]/sqrt(m)');
hc=colorbar;
ylabel(hc,'log10(Height) [nm]/sqrt(m)');

cmap=FT_pos_colormap();
colormap(cmap);
c_max=max(max(plotlsd));
c_min=min(min(plotlsd));
set(gca, 'CLim', [c_min c_max]);
set(gca, 'ZLim', [c_min c_max]);


%% plot 1D psd
figure(112)
xvec=mygrid.fft_X(1,centerp+1:points);
vpoints=length(xvec);
vstart=min(find(xvec>0));
loglog(xvec(vstart:vpoints),sqrt(psd1d(vstart:vpoints)));
xlabel('kx ');
set(gca,'XLim', [xvec(vstart) xvec(vpoints)]);
ylabel('log10(Height) [nm]/sqrt(m)');
grid on;



