function mut(varargin)

if sum(double(size(varargin) ~= [1 4]))
    disp(' ')
    disp('COMMAND:')
    disp('    mut(data, initialparameters, ''id'', regression)')
    disp(' ')
    disp('ARGUMENTS:')
    disp('    data(:,1) = time                        (s)')
    disp('    data(:,2) = rotation                    (deg)')
    disp(' ')
    disp('    initialparameters(1) = alfa0            (deg g-1 cm3 dm-1)')
    disp('    initialparameters(2) = alfaeq           (deg)')
    disp('    initialparameters(3) = psi(+)+psi(-)    (s-1)')
    disp(' ')
    disp('    id = identifying text string')
    disp(' ')
    disp('    regression = 0 => only plots')
    disp('    regression = 1 => plots and regression')
    disp(' ')
    return
end

% renomeação de variáveis, ordenamento dos dados em X(:,1:2) por temperaturas cescentes
dados = varargin{1}; beta0 = varargin{2}; id = varargin{3}; regression = varargin{4};
[X(:,1) IX] = sort(dados(:,1));
X(:,2) = dados(IX,2);

if regression == 1
    % com regressão
    options = statset('MaxIter',100,'TolX',1e-8,'TolFun',1e-8);         % opções de regressão não linear
    [beta,r,J,COVB,mse] = nlinfit(X,X(:,2)',@muting,beta0,options);    % regressão não linear
    errs95 = nlparci(real(beta),real(r),'jacobian',real(J));            % estimativa de erros dos parâmetros, intervalo de confiança de 95%
    betas = beta;
    xerrs = num2str([(errs95(:,2)-errs95(:,1))/2]);
else
    % sem regressão
    betas = beta0;
    xerrs = num2str(zeros(3,1));
end

% preparação dos resultados para inserção no gráfico dos dados
[yfit tempog alfag] = muting(betas,X);                               % cálculo de absorvâncias e traçado do gráfico absorvância = f(temperatura)
data = num2str(real(betas)');
lege1 = ['        alfa0';'       alfaeq';'psi(+)+psi(-)'];
A = [lege1 ['  ';'  ';'  '] data ...
    ['  ';'  ';'  '] xerrs ...
    ['   deg g-1 cm3 dm-1    ';'   deg                 ';'   s-1                 ']];
B = '                         conf. intrv. 95%';
subplot(4,2,1:4)
title([id '  ' datestr(now)],'Interpreter','none');
% text(min(X(:,1)),max(X(:,2)),strvcat(B,A), ...
%     'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');      % tabela de resultados a inserir nos gráficos
disp(' ');
disp(strvcat(B,A));

% preparação dos resultados para inserção no gráfico dos resíduos
r = yfit' - X(:,2);                                            % cálculo dos resíduos
% cc = real(corrcoef([yfit' X(:,2)]));                         % coeficiente de correlação NÃO É adequado para regressão não-linear 
[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 das estatísticas dos resíduos
D = ['residue average                  ' num2str(mean(r)) ' (deg)'];
E = ['residue module average           ' num2str(mean(abs(r))) ' (deg)'];
F = ['standard error of regression     ' num2str(sqrt(sum(r.*r)/(length(r)-1))) ' (deg)'];      % 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(' ');
disp(textor);
disp(' ');

% gráficos de análise dos resíduos
subplot(4,2,6), normplot(r)
xlabel('residues (deg)');
ylabel('probability');
subplot(4,2,7), plot(r(1:length(r)-1),r(2:length(r)),'-o');
title('Sequencial Self-correlation');
xlabel('residues (deg)');
ylabel('residues (deg)');
subplot(4,2,8), plot(X(:,1)-273.15,r,'-o');
% text(mean(X(:,1))-273.15,mean(r),textor, ...
%     'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');      % insere tabela de resultados estatísticos no gráfico dos resíduos
title('Residue Sample');
xlabel('time (s)');
ylabel('residues (deg)');




function [alfac tempog alfag] = muting(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
tempox = X(:,1);      % tempos experimentais (s)
alfax = X(:,2);       % rotações experimentais (graus)

% parâmetros ajustáveis
alfa0 = par(1);
alfaeq = par(2);
k1maisk2 = par(3);

% cálculo de rotações estimadas (alfac) para os valores de parâmetros
% correntes (alfa0, alfaeq, k1+k2) e para os tempos experimentais.
alfac = [];
for t = tempox'
    alfa = alfaeq + (alfa0 - alfaeq)*exp(-(k1maisk2)*t);
    alfac = [alfac; alfa];
end
    alfac = alfac';

% produção de linha de ajuste ligando 50 pontos
passo = (max(tempox) - min(tempox))/50;
tempog = min(tempox):passo:max(tempox);
alfag = [];
for t = tempog
    alfa = alfaeq + (alfa0 - alfaeq)*exp(-(k1maisk2)*t);
    alfag = [alfag; alfa];
end

% gráficos
figure(1); clf;

% gráfico dos dados experimentais e da linha de ajuste
subplot(4,2,1:4)
plot(tempox,alfax,'or','MarkerFaceColor','r','MarkerSize',6);
hold on
plot(tempog,alfag,'-k','LineWidth',2);
xlabel('time (s)');
ylabel('optical rotation (deg)');
hold off

1==1;


