function melt(varargin)

if sum(double(size(varargin) ~= [1 4]))
    disp(' ')
    disp('COMMAND:')
    disp('    melt(data, initialparameters, ''id'', regression)')
    disp(' ')
    disp('ARGUMENTS:')
    disp('    data(:,1) = temperature              (ºC)')
    disp('    data(:,2) = absorbance')
    disp('    data(1,3) = [oligonucleotide]        (mol/dm3, total concentration of strands)')
    disp('    data(2,3) = solution path length     (cm)')
    disp('    data(3,3) = [NaCl]                   (mol/dm3)')
    disp(' ')
    disp('    initialparameters(1) = deltaH0       (J/mol)')
    disp('    initialparameters(2) = deltaS0       (J/K/mol)')
    disp('    initialparameters(3) = epsilon rate of change with temperatura, paired strands')
    disp('    initialparameters(4) = epsilon0, paired strands')
    disp('    initialparameters(5) = epsilon rate of change with temperatura, unpaired strands')
    disp('    initialparameters(6) = epsilon0, unpaired strands')
    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);
X(:,1) = X(:,1) + 273.15;
X(1:3,3) = dados(1:3,3);
X(1,4) = 0;                                                           % sinal para desenhar diagrama de distribuição

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)',@melting,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;0;0;0]);
else
    % sem regressão
    betas = beta0;
    xerrs = num2str(zeros(9,1));
end

% preparação dos resultados para inserção no gráfico dos dados
X(1,4) = 1;
[yfit Tg Ag Tmeltx] = melting(betas,X);                               % cálculo de absorvâncias e traçado do gráfico absorvância = f(temperatura)
data = num2str([real(betas)';X(1,3);X(2,3);X(3,3)]);
lege1 = ['  deltaH0';'  deltaS0';'     r-ds';'    E0-ds';'     r-ss';'    E0-ss';'  [oligo]';'        L';'   [NaCl]'];
A = [lege1 ['  ';'  ';'  ';'  ';'  ';'  ';'  ';'  ';'  '] data ...
    ['  ';'  ';'  ';'  ';'  ';'  ';'  ';'  ';'  '] xerrs ...
    ['       J/mol';'     J/K/mol';'  1/(M cm K)';'    1/(M cm)';'  1/(M cm K)';'    1/(M cm)';'           M';'          cm';'           M']];
B = '                         confidence 95%';
disp(' ');
disp(strvcat(B,A));
disp(' ');
disp(Tmeltx);
disp(' ');

% 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)) ' (A)'];
E = ['residue module average           ' num2str(mean(abs(r))) ' (A)'];
F = ['standard error of regression     ' num2str(sqrt(sum(r.*r)/(length(r)-1))) ' (A)'];      % 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(' ');

subplot(2,2,[1 3])
title([id '  ' datestr(now)],'Interpreter','none');
xlabel('Temperature (ºC)');
ylabel('absorbance');

figure(2), clf;

% gráficos de análise dos resíduos
subplot(2,2,1), normplot(r)
xlabel('residues (A)');
ylabel('probability');
subplot(2,2,2), plot(r(1:length(r)-1),r(2:length(r)),'-o')
title('Sequencial Self-correlation');
xlabel('residues (A)');
ylabel('residues (A)');
subplot(2,2,3), plot(X(:,1)-273.15,r,'-o')
title('Residue Sequence');
xlabel('temperature (ºC)');
ylabel('residues (A)');

subplot(2,2,4)
hist(r);
title('Histogram');
ylabel('frequency');
xlabel('residues (A)');
1==1;




function [Ac Tg Ag texto] = melting(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
Ax = X(:,2);     % absorvâncias experimentais
T = X(:,1);      % temperaturas experimentais (Kelvin)
CT = X(1,3);     % concentração do oligómero (molar)
L = X(2,3);      % espessura da cuvete em centímetros            
R = 8.314;
dd = X(1,4);

% parâmetros ajustáveis
H0 = par(1);
S0 = par(2);
pemp1 = par(3);
pemp2 = par(4);
pdesemp1 = par(5);
pdesemp2 = par(6);


% cálculo de absorvâncias estimadas (Ac) para os valores de parâmetros
% correntes (H0, S0, pemp1, pemp2, pdesemp1, pdesemp2 e C) e para
% as temperaturas experimentais.
Ac = [];
t = T;
for t = T'
    
    K = exp(-(H0-t*S0)/(R*t));
    
    a = 2*K;
    b = 1;
    c = -CT;
    C = (-b+sqrt(b^2-4*a*c))/(2*a);
    
    a = 4;
    b = -(4*CT+(1/K));
    c = CT^2;
    C2 = (-b-sqrt(b^2-4*a*c))/(2*a);
    
    Eemp = pemp1*t + pemp2;
    Edesemp = pdesemp1*t + pdesemp2;
    Ac = [Ac; (2*Eemp*C2 + Edesemp*C)*L];

end
    Ac = Ac';

% produção de linha de ajuste ligando 50 pontos
passo = (max(T) - min(T))/50;
Tg = min(T):passo:max(T);
Ag = []; Cg = []; C2g = [];
for t = Tg
    
    K = exp(-(H0-t*S0)/(R*t));
    
    a = 2*K;
    b = 1;
    c = -CT;
    C = (-b+sqrt(b^2-4*a*c))/(2*a);
    if dd; Cg = [Cg; C]; end;
    
    a = 4;
    b = -(4*CT+(1/K));
    c = CT^2;
    C2 = (-b-sqrt(b^2-4*a*c))/(2*a);
    if dd; C2g = [C2g; C2]; end;

    Eemp = pemp1*t + pemp2;
    Edesemp = pdesemp1*t + pdesemp2;
    Ag = [Ag; (2*Eemp*C2 + Edesemp*C)*L];
    
end

% gráficos
figure(1); clf;

% gráfico dos dados experimentais e da linha de ajuste
subplot(2,2,[1 3])
plot(T-273.15,Ax,'or','MarkerFaceColor','r','MarkerSize',6);
hold on
plot(Tg-273.15,Ag,'-k','LineWidth',2);
xlabel('temperature (ºC)');
ylabel('A (absorbance)');

% gráfico primeira derivada dA/dT
subplot(2,2,4)
Tmelt = [];                  % temperatura no ponto de inflexão de A = f(T)
passo = Tg(2) - Tg(1);
xprimderiv = Tg(1:length(Tg)-1) + passo/2;
xsegnderiv = xprimderiv(1:length(xprimderiv)-1) + passo/2;
yprimderiv = diff(Ag)/passo;                   % primeira derivada
ysegnderiv = diff(Ag,2)/passo^2;               %segunda derivada
for iyprim = 3:length(yprimderiv)
    if yprimderiv(iyprim-2) < yprimderiv(iyprim-1) && ...
            yprimderiv(iyprim-1) > yprimderiv(iyprim);
        Tmelt = xprimderiv(iyprim-1);
    end
end
plot(xsegnderiv-273.15,ysegnderiv);
if ~isempty(Tmelt)
    texto = ['Melting temperature: '  num2str(Tmelt-273.15) ' ºC'];
    % text(mean(xprimderiv)-273.15,0,texto,'FontName','Courier New','FontSize',10,'HorizontalAlignment','left');
end
xlabel('temperature (ºC)');
ylabel('ddA / ddT');
title('Second Derivative');

% gráfico do diagrama de distribuição
if dd
    subplot(2,2,2)
    plot(Tg-273.15,Cg,Tg-273.15,C2g);
    hdd = legend('single strands','double strands');
    set(hdd,'Location','NorthWest');
    xlabel('temperature (ºC)');
    ylabel('concentration (M)');
    title('Distribution Diagram');
    1==1;
end

hold off

1==1;


