%% Sample script to test sidf_io
%% JH Taylor - 17 July 2002

clear all;
% define the problem
sim_model = 'lim_filt2';
ns = 4;
mtol = 2; ptol = 5; progress = 1;
% create the vector of amplitudes
for jj=1:3
   a = 4^(jj-2); %% a = .25, 1, 4
   avec(jj) = a;
end;
% create the frequency vector
for ii=1:20
   w = 3*10^((ii-16)/15);  %% w_min = 0.3, w_max = ~5.5
   wvec(ii) = w;
end;
% invoke sidf_io
timer = clock;
[mag,phase] = sidf_io(sim_model,ns,avec,wvec,mtol,ptol,progress);
run_time = etime(clock,timer)
%% plotting wired for #a = 3
subplot(211),loglog(wvec,mag(:,1),wvec,mag(:,2),wvec,mag(:,3));
subplot(211),axis([.2 10 .01 10]);
subplot(211),title('|G(jw,a)| for lim\_filt2 obtained by simulation with Fourier analysis')
subplot(211),text(1.0*wvec(1),0.04,date)
subplot(211),text(1.4,3.10,['a = ',num2str(avec(1))]) 
subplot(211),text(0.8,0.65,['a = ',num2str(avec(2))]) 
subplot(211),text(0.8,0.18,['a = ',num2str(avec(3))]) 
subplot(211),xlabel('frequency (rad/sec)');
subplot(211),ylabel('magnitude');
subplot(212),semilogx(wvec,phase(:,1),wvec,phase(:,2),wvec,phase(:,3));
subplot(212),title('<G(jw,a) for lim\_filt2')
subplot(212),text(1.4,-45.,['a = ',num2str(avec(1)),', ',num2str(avec(2)),', ',num2str(avec(3))])
subplot(212),xlabel('frequency (rad/sec)');
subplot(212),ylabel('phase (deg)');
subplot(212),axis([.2 10 -200 0]);
%% print -dps rf2_out
%% save rf2_out sim_model mtol ptol avec wvec mag phase 
