%====================================================================
% Script per il Metodo dei 3/8
% Necessita della Function treottavi.m
%====================================================================

clear all
%clc
close all

file = fopen('tabella_38.txt','w');

disp('=======================================================')
disp('INTEGRAZIONE DEFINITA DI UNA FUNZIONE');
%exprf = input('Inserisci la funzione (senza apici) = ','s');
exprf = 'sin(x.^2)';
% trasforma la stringa in anonymous function
fs = ['@(x)',exprf];
f = str2func(fs);
fprintf('Integrazione Numerica - Metodo dei 3/8 \n');
fprintf(file,'Integrazione Numerica - Metodo dei 3/8 \n');
fprintf('Funzione integranda = %s \n', exprf);
fprintf(file,'Funzione integranda = %s \n', exprf);
a = 0;
b = pi/4;
fprintf('Intervallo di integrazione = [ %d , %18.16f ] \n', a,b);
fprintf(file,'Intervallo di integrazione = [ %d , %18.16f ] \n', a,b);
toll = 1.e-15;
fprintf('Tolleranza toll = %2.0e \n', toll);
fprintf(file,'Tolleranza toll = %2.0e \n', toll);
m = 3;
fprintf('Suddivisione iniziale m = %d \n',m);
fprintf(file,'Suddivisione iniziale m = %d \n',m);
h = (b-a)/3;
fprintf('Passo iniziale h = %10.6f \n', h);
fprintf(file,'Passo iniziale h = %10.6f \n', h);
nmax = 15;
fprintf('Numero massimo di iterazioni nmax = %d \n', nmax);
fprintf(file,'Numero massimo di iterazioni nmax = %d \n', nmax);
disp('=======================================================')
% Parte esecutiva

% Inizializzazioni per ciclo while (prima iterazione) 
intvec(1) = treottavi_LIN (f,a,b,m);
fprintf('%s \t %s \t\t\t %s \t\t %s \n\n','   m','h','  Integrale',' diff');
fprintf(file,'\n %s \t %s \t\t\t %s \t\t %s \n\n','   m','h','  Integrale',' diff');
fprintf('%4.0f %16.5e %19.16f \n',[m,h,intvec(1)]);
fprintf(file,'%4.0f %16.5e %19.16f \n',[m,h,intvec(1)]);
inold=intvec(1);
diff = toll+1;
n = 1;
% ciclo while
while (diff > toll) && (n<nmax)
    n = n+1;
    m = m*2;
    h = h/2;
    intvec(n) = treottavi_LIN (f,a,b,m);
    diff = abs(inold-intvec(n));
    fprintf('%4.0f %16.5e %19.16f %11.4e \n',[m,h,intvec(n),diff]);
    fprintf(file,'%4.0f %16.5e %19.16f %11.4e \n',[m,h,intvec(n),diff]);
    inold = intvec(n);
end   
if n==nmax
    disp(' ')
    disp('*** ATTENZIONE RAGGIUNTO IL NUMERO MASSIMO DI APPROSSIMAZIONI ***')
end

% Calcola la soluzione di riferimento
q = integral(f,a,b);

fprintf('\nValore integrale di riferimento q = %20.16f \n',q);
fprintf(file,'\nValore integrale di riferimento q = %20.16f \n',q);
fclose(file);

% Grafici

% Crea la figura che contiene la funzione errore
figure(2);
semilogy([1:length(intvec)],abs(q-intvec),'r');
% inserisce il titolo della figura 
% PRE 2020b
title({'Errori assoluti rispetto alla soluzione di riferimento'; ...
     ['Funzione = ', exprf,' - Metodo 3/8']});
% POST 2020b
%title('Errori assoluti rispetto alla soluzione di riferimento')
%subtitle(['Funzione = ', exprf,' -, Metodo 3/8']);
xlabel('iterazioni');
ylabel('errore assoluto');

%gtext('Nome e Cognome')