%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to analyse the spectrum   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff
global corr;

if exist('corr')==0 
 load corr.out;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%
% TIME STEP = 0.0004 PS %
%%%%%%%%%%%%%%%%%%%%%%%%%
TIME_STEP = 0.0004;
a = size(corr);
STEPS = a(1);
FREQ  = 2*pi/TIME_STEP/STEPS*1E+12;
WAVENUMBER = FREQ/2/pi/3E+10;
X = 1:STEPS;

a = abs(fft(corr));
plot(X,a);

%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT OUT THE SPECTRUM %
%%%%%%%%%%%%%%%%%%%%%%%%%
ss = input('Change scope (1:1000 etc): ','s');
while isstr(ss)
 eval(['X = ' ss ';']);
 plot(X,a(X));
 ss = input('Change scope: ','s');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAYBE YOU WOULD LIKE TO FILTER %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ss = input('Want to filter? ','s');
if isstr(ss) 
 b = a;
 while(isstr(ss))
  b = filtfilt([0.90, 0.10],[1],b);
  plot(X,b(X));
  ss = input('Want to filter more? ','s');
 end;
 ss = input('Want the filted? ','s');
 if isstr(ss) a=b; end;
 plot(X,a(X));
end;

%%%%%%%%%%%%%%
% THE TITLE  %
%%%%%%%%%%%%%%
ss =input('Title: ','s');
if isstr(ss) 
 t=ss; 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT OUT WITH UNITS OF RAD/S %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(X*FREQ,a(X));
v = axis;
axis([min(X*FREQ),max(X*FREQ),v(3),v(4)]);
title(t);
xlabel('Rad/s');
PRINT;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT OUT WITH UNITS OF 1/CM  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(X*WAVENUMBER,a(X));
v = axis;
axis([min(X*WAVENUMBER),max(X*WAVENUMBER),v(3),v(4)]);
title(t);
xlabel('1/cm');
PRINT;

end;









