model.knots = [365 730 1095 1461 1826 2191 2556 2922 3287 3652 4017 4383 4748]';

frac4=model.knots./1461;
complete4= floor(frac4);
mod4=frac4-complete4;
daysleft4=model.knots - 1461*complete4;

for i=1:length(model.knots)
if (mod4(i) >= 1096/1461)    % Test for leap year
     knot_year(i) = 2000 + 4*complete4(i) + 3 + (daysleft4(i)-1095)/366;  % Leap year
else
    knot_year(i)= 2000 + 4*complete4(i) + daysleft4(i)/365;   % Non-leap year
end
end

knot_year

function decimal_year = mjd2000_to_decimal_year_v2(mjd2000)
%
% decimal_year = mjd2000_to_decimal_year_v2(mjd2000)
%
% the conversion is valid for MJD2000 from 36525 to -36159
% (i.e., decimal years between 1901.0 and 2100.0)
%
% Convert MJD2000 to decimal years.
%
    frac4 = mjd2000 ./ 1461;
    complete4 = floor(frac4);
    mod4 = frac4 - complete4;
    daysleft4 = mjd2000 - 1461*complete4;

    mask_leap_years = mod4 <= 366/1461;
    mask_regular_years = ~mask_leap_years;
    
    decimal_year = 2000 + 4*complete4;
    decimal_year(mask_leap_years) += daysleft4(mask_leap_years)/366;
    decimal_year(mask_regular_years) += 1 + (daysleft4(mask_regular_years) - 366)/365;

end % mjd2000_to_decimal_year_2

model.knots = [0 366 731 1096 1461 1827 2192 2557 2922 3288 3653 4018 4383 4749 36525 -36159];
mjd2000_to_decimal_year_v2(model.knots)


model.knots = [0 366 731 1096 1461 1827 2192 2557 2922 3288 3653 4018 4383 4749 36525 -36159];
knot_year = mjd2000_to_decimal_year(model.knots)
knot_day = decimal_year_to_mjd2000(knot_year)

N = 20000
mjd2000_min = decimal_year_to_mjd2000(1901)
mjd2000_max = decimal_year_to_mjd2000(2100)

mjd2000 = (mjd2000_max - mjd2000_min)*rand(1, N) + mjd2000_min;

decimal_year = mjd2000_to_decimal_year_v2(mjd2000);
decimal_year_ref = mjd2000_to_decimal_year(mjd2000);

# the values should be equal and the deffierence close to zero
max_diff = max(abs(decimal_year .- decimal_year_ref))