ncfile  = '../data/volumevar.nc';
txtfile = '../python/nourishments.txt';

fid = fopen(txtfile);
while ~feof(fid)
    line = fgetl(fid);
    transect = regexp(line, '^\d+', 'match');
    years = regexp(line, '(?<= )\d{4}', 'match');
    years = cellfun(@str2double, years);
    fprintf('%s\n', transect{1})
    nour.(['id' transect{1}]) = years;
end
fclose(fid);

id = nc_varget(ncfile, 'id');
volC = nc_varget(ncfile, 'volumeC');
%mvolC = nc_varget(ncfile, 'mean_volumeC');
dv = datevec(nc_cf_time(ncfile));

idxs = ismember(id(:)', [5004898 6002011 7003975 8008950 9011450]);
fid = fopen('distributions.txt', 'w');

fprintf(fid, '%10s\t%s\t%s\t%s\t%s\t%s\n', 'transect', 'histvar', 'a b', 'dV', 'a b', 'factor');
for idx = find(idxs)
%    idx = id==idxs(i);
    nouryears = nour.(sprintf('id%i',id(idx)));
    nouridx = ismember(dv(:,1), nouryears);
    histvar = volC(:,idx) - nanmean(volC(:,idx));
    maxhistvar = nanmax(abs(histvar));
    [a, b, c] = DetermineParametersVolumeVariation(...
        'TargetVolume', maxhistvar);
    
    annvar = diff(volC(:,idx) - nanmean(volC(:,idx)));
    if id(idx) == 6002011
        annvar(annvar>200) = nan;
    else
        annvar(annvar>1000) = nan;
    end
    maxannvar = nanmax(abs(annvar));
    [a2, b2, c2] = DetermineParametersVolumeVariation(...
        'TargetVolume', maxannvar);
    figure
    plot(dv(:,1), histvar,...
        dv(1:end-1,1), annvar, '-x',...
        dv(nouridx)-1, annvar(nouridx(2:end)), 'o')
    title(sprintf('%i', id(idx)))
    fprintf(fid, '%10i\t%.0f\t%.0f\t%.0f\t%.0f\t%.1f\n', id(idx), maxhistvar, b, maxannvar, b2, maxhistvar/maxannvar);
    
    if idx==find(idxs, 1, 'last')
        
        tidx = dv(:,1)<2008;
        histvar = volC(tidx,idx) - nanmean(volC(:,idx));
        maxhistvar = nanmax(abs(histvar));
        [a, b, c] = DetermineParametersVolumeVariation(...
            'TargetVolume', maxhistvar);
        
        annvar = diff(volC(tidx,idx) - nanmean(volC(:,idx)));
        maxannvar = nanmax(abs(annvar));
        [a2, b2, c2] = DetermineParametersVolumeVariation(...
            'TargetVolume', maxannvar);
        
        fprintf(fid, '%10i\t%.0f\t%.0f\t%.0f\t%.0f\t%.1f\n', id(idx), maxhistvar, b, maxannvar, b2, maxhistvar/maxannvar);
    end
end
fclose(fid);