function pot(varargin)

if sum(double(size(varargin) ~= [1 4]))
    disp(' ')
    disp('    COMMAND: pot(data, initialparameters, ''id'', regression)')
    disp(' ')
    disp('    data(:,1) = volume_NaOH                 (cm3)')
    disp('    data(:,2) = pH')
    disp('    data(1,3) = vol_analyte                 (cm3)')
    disp('    data(2,3) = [NaOH]                      (mol/dm3)')
    disp('    data(3,3) = pKw')
    disp('    data(4,3) = temperature                 (ºC)')
    disp(' ')
    disp('    initialparameters(1) = [H2A]            (mol/dm3)')
    disp('    initialparameters(2) = pKa1')
    disp('    initialparameters(3) = pKa2')
    disp('    initialparameters(4) = excess_vol_NaOH  (cm3)')
    disp(' ')
    disp('    id = identifying text string')
    disp(' ')
    disp('    regression = 0 => only plots')
    disp('    regression = 1 => regression and plots')
    disp(' ')
    return
end

X = varargin{1}; beta0 = varargin{2}; id = varargin{3}; regression = varargin{4};
beta0 = [beta0(1) 10^-beta0(2) 10^-beta0(3) beta0(4)];
X(1,4) = 0;                                                           % sinal para desenhar diagrama de distribuição

if regression == 1                                                    % regressão não linear
    options = statset('MaxIter',100,'TolX',1e-8,'TolFun',1e-8);
    [beta,r,J,COVB,mse] = nlinfit(X,X(:,2)',@poting,beta0,options);
    errs95 = nlparci(real(beta),real(r),'jacobian',real(J));          % intervalos de incerteza dos parâmetros com 95% de confiança
    X(1,4) = 1;
    [yfit Vline pHline Eqvolx] = poting(beta,X);
else
    beta = beta0;
    errs95 = zeros(4,2);
    X(1,4) = 1;
    [yfit Vline pHline Eqvolx] = poting(beta0,X);                            % cálculo da função com valores inicias de parâmetros
end

% estatísticas dos resultados do ajuste
% cc = real(corrcoef([yfit' X(:,2)]));            % coeficiente de correlação NÃO É adequado para regressão não-linear
r = yfit' - X(:,2);                               % resíduos
errs95 = [errs95; zeros(4,2)];                    % intervalos de incerteza dos parâmetros com 95% de confiança
errs95(2:3,:) = -log10(errs95(2:3,:));
errs95 = abs((errs95(:,2)-errs95(:,1))/2);
[hX2 pX2 statsX] = chi2gof(r,'cdf',{@normcdf,0,std(r)});       % estatística de teste chi-quadrado dos resíduos, nível de significância 5%
[htS ptS ci statst] = ttest(r,0);                              % estatística de teste t-Student dos resíduos, nível de significância 5%

% texto dos resultados, parâmetros
espa = ['  ';'  ';'  ';'  ';'  ';'  ';'  ';'  '];
col1 = ['   [H2A]analyte'; ...
        '           pKa1'; ...
        '           pKa2'; ...
        '              v'; ...
        '       Vanalyte'; ...
        '  [NaOH]titrant'; ...
        '            pKw'; ...
        '    Temperature'];
col2 = num2str([beta(1);-log10(beta(2));-log10(beta(3));beta(4);X(1,3);X(2,3);X(3,3);X(4,3)]);
col3 = num2str(errs95);
col4 = ['M      ';'       ';'       ';'ml     ';'ml     ';'M      ';'       ';'ºC     '];
A = '                            conf. interv. 95%';
B = [col1 espa col2 espa col3 espa col4];
texto = strvcat(A,' ',B);
disp(' ');
disp(texto);
disp(' ');
disp(Eqvolx);
disp(' ');

% texto dos resultados, resíduos
D = ['residue average                  ' num2str(mean(r)) ' (pH)'];
E = ['residue module average           ' num2str(mean(abs(r))) ' (pH)'];
F = ['standard error of regression     ' num2str(sqrt(sum(r.*r)/(length(r)-1))) ' (pH)'];      % para o modelo ajustado é o mesmo que std(r) 
G = ['% positive residues:             ' num2str(100*sum(r > 0)/length(r))];
H = ['% negative residues:             ' num2str(100*sum(r < 0)/length(r))];
I = ['t-statistic for residue sample   ' num2str(statst.tstat)];
J = ['p-value for t-statistic          ' num2str(ptS)];
% K = ['chi quadrado para resíduos       ' num2str(statsX.chi2stat)];
% L = ['valor p para chi quadrado        ' num2str(pX2)];
% textor = strvcat(D,E,F,G,H,I,J,K,L);
textor = strvcat(D,E,F,G,H,I,J);
disp(textor);
disp(' ');

% inserir texto em gráficos das curvas titulação
subplot(4,2,1:4)
title([id '  ' datestr(now)],'Interpreter','none');
% text(mean(X(:,1))-1,min(X(:,2)+2),texto,'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');
xlabel('NaOH volume (ml)');
ylabel('pH');

% gráficos dos resíduos da regressão e inserção de texto
subplot(4,2,6), normplot(r)
xlabel('residues (pH)');
ylabel('probability');
subplot(4,2,7), plot(r(1:length(r)-1),r(2:length(r)),'-o')
title('Sequencial Self-correlation')
xlabel('residues (pH)');
ylabel('residues (pH)');
subplot(4,2,8), plot(X(:,1),r,'-o')
title('Residue Sample');
xlabel('NaOH volume (ml)');
ylabel('residues (pH)');
% text(mean(X(:,1)),mean(r),textor,'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');


1==1;


function [pH Vline pHline texto] = poting(par,X)

% CONTÉM FUNÇÃO A AJUSTAR. PRODUZ GRÁFICO COM PONTOS EXPERIMENTAIS E LINHA
% DO AJUSTE PARA VALORES CORRENTES DOS PARÂMETROS

% renomeação de variáveis
V = X(:,1);                                       % volume de NaOH (ml)
pHexp = X(:,2);                                   % pH experimental
VHAoo = X(1,3);                                   % volume inicial da solução de glicina
OHoo = X(2,3);                                    % concentração do NaOH titulante
Kw = 10^(-X(3,3));                                % constante de autoptotólise da água
a = par(1); b = par(2); c = par(3); d = par(4);   % parâmetros
dd = X(1,4);

% cálculo de pH estimado para os valores de parâmetros par() e volume V de NaOH adicionado.
% Cálculo da função de ajuste com 200 pontos (pHline).
passo = (max(V) - min(V))/200;
Vline = min(V):passo:max(V);
pH = [];
pHline = []; kapa = []; dabliu = []; ipsilon = []; h2a = [];
Vs = [V; Vline'];
for ix = 1:length(Vs)
%     OHo = OHoo * (Vs(ix) - d+(a*VHAoo/OHoo)) / (VHAoo + Vs(ix));
    OHo = OHoo * (Vs(ix) - d) / (VHAoo + Vs(ix));
    HAo = a * VHAoo / (VHAoo + Vs(ix));
    c4 = -1;
    c3 = -b -OHo;
    c2 = -b*c +b*HAo +Kw -OHo*b;
    c1 = 2*b*c*HAo +b*Kw -OHo*b*c;
    c0 = b*c*Kw;
    y1 = roots([c4 c3 c2 c1 c0]);
    y1 = y1(y1 > 0);
    if length(y1) > 1
        beep
        disp(' ')
        disp(' ')
        error('More than one positive root for quartic polynomial! Try different initial parameter values.');
    end
    if ix > length(V)
        pHline = [pHline; -log10(y1)];
        if dd
            ze = y1;
            kapa1 = OHo-Kw/ze;
            dabliu1 = b*c*HAo/(b*ze+b*c+ze^2);    %c*(ze+kapa1)/(ze+2*c);
            ipsilon1 = b*HAo*ze/(ze^2+b*c+b*ze);  %dabliu1*ze/c;
            h2a1 = HAo-ipsilon1-dabliu1;            %ipsilon1*ze/b;
            dabliu1 = dabliu1 * (VHAoo + Vs(ix));
            ipsilon1 = ipsilon1 * (VHAoo + Vs(ix));
            h2a1 = h2a1 * (VHAoo + Vs(ix));
            kapa = [kapa; kapa1];
            dabliu = [dabliu; dabliu1];
            ipsilon = [ipsilon; ipsilon1];
            h2a = [h2a; h2a1];
        end
    else
        pH = [pH; -log10(y1)];
    end

end
pH = pH';

ntotal = ipsilon + dabliu + h2a;
ntotal = mean(ntotal);
ipsilon = ipsilon/ntotal; dabliu = dabliu/ntotal; h2a = h2a/ntotal; 

if dd
    figure(2); clf;
    plot(pHline,h2a,pHline,ipsilon,pHline,dabliu);
    title('Glycine Species Mole Fractions');
    xlabel('pH');
    ylabel('mole fraction');
    hdd = legend('H_2G^+','HG','G^-');
    set(hdd,'Location','North');
end

    
% gráficos
figure(1); clf;

% gráfico dos dados experimentais e da função a ajustar, para valores de parâmetros correntes
subplot(4,2,1:4)
plot(V,pHexp,'or','MarkerFaceColor','r','MarkerSize',6);
hold on
plot(Vline,pHline,'-k','LineWidth',2);
axis([min(Vline),max(Vline),min(pHline)-0.5,max(pHline)+0.5]);

% gráfico da segunda derivada da função a ajustar, para os valores correntes dos parâmetros
subplot(4,2,5)
passo = Vline(2) - Vline(1);
xprimderiv = Vline(1:length(Vline)-1) + passo/2;
xsegnderiv = xprimderiv(1:length(xprimderiv)-1) + passo/2;
yprimderiv = diff(pHline)/passo;
ysegnderiv = diff(pHline,2)/passo^2;
veq = [];
for iyprim = 3:length(yprimderiv)
    if yprimderiv(iyprim-2) < yprimderiv(iyprim-1) && ...
            yprimderiv(iyprim-1) > yprimderiv(iyprim);
        veq = [veq Vline(iyprim-1)];
    end
end
% plot(xprimderiv,log10(yprimderiv));
plot(xsegnderiv,ysegnderiv);
axis([min(xprimderiv),max(xprimderiv),min(log10(yprimderiv)),max(log10(yprimderiv))]);
espa = ['  ';'  '];
col1 = ['           Vendpoint1';'           Vendpoint2'];
if length(veq) == 1;
    veq(2) = veq; veq(1) = 0;
end
col2 = num2str([veq(1);veq(2)]);
col4 = ['ml     ';'ml     '];
texto = [col1 espa col2 espa col4];
% text(min(xprimderiv),0,texto,'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');
xlabel('NaOH volume (ml)');
ylabel('log10( d(d(pH))/d(d(vol. NaOH)) )');
title('Second Derivative');

1==1;




