test  = { 'square', 'smooth' } ;

nflux = {'ToroBillet'}  ;

nflux1 = { 'LaxFriedrichs', ...
           'force', ...
           'Richtmyer', ...
           'slic', ...
           'GodunovCentered', ...
           'ToroBillet', ...
           'Fromm', ...
           'WarmingBeam', ...
           'Waf' } ;

close all ;

for j=1:length(nflux)
    

  ok = input( ['Start numerical flux: ' nflux{j} ' [q=quit] ' ], 's' ) ;
  if ok == 'q'
    break ;
  end ;

  handleFig = figure('Position',[10 10 1000 600],'color','w') ;

  for i=1:length(test)

    LinearAdvectionModel( 'setup',    -0.4 ) ;
    LinearAdvectionModel( 'set_test', test{i} ) ;
    LinearAdvectionModel( 'set_numerical_flux', nflux{j} ) ;
    LinearAdvectionModel( 'set_limiter', 'superbee' ) ;

    PRB . CFL       = 0.9*LinearAdvectionModel( 'get_cfl' ) ;
    PRB . FinalTime = LinearAdvectionModel( 'get_final_time' ) ;
    PRB . Ncells    = 100  ; % select 100 cell inteval subdivision
    PRB . StepTime  = PRB . FinalTime / 100 ;

    SlopeLimiterModel( 'setup', 'superbee' ) ;
    % others limiters:
    % 'minbee', 'superbee', 'vanleer', 'vanalbada', 'minmod', 'warmingbeam', 'laxwendroff', 'fromm'

    % call the solver
    [ UFRAMES2, TFRAMES ] = Solver(PRB,@LinearAdvectionModel) ;

    PRB . CFL = 0.9 ; % change CFL for MUSCL
    [ UFRAMES1, TFRAMES ] = MUSCLSolver(PRB, @LinearAdvectionModel, @SlopeLimiterModel) ;
    [ Xe, Ue ] = LinearAdvectionModel( 'exact', PRB . FinalTime, 1000 ) ; 

    U1 = UFRAMES1(:,:,end) ;
    U2 = UFRAMES2(:,:,end) ;
    X  = LinearAdvectionModel( 'get_x' ) ;
  
    subplot(length(test),1,i);
    plot(Xe, Ue,'-b', 'LineWidth',2, 'Color', 'green') ; hold on ;
    plot(X,  U1,'-ok','MarkerFaceColor','r');
    plot(X,  U2,'-ob','MarkerFaceColor','b');
    title(['test case: ' test{i} '  solver: ' nflux{j} ' and MUSCL'], ...
          'FontSize', 24, 'FontName', 'Arial' );
    axis([X(1) X(end) -0.2 1.2]) ;
  end
end
