clear all
clc

pkg load optim

% Initialize some other values
BBR_S = [ ...
    120.272     93.782     69.465     49.861     34.546     22.970; ...
    276.088    230.073    185.954    145.986    111.780     83.573; ...
    550.313    483.285    414.183    347.850    286.163    228.998; ...
    720.538    671.688    605.032    547.171    480.848    415.814; ...
   1006.866    947.639    886.043    822.754    739.724    667.029];

BBR_D = 1./BBR_S;

BBR_tr = [ ...
434631.32360 814933.73175 1629867.46349 3259734.92699 6519469.85398 13038939.70796; ...
 29101.54173  54565.39075  109130.78150  218261.56301  436523.12602   873046.25204; ...
  1565.26858   2934.87860    5869.75719   11739.51438   23479.02877    46958.05753; ...
   149.31155    279.95915     559.91830    1119.83659    2239.67319     4479.34638; ...
     8.00000     15.00000      30.00000      60.00000     120.00000      240.00000];

% initialise the GPL parameters
GPL_x = [0; 0; 0];  % in order D0 D1 m

% the first estimate of model parameters
% in the GPL model D(t) = D0 + D1*(t)^m
GPL_D0 = 0.0005; 
GPL_D1 = 0.0002; 
GPL_m  = 0.4;      

% D_test = fn_GPL_model (GPL_D0, GPL_D1, GPL_m, BBR_tr)
  
GPL_D  = BBR_D(1,:);
GPL_tr = BBR_tr(1,:);
for ii = 2 : rows(BBR_D)
  GPL_D  = [GPL_D, BBR_D(ii,:)];
  GPL_tr = [GPL_tr, BBR_tr(ii,:)];
endfor

fun_GPL = @(GPL_x)fn_BBR_GPL_RMSE (GPL_x, GPL_D, GPL_tr);
GPL_x0 = [GPL_D0; GPL_D1; GPL_m];
GPL_options = optimset ("fsolve");
GPL_options.MaxIter = 1000;
%GPL_options.Jacobian = 'on';
[GPL_x, GPL_RMSE, GPL_cvg, GPL_outp] =...
 fsolve(fun_GPL, GPL_x0, GPL_options) % alternative fsolve
GPL_D0 = GPL_x(1);
GPL_D1 = GPL_x(2);
GPL_m  = GPL_x(3);

disp(' ');
fprintf('GPL parameters:\n');
fprintf('D0 = %10.9f \n',GPL_D0);
fprintf('D1 = %10.9f \n',GPL_D1);
fprintf('m  = %7.6f \n',GPL_m);
fprintf('GPL_RMSE  = %7.6f \n',GPL_RMSE);
fprintf('fsolve outcome: GPL_cvg = %i:  ',GPL_cvg);
switch GPL_cvg
  case 0
    fprintf('maximum number of iterations exceeded');
    fprintf('   Iterations = %i \n',GPL_outp.iterations);
  case 1 
    fprintf('Converged to a solution point');
  case 2 
    fprintf('parameter change less than specified precision in two consecutive iterations');
  case 3 
    fprintf('improvement in objective function less than specified');
  case -1 
    fprintf('algorithm aborted by a user function');
  case -3 
    fprintf('The trust region radius became excessively small');
  otherwise
    fprintf('unknown outcome');
endswitch

disp(' ');
disp('DONE');