% simulink model to generate FRFs from numerical simulation


% define parameters for data acquisition

dt=0.01;
N = 2048;
T = N*dt;
fs = 1/dt;
% define filter

[b,a] = ellip(16,.1,60,2*pi*fs/2.1,'s'); % input bandlimiter

F = 1;
seed = 1;
Tburst = T/2;    % lenght of burst for random burst input

sim('test1',[0,T])

% find the FRFs

%window the data
w=(sin(pi*(0:N-1)/(N-1)).^2)';
%w=ones(N,1);
x=yout(1:N,1);
f=force(1:N,1);
xw=x.*w;
fw=f.*w;
FX=fft(xw);
FF=fft(fw);
SXF=FF.*conj(FX);
SXX=FX.*conj(FX);
SFF=FF.*conj(FF);
SFX=FX.*conj(FF);
TXF=SXX./SXF;
TXF2=SFX./SFF;
lTXF2=length(TXF)/2;
freq=[0 ((1:lTXF2-1)/lTXF2/2/dt)]';
TXF_p=TXF(1:lTXF2);
mag=abs(TXF_p);
ang=angle(TXF_p)*180/pi;

%figure


% mfile to pick off peaks
set(0,'DefaultUIcontrolUnits','normal')



fmin = min(freq);
fmax = max(freq);
cf = (fmax-fmin)/2;
f1 = cf;
f1_t=cf;
df = 1/T;
numbins = (fmax-fmin)/df;
lnst = ['- ';'- '];
lnco = ['y';'r'];
 

figNumber = figure(1);

% define location for the figure

hand(1) = axes('position',[.1 .2 .8 .7],'visible','on');

% plot the magnitude
semilogy(freq,mag)
title('Txf - Transfer function magnitude'), ...
xlabel('Frequency (Hz)')


ylim = get(hand(1),'ylim');
flim = get(hand(1),'xlim');
yminm=ylim(1);
ymaxm=ylim(2);

% define slider
hand(2) = uicontrol(figNumber, ...
        'style','slider', ...
        'pos',[.15 .01 .2 0.04], ...
        'back','w', ...
        'max',fmax,'min',fmin,'sliderstep',[df/(fmax-fmin) .1], ...
        'userdata',freq, ...
  'call',['f1_t=get(hand(2),''value'');',...
     '[f1,hand(5:6)] = slidfreq(f1_t,hand(2),',...
     'hand(5:6),hand(1),ymaxm,yminm,lnst(2,:),lnco(2,:));',...    
     'set(hand(4),''string'',num2str(f1));', ...
     'set(hand(8),''string'',num2str(mag(find(f1-freq==0))))']);

%define text and edit box for frequency
hand(3)=uicontrol(figNumber, ...
        'style','text', ...
	  'pos',[.15 .05 .1 .04], ...
     'string','Freq:');

hand(4) = uicontrol(figNumber, ...
	'style','edit', ...
	'pos',[.25 .05 .1 .04], ...
	'back','w', ...
	'callback',['f1_t=str2num(get(hand(4),''string''));', ...
		 '[f1,hand(5:6)] = slidfreq(f1_t,hand(2),',...
		    'hand(5:6),hand(1),ymaxm,yminm,lnst(2,:),lnco(2,:));',...    
          'set(hand(4),''string'',num2str(f1));', ...
          'set(hand(8),''string'',num2str(mag(find(f1-freq==0))))']);

% define some elements need for plotting the peakpick line 
hand(5) = line([0 1],[0 1]);
hand(6) = text(0,0,' ');

% define "use mouse" control
hand(7)=uicontrol(figNumber, ...
        'style','pushbutton', ...
        'pos',[.15 .10 .2 .04], ...
        'string','Use Mouse', ...
	'callback',['[f1_t,temp] = ginput(1);',...
      '[f1,hand(5:6)] = slidfreq(f1_t,hand(2),',...
		    'hand(5:6),hand(1),ymaxm,yminm,lnst(2,:),lnco(2,:));',...    
          'set(hand(4),''string'',num2str(f1));', ...
          'set(hand(8),''string'',num2str(mag(find(f1-freq==0))))']);
    
% define box that has the magnitude value    
hand(8)=uicontrol(figNumber, ...
   'style','text', ...
   'back','w', ...
   'pos',[.7 .1 .2 .04], ...
   'string',num2str(mag(find(f1-freq==0))));

hand(9)=uicontrol(figNumber, ...
        'style','text', ...
	  'pos',[.6 .1 .1 .04], ...
     'string','Mag:');
     
% draw the initial vertical line    
[f1,hand(5:6)]=slidfreq(f1_t,hand(2), ...
       hand(5:6),hand(1),ymaxm,yminm,lnst(2,:),lnco(2,:));