% Simulation of loudspeaker A + filters
% 100 realizations, 25 scales

%   CBVK Toolbox
%   JHMG, PMSB 2016

clear

arqin='inputs_r100';
%arqout='../../../mat/LoudspeakerA/outputs_AfiltA_r100_s25';
arqout='outputs_AfiltA_r100_s25';
diff='diffloudApfiltA';
                       
scales=10.^([-20:1.5:10 12.5:2.5:20]/20);
fs=5000;

R=100;                  % Number of realizations
ini=1;                  % Initial realization
x0 = zeros(7,1);        % Initial state
                        
%---------------------------------------

[f,g,h,sp]=dloudApfiltA_com;     % Regenerate model to be sure

if matlabpool('size')==0
  matlabpool
end

S=length(scales);

outputs_r = cell(S,1);
outputs = cell(R,1);

for s=1:S
    PP{s}.s = 0;
end

% Parameters and variables for simulation
P0.fs = fs;
P0.r = sp.r;
P0.p = sp.p;


load(arqin);

Ts=1/fs;
% Realization loop 
for r=ini:ini+R-1
    
    u0 = inputs{r};    
    N=length(u0);             % Realization lengths may differ
    tspan=(0:N-1)*Ts;  
    
    % For computation of output of input filter excited with impulses
    P0.e = cumsum(diag(u0)*exp(-Ts*(0:N-1)' * transp(sp.p)));
        
    % Scale loop
    parfor s=1:S
        % Scale and 1/fs factor
        PP{s}.s = scales(s)*Ts;
        
        y = odedev(diff,P0,PP{s},h,tspan,x0);    
        outputs_r{s} = y;        
        
        disp(['Run = ' num2str(r) ' : Scale ' num2str(s) ' computed.']);
    end;
    
    outputs{r} = outputs_r;
    
    save(arqout,'outputs','scales' );
    %disp(['>>>>>> Run ' num2str(tmc) ' took ' num2str(toc/60) ...
     % ' minutes <<<<<<<'])
end;    


