base_functions.tar

./base_functions/0000700442630740010010000000000014243756265014703 5ustar  adelmanDomain Users./base_functions/bandpass_err_fcn.m0000700442630740010010000004152513722027514020351 0ustar  adelmanDomain Usersfunction [bandpassed_array,trend_array,nan_ind] = bandpass_err_fcn(input_array,n_dim,delta_dim,low_bound,high_bound,steepness_factor,trend_handling_opt,nan_handling_opt,edge_handling_opt,uneven_edge_handling_opt)

% delta_dim: spacing of grid in dimension of bandpassing
% low_bound: low frequency-wavenumber bound
% high_bound: high frequency-wavenumber bound
% n_dim: number of dimension to bandpass along
% steepness_factor: how rapidly the bandpass filter cuts off at the
% boundaries (must be positive; higher values have a steeper cutoff, and
% possibly more spectral ringing)
% trend_handling_opt: 0 = do not add trend back to bandpassed array, 1 = add trend
% back to bandpassed array
% nan_handling_opt: 0 = do not put NaNs in array, 1 = include NaNs in array
% (trend array and NaN indices are included as outputs regardless)
% edge_handling_opt: 0 = do not replace edge values with NaNs, 1 = replace
% edge values with NaNs
% uneven_edge_handling_opt: 0 = treat edges as even, with zeros for NaNs;
% 1 = adjust padding positions for uneven edges

dim_length = size(input_array,n_dim);

f_vec = (1/(delta_dim*(2*dim_length)))*(mod((0:1:((2*dim_length) - 1))' + dim_length,(2*dim_length)) - dim_length);

% n_avg_ind = round(((1/8)*exp(-mean(log([low_bound high_bound]))))/delta_dim);  
n_avg_ind = round(((1/8)*exp(-mean(log([max([low_bound (1/(delta_dim*dim_length))]) min([high_bound (1/(2*delta_dim))])]))))/delta_dim);  

% bandpass_filter = 2*0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));
bandpass_filter = 0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));

% G = [ones(dim_length,1) (1:1:dim_length)'];
% reg_operator_trend = G*(((G')*G)^(-1))*(G');
padding_ind_inside = floor(dim_length/2) + (1:1:dim_length)';
if abs(uneven_edge_handling_opt) < 1e-5
    padding_ind_begin = (1:1:floor(dim_length/2))';
    padding_ind_end = ((floor(dim_length/2) + dim_length + 1):1:(2*dim_length))';
end

size_array = size(input_array);
prod_size_not_inc_bp_dim = prod(size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));


input_array_permuted = reshape(permute(input_array,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_bp_dim]);
nan_ind_permuted = union(union(find(isnan(input_array_permuted) == 1),find(isinf(input_array_permuted) == 1)),find(abs(input_array_permuted) < (1e-15)*max(abs(input_array_permuted(isinf(input_array_permuted) == 0)))));
clear input_array
input_array_permuted(nan_ind_permuted) = 0;
nan_mask_input_permuted = ones(size(input_array_permuted));
nan_mask_input_permuted(nan_ind_permuted) = 0;
% trend_array_permuted = reshape(reg_operator_trend*(reshape(input_array_permuted,[dim_length prod_size_not_inc_bp_dim])),size_array([n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));

% remove means and trends
index_array_nomean = repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]) - repmat((sum(nan_mask_input_permuted.*repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]),1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
input_array_permuted_nomean = input_array_permuted - repmat((sum(nan_mask_input_permuted.*input_array_permuted,1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
trend_array_permuted = (index_array_nomean.*repmat((1./sum((nan_mask_input_permuted.*index_array_nomean).^2,1)).*sum(nan_mask_input_permuted.*index_array_nomean.*input_array_permuted_nomean,1),[size(input_array_permuted,1) 1]));
trend_array_permuted = trend_array_permuted + (input_array_permuted - input_array_permuted_nomean);

input_padded = zeros([(2*dim_length) prod_size_not_inc_bp_dim]);
input_array_permuted_minus_trend = input_array_permuted - trend_array_permuted;
clear input_array_permuted

input_array_permuted_minus_trend(nan_ind_permuted) = 0;


% input_array_variance = sum(nan_mask_input_permuted.*(input_array_permuted_minus_trend.^2),1)./sum(nan_mask_input_permuted,1);


if abs(edge_handling_opt - 1) < 1e-5
    test_lags = [1; 2; 3; round(exp((1.2:0.3:log(dim_length))'))];
    N = length(test_lags);
    autocorr_at_lags = NaN([N prod_size_not_inc_bp_dim]);
    integral_timescale = NaN([1 prod_size_not_inc_bp_dim]);
    n_lag_step = 1;
    complete_flag = 0;
    while ((complete_flag == 0) && (n_lag_step <= N))
        lag_ind = test_lags(n_lag_step);
        autocorr_at_lags(n_lag_step,:) = sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*input_array_permuted_minus_trend(1:(dim_length - lag_ind),:).*input_array_permuted_minus_trend((1 + lag_ind):dim_length,:),1)./((sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*(input_array_permuted_minus_trend(1:(dim_length - lag_ind),:).^2),1).^(1/2)).*(sum((nan_mask_input_permuted(1:(dim_length - lag_ind),:) & nan_mask_input_permuted((1 + lag_ind):dim_length,:)).*(input_array_permuted_minus_trend((1 + lag_ind):dim_length,:).^2),1).^(1/2)));

        if n_lag_step == N
            first_neg_ind = find(isnan(integral_timescale) == 1);                
        else
            first_neg_ind = find((autocorr_at_lags(n_lag_step,:) < 0) & (isnan(integral_timescale) == 1));                
        end
        if isempty(first_neg_ind) == 0
            first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_bp_dim]);
            first_neg_ind_mask_2D(first_neg_ind) = 1;
            first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
        %     first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
            curr_integral_timescale = delta_dim*(first_neg_ind_mask(1,:) + (2*sum(repmat(0.5*(test_lags([(2:1:n_lag_step)'; n_lag_step]) - [0; test_lags(1:(n_lag_step - 1))]),[1 prod_size_not_inc_bp_dim]).*first_neg_ind_mask.*autocorr_at_lags(1:n_lag_step,:),1)));
            integral_timescale(first_neg_ind) = curr_integral_timescale(first_neg_ind);
        end
        if sum(sum(isnan(integral_timescale))) == 0
            complete_flag = 1;
        end
        
        n_lag_step = n_lag_step + 1;
    end
    % input_array_permuted_lag1_autocorr = sum((nan_mask_input_permuted(1:(dim_length - 1),:).*input_array_permuted_minus_trend(1:(dim_length - 1),:)).*(nan_mask_input_permuted(2:dim_length,:).*input_array_permuted_minus_trend(2:dim_length,:)),1)./(((sum(nan_mask_input_permuted(1:(dim_length - 1),:).*(input_array_permuted_minus_trend(1:(dim_length - 1),:).^2),1)).^(1/2)).*((sum(nan_mask_input_permuted(2:dim_length,:).*(input_array_permuted_minus_trend(2:dim_length,:).^2),1)).^(1/2)));
    % freq_dominant_approx = 1./(4*delta_dim*((1/(1 - input_array_permuted_lag1_autocorr)) - 1));
    freq_dominant_approx = 1./(pi*integral_timescale);
end

input_padded(padding_ind_inside,:) = input_array_permuted_minus_trend;
clear autocorr_at_lags input_array_permuted_minus_trend

% pad with erf slopes on either end of the range
if abs(uneven_edge_handling_opt - 1) < 1e-5
    first_good_ind = ones([1 prod_size_not_inc_bp_dim]);
    last_good_ind = dim_length*ones([1 prod_size_not_inc_bp_dim]);
    cum_forward_nan_mask_input_permuted = cumsum(nan_mask_input_permuted,1);
    cum_reverse_nan_mask_input_permuted = flipdim(cumsum(flipdim(nan_mask_input_permuted,1)),1);
    first_good_ind(abs(cum_forward_nan_mask_input_permuted(dim_length,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - 1) < 1e-5) & (abs(cum_forward_nan_mask_input_permuted - 1) < 1e-5)),find((abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(1:(dim_length - 1),:)]) < 1e-5) & (abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(2:dim_length,:)] - 1) < 1e-5))) - 1,dim_length) + 1)';
    last_good_ind(abs(cum_reverse_nan_mask_input_permuted(1,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - dim_length) < 1e-5) & (abs(cum_reverse_nan_mask_input_permuted - 1) < 1e-5)),find((abs([cum_reverse_nan_mask_input_permuted(2:dim_length,:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))]) < 1e-5) & (abs([cum_reverse_nan_mask_input_permuted(1:(dim_length - 1),:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))] - 1) < 1e-5))) - 1,dim_length) + 1)';
end
if n_avg_ind > 0
    if abs(uneven_edge_handling_opt) < 1e-5
        pad_slope_begin = (mean(input_padded(floor(dim_length/2) + ((n_avg_ind + 1):(2*n_avg_ind)),:),1) - mean(input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1))/n_avg_ind;
        pad_val_begin = mean(input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1) - (((n_avg_ind - 1)/2)*pad_slope_begin);
        pad_slope_end = (mean(input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) - mean(input_padded(floor(dim_length/2) + ((dim_length - (2*n_avg_ind) + 1):(dim_length - n_avg_ind)),:),1))/n_avg_ind;
        pad_val_end = mean(input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) + (((n_avg_ind - 1)/2)*pad_slope_end);
    else
        along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
        avg_range_mask_1 = zeros(size(input_padded));
        avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= n_avg_ind)) = 1;
        avg_range_mask_2 = zeros(size(input_padded));
        avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= n_avg_ind + 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= (2*n_avg_ind))) = 1;
        pad_slope_begin = ((sum(avg_range_mask_2.*input_padded,1)./sum(avg_range_mask_2,1)) - (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
        pad_val_begin = (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)) - (((n_avg_ind - 1)/2)*pad_slope_begin);
        avg_range_mask_1 = zeros(size(input_padded));
        avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -1) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -n_avg_ind)) = 1;
        avg_range_mask_2 = zeros(size(input_padded));
        avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -(n_avg_ind + 1)) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -(2*n_avg_ind))) = 1;
        pad_slope_end = (-(sum(avg_range_mask_2.*input_padded,1)./sum(avg_range_mask_2,1)) + (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
        pad_val_end = (sum(avg_range_mask_1.*input_padded,1)./sum(avg_range_mask_1,1)) + (((n_avg_ind - 1)/2)*pad_slope_end);
        clear avg_range_mask_*
    end
else
    if abs(uneven_edge_handling_opt) < 1e-5
        pad_slope_begin = input_padded(floor(dim_length/2) + 1,:)/0.5;
        pad_val_begin = input_padded(floor(dim_length/2) + 1,:);
        pad_slope_end = -input_padded(floor(dim_length/2) + dim_length,:)/0.5;
        pad_val_end = input_padded(floor(dim_length/2) + dim_length,:);
    else
        along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
        avg_range_mask = zeros(size(input_padded));
        avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
        pad_slope_begin = (sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1))/0.5;
        pad_val_begin = sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1);
        avg_range_mask = zeros(size(input_padded));
        avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
        pad_slope_end = (sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1))/0.5;
        pad_val_end = sum(avg_range_mask.*input_padded,1)./sum(avg_range_mask,1);
        clear avg_range_mask
    end
end
if abs(uneven_edge_handling_opt) < 1e-5
    input_padded(padding_ind_begin,:) = repmat(pad_val_begin,[length(padding_ind_begin) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[length(padding_ind_begin) 1]).*(repmat(reshape(padding_ind_begin,[length(padding_ind_begin) 1]),[1 prod_size_not_inc_bp_dim]) - min(padding_ind_inside))));
    input_padded(padding_ind_end,:) = repmat(pad_val_end,[length(padding_ind_end) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[length(padding_ind_end) 1]).*(repmat(reshape(padding_ind_end,[length(padding_ind_end) 1]),[1 prod_size_not_inc_bp_dim]) - max(padding_ind_inside))));
else
    begin_padding_array = repmat(pad_val_begin,[(2*dim_length) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[(2*dim_length) 1]).*(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]))));
    input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0) = begin_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0);
    end_padding_array = repmat(pad_val_end,[(2*dim_length) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[(2*dim_length) 1]).*(repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]) - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]))));
    input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0) = end_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0);
    clear along_dim_ind_array
end
input_padded((isnan(input_padded) == 1) | (isinf(input_padded) == 1)) = 0;


% padded_array_variance = sum(input_padded.^2,1)/size(input_padded,1);
bandpassed_array_padded = ifft(repmat(bandpass_filter,[1 prod_size_not_inc_bp_dim]).*fft(input_padded,[],1),[],1);

mean_input_padded = mean(input_padded,1);
clear input_padded

% % make variance adjustment
% bandpassed_array_padded = repmat((input_array_variance./padded_array_variance).^(1/2),[(2*dim_length) 1]).*bandpassed_array_padded;


% add trend back to array, depending on option given
if abs(trend_handling_opt - 0) < 1e-5
    bandpassed_array_not_padded = bandpassed_array_padded(padding_ind_inside,:);
elseif abs(trend_handling_opt - 1) < 1e-5
    bandpassed_array_not_padded = bandpassed_array_padded(padding_ind_inside,:) + reshape(trend_array_permuted,[dim_length prod_size_not_inc_bp_dim]) + repmat(mean_input_padded,[dim_length 1]);
end
clear bandpassed_array_padded

bandpassed_array_not_padded(nan_ind_permuted) = NaN;


% % pad edges with NaNs where there is low precision about bandpassed values

% n_edge_err = floor(((1/high_bound)/3)/delta_dim);
if abs(edge_handling_opt - 1) < 1e-5
    denom_edge_scale = (delta_dim.*high_bound).*((1 - ((low_bound./high_bound).^0.5)).^1.25);
    n_edge_err = round(max([min([((((0.9*(low_bound./((pi/2)*freq_dominant_approx))) - (3/40))./(delta_dim.*((((pi/2)*freq_dominant_approx).*low_bound).^(1/2)))) + ((0.45 - (1.8*(low_bound./((pi/2)*freq_dominant_approx))))./denom_edge_scale)); (0.15./(delta_dim.*((((pi/2)*freq_dominant_approx).*low_bound).^(1/2))))],[],1); (0.3./repmat(denom_edge_scale,[1 size(freq_dominant_approx,2)]))],[],1));
    
%     bandpassed_array_not_padded(1:n_edge_err,:) = NaN;
%     bandpassed_array_not_padded((dim_length - n_edge_err + 1):dim_length,:) = NaN;
    if abs(uneven_edge_handling_opt) < 1e-5
        bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(1 + n_edge_err,[dim_length 1]) < 0) = NaN;
        bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(dim_length - n_edge_err,[dim_length 1]) > 0) = NaN;
    else
        bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(first_good_ind + n_edge_err,[dim_length 1]) < 0) = NaN;
        bandpassed_array_not_padded(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - repmat(last_good_ind - n_edge_err,[dim_length 1]) > 0) = NaN;
    end
end

bandpassed_array = real(permute(reshape(bandpassed_array_not_padded,[dim_length size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]));
trend_array = permute(reshape(trend_array_permuted,[dim_length size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
clear bandpassed_array_not_padded trend_array_permuted

nan_ind = find(isnan(bandpassed_array) == 1);

% remove NaNs from array, depending on option given
if abs(nan_handling_opt - 0) < 1e-5
    bandpassed_array(nan_ind) = 0;
end./base_functions/bandpass_err_fcn_gain_coeffs.m0000700442630740010010000002765213722056464022707 0ustar  adelmanDomain Usersfunction [bandpass_gain_coeff_array] = bandpass_err_fcn_gain_coeffs(input_array,n_dim,delta_dim,low_bound,high_bound,steepness_factor,mean_handling_opt,uneven_edge_handling_opt)

% delta_dim: spacing of grid in dimension of bandpassing
% low_bound: low frequency-wavenumber bound
% high_bound: high frequency-wavenumber bound
% n_dim: number of dimension to bandpass along
% steepness_factor: how rapidly the bandpass filter cuts off at the
% boundaries (must be positive; higher values have a steeper cutoff, and
% possibly more spectral ringing)
% mean_handling_opt: 0 = gain coefficients do not include mean; 1 = gain coeffients include mean (trend is not removed in this function regardless)
% uneven_edge_handling_opt: 0 = treat edges as even, with zeros for NaNs;
% 1 = adjust padding positions for uneven edges

dim_length = size(input_array,n_dim);

f_vec = (1/(delta_dim*(2*dim_length)))*(mod((0:1:((2*dim_length) - 1))' + dim_length,(2*dim_length)) - dim_length);

% n_avg_ind = round(((1/8)*exp(-mean(log([low_bound high_bound]))))/delta_dim);  
n_avg_ind = round(((1/8)*exp(-mean(log([max([low_bound (1/(delta_dim*dim_length))]) min([high_bound (1/(2*delta_dim))])]))))/delta_dim);  

% bandpass_filter = 2*0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));
bandpass_filter = 0.5*(erf(steepness_factor*(log(abs(f_vec)) - log(low_bound))) - erf(steepness_factor*(log(abs(f_vec)) - log(high_bound))));
if abs(mean_handling_opt - 1) < 1e-5
    bandpass_filter(1) = 1;
end

% G = [ones(dim_length,1) (1:1:dim_length)'];
% reg_operator_trend = G*(((G')*G)^(-1))*(G');
padding_ind_inside = floor(dim_length/2) + (1:1:dim_length)';
if abs(uneven_edge_handling_opt) < 1e-5
    padding_ind_begin = (1:1:floor(dim_length/2))';
    padding_ind_end = ((floor(dim_length/2) + dim_length + 1):1:(2*dim_length))';
end

size_array = size(input_array);
prod_size_not_inc_bp_dim = prod(size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));


input_array_permuted = reshape(permute(input_array,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_bp_dim]);
nan_ind_permuted = union(union(find(isnan(input_array_permuted) == 1),find(isinf(input_array_permuted) == 1)),find(abs(input_array_permuted) < (1e-15)*max(abs(input_array_permuted(isinf(input_array_permuted) == 0)))));
clear input_array
input_array_permuted(nan_ind_permuted) = 0;
nan_mask_input_permuted = ones(size(input_array_permuted));
nan_mask_input_permuted(nan_ind_permuted) = 0;
% % trend_array_permuted = reshape(reg_operator_trend*(reshape(input_array_permuted,[dim_length prod_size_not_inc_bp_dim])),size_array([n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));
% 
% % remove means and trends
% index_array_nomean = repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]) - repmat((sum(nan_mask_input_permuted.*repmat((1:1:size(input_array_permuted,1))',[1 prod_size_not_inc_bp_dim]),1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
% input_array_permuted_nomean = input_array_permuted - repmat((sum(nan_mask_input_permuted.*input_array_permuted,1)./sum(nan_mask_input_permuted,1)),[size(input_array_permuted,1) 1]);
% trend_array_permuted = (index_array_nomean.*repmat((1./sum((nan_mask_input_permuted.*index_array_nomean).^2,1)).*sum(nan_mask_input_permuted.*index_array_nomean.*input_array_permuted_nomean,1),[size(input_array_permuted,1) 1]));
% trend_array_permuted = trend_array_permuted + (input_array_permuted - input_array_permuted_nomean);
% 
% input_array_permuted_minus_trend = input_array_permuted - trend_array_permuted;
% clear input_array_permuted
% 
% input_array_permuted_minus_trend(nan_ind_permuted) = 0;
% 
% 
% % input_array_variance = sum(nan_mask_input_permuted.*(input_array_permuted_minus_trend.^2),1)./sum(nan_mask_input_permuted,1);


bandpass_gain_coeff_permuted = NaN([dim_length prod_size_not_inc_bp_dim dim_length]);
for curr_ind = 1:dim_length
    curr_input_padded = zeros([(2*dim_length) prod_size_not_inc_bp_dim]);
    curr_input_padded(padding_ind_inside(1) - 1 + curr_ind,:) = 1;

    % pad with erf slopes on either end of the range
    if abs(uneven_edge_handling_opt - 1) < 1e-5
        first_good_ind = ones([1 prod_size_not_inc_bp_dim]);
        last_good_ind = dim_length*ones([1 prod_size_not_inc_bp_dim]);
        cum_forward_nan_mask_input_permuted = cumsum(nan_mask_input_permuted,1);
        cum_reverse_nan_mask_input_permuted = flipdim(cumsum(flipdim(nan_mask_input_permuted,1)),1);
        first_good_ind(abs(cum_forward_nan_mask_input_permuted(dim_length,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - 1) < 1e-5) & (abs(cum_forward_nan_mask_input_permuted - 1) < 1e-5)),find((abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(1:(dim_length - 1),:)]) < 1e-5) & (abs([((-1)*ones([1 prod_size_not_inc_bp_dim])); cum_forward_nan_mask_input_permuted(2:dim_length,:)] - 1) < 1e-5))) - 1,dim_length) + 1)';
        last_good_ind(abs(cum_reverse_nan_mask_input_permuted(1,:)) > 1e-5) = (mod(union(find((abs(repmat((1:1:dim_length)',[1 prod_size_not_inc_bp_dim]) - dim_length) < 1e-5) & (abs(cum_reverse_nan_mask_input_permuted - 1) < 1e-5)),find((abs([cum_reverse_nan_mask_input_permuted(2:dim_length,:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))]) < 1e-5) & (abs([cum_reverse_nan_mask_input_permuted(1:(dim_length - 1),:); ((-1)*ones([1 prod_size_not_inc_bp_dim]))] - 1) < 1e-5))) - 1,dim_length) + 1)';
    end
    if n_avg_ind > 0
        if abs(uneven_edge_handling_opt) < 1e-5
            pad_slope_begin = (mean(curr_input_padded(floor(dim_length/2) + ((n_avg_ind + 1):(2*n_avg_ind)),:),1) - mean(curr_input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1))/n_avg_ind;
            pad_val_begin = mean(curr_input_padded(floor(dim_length/2) + (1:n_avg_ind),:),1) - (((n_avg_ind - 1)/2)*pad_slope_begin);
            pad_slope_end = (mean(curr_input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) - mean(curr_input_padded(floor(dim_length/2) + ((dim_length - (2*n_avg_ind) + 1):(dim_length - n_avg_ind)),:),1))/n_avg_ind;
            pad_val_end = mean(curr_input_padded(floor(dim_length/2) + ((dim_length - n_avg_ind + 1):dim_length),:),1) + (((n_avg_ind - 1)/2)*pad_slope_end);
        else
            along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
            avg_range_mask_1 = zeros(size(curr_input_padded));
            avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= n_avg_ind)) = 1;
            avg_range_mask_2 = zeros(size(curr_input_padded));
            avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 >= n_avg_ind + 1) & (along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) + 1 <= (2*n_avg_ind))) = 1;
            pad_slope_begin = ((sum(avg_range_mask_2.*curr_input_padded,1)./sum(avg_range_mask_2,1)) - (sum(avg_range_mask_1.*curr_input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
            pad_val_begin = (sum(avg_range_mask_1.*curr_input_padded,1)./sum(avg_range_mask_1,1)) - (((n_avg_ind - 1)/2)*pad_slope_begin);
            avg_range_mask_1 = zeros(size(curr_input_padded));
            avg_range_mask_1((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -1) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -n_avg_ind)) = 1;
            avg_range_mask_2 = zeros(size(curr_input_padded));
            avg_range_mask_2((along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 <= -(n_avg_ind + 1)) & (along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) - 1 >= -(2*n_avg_ind))) = 1;
            pad_slope_end = (-(sum(avg_range_mask_2.*curr_input_padded,1)./sum(avg_range_mask_2,1)) + (sum(avg_range_mask_1.*curr_input_padded,1)./sum(avg_range_mask_1,1)))/n_avg_ind;
            pad_val_end = (sum(avg_range_mask_1.*curr_input_padded,1)./sum(avg_range_mask_1,1)) + (((n_avg_ind - 1)/2)*pad_slope_end);
            clear avg_range_mask_*
        end
    else
        if abs(uneven_edge_handling_opt) < 1e-5
            pad_slope_begin = curr_input_padded(floor(dim_length/2) + 1,:)/0.5;
            pad_val_begin = curr_input_padded(floor(dim_length/2) + 1,:);
            pad_slope_end = -curr_input_padded(floor(dim_length/2) + dim_length,:)/0.5;
            pad_val_end = curr_input_padded(floor(dim_length/2) + dim_length,:);
        else
            along_dim_ind_array = repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]);
            avg_range_mask = zeros(size(curr_input_padded));
            avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
            pad_slope_begin = (sum(avg_range_mask.*curr_input_padded,1)./sum(avg_range_mask,1))/0.5;
            pad_val_begin = sum(avg_range_mask.*curr_input_padded,1)./sum(avg_range_mask,1);
            avg_range_mask = zeros(size(curr_input_padded));
            avg_range_mask(abs(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1])) < 1e-5) = 1;
            pad_slope_end = (sum(avg_range_mask.*curr_input_padded,1)./sum(avg_range_mask,1))/0.5;
            pad_val_end = sum(avg_range_mask.*curr_input_padded,1)./sum(avg_range_mask,1);
            clear avg_range_mask
        end
    end
    if abs(uneven_edge_handling_opt) < 1e-5
        curr_input_padded(padding_ind_begin,:) = repmat(pad_val_begin,[length(padding_ind_begin) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[length(padding_ind_begin) 1]).*(repmat(reshape(padding_ind_begin,[length(padding_ind_begin) 1]),[1 prod_size_not_inc_bp_dim]) - min(padding_ind_inside))));
        curr_input_padded(padding_ind_end,:) = repmat(pad_val_end,[length(padding_ind_end) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[length(padding_ind_end) 1]).*(repmat(reshape(padding_ind_end,[length(padding_ind_end) 1]),[1 prod_size_not_inc_bp_dim]) - max(padding_ind_inside))));
    else
        begin_padding_array = repmat(pad_val_begin,[(2*dim_length) 1]).*(1 + erf(((pi^(1/2))/2)*repmat(abs(pad_slope_begin./pad_val_begin),[(2*dim_length) 1]).*(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]))));
        curr_input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0) = begin_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + first_good_ind,[(2*dim_length) 1]) < 0);
        end_padding_array = repmat(pad_val_end,[(2*dim_length) 1]).*(1 - erf(((pi^(1/2))/2)*repmat(abs(pad_slope_end./pad_val_end),[(2*dim_length) 1]).*(repmat((1:1:(2*dim_length))',[1 prod_size_not_inc_bp_dim]) - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]))));
        curr_input_padded(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0) = end_padding_array(along_dim_ind_array - repmat(floor(dim_length/2) + last_good_ind,[(2*dim_length) 1]) > 0);
        clear along_dim_ind_array
    end
    curr_input_padded((isnan(curr_input_padded) == 1) | (isinf(curr_input_padded) == 1)) = 0;


    curr_bandpass_gain_coeff_padded = ifft(repmat(bandpass_filter,[1 prod_size_not_inc_bp_dim]).*fft(curr_input_padded,[],1),[],1);
    bandpass_gain_coeff_permuted(:,:,curr_ind) = curr_bandpass_gain_coeff_padded(padding_ind_inside,:);
    
end

bandpass_gain_coeff_array = permute(reshape(bandpass_gain_coeff_permuted,[dim_length size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]) dim_length]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array)) (length(size_array) + 1)]);
clear bandpass_gain_coeff_permuted
./base_functions/bcyr.m0000700442630740010010000000605512741513710016014 0ustar  adelmanDomain Usersfunction [output_colorbar] = bcyr(N)

% generate "blue-cyan-yellow-red" colormap (more
% colorblind-friendly)



% generate initial "jet" colormap with N colors

% if mod(N,2) == 0
%     starting_colorbar = jet(N);
% else
%     starting_colorbar = jet(N + 1);
%     starting_colorbar = starting_colorbar(2:size(starting_colorbar,1),:);
% end

if ismember(mod(N,8),[4 5 7]) == 1
    starting_colorbar = jet(N);
elseif ismember(mod(N,8),[0 6]) == 1
    starting_colorbar = jet(N + 3);
    starting_colorbar = starting_colorbar(3:(size(starting_colorbar,1) - 1),:);
elseif mod(N,8) == 1
    starting_colorbar = jet(N + 4);
    starting_colorbar = starting_colorbar(3:(size(starting_colorbar,1) - 2),:);
elseif mod(N,8) == 2
    starting_colorbar = jet(N + 1);
    starting_colorbar = starting_colorbar(2:size(starting_colorbar,1),:);
elseif mod(N,8) == 3
    starting_colorbar = jet(N + 2);
    starting_colorbar = starting_colorbar(2:(size(starting_colorbar,1) - 1),:);
end
red = double(starting_colorbar(:,1));
green = double(starting_colorbar(:,2));
blue = double(starting_colorbar(:,3));


% scale colormap to convert from "jet" to "blue-white-red" colormap

contain_blue_ind = find(blue > 0);
red(contain_blue_ind) = min([(2*red(contain_blue_ind)) (1*ones(length(contain_blue_ind),1))],[],2);
contain_red_ind = find(red > 0);
blue(contain_red_ind) = min([(2*blue(contain_red_ind)) (1*ones(length(contain_red_ind),1))],[],2);
% blue_end_ind = find((red < 1) & (blue > 10));
blue_end_ind = find((blue - red) > 0);
green_weighting_over_red = 10.0;
sum_red_green = red(blue_end_ind) + green(blue_end_ind);
brighter_ind = find(sum_red_green > 1);
red(blue_end_ind(brighter_ind)) = 1 - ((green_weighting_over_red/(green_weighting_over_red + 1))*(2 - sum_red_green(brighter_ind)));
dimmer_ind = find(sum_red_green <= 1);
red(blue_end_ind(dimmer_ind)) = (1/(green_weighting_over_red + 1))*sum_red_green(dimmer_ind);
green(blue_end_ind) = sum_red_green - red(blue_end_ind);
% red(blue_end_ind) = min([(0.7*mean([red(blue_end_ind) green(blue_end_ind)],2)) (1*ones(length(blue_end_ind),1))],[],2);
% green(blue_end_ind) = min([(1.3*mean([red(blue_end_ind) green(blue_end_ind)],2)) (1*ones(length(blue_end_ind),1))],[],2);
% red_end_ind = find((blue < 1) & (red > 10));
red_end_ind = find((red - blue) > 0);
green_weighting_over_blue = 10.0;
sum_blue_green = blue(red_end_ind) + green(red_end_ind);
brighter_ind = find(sum_blue_green > 1);
blue(red_end_ind(brighter_ind)) = 1 - ((green_weighting_over_blue/(green_weighting_over_blue + 1))*(2 - sum_blue_green(brighter_ind)));
dimmer_ind = find(sum_blue_green <= 1);
blue(red_end_ind(dimmer_ind)) = (1/(green_weighting_over_blue + 1))*sum_blue_green(dimmer_ind);
green(red_end_ind) = sum_blue_green - blue(red_end_ind);
% green(red_end_ind) = mean([green(red_end_ind) blue(red_end_ind)],2);
% blue(red_end_ind) = mean([green(red_end_ind) blue(red_end_ind)],2);


% create new color array

output_colorbar = zeros(N,3);
output_colorbar(:,1) = red;
output_colorbar(:,2) = green;
output_colorbar(:,3) = blue;
./base_functions/clevels_define_bcyr.m0000700442630740010010000000315213674531216021045 0ustar  adelmanDomain Usersfunction [c_levels,cmap] = clevels_define(curr_array,length_cbar,end_ref_vec,range_ref_frac)

% function to define contour levels automatically based on distribution of values in curr_array
% length_cbar = number of colors in color bar
% end_ref_vec = acceptable "end" reference numbers for color bar (will be automatically scaled by appropriate power of 10)
% range_ref_frac = fraction/percentile value to use in defining the end reference values


curr_array_sorted = sort(reshape(curr_array,[numel(curr_array) 1]),1,'ascend');
curr_array_sorted = curr_array_sorted((isnan(curr_array_sorted) == 0) & (abs(curr_array_sorted) > 1e-15));
curr_array_lowrange_ref = curr_array_sorted(round(range_ref_frac*numel(curr_array_sorted)));
curr_array_highrange_ref = curr_array_sorted(round((1 - range_ref_frac)*numel(curr_array_sorted)));
curr_array_range_ref = max(abs([curr_array_lowrange_ref curr_array_highrange_ref]));
curr_array_log_ref = log(curr_array_range_ref)/log(10);
log_ref_vec = log(end_ref_vec)/log(10);
[~,closest_log_ref_ind] = min(abs(mod(log_ref_vec - curr_array_log_ref + (1/2),1) - (1/2)));
closest_exp = round(curr_array_log_ref - log_ref_vec(closest_log_ref_ind));
curr_ref = (end_ref_vec(closest_log_ref_ind))*(10^closest_exp);

c_levels = [min([(curr_array_sorted(1) - 1e-5) (-curr_ref - 1e-5)]); ((-curr_ref):(curr_ref/((length_cbar - 2)/2)):curr_ref)'; max([(curr_array_sorted(numel(curr_array_sorted)) + 1e-5) (curr_ref + 1e-5)])];
cmap = colormap(0.85*bcyr(length(c_levels) - 3));    % colormap for contours
cmap = [((2*cmap(1,:)) + ((-1)*cmap(2,:))); cmap; (((-1)*cmap(size(cmap,1) - 1,:)) + (2*cmap(size(cmap,1),:)))];
./base_functions/contourm_overlay.m0000700442630740010010000000146713317611112020461 0ustar  adelmanDomain Usersfunction [] = contourm_overlay(x_array,y_array,z_array,contour_values,line_color,line_style,line_width)

contour_matrix = contourc(x_array,y_array,z_array,contour_values);
pt_ind = 1;
while pt_ind <= size(contour_matrix,2)
    n_pairs_inseg = contour_matrix(2,pt_ind);
%     x_pts_inseg = interp1((1:1:size(x_array,2)),x_array,contour_matrix(1,(pt_ind + 1):(pt_ind + n_pairs_inseg)));
%     y_pts_inseg = interp1((1:1:size(y_array,1)),y_array,contour_matrix(2,(pt_ind + 1):(pt_ind + n_pairs_inseg)));
    x_pts_inseg = contour_matrix(1,(pt_ind + 1):(pt_ind + n_pairs_inseg));
    y_pts_inseg = contour_matrix(2,(pt_ind + 1):(pt_ind + n_pairs_inseg));
    linem(y_pts_inseg,x_pts_inseg,ones(1,n_pairs_inseg),'Color',line_color,'LineStyle',line_style,'LineWidth',line_width)
    
    pt_ind = pt_ind + n_pairs_inseg + 1;
end
./base_functions/contourm_plot_with_landmask.m0000700442630740010010000001256414156521012022664 0ustar  adelmanDomain Usersfunction [h,axm,cbar] = contourm_plot_with_landmask(fig_handle,longitude,latitude,value_array,c_levels,cmap,cbar_location,longwest,longeast,latsouth,latnorth,projection,origin)

% sets up a contour plot that has land masked out



% create first "dummy" contour plot

% lon_offset = mean(longitude([1 length(longitude)])) - mean([longwest longeast]);
if max(max(longitude)) - min(min(longitude)) > 355
    lon_offset = mean([min(min(longitude)) max(max(longitude))]) - mean([longwest longeast]);
else
    lon_offset = 0;
end
longwest_offset = longwest + lon_offset;
longeast_offset = longeast + lon_offset;


fig1000 = figure(1000);
colormap(cmap)
fig_paper_pos = get(fig1000,'PaperPosition');
fig_paper_pos(4) = ((latnorth - latsouth)/((longeast - longwest)*(cosd(mean([latsouth latnorth])))))*fig_paper_pos(3);
fig_paper_pos(3:4) = 4*fig_paper_pos(3:4);
fig_pos = get(fig1000,'Position');
fig_pos(4) = ((latnorth - latsouth)/((longeast - longwest)*(cosd(mean([latsouth latnorth])))))*fig_pos(3);
fig_pos(3:4) = 4*fig_pos(3:4);
set(fig1000,'PaperPosition',fig_paper_pos,'Position',fig_pos,'PaperSize',[44 34])
[C,g] = contourf(longitude,latitude,value_array,c_levels);
caxis([min(c_levels) max(c_levels)])
set(g,'edgecolor','none')
set(gca,'xlim',[longwest_offset longeast_offset],'ylim',[latsouth latnorth])
daspect([1 cosd(mean([latsouth latnorth])) 1])
set(gca,'position',[0 0 1 1],'units','normalized','Visible','off')
pause(20)
fig_pos = get(fig1000,'Position');
fig_pos(4) = ((latnorth - latsouth)/((longeast - longwest)*(cosd(mean([latsouth latnorth])))))*fig_pos(3);
set(fig1000,'Position',fig_pos)
print(fig1000,'contour_plot_temporary.png','-dpng','-r300')

% get(fig1000,'Position')
close(fig1000)


if longwest < -360
    lon_offset_landfind = lon_offset + (360*(floor((longwest - (-360))/360)));
elseif longeast > 360
    lon_offset_landfind = lon_offset - (360*(ceil((longeast - 360)/360)));
end
longwest_offset = longwest + lon_offset_landfind;
longeast_offset = longeast + lon_offset_landfind;


% get contour plot and overlay land mask

rgb_array = imread('contour_plot_temporary.png');
rgb_array = flip(permute(rgb_array,[2 1 3]),2);

delete('contour_plot_temporary.png')


size_rgb_array = size(rgb_array);
frac_first_part = (min([longeast_offset 0]) - max([(-360) longwest_offset]))/(min([longeast_offset 360]) - max([(-360) longwest_offset]));
if frac_first_part <= 0
    landmask_ind = landfind_indices(longwest_offset,longeast_offset,latsouth,latnorth,size_rgb_array(1:2));
elseif frac_first_part >= 1
    landmask_ind = landfind_indices(longwest_offset + 360,longeast_offset + 360,latsouth,latnorth,size_rgb_array(1:2));
else
    landmask_i_ind = [(mod(landfind_indices(longwest_offset + 360,(longwest_offset*(1 - (floor(frac_first_part*(size_rgb_array(1)))/(frac_first_part*(size_rgb_array(1)))))) + 360,latsouth,latnorth,floor([frac_first_part 1].*(size_rgb_array(1:2)))) - 1,floor(frac_first_part*(size_rgb_array(1)))) + 1); (mod(landfind_indices(longwest_offset*(1 - (floor(frac_first_part*(size_rgb_array(1)))/(frac_first_part*(size_rgb_array(1))))),longeast_offset,latsouth,latnorth,size_rgb_array(1:2) - floor([frac_first_part 0].*(size_rgb_array(1:2)))) - 1,size_rgb_array(1) - floor(frac_first_part*(size_rgb_array(1)))) + 1 + (floor(frac_first_part*(size_rgb_array(1)))))];
    landmask_j_ind = [(ceil(landfind_indices(longwest_offset + 360,(longwest_offset*(1 - (floor(frac_first_part*(size_rgb_array(1)))/(frac_first_part*(size_rgb_array(1)))))) + 360,latsouth,latnorth,floor([frac_first_part 1].*(size_rgb_array(1:2))))/(floor(frac_first_part*(size_rgb_array(1)))))); (ceil(landfind_indices(longwest_offset*(1 - (floor(frac_first_part*(size_rgb_array(1)))/(frac_first_part*(size_rgb_array(1))))),longeast_offset,latsouth,latnorth,size_rgb_array(1:2) - floor([frac_first_part 0].*(size_rgb_array(1:2))))/(size_rgb_array(1) - floor(frac_first_part*(size_rgb_array(1))))))];
    landmask_ind = double(landmask_i_ind) + ((size_rgb_array(1))*(double(landmask_j_ind) - 1));
end
rgb_array_reshaped = reshape(rgb_array,[prod(size_rgb_array(1:2)) 3]);

% convert black undefined areas in SSH field to white shading
black_shaded_ind = find((rgb_array_reshaped(:,1) == 0) & (rgb_array_reshaped(:,2) == 0) & (rgb_array_reshaped(:,3) == 0));
rgb_array_reshaped(black_shaded_ind,:) = 255*ones(length(black_shaded_ind),3);

% put black mask on land areas
rgb_array_reshaped(landmask_ind,:) = zeros(length(landmask_ind),3);

rgb_array_masked = reshape(rgb_array_reshaped,size_rgb_array);


% plot contour map with land mask

fig_handle;
axm = axesm(projection,'origin',origin);
% R = georefcells([latsouth latnorth],[longwest longeast],(size_rgb_array(1:2))');
h = geoshow(repmat((latsouth:((latnorth - latsouth)/(size(rgb_array,2) - 1)):latnorth)',[1 size(rgb_array,1)]),repmat(longwest_offset:((longeast_offset - longwest_offset)/(size(rgb_array,1) - 1)):longeast_offset,[size(rgb_array,2) 1]),permute(rgb_array_masked,[2 1 3]));

% plot color bar
colormap(cmap)
cbar = colorbar('location',cbar_location);
cbar_labels = cell(1,length(c_levels));
if length(c_levels) < 10
    for n_label = 1:length(c_levels)
        cbar_labels{n_label} = num2str(c_levels(n_label));
    end
else
    for n_label = 1:length(c_levels)
        if mod(n_label,2) == 1
            cbar_labels{n_label} = num2str(c_levels(n_label));
        else
            cbar_labels{n_label} = '';
        end
    end
end
set(cbar,'xtick',(0:1:length(cmap)) - 0.5,'xticklabel',cbar_labels,'FontSize',16)./base_functions/contour_overlay.m0000700442630740010010000000146512741513710020307 0ustar  adelmanDomain Usersfunction [] = contour_overlay(x_array,y_array,z_array,contour_values,line_color,line_style,line_width)

contour_matrix = contourc(x_array,y_array,z_array,contour_values);
pt_ind = 1;
while pt_ind <= size(contour_matrix,2)
    n_pairs_inseg = contour_matrix(2,pt_ind);
%     x_pts_inseg = interp1((1:1:size(x_array,2)),x_array,contour_matrix(1,(pt_ind + 1):(pt_ind + n_pairs_inseg)));
%     y_pts_inseg = interp1((1:1:size(y_array,1)),y_array,contour_matrix(2,(pt_ind + 1):(pt_ind + n_pairs_inseg)));
    x_pts_inseg = contour_matrix(1,(pt_ind + 1):(pt_ind + n_pairs_inseg));
    y_pts_inseg = contour_matrix(2,(pt_ind + 1):(pt_ind + n_pairs_inseg));
    line(x_pts_inseg,y_pts_inseg,ones(1,n_pairs_inseg),'Color',line_color,'LineStyle',line_style,'LineWidth',line_width)
    
    pt_ind = pt_ind + n_pairs_inseg + 1;
end
./base_functions/contour_plot_with_landmask.m0000700442630740010010000000573013473543334022520 0ustar  adelmanDomain Usersfunction [h,cbar] = contour_plot_with_landmask(fig_handle,longitude,latitude,value_array,c_levels,cmap,cbar_location,longwest,longeast,latsouth,latnorth)

% sets up a contour plot that has land masked out


% adjust for continuity of longitude values
if longeast <= longwest
    longeast = longeast + (360*(floor((longwest - longeast)/360) + 1));
end


% create first "dummy" contour plot

fig1000 = figure(1000);
colormap(cmap)
fig_paper_pos = get(fig1000,'PaperPosition');
fig_paper_pos(4) = ((latnorth - latsouth)/(cosd(mean([latsouth latnorth]))*(longeast - longwest)))*fig_paper_pos(3);
fig_paper_pos(3:4) = 2*fig_paper_pos(3:4);
set(fig1000,'PaperPosition',fig_paper_pos,'PaperSize',[22 17])
[C,g] = contourf(longitude,latitude,value_array,c_levels);
caxis([min(c_levels) max(c_levels)])
set(g,'edgecolor','none')
set(gca,'xlim',[longwest longeast],'ylim',[latsouth latnorth])
daspect([1 cosd(mean([latsouth latnorth])) 1])
set(gca,'position',[0 0 1 1],'units','normalized','Visible','off')
pause(5)
plotboxaspectratio = get(gca,'PlotBoxAspectRatio');
fig_pos = get(gcf,'Position');
ratio_plotbox = plotboxaspectratio(2)/plotboxaspectratio(1);
ratio_fig = fig_pos(4)/fig_pos(3);
if ratio_fig > ratio_plotbox
    fig_pos(4) = ratio_plotbox*fig_pos(3);
else
    fig_pos(3) = fig_pos(4)/ratio_plotbox;
end
set(gcf,'Position',fig_pos)

print(gcf,'contour_plot_temporary.png','-dpng','-r300')
close(fig1000)



% get contour plot and overlay land mask

rgb_array = imread('contour_plot_temporary.png');
rgb_array = flip(permute(rgb_array,[2 1 3]),2);

delete('contour_plot_temporary.png')


size_rgb_array = size(rgb_array);
landmask_ind = landfind_indices(longwest,longeast,latsouth,latnorth,size_rgb_array(1:2));
rgb_array_reshaped = reshape(rgb_array,[prod(size_rgb_array(1:2)) 3]);

% convert black undefined areas in SSH field to white shading
black_shaded_ind = find((rgb_array_reshaped(:,1) == 0) & (rgb_array_reshaped(:,2) == 0) & (rgb_array_reshaped(:,3) == 0));
rgb_array_reshaped(black_shaded_ind,:) = 255*ones(length(black_shaded_ind),3);

% put black mask on land areas
rgb_array_reshaped(landmask_ind,:) = zeros(length(landmask_ind),3);

rgb_array_masked = reshape(rgb_array_reshaped,size_rgb_array);


% plot contour map with land mask

fig_handle;
h = image(longwest:((longeast - longwest)/(size(rgb_array,1) - 1)):longeast,(latsouth:((latnorth - latsouth)/(size(rgb_array,2) - 1)):latnorth)',permute(rgb_array_masked,[2 1 3]));
daspect([1 cosd(mean([latsouth latnorth])) 1])

% plot color bar
colormap(cmap)
cbar = colorbar('location',cbar_location);
cbar_labels = cell(1,length(c_levels));
if length(c_levels) < 10
    for n_label = 1:length(c_levels)
        cbar_labels{n_label} = num2str(c_levels(n_label));
    end
else
    for n_label = 1:length(c_levels)
        if mod(n_label,2) == 1
            cbar_labels{n_label} = num2str(c_levels(n_label));
        else
            cbar_labels{n_label} = '';
        end
    end
end
set(cbar,'xtick',(0:1:length(cmap)) - 0.5,'xticklabel',cbar_labels,'FontSize',16)
./base_functions/correlation_scalar_scalar_uncert_bounds.m0000700442630740010010000003152013560072520025174 0ustar  adelmanDomain Usersfunction [correlation_array_at_lags,correlation_dof_array_zero_lag,integral_scale_1,integral_scale_2,correlation_array_low_mag_bound,correlation_array_high_mag_bound,lags] = correlation_scalar_scalar_uncert_bounds(input_array_1,input_array_2,n_dim,delta_dim,delta_lag,lag_range_to_test,confidence_level)


dim_length = size(input_array_1,n_dim);

size_array = size(input_array_1);
prod_size_not_inc_corr_dim = prod(size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));

input_array_1_permuted = reshape(permute(input_array_1,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_corr_dim]);
input_array_2_permuted = reshape(permute(input_array_2,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_corr_dim]);

nan_mask_input_1 = ones(size(input_array_1_permuted));
nan_mask_input_1(union(union(find(isnan(input_array_1_permuted) == 1),find(isinf(input_array_1_permuted) == 1)),find(abs(input_array_1_permuted) < (1e-15)*max(max(abs(input_array_1_permuted)))))) = 0;
nan_mask_input_2 = ones(size(input_array_2_permuted));
nan_mask_input_2(union(union(find(isnan(input_array_2_permuted) == 1),find(isinf(input_array_2_permuted) == 1)),find(abs(input_array_2_permuted) < (1e-15)*max(max(abs(input_array_2_permuted)))))) = 0;
% nan_mask_input = nan_mask_input_1 & nan_mask_input_2;


input_array_1_zeronans = input_array_1_permuted;
input_array_1_zeronans(abs(nan_mask_input_1) < 1e-5) = 0;
input_array_1_zeronans = input_array_1_zeronans - repmat(sum(nan_mask_input_1.*input_array_1_zeronans,1)./sum(nan_mask_input_1,1),[dim_length 1]);
input_array_2_zeronans = input_array_2_permuted;
input_array_2_zeronans(abs(nan_mask_input_2) < 1e-5) = 0;
input_array_2_zeronans = input_array_2_zeronans - repmat(sum(nan_mask_input_2.*input_array_2_zeronans,1)./sum(nan_mask_input_2,1),[dim_length 1]);

% remove means and trends
index_array_nomean = repmat((1:1:size(input_array_1_zeronans,1))',[1 prod_size_not_inc_corr_dim]) - repmat((sum(nan_mask_input_1.*repmat((1:1:size(input_array_1_zeronans,1))',[1 prod_size_not_inc_corr_dim]),1)./sum(nan_mask_input_1,1)),[size(input_array_1_zeronans,1) 1]);
input_array_1_zeronans_nomean = input_array_1_zeronans - repmat((sum(nan_mask_input_1.*input_array_1_zeronans,1)./sum(nan_mask_input_1,1)),[size(input_array_1_zeronans,1) 1]);
input_array_1_zeronans = input_array_1_zeronans_nomean - (index_array_nomean.*repmat((1./sum((nan_mask_input_1.*index_array_nomean).^2,1)).*sum(nan_mask_input_1.*index_array_nomean.*input_array_1_zeronans_nomean,1),[size(input_array_1_zeronans,1) 1]));
index_array_nomean = repmat((1:1:size(input_array_2_zeronans,1))',[1 prod_size_not_inc_corr_dim]) - repmat((sum(nan_mask_input_2.*repmat((1:1:size(input_array_2_zeronans,1))',[1 prod_size_not_inc_corr_dim]),1)./sum(nan_mask_input_2,1)),[size(input_array_2_zeronans,1) 1]);
input_array_2_zeronans_nomean = input_array_2_zeronans - repmat((sum(nan_mask_input_2.*input_array_2_zeronans,1)./sum(nan_mask_input_2,1)),[size(input_array_2_zeronans,1) 1]);
input_array_2_zeronans = input_array_2_zeronans_nomean - (index_array_nomean.*repmat((1./sum((nan_mask_input_2.*index_array_nomean).^2,1)).*sum(nan_mask_input_2.*index_array_nomean.*input_array_2_zeronans_nomean,1),[size(input_array_2_zeronans,1) 1]));


% compute degrees of freedom

delta_lag_ind = round(delta_lag/delta_dim);
lags = round((delta_lag*(ceil(lag_range_to_test(1)/delta_lag):1:floor(lag_range_to_test(2)/delta_lag))')/delta_dim)*delta_dim;

autocorr_lags = ((delta_dim*delta_lag_ind):(delta_dim*delta_lag_ind):max(abs(lags)))';
N = length(autocorr_lags);


lag_step = delta_lag_ind;    % number of time points between lags to test for dof

input_1_autocorr_at_lags = NaN([N prod_size_not_inc_corr_dim]);
input_2_autocorr_at_lags = NaN([N prod_size_not_inc_corr_dim]);
input_1_integral_scale = NaN([1 prod_size_not_inc_corr_dim]);
input_2_integral_scale = NaN([1 prod_size_not_inc_corr_dim]);
n_lag_step = 1;
lag_ind = lag_step;
complete_flag = 0;
while ((complete_flag == 0) && (n_lag_step <= N))

    part_1_ind = (1:1:(dim_length - lag_ind))';
    part_2_ind = ((1 + lag_ind):1:dim_length)';
    nan_mask_combined = nan_mask_input_1(part_1_ind,:) & nan_mask_input_1(part_2_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_1_zeronans(part_1_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_1_zeronans(part_2_ind,:),1)./sum(nan_mask_combined,1);
    input_1_autocorr_at_lags(n_lag_step,:) = sum(nan_mask_combined.*(input_array_1_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).*(input_array_1_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_1_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_1_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])).^2),1)).^(1/2)));
    
    nan_mask_combined = nan_mask_input_2(part_1_ind,:) & nan_mask_input_2(part_2_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_2_zeronans(part_1_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_2_zeronans(part_2_ind,:),1)./sum(nan_mask_combined,1);
    input_2_autocorr_at_lags(n_lag_step,:) = sum(nan_mask_combined.*(input_array_2_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).*(input_array_2_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_2_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_2_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])).^2),1)).^(1/2)));
    
    if n_lag_step == N
        first_neg_ind_1 = find((isnan(input_1_integral_scale) == 1) & (isnan(input_1_autocorr_at_lags(1,:)) == 0));
        first_neg_ind_2 = find((isnan(input_2_integral_scale) == 1) & (isnan(input_2_autocorr_at_lags(1,:)) == 0));
    else
        first_neg_ind_1 = find(((input_1_autocorr_at_lags(n_lag_step,:) < 0) | ((isnan(input_1_autocorr_at_lags(n_lag_step,:)) == 1) & (isnan(input_1_autocorr_at_lags(1,:)) == 0))) & (isnan(input_1_integral_scale) == 1));
        first_neg_ind_2 = find(((input_2_autocorr_at_lags(n_lag_step,:) < 0) | ((isnan(input_2_autocorr_at_lags(n_lag_step,:)) == 1) & (isnan(input_2_autocorr_at_lags(1,:)) == 0))) & (isnan(input_2_integral_scale) == 1));
    end
    
    first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_corr_dim]);
    first_neg_ind_mask_2D(first_neg_ind_1) = 1;
    first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
    first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
    input_1_autocorr_at_lags_zeronans = input_1_autocorr_at_lags;
    input_1_autocorr_at_lags_zeronans(isnan(input_1_autocorr_at_lags) == 1) = 0;
    curr_input_1_integral_scale = delta_dim*lag_step*(ceil(first_neg_ind_mask(1,:)) + max([zeros([1 prod_size_not_inc_corr_dim]); (2*sum(first_neg_ind_mask.*input_1_autocorr_at_lags_zeronans(1:n_lag_step,:),1))],[],1));
    clear input_1_autocorr_at_lags_zeronans
    if n_lag_step > 1
        nan_case_integral_scale = min([repmat((delta_dim*lag_ind),[1 prod_size_not_inc_corr_dim]); (delta_dim*lag_step*(first_neg_ind_mask(1,:) + (2*sum((first_neg_ind_mask(1:(n_lag_step - 1),:)).*input_1_autocorr_at_lags(1:(n_lag_step - 1),:),1))))],[],1);
        curr_input_1_integral_scale(first_neg_ind_1(isnan(curr_input_1_integral_scale(first_neg_ind_1)) == 1)) = nan_case_integral_scale(first_neg_ind_1(isnan(curr_input_1_integral_scale(first_neg_ind_1)) == 1));
    end
    input_1_integral_scale(first_neg_ind_1) = curr_input_1_integral_scale(first_neg_ind_1);

    first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_corr_dim]);
    first_neg_ind_mask_2D(first_neg_ind_2) = 1;
    first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
    first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
    input_2_autocorr_at_lags_zeronans = input_2_autocorr_at_lags;
    input_2_autocorr_at_lags_zeronans(isnan(input_2_autocorr_at_lags) == 1) = 0;
    curr_input_2_integral_scale = delta_dim*lag_step*(ceil(first_neg_ind_mask(1,:)) + max([zeros([1 prod_size_not_inc_corr_dim]); (2*sum(first_neg_ind_mask.*input_2_autocorr_at_lags_zeronans(1:n_lag_step,:),1))],[],1));
    clear input_2_autocorr_at_lags_zeronans
    if n_lag_step > 1
        nan_case_integral_scale = min([repmat((delta_dim*lag_ind),[1 prod_size_not_inc_corr_dim]); (delta_dim*lag_step*(first_neg_ind_mask(1,:) + (2*sum((first_neg_ind_mask(1:(n_lag_step - 1),:)).*input_2_autocorr_at_lags(1:(n_lag_step - 1),:),1))))],[],1);
        curr_input_2_integral_scale(first_neg_ind_2(isnan(curr_input_2_integral_scale(first_neg_ind_2)) == 1)) = nan_case_integral_scale(first_neg_ind_2(isnan(curr_input_2_integral_scale(first_neg_ind_2)) == 1));
    end
    input_2_integral_scale(first_neg_ind_2) = curr_input_2_integral_scale(first_neg_ind_2);
    
    if ((sum(sum(isnan(input_1_integral_scale))) == 0) && (sum(sum(isnan(input_2_integral_scale))) == 0))
        complete_flag = 1;
    end
    
%     disp(['lag_ind = ',num2str(lag_ind),' (out of ',num2str(N*lag_step),')'])
    
    n_lag_step = n_lag_step + 1;
    lag_ind = lag_ind + lag_step;

end

input_corr_integral_scale = min([input_1_integral_scale; input_2_integral_scale],[],1);
corr_dof_array = (delta_dim*sum((nan_mask_input_1 & nan_mask_input_2),1))./input_corr_integral_scale;

integral_scale_1 = permute(reshape(input_1_integral_scale,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
integral_scale_2 = permute(reshape(input_2_integral_scale,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
correlation_dof_array_zero_lag = permute(reshape(corr_dof_array,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);

% compute pointwise correlations

corrcoeff_array = NaN([length(lags) prod_size_not_inc_corr_dim]);
corrcoeff_low_mag = NaN([length(lags) prod_size_not_inc_corr_dim]);
corrcoeff_high_mag = NaN([length(lags) prod_size_not_inc_corr_dim]);
for n_lag_ind = 1:length(lags)
    curr_lag = lags(n_lag_ind);
    lag_ind = round(curr_lag/delta_dim);
    input_1_corr_ind = (max([1 (1 - lag_ind)]):min([dim_length (dim_length - lag_ind)]))';
    input_2_corr_ind = (max([1 (1 + lag_ind)]):min([dim_length (dim_length + lag_ind)]))';
    
    nan_mask_combined = nan_mask_input_1(input_1_corr_ind,:) & nan_mask_input_2(input_2_corr_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_1_zeronans(input_1_corr_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_2_zeronans(input_2_corr_ind,:),1)./sum(nan_mask_combined,1);
    corrcoeff_array(n_lag_ind,:) = sum(nan_mask_combined.*(input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).*(input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])).^2),1)).^(1/2)));
    
    curr_corr_dof_array = (delta_dim*sum(nan_mask_combined,1))./input_corr_integral_scale;
    
    stderror_Z = (curr_corr_dof_array - 3).^(-1/2);

    stderror_Z(curr_corr_dof_array < 3) = Inf;

    normalized_Z_value = 0.5*(log(1 + abs(corrcoeff_array(n_lag_ind,:))) - log(1 - abs(corrcoeff_array(n_lag_ind,:))));
    Z_confidence_lower_mag = normalized_Z_value - (stderror_Z*(erfinv(confidence_level))*(2^(1/2)));
    r_confidence_lower_mag = 1 - (2./(exp(2*Z_confidence_lower_mag) + 1));
    Z_confidence_higher_mag = normalized_Z_value + (stderror_Z*(erfinv(confidence_level))*(2^(1/2)));
    r_confidence_higher_mag = 1 - (2./(exp(2*Z_confidence_higher_mag) + 1));
    
    corrcoeff_low_mag(n_lag_ind,:) = (sign(corrcoeff_array(n_lag_ind,:))).*r_confidence_lower_mag;
    corrcoeff_high_mag(n_lag_ind,:) = (sign(corrcoeff_array(n_lag_ind,:))).*r_confidence_higher_mag;
end


correlation_array_at_lags = permute(reshape(corrcoeff_array,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
correlation_array_low_mag_bound = permute(reshape(corrcoeff_low_mag,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
correlation_array_high_mag_bound = permute(reshape(corrcoeff_high_mag,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);./base_functions/interp2_fft.m0000700442630740010010000001005314243217252017271 0ustar  adelmanDomain Usersfunction [output_field] = interp2_fft(input_i_vec,input_j_vec,input_field,output_i,output_j)

% interpolate from 2-D grid using FFT (works best when input grid given by input_i_vec and input_j_vec has consistent spacing
% Note: (output_i, output_j) are vectors containing coordinates for individual points, NOT basis coordinates for arrays like the input vectors.  If output is desired in array format, use reshape accordingly.

warning('off','MATLAB:nearlySingularMatrix')
warning('off','MATLAB:singularMatrix')

size_input_field = size(input_field);

interp_coeffs = fft(fft(input_field,[],1)/(size(input_field,1)),[],2)/(size(input_field,2));

interp_i_ind = interp1(input_i_vec,0:1:(length(input_i_vec) - 1),reshape(output_i,[numel(output_i) 1]));
interp_j_ind = interp1(input_j_vec,0:1:(length(input_j_vec) - 1),reshape(output_j,[numel(output_j) 1]));
k_vec = (mod((0:1:(length(input_i_vec) - 1))/length(input_i_vec) + (1/2),1) - (1/2))';
l_vec = (mod((0:1:(length(input_j_vec) - 1))/length(input_j_vec) + (1/2),1) - (1/2))';


if (length(k_vec))*(length(l_vec)) >= length(interp_i_ind)
    numel_array_limit = 1e8;
    largest_array_size = prod([size_input_field length(output_i)]);
    n_subsets = ceil(largest_array_size/numel_array_limit);
    length_subset = ceil(length(output_i)/n_subsets);
    start_ind_subset = (1:length_subset:length(output_i))';
    for curr_subset = 1:length(start_ind_subset)
        curr_ind_subset = (start_ind_subset(curr_subset):1:min([(start_ind_subset(curr_subset) + length_subset) length(output_i)]))';
        G = reshape(exp(1i*(((2*pi)*repmat(repmat(reshape(k_vec,[1 length(k_vec) 1]),[length(curr_ind_subset) 1]).*repmat(interp_i_ind(curr_ind_subset),[1 length(k_vec)]),[1 1 length(l_vec)])) + ((2*pi)*repmat(repmat(reshape(l_vec,[1 1 length(l_vec)]),[length(curr_ind_subset) 1 1]).*repmat(interp_j_ind(curr_ind_subset),[1 1 length(l_vec)]),[1 length(k_vec) 1])))),[length(curr_ind_subset) (length(k_vec)*length(l_vec))]);
        if length(size_input_field) > 2
            if curr_subset == 1
                output_field = NaN([size(output_i) size_input_field(3:length(size_input_field))]);
            end
            output_ind_str = ['(curr_ind_subset,1',repmat(',:',[1 (length(size_input_field) - 2)]),')'];
            eval(['output_field',output_ind_str,' = reshape(real(G*reshape(interp_coeffs,[prod(size_input_field(1:2)) prod(size_input_field(3:length(size_input_field)))])),[length(curr_ind_subset) 1 size_input_field(3:length(size_input_field))]);'])
        else
            if curr_subset == 1
                output_field = NaN(size(output_i));
            end
            output_field(curr_ind_subset) = reshape(real(G*reshape(interp_coeffs,[numel(interp_coeffs) 1])),[length(curr_ind_subset) 1]);
        end
    end
else
    if length(size_input_field) > 2
        m = reshape(interp_coeffs,[prod(size_input_field(1:2)) size_input_field(3:(length(size_input_field)))]);
        output_field = NaN([size(output_i) size_input_field(3:length(size_input_field))]);
        for interp_i = 1:length(interp_i_ind)
            curr_i = interp_i_ind(interp_i);
            curr_j = interp_j_ind(interp_i);
            output_ind_str = ['(interp_i,1',repmat(',:',[1 (length(size_input_field) - 2)]),')'];
            eval(['output_field',output_ind_str,' = reshape(real((exp(1i*(((2*pi)*repmat(reshape(k_vec,[1 length(k_vec)])*curr_i,[1 length(l_vec)])) + ((2*pi)*reshape(repmat(reshape(l_vec,[1 length(l_vec)]),[length(k_vec) 1]),[1 ((length(k_vec))*(length(l_vec)))])*curr_j))))*m),[1 1 size_input_field(3:length(size_input_field))]);'])
        end
    else
        m = reshape(interp_coeffs,[numel(interp_coeffs) 1]);
        output_field = NaN(size(output_i));
        for interp_i = 1:length(interp_i_ind)
            curr_i = interp_i_ind(interp_i);
            curr_j = interp_j_ind(interp_i);
            output_field(interp_i) = real((exp(1i*(((2*pi)*repmat(reshape(k_vec,[1 length(k_vec)])*curr_i,[1 length(l_vec)])) + ((2*pi)*reshape(repmat(reshape(l_vec,[1 length(l_vec)]),[length(k_vec) 1]),[1 ((length(k_vec))*(length(l_vec)))])*curr_j))))*m);
        end
    end
end./base_functions/irregular_clevels_plot.m0000700442630740010010000000144413473615506021632 0ustar  adelmanDomain Usersfunction [output_array,output_clabels] = irregular_clevels_plot(input_array,c_levels);

% contour plotting with irregularly spaced contours

output_array = input_array;

for n_c_level = 1:(length(c_levels) - 1)
    lower_c_level = c_levels(n_c_level);
    upper_c_level = c_levels(n_c_level + 1);
    
    in_c_level_ind = find((input_array >= lower_c_level) & (input_array < upper_c_level));
    
%     output_array(in_c_level_ind) = n_c_level - 0.5;
    output_array(in_c_level_ind) = n_c_level - 1 + ((input_array(in_c_level_ind) - lower_c_level)/(upper_c_level - lower_c_level));
    
end

% output_clabels = mat2cell(num2str(c_levels),[1 length(c_levels)]);
output_clabels = cell(length(c_levels),1);
for n_clabel = 1:length(c_levels)
    output_clabels{n_clabel} = num2str(c_levels(n_clabel));
end
./base_functions/landfind_indices.m0000700442630740010010000001320113473542412020325 0ustar  adelmanDomain Usersfunction [landmask_ind] = landfind_indices(longwest,longeast,latsouth,latnorth,pixel_dim)

% function that outputs land mask indices, given four coordinate bounds (MUST be the same as the min and max of the image being overlaid) and the pixel dimensions


%begin bathymetry calculations

ncid_topo = netcdf.open('topo2.grd','NC_NOWRITE');
% x_range = netcdf.getVar(ncid_topo,0);
% y_range = netcdf.getVar(ncid_topo,1);
% z_range = netcdf.getVar(ncid_topo,2);
% spacing = netcdf.getVar(ncid_topo,3);
% dimension = netcdf.getVar(ncid_topo,4);
z = netcdf.getVar(ncid_topo,5);

spacing = 1/30;     %spacing of bathymetry data, in degrees

%masking bathymetry threshold (in meters elevation)

mask_threshold = 0;

% %isolate bathymetry data
% 
% bathy_ns_indices = (floor((90-latnorth)/spacing)+1:ceil((90-latsouth)/spacing))';
% 
% lat_bathy = 90 - (double(bathy_ns_indices-1)*spacing + spacing/2);
% 
% bathy_ew_indices = (floor(longwest/spacing)+1:ceil(longeast/spacing))';
% 
% long_bathy = double(bathy_ew_indices-1)*spacing + spacing/2;

lat_bathy = 90 - ((0.5*spacing):spacing:180)';

long_bathy = ((0.5*spacing):spacing:360)';


z_reshaped = reshape(z,360/spacing,180/spacing);

clear z


% z_range = z_reshaped(bathy_ew_indices,bathy_ns_indices);
% 
% clear z_reshaped
% 
% z_range = flipdim(z_range,2);

z_reshaped = flipdim(z_reshaped,2);

lat_bathy = flipdim(lat_bathy,1);

mask_ind = find(z_reshaped >= mask_threshold);

zero_padding_dim = [5 5];

z_land_mask = zeros(size(z_reshaped));
% mask_ind_padded = mask_ind + (zero_padding_dim(2)*(size(z_reshaped,1) + (2*zero_padding_dim(1)))) + zero_padding(1);
z_land_mask(mask_ind) = 1;

% resolutions of imported image

% x_res_image = (longeast - longwest)/(pixel_dim(1) - 1);
% y_res_image = (latnorth - latsouth)/(pixel_dim(2) - 1);
x_res_image = (mod(longeast - longwest - 1e-10,360) + 1e-10)/(pixel_dim(1));
y_res_image = (latnorth - latsouth)/(pixel_dim(2));

x_res_ratio = x_res_image/spacing;
y_res_ratio = y_res_image/spacing;


% apply smoothing filter
x_smooth_filter_vec = zeros(size(z_land_mask,1),1);
y_smooth_filter_vec = zeros(size(z_land_mask,2),1);

edge_full_ind = floor((x_res_ratio + 1)/2);
x_smooth_filter_vec(1:edge_full_ind) = 1;
x_smooth_filter_vec(edge_full_ind + 1) = ((x_res_ratio + 1)/2) - edge_full_ind;
x_smooth_filter_vec((length(x_smooth_filter_vec) - edge_full_ind + 1):length(x_smooth_filter_vec)) = flipdim(x_smooth_filter_vec(2:(edge_full_ind + 1)),1);

y_smooth_filter_vec(1:edge_full_ind) = 1;
y_smooth_filter_vec(edge_full_ind + 1) = ((y_res_ratio + 1)/2) - edge_full_ind;
y_smooth_filter_vec((length(y_smooth_filter_vec) - edge_full_ind + 1):length(y_smooth_filter_vec)) = flipdim(y_smooth_filter_vec(2:(edge_full_ind + 1)),1);

x_smooth_filter_array = repmat(x_smooth_filter_vec,[1 size(z_land_mask,2)]);
y_smooth_filter_array = repmat(y_smooth_filter_vec',[size(z_land_mask,1) 1]);

smooth_filter_array = x_smooth_filter_array.*y_smooth_filter_array;
smooth_filter_array = smooth_filter_array/(sum(sum(smooth_filter_array)));

fft_smooth_filter_array = fft(fft(smooth_filter_array,[],1),[],2);

z_land_mask_smoothed = ifft(ifft(fft_smooth_filter_array.*fft(fft(z_land_mask,[],1),[],2),[],2),[],1);


% "interpolate" to the same grid as the image

if mod(longwest,360) >= mod(longeast,360)
    in_bathy_lon_range_ind_1 = find(long_bathy >= mod(longwest,360));
    in_bathy_lon_range_ind_2 = find(long_bathy <= mod(longeast,360));
    in_bathy_lon_range_ind = [in_bathy_lon_range_ind_1; in_bathy_lon_range_ind_2];
else
    in_bathy_lon_range_ind = find((long_bathy >= mod(longwest,360)) & (long_bathy <= mod(longeast,360)));
end
in_bathy_lat_range_ind = find((lat_bathy >= latsouth) & (lat_bathy <= latnorth));
% SW_corner_x_weight_inner = (longwest - long_bathy(in_bathy_lon_range_ind(1) - 1))/spacing;
% SW_corner_y_weight_inner = (latsouth - lat_bathy(min(in_bathy_lat_range_ind) - 1))/spacing;
SW_corner_x_weight_inner = ((mod(longwest - long_bathy(in_bathy_lon_range_ind(1)) + 180,360) - 180)/spacing) + 1;
SW_corner_y_weight_inner = ((latsouth - lat_bathy(min(in_bathy_lat_range_ind)))/spacing) + 1;

x_rel_placing_vec = SW_corner_x_weight_inner + (x_res_ratio*((0:1:(pixel_dim(1) - 1))'));

land_mask_x_interp = NaN(pixel_dim(1),size(z_land_mask_smoothed,2));
for x_ind = 1:length(x_rel_placing_vec)
    curr_x_rel_placing = x_rel_placing_vec(x_ind);
    
    x_bathy_lower_ind = mod(in_bathy_lon_range_ind(1) - 1 + floor(curr_x_rel_placing) - 1,length(long_bathy)) + 1;
    x_bathy_upper_ind = mod(x_bathy_lower_ind + 1 - 1,length(long_bathy)) + 1;
    
    x_upper_weight = curr_x_rel_placing - floor(curr_x_rel_placing);
    x_lower_weight = 1 - x_upper_weight;
    
    land_mask_x_interp(x_ind,:) = (x_lower_weight*z_land_mask_smoothed(x_bathy_lower_ind,:)) + (x_upper_weight*z_land_mask_smoothed(x_bathy_upper_ind,:)); 
    
end

y_rel_placing_vec = SW_corner_y_weight_inner + (y_res_ratio*((0:1:(pixel_dim(2) - 1))'));

land_mask = NaN(pixel_dim);
for y_ind = 1:length(y_rel_placing_vec)
    curr_y_rel_placing = y_rel_placing_vec(y_ind);
    
    y_bathy_lower_ind = max([1 (min(in_bathy_lat_range_ind) - 1 + floor(curr_y_rel_placing))]);
    y_bathy_upper_ind = min([length(lat_bathy) (y_bathy_lower_ind + 1)]);
    
    y_upper_weight = curr_y_rel_placing - floor(curr_y_rel_placing);
    y_lower_weight = 1 - y_upper_weight;
    
    land_mask(:,y_ind) = (y_lower_weight*land_mask_x_interp(:,y_bathy_lower_ind)) + (y_upper_weight*land_mask_x_interp(:,y_bathy_upper_ind)); 
    
end


% find indices of land points
landmask_ind = find(land_mask > 0.5);


% close netCDF file
netcdf.close(ncid_topo);
./base_functions/pcolor_plot_with_landmask.m0000700442630740010010000000567013735021450022317 0ustar  adelmanDomain Usersfunction [h,cbar] = pcolor_plot_with_landmask(fig_handle,longitude,latitude,value_array,c_levels,cmap,cbar_location,longwest,longeast,latsouth,latnorth)

% sets up a contour plot that has land masked out


% adjust for continuity of longitude values
if longeast <= longwest
    longeast = longeast + (360*(floor((longwest - longeast)/360) + 1));
end


% create first "dummy" contour plot

fig1000 = figure(1000);
colormap(cmap)
fig_paper_pos = get(fig1000,'PaperPosition');
fig_paper_pos(4) = ((latnorth - latsouth)/(cosd(mean([latsouth latnorth]))*(longeast - longwest)))*fig_paper_pos(3);
fig_paper_pos(3:4) = 2*fig_paper_pos(3:4);
set(fig1000,'PaperPosition',fig_paper_pos,'PaperSize',[22 17])
pcolor(longitude,latitude,value_array);
shading flat
caxis([min(c_levels) max(c_levels)])
set(gca,'xlim',[longwest longeast],'ylim',[latsouth latnorth])
daspect([1 cosd(mean([latsouth latnorth])) 1])
set(gca,'position',[0 0 1 1],'units','normalized','Visible','off')
pause(10)
plotboxaspectratio = get(gca,'PlotBoxAspectRatio');
fig_pos = get(gcf,'Position');
ratio_plotbox = plotboxaspectratio(2)/plotboxaspectratio(1);
ratio_fig = fig_pos(4)/fig_pos(3);
if ratio_fig > ratio_plotbox
    fig_pos(4) = ratio_plotbox*fig_pos(3);
else
    fig_pos(3) = fig_pos(4)/ratio_plotbox;
end
set(gcf,'Position',fig_pos)

print(gcf,'contour_plot_temporary.png','-dpng','-r300')
close(fig1000)



% get contour plot and overlay land mask

rgb_array = imread('contour_plot_temporary.png');
rgb_array = flip(permute(rgb_array,[2 1 3]),2);

delete('contour_plot_temporary.png')


size_rgb_array = size(rgb_array);
landmask_ind = landfind_indices(longwest,longeast,latsouth,latnorth,size_rgb_array(1:2));
rgb_array_reshaped = reshape(rgb_array,[prod(size_rgb_array(1:2)) 3]);

% convert black undefined areas in SSH field to white shading
black_shaded_ind = find((rgb_array_reshaped(:,1) == 0) & (rgb_array_reshaped(:,2) == 0) & (rgb_array_reshaped(:,3) == 0));
rgb_array_reshaped(black_shaded_ind,:) = 255*ones(length(black_shaded_ind),3);

% put black mask on land areas
rgb_array_reshaped(landmask_ind,:) = zeros(length(landmask_ind),3);

rgb_array_masked = reshape(rgb_array_reshaped,size_rgb_array);


% plot contour map with land mask

fig_handle;
h = image(longwest:((longeast - longwest)/(size(rgb_array,1) - 1)):longeast,(latsouth:((latnorth - latsouth)/(size(rgb_array,2) - 1)):latnorth)',permute(rgb_array_masked,[2 1 3]));
daspect([1 cosd(mean([latsouth latnorth])) 1])

% plot color bar
colormap(cmap)
cbar = colorbar('location',cbar_location);
cbar_labels = cell(1,length(c_levels));
if length(c_levels) < 10
    for n_label = 1:length(c_levels)
        cbar_labels{n_label} = num2str(c_levels(n_label));
    end
else
    for n_label = 1:length(c_levels)
        if mod(n_label,2) == 1
            cbar_labels{n_label} = num2str(c_levels(n_label));
        else
            cbar_labels{n_label} = '';
        end
    end
end
set(cbar,'xtick',(0:1:length(cmap)) - 0.5,'xticklabel',cbar_labels,'FontSize',16)
./base_functions/regression_linear_scalar_scalar.m0000700442630740010010000003654513540461302023445 0ustar  adelmanDomain Usersfunction [regression_array_at_lags,regression_dof_array_zero_lag,regression_array_low_mag_bound,regression_array_high_mag_bound,mean_array_1,mean_array_2,std_dev_estimate_1,std_dev_estimate_2,lags] = regression_linear_scalar_scalar(input_array_1,input_array_2,n_dim,delta_dim,delta_lag,lag_range_to_test,confidence_level,trend_handling_opt)

% trend_handling_opt: 0 = carry out regression (and compute standard deviation) with trend removed, 1 = retain trend in regression and standard deviation calculations

dim_length = size(input_array_1,n_dim);

size_array = size(input_array_1);
prod_size_not_inc_corr_dim = prod(size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]));

input_array_1_permuted = reshape(permute(input_array_1,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_corr_dim]);
input_array_2_permuted = reshape(permute(input_array_2,[n_dim (1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))]),[dim_length prod_size_not_inc_corr_dim]);

nan_mask_input_1 = ones(size(input_array_1_permuted));
nan_mask_input_1(union(union(find(isnan(input_array_1_permuted) == 1),find(isinf(input_array_1_permuted) == 1)),find(abs(input_array_1_permuted) < (1e-15)*max(max(abs(input_array_1_permuted)))))) = 0;
nan_mask_input_2 = ones(size(input_array_2_permuted));
nan_mask_input_2(union(union(find(isnan(input_array_2_permuted) == 1),find(isinf(input_array_2_permuted) == 1)),find(abs(input_array_2_permuted) < (1e-15)*max(max(abs(input_array_2_permuted)))))) = 0;
% nan_mask_input = nan_mask_input_1 & nan_mask_input_2;


input_array_1_zeronans = input_array_1_permuted;
input_array_1_zeronans(abs(nan_mask_input_1) < 1e-5) = 0;
input_array_2_zeronans = input_array_2_permuted;
input_array_2_zeronans(abs(nan_mask_input_2) < 1e-5) = 0;

% remove means and trends
index_array_nomean = repmat((1:1:dim_length)',[1 prod_size_not_inc_corr_dim]) - repmat((sum(nan_mask_input_1.*repmat((1:1:dim_length)',[1 prod_size_not_inc_corr_dim]),1)./sum(nan_mask_input_1,1)),[dim_length 1]);
input_array_1_zeronans_nomean = input_array_1_zeronans - repmat(sum(nan_mask_input_1.*input_array_1_zeronans,1)./sum(nan_mask_input_1,1),[dim_length 1]);
if abs(trend_handling_opt - 0) < 1e-5
    trend_array_1 = index_array_nomean.*repmat((1./sum((nan_mask_input_1.*index_array_nomean).^2,1)).*sum(nan_mask_input_1.*index_array_nomean.*input_array_1_zeronans_nomean,1),[dim_length 1]);
    trend_array_1 = trend_array_1 + (input_array_1_zeronans - input_array_1_zeronans_nomean);
else
    trend_array_1 = input_array_1_zeronans - input_array_1_zeronans_nomean;
end
input_array_1_zeronans = input_array_1_zeronans - trend_array_1;
index_array_nomean = repmat((1:1:dim_length)',[1 prod_size_not_inc_corr_dim]) - repmat((sum(nan_mask_input_2.*repmat((1:1:dim_length)',[1 prod_size_not_inc_corr_dim]),1)./sum(nan_mask_input_2,1)),[dim_length 1]);
input_array_2_zeronans_nomean = input_array_2_zeronans - repmat(sum(nan_mask_input_2.*input_array_2_zeronans,1)./sum(nan_mask_input_2,1),[dim_length 1]);
if abs(trend_handling_opt - 0) < 1e-5
    trend_array_2 = index_array_nomean.*repmat((1./sum((nan_mask_input_2.*index_array_nomean).^2,1)).*sum(nan_mask_input_2.*index_array_nomean.*input_array_2_zeronans_nomean,1),[dim_length 1]);
    trend_array_2 = trend_array_2 + (input_array_2_zeronans - input_array_2_zeronans_nomean);
else
    trend_array_2 = input_array_2_zeronans - input_array_2_zeronans_nomean;
end
input_array_2_zeronans = input_array_2_zeronans - trend_array_2;


% compute degrees of freedom

delta_lag_ind = round(delta_lag/delta_dim);
lags = round((delta_lag*(ceil(lag_range_to_test(1)/delta_lag):1:floor(lag_range_to_test(2)/delta_lag))')/delta_dim)*delta_dim;

autocorr_lags = ((delta_dim*delta_lag_ind):(delta_dim*delta_lag_ind):max(abs(lags)))';
N = length(autocorr_lags);


lag_step = delta_lag_ind;    % number of time points between lags to test for dof

input_1_autocorr_at_lags = NaN([N prod_size_not_inc_corr_dim]);
input_2_autocorr_at_lags = NaN([N prod_size_not_inc_corr_dim]);
input_1_integral_scale = NaN([1 prod_size_not_inc_corr_dim]);
input_2_integral_scale = NaN([1 prod_size_not_inc_corr_dim]);
n_lag_step = 1;
lag_ind = lag_step;
complete_flag = 0;
while ((complete_flag == 0) && (n_lag_step <= N))

    part_1_ind = (1:1:(dim_length - lag_ind))';
    part_2_ind = ((1 + lag_ind):1:dim_length)';
    nan_mask_combined = nan_mask_input_1(part_1_ind,:) & nan_mask_input_1(part_2_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_1_zeronans(part_1_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_1_zeronans(part_2_ind,:),1)./sum(nan_mask_combined,1);
    input_1_autocorr_at_lags(n_lag_step,:) = sum(nan_mask_combined.*(input_array_1_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).*(input_array_1_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_1_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_1_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])).^2),1)).^(1/2)));
    
    nan_mask_combined = nan_mask_input_2(part_1_ind,:) & nan_mask_input_2(part_2_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_2_zeronans(part_1_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_2_zeronans(part_2_ind,:),1)./sum(nan_mask_combined,1);
    input_2_autocorr_at_lags(n_lag_step,:) = sum(nan_mask_combined.*(input_array_2_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).*(input_array_2_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_2_zeronans(part_1_ind,:) - repmat(part_1_mean,[length(part_1_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_2_zeronans(part_2_ind,:) - repmat(part_2_mean,[length(part_2_ind) 1])).^2),1)).^(1/2)));

    if n_lag_step == N
        first_neg_ind_1 = find((isnan(input_1_integral_scale) == 1) & (isnan(input_1_autocorr_at_lags(1,:)) == 0));
        first_neg_ind_2 = find((isnan(input_2_integral_scale) == 1) & (isnan(input_2_autocorr_at_lags(1,:)) == 0));
    else
        first_neg_ind_1 = find(((input_1_autocorr_at_lags(n_lag_step,:) < 0) | ((isnan(input_1_autocorr_at_lags(n_lag_step,:)) == 1) & (isnan(input_1_autocorr_at_lags(1,:)) == 0))) & (isnan(input_1_integral_scale) == 1));
        first_neg_ind_2 = find(((input_2_autocorr_at_lags(n_lag_step,:) < 0) | ((isnan(input_2_autocorr_at_lags(n_lag_step,:)) == 1) & (isnan(input_2_autocorr_at_lags(1,:)) == 0))) & (isnan(input_2_integral_scale) == 1));
    end
    
    first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_corr_dim]);
    first_neg_ind_mask_2D(first_neg_ind_1) = 1;
    first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
    first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
    curr_input_1_integral_scale = delta_dim*lag_step*(ceil(first_neg_ind_mask(1,:)) + max([zeros([1 prod_size_not_inc_corr_dim]); (2*sum(first_neg_ind_mask.*input_1_autocorr_at_lags(1:n_lag_step,:),1))],[],1));
    if n_lag_step > 1
        nan_case_integral_scale = min([repmat((delta_dim*lag_ind),[1 prod_size_not_inc_corr_dim]); (delta_dim*lag_step*(first_neg_ind_mask(1,:) + (2*sum((first_neg_ind_mask(1:(n_lag_step - 1),:)).*input_1_autocorr_at_lags(1:(n_lag_step - 1),:),1))))],[],1);
        curr_input_1_integral_scale(first_neg_ind_1(isnan(curr_input_1_integral_scale(first_neg_ind_1)) == 1)) = nan_case_integral_scale(first_neg_ind_1(isnan(curr_input_1_integral_scale(first_neg_ind_1)) == 1));
    end
    input_1_integral_scale(first_neg_ind_1) = curr_input_1_integral_scale(first_neg_ind_1);

    first_neg_ind_mask_2D = zeros([1 prod_size_not_inc_corr_dim]);
    first_neg_ind_mask_2D(first_neg_ind_2) = 1;
    first_neg_ind_mask = repmat(first_neg_ind_mask_2D,[n_lag_step 1]);
    first_neg_ind_mask(n_lag_step,:) = 0.5*first_neg_ind_mask(n_lag_step,:);
    curr_input_2_integral_scale = delta_dim*lag_step*(ceil(first_neg_ind_mask(1,:)) + max([zeros([1 prod_size_not_inc_corr_dim]); (2*sum(first_neg_ind_mask.*input_2_autocorr_at_lags(1:n_lag_step,:),1))],[],1));
    if n_lag_step > 1
        nan_case_integral_scale = min([repmat((delta_dim*lag_ind),[1 prod_size_not_inc_corr_dim]); (delta_dim*lag_step*(first_neg_ind_mask(1,:) + (2*sum((first_neg_ind_mask(1:(n_lag_step - 1),:)).*input_2_autocorr_at_lags(1:(n_lag_step - 1),:),1))))],[],1);
        curr_input_2_integral_scale(first_neg_ind_2(isnan(curr_input_2_integral_scale(first_neg_ind_2)) == 1)) = nan_case_integral_scale(first_neg_ind_2(isnan(curr_input_2_integral_scale(first_neg_ind_2)) == 1));
    end
    input_2_integral_scale(first_neg_ind_2) = curr_input_2_integral_scale(first_neg_ind_2);
    
    if ((sum(sum(isnan(input_1_integral_scale))) == 0) && (sum(sum(isnan(input_2_integral_scale))) == 0))
        complete_flag = 1;
    end

%     disp(['lag_ind = ',num2str(lag_ind),' (out of ',num2str(N*lag_step),')'])
    
    n_lag_step = n_lag_step + 1;
    lag_ind = lag_ind + lag_step;

end

input_corr_integral_scale = min([input_1_integral_scale; input_2_integral_scale],[],1);
reg_dof_array = (delta_dim*sum((nan_mask_input_1 & nan_mask_input_2),1))./input_corr_integral_scale;

% integral_scale_1 = permute(reshape(input_1_integral_scale,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
% integral_scale_2 = permute(reshape(input_2_integral_scale,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
regression_dof_array_zero_lag = permute(reshape(reg_dof_array,[1 size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);


% compute pointwise correlations

regression_array = NaN([length(lags) prod_size_not_inc_corr_dim]);
regression_low_mag = NaN([length(lags) prod_size_not_inc_corr_dim]);
regression_high_mag = NaN([length(lags) prod_size_not_inc_corr_dim]);
mean_1 = NaN([length(lags) prod_size_not_inc_corr_dim]);
mean_2 = NaN([length(lags) prod_size_not_inc_corr_dim]);
std_dev_est_1 = NaN([length(lags) prod_size_not_inc_corr_dim]);
std_dev_est_2 = NaN([length(lags) prod_size_not_inc_corr_dim]);
for n_lag_ind = 1:length(lags)
    curr_lag = lags(n_lag_ind);
    lag_ind = round(curr_lag/delta_dim);
    input_1_corr_ind = (max([1 (1 - lag_ind)]):min([dim_length (dim_length - lag_ind)]))';
    input_2_corr_ind = (max([1 (1 + lag_ind)]):min([dim_length (dim_length + lag_ind)]))';
    
    curr_dof_1 = repmat((delta_dim*(length(input_1_corr_ind))),size(input_1_integral_scale))./input_1_integral_scale;
    curr_dof_2 = repmat((delta_dim*(length(input_2_corr_ind))),size(input_2_integral_scale))./input_2_integral_scale;
%     curr_error_perc_std_dev_1 = (erfinv(confidence_level)/(2^(1/2)))*((2*(curr_dof_1 - 1)).^(-1/2));
    curr_error_perc_std_dev_2 = (erfinv(confidence_level)/(2^(1/2)))*((2*(curr_dof_2 - 1)).^(-1/2));
    
    nan_mask_combined = nan_mask_input_1(input_1_corr_ind,:) & nan_mask_input_2(input_2_corr_ind,:);
    part_1_mean = sum(nan_mask_combined.*input_array_1_zeronans(input_1_corr_ind,:),1)./sum(nan_mask_combined,1);
    part_2_mean = sum(nan_mask_combined.*input_array_2_zeronans(input_2_corr_ind,:),1)./sum(nan_mask_combined,1);
    curr_corrcoeff = sum(nan_mask_combined.*(input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).*(input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])),1)./(((sum(nan_mask_combined.*((input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).^2),1)).^(1/2)).*((sum(nan_mask_combined.*((input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])).^2),1)).^(1/2)));
    regression_array(n_lag_ind,:) = sum(nan_mask_combined.*(input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).*(input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])),1)./(sum(nan_mask_combined.*((input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).^2),1));
    
    curr_reg_dof_array = (delta_dim*sum(nan_mask_combined,1))./input_corr_integral_scale;
    
    stderror_Z = (curr_reg_dof_array - 3).^(-1/2);

    stderror_Z(curr_reg_dof_array < 3) = Inf;

    normalized_Z_value = 0.5*(log(1 + abs(curr_corrcoeff)) - log(1 - abs(curr_corrcoeff)));
    Z_confidence_lower_mag = normalized_Z_value - ((((stderror_Z*(erfinv(confidence_level))*(2^(1/2))).^2) + ((normalized_Z_value.*curr_error_perc_std_dev_2).^2)).^(1/2));
    r_confidence_lower_mag = 1 - (2./(exp(2*Z_confidence_lower_mag) + 1));
    Z_confidence_higher_mag = normalized_Z_value + ((((stderror_Z*(erfinv(confidence_level))*(2^(1/2))).^2) + ((normalized_Z_value.*curr_error_perc_std_dev_2).^2)).^(1/2));
    r_confidence_higher_mag = 1 - (2./(exp(2*Z_confidence_higher_mag) + 1));
    
    mean_1(n_lag_ind,:) = part_1_mean + (sum(nan_mask_combined.*trend_array_1(input_1_corr_ind,:),1)./sum(nan_mask_combined,1));
    mean_2(n_lag_ind,:) = part_2_mean + (sum(nan_mask_combined.*trend_array_2(input_2_corr_ind,:),1)./sum(nan_mask_combined,1));
    std_dev_est_1(n_lag_ind,:) = (((sum(nan_mask_combined.*((input_array_1_zeronans(input_1_corr_ind,:) - repmat(part_1_mean,[length(input_1_corr_ind) 1])).^2),1))./(sum(nan_mask_combined,1))).^(1/2)).*((curr_dof_1./(curr_dof_1 - 1)).^(1/2));
    std_dev_est_2(n_lag_ind,:) = (((sum(nan_mask_combined.*((input_array_2_zeronans(input_2_corr_ind,:) - repmat(part_2_mean,[length(input_2_corr_ind) 1])).^2),1))./(sum(nan_mask_combined,1))).^(1/2)).*((curr_dof_2./(curr_dof_2 - 1)).^(1/2));
    regression_low_mag(n_lag_ind,:) = (sign(regression_array(n_lag_ind,:))).*r_confidence_lower_mag.*((std_dev_est_2(n_lag_ind,:))./(std_dev_est_1(n_lag_ind,:)));
    regression_high_mag(n_lag_ind,:) = (sign(regression_array(n_lag_ind,:))).*r_confidence_higher_mag.*((std_dev_est_2(n_lag_ind,:))./(std_dev_est_1(n_lag_ind,:)));
    
end

regression_array_at_lags = permute(reshape(regression_array,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
regression_array_low_mag_bound = permute(reshape(regression_low_mag,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
regression_array_high_mag_bound = permute(reshape(regression_high_mag,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
mean_array_1 = permute(reshape(mean_1,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
mean_array_2 = permute(reshape(mean_2,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
std_dev_estimate_1 = permute(reshape(std_dev_est_1,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);
std_dev_estimate_2 = permute(reshape(std_dev_est_2,[length(lags) size_array([(1:1:(n_dim - 1)) ((n_dim + 1):1:length(size_array))])]),[(2:1:n_dim) 1 ((n_dim + 1):1:length(size_array))]);./base_functions/taylordiag.m0000700442630740010010000004443614174631404017224 0ustar  adelmanDomain Users% TAYLORDIAG Plot a Taylor Diagram
%
% [hp ht axl] = taylordiag(STDs,RMSs,CORs,['option',value])
%
% Plot a Taylor diagram from statistics of different series.
%
% INPUTS:
%	STDs: Standard deviations
%	RMSs: Centered Root Mean Square Difference 
%	CORs: Correlation
%
%	Each of these inputs are one dimensional with same length. First
%	indice corresponds to the reference serie for the diagram. For exemple
%	STDs(1) is the standard deviation of the reference serie and STDs(2:N) 
%	are the standard deviations of the other series.
% 
%	Note that by definition the following relation must be true for all series i:
%	RMSs(i) - sqrt(STDs(i).^2 + STDs(1)^2 - 2*STDs(i)*STDs(1).*CORs(i)) = 0
%	This relation is checked and if not verified an error message is sent. Please see
%	Taylor's JGR article for more informations about this.	
%	You can use the ALLSTATS function to avoid this to happen, I guess ;-). You can get
%	it somewhere from: http://codes.guillaumemaze.org/matlab
%
% OUTPUTS:
% 	hp: returns handles of plotted points
%	ht: returns handles of the text legend of points
%	axl: returns a structure of handles of axis labels
%
% LIST OF OPTIONS:
%	For an exhaustive list of options to customize your diagram, please call the function
%	without arguments:
%		>> taylordiag
%
% SHORT TUTORIAL (see taylordiag_test.m for more informations):
%	 An easy way to get compute inputs is to use the ALLSTATS function you can get from:
%	 	http://codes.guillaumemaze.org/matlab
%	 Let's say you gathered all the series you want to put in the Taylor diagram in a 
%	 single matrix BUOY(N,nt) with N the number of series and nt their (similar) length.
%	 If BUOY(1,:) is the serie of reference for the diagram:
%		 for iserie = 2 : size(BUOY,1)
%		    S = allstats(BUOY(1,:),BUOY(iserie,:));
%		    MYSTATS(iserie,:) = S(:,2); % We get stats versus reference
%		 end%for iserie
%		 MYSTATS(1,:) = S(:,1); % We assign reference stats to the first row
%	 Note that the ALLSTATS function can handle NaNs, so be careful to compute statistics
%	 with enough points !
%	 Then you're ready to simply run:
%		taylordiag(MYSTATS(:,2),MYSTATS(:,3),MYSTATS(:,4));
%	
% REF: 	K. Taylor 
%		Summarizing multiple aspects of model performance in a single diagram
%		Journal of Geophysical Research-Atmospheres, 2001, V106, D7.
%
% Rev. by Guillaume Maze on 2010-02-10: Help more helpful ! Options now displayed by call.
% Copyright (c) 2008 Guillaume Maze. 
% http://codes.guillaumemaze.org
% All rights reserved.

% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% 	* Redistributions of source code must retain the above copyright notice, this list of 
% 	conditions and the following disclaimer.
% 	* Redistributions in binary form must reproduce the above copyright notice, this list 
% 	of conditions and the following disclaimer in the documentation and/or other materials 
% 	provided with the distribution.
% 	* Neither the name of the Laboratoire de Physique des Oceans nor the names of its contributors may be used 
%	to endorse or promote products derived from this software without specific prior 
%	written permission.
%
% THIS SOFTWARE IS PROVIDED BY Guillaume Maze ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Guillaume Maze BE LIABLE FOR ANY 
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
% BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 
% STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%




function varargout = taylordiag(varargin)

%%
if nargin == 0
	disp_optionslist;
	return
else
	narg = nargin - 3;
	if mod(narg,2) ~=0 
		error('taylordiag.m : Wrong number of arguments')
	end
end

STDs = varargin{1};
RMSs = varargin{2};
CORs = varargin{3};

%% CHECK THE INPUT FIELDS:
apro = 100;
di   = fix(RMSs*apro)/apro - fix(sqrt(STDs.^2 + STDs(1)^2 - 2*STDs*STDs(1).*CORs)*apro)/apro;
if find(di~=0)
%	help taylordiag.m
	ii = find(di~=0);
	if length(ii) == length(di)
		error(sprintf('taylordiag.m : Something''s wrong with ALL the datas\nYou must have:\nRMSs - sqrt(STDs.^2 + STDs(1)^2 - 2*STDs*STDs(1).*CORs) = 0 !'))
	else
		error(sprintf('taylordiag.m : Something''s wrong with data indice(s): [%i]\nYou must have:\nRMSs - sqrt(STDs.^2 + STDs(1)^2 - 2*STDs*STDs(1).*CORs) = 0 !',ii))		
	end
end
		
%% IN POLAR COORDINATES:
rho   = STDs;
theta = real(acos(CORs));
dx    = rho(1);	% Observed STD

%%



%% BEGIN THE PLOT HERE TO GET AXIS VALUES:
hold off
cax = gca;
tc = get(cax,'xcolor');
%ls = get(cax,'gridlinestyle');
ls = '-'; % DEFINE HERE THE GRID STYLE
next = lower(get(cax,'NextPlot'));

%% LOAD CUSTOM OPTION OF AXE LIMIT:
nopt = narg/2; foundrmax = 0;
for iopt = 4 : 2 : narg+3
	optvalue = varargin{iopt+1};
	switch lower(varargin{iopt}), case 'limstd', rmax = optvalue; foundrmax=1; end
end

% make a radial grid
hold(cax,'on');
if foundrmax==0
	maxrho = max(abs(rho(:)));
else
	maxrho = rmax;
end
hhh = line([-maxrho -maxrho maxrho maxrho],[-maxrho maxrho maxrho -maxrho],'parent',cax);
set(cax,'dataaspectratio',[1 1 1],'plotboxaspectratiomode','auto')
v = [get(cax,'xlim') get(cax,'ylim')];
ticks = sum(get(cax,'ytick')>=0);
delete(hhh);


% check radial limits and ticks
rmin = 0; 
if foundrmax == 0;
	rmax = v(4); 
end
rticks = max(ticks-1,2);
if rticks > 5   % see if we can reduce the number
    if rem(rticks,2) == 0
        rticks = rticks/2;
    elseif rem(rticks,3) == 0
        rticks = rticks/3;
    end
end
rinc  = (rmax-rmin)/rticks;
tick  = (rmin+rinc):rinc:rmax;

%% LOAD DEFAULT PARAMETERS:
if find(CORs<0)
	Npan = 2; % double panel
else
	Npan = 1;
end
tickRMSangle  = 135;	
showlabelsRMS = 1;
showlabelsSTD = 1;
showlabelsCOR = 1;
colSTD = [0 0 0];
colRMS = [0 .6 0];
colCOR = [0 0 1];
tickCOR(1).val = [1 .99 .95 .9:-.1:0];
tickCOR(2).val = [1 .99 .95 .9:-.1:0 -.1:-.1:-.9 -.95 -.99 -1];
widthCOR = .8;
widthRMS = .8;
widthSTD = .8;
styleCOR = '-.';
styleRMS = '--';
styleSTD = ':';
titleRMS = 1;
titleCOR = 1;
titleSTD = 1;
tickRMS = tick; rincRMS = rinc;
tickSTD = tick; rincSTD = rinc;
markersize = 12;
markerstyle = 'o';
markercolor = [0 0 1];


%% LOAD CUSTOM OPTIONS:
nopt = narg/2;
for iopt = 4 : 2 : narg+3
	optname  = varargin{iopt};
	optvalue = varargin{iopt+1};
	switch lower(optname)
		
		case 'tickrms'
			tickRMS = sort(optvalue);
			rincRMS = (max(tickRMS)-min(tickRMS))/length(tickRMS);				
		case 'showlabelsrms'
			showlabelsRMS = optvalue;
		case 'tickrmsangle'
			tickRMSangle = optvalue;
		case 'colrms'
			colRMS = optvalue;
		case 'widthrms'
			widthRMS = optvalue;
		case 'stylerms'
			styleRMS = optvalue;
		case 'titlerms'
			titleRMS = optvalue;

		case 'tickstd'
			tickSTD = sort(optvalue);
			rincSTD = (max(tickSTD)-min(tickSTD))/length(tickSTD);	
		case 'showlabelsstd'
			showlabelsSTD = optvalue;
		case 'colstd'
			colstd = optvalue;
		case 'widthstd'
			widthSTD = optvalue;
		case 'stylestd'
			styleSTD = optvalue;
		case 'titlestd'
			titleSTD = optvalue;
		case 'npan'
			Npan = optvalue;	
			
		case 'tickcor'
			tickCOR(Npan).val = optvalue;
		case 'colcor'
			colCOR = optvalue;
		case 'widthcor'
			widthCOR = optvalue;
		case 'stylecor'
			styleCOR = optvalue;
		case 'titlecor'
			titleCOR = optvalue;
		case 'showlabelscor'
			showlabelsCOR = optvalue;
                case 'markersize'
                        markersize = optvalue;
                case 'markerstyle'
                        markerstyle = optvalue;
                case 'markercolor'
                        markercolor = optvalue;
	end
end
		
	
%% CONTINUE THE PLOT WITH UPDATED OPTIONS:
		
% define a circle
    th = 0:pi/150:2*pi;
    xunit = cos(th);
    yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
    inds = 1:(length(th)-1)/4:length(th);
    xunit(inds(2:2:4)) = zeros(2,1);
    yunit(inds(1:2:5)) = zeros(3,1);
% plot background if necessary
   if ~ischar(get(cax,'color')),
%		ig = find(th>=0 & th<=pi);
		ig = 1:length(th);
      patch('xdata',xunit(ig)*rmax,'ydata',yunit(ig)*rmax, ...
            'edgecolor',tc,'facecolor',get(cax,'color'),...
            'handlevisibility','off','parent',cax);
   end

% DRAW RMS CIRCLES:
	% ANGLE OF THE TICK LABELS
    c82 = cos(tickRMSangle*pi/180); 
	s82 = sin(tickRMSangle*pi/180); 
    for ic = 1 : length(tickRMS)
		i = tickRMS(ic);
		iphic = find( sqrt(dx^2+rmax^2-2*dx*rmax*xunit) >= i ,1);
		ig = find(i*cos(th)+dx <= rmax*cos(th(iphic)));
		hhh = line(xunit(ig)*i+dx,yunit(ig)*i,'linestyle',styleRMS,'color',colRMS,'linewidth',widthRMS,...
                   'handlevisibility','off','parent',cax);	
		if showlabelsRMS
	        text((i+rincRMS/20)*c82+dx,(i+rincRMS/20)*s82, ...
	            ['  ' num2str(i)],'verticalalignment','bottom',...
	            'handlevisibility','off','parent',cax,'color',colRMS,'rotation',tickRMSangle-90)
		end
    end
	
% DRAW DIFFERENTLY THE CIRCLE CORRESPONDING TO THE OBSERVED VALUE
%      hhh = line((cos(th)*dx),sin(th)*dx,'linestyle','--','color',colSTD,'linewidth',1,...
%                   'handlevisibility','off','parent',cax);


% DRAW STD CIRCLES:
	% draw radial circles
    for ic = 1 : length(tickSTD)
		i = tickSTD(ic);
        hhh = line(xunit*i,yunit*i,'linestyle',styleSTD,'color',colSTD,'linewidth',widthSTD,...
                   'handlevisibility','off','parent',cax);
		if showlabelsSTD
			if Npan == 2
				if length(find(tickSTD==0)) == 0
					text(0,-rinc/20,'0','verticalalignment','top','horizontalAlignment','center',...
		            'handlevisibility','off','parent',cax,'color',colSTD);
				end
		        text(i,-rinc/20, ...
		             num2str(i),'verticalalignment','top','horizontalAlignment','center',...
		            'handlevisibility','off','parent',cax,'color',colSTD)
			else
				if length(find(tickSTD==0)) == 0
					text(-rinc/20,rinc/20,'0','verticalalignment','middle','horizontalAlignment','right',...
		            'handlevisibility','off','parent',cax,'color',colSTD);
				end
		        text(-rinc/20,i, ...
		             num2str(i),'verticalalignment','middle','horizontalAlignment','right',...
		            'handlevisibility','off','parent',cax,'color',colSTD)
			end
		end
    end
    set(hhh,'linestyle','-') % Make outer circle solid

% DRAW CORRELATIONS LINES EMANATING FROM THE ORIGIN:
	corr = tickCOR(Npan).val;
	th  = acos(corr);
    cst = cos(th); snt = sin(th);
    cs = [-cst; cst];
    sn = [-snt; snt];
    line(rmax*cs,rmax*sn,'linestyle',styleCOR,'color',colCOR,'linewidth',widthCOR,...
         'handlevisibility','off','parent',cax)

	% annotate them in correlation coef
	if showlabelsCOR
		rt = 1.05*rmax;		
	    for i = 1:length(corr)
	        text(rt*cst(i),rt*snt(i),num2str(corr(i)),...
	             'horizontalalignment','center',...
	             'handlevisibility','off','parent',cax,'color',colCOR);
	        if i == length(corr)
	            loc = int2str(0);
				loc = '1';
	        else
	            loc = int2str(180+i*30);
				loc = '-1';
	        end
	    end
	end

% AXIS TITLES
	axlabweight = 'bold';
	ix = 0;
	if Npan == 1
		if titleSTD
			ix = ix + 1;
			ax(ix).handle = ylabel('Standard deviation','color',colSTD,'fontweight',axlabweight);
		end
		
		if titleCOR
			ix = ix + 1;
			clear ttt
			pos1 = 45;	DA = 6;
			lab = 'Correlation';
% 			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.15*rmax;	ii = 0;
% 			for ic = 1 : length(c)
% 				ith = c(ic);
% 				ii = ii + 1;
% 				ttt(ii)=text(dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
% 				set(ttt(ii),'rotation',ith-90,'color',colCOR,'horizontalalignment','center',...
% 					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
% 			end
                        ttt = text(dd*cos((pos1+DA)*pi/180),dd*sin((pos1+DA)*pi/180),lab);
                        set(ttt,'rotation',pos1 - 90,'color',colCOR,'fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight)
			ax(ix).handle = ttt;
		end
		
		if titleRMS
			ix = ix + 1;
			clear ttt
			pos1 = tickRMSangle+(180-tickRMSangle)/2; DA = 5; pos1 = 85;
			lab = 'RMSD';
% 			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.05*tickRMS(1); 
			dd = .95*tickRMS(2);
% 			ii = 0;
% 			for ic = 1 : length(c)
% 				ith = c(ic);
% 				ii = ii + 1;
% 				ttt(ii)=text(dx+dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
% 				set(ttt(ii),'rotation',ith-90,'color',colRMS,'horizontalalignment','center',...
% 					'verticalalignment','top','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
% 			end
			ttt = text(dx+dd*cos((pos1+DA)*pi/180),dd*sin((pos1+DA)*pi/180),lab);
                        set(ttt,'rotation',pos1 - 90,'color',colRMS,'fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight)
			ax(ix).handle = ttt;
		end
		
		
	else
		if titleSTD
			ix = ix + 1;
			ax(ix).handle =xlabel('Standard deviation','fontweight',axlabweight,'color',colSTD);
		end
		
		if titleCOR
			ix = ix + 1;
			clear ttt
			pos1 = 90;	DA = 5;
			lab = 'Correlation';
			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.1*rmax;	ii = 0;
			for ic = 1 : length(c)
				ith = c(ic);
				ii = ii + 1;
				ttt(ii)=text(dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
				set(ttt(ii),'rotation',ith-90,'color',colCOR,'horizontalalignment','center',...
					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
			end
			ax(ix).handle = ttt;
		end
		
		if titleRMS
			ix = ix + 1;
			clear ttt
			pos1 = 160; DA = 10;
			lab = 'RMSD';
			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.05*tickRMS(1); ii = 0;
			for ic = 1 : length(c)
				ith = c(ic);
				ii = ii + 1;
				ttt(ii)=text(dx+dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
				set(ttt(ii),'rotation',ith-90,'color',colRMS,'horizontalalignment','center',...
					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
			end
			ax(ix).handle = ttt;
		end
	end
	

% VARIOUS ADJUSTMENTS TO THE PLOT:
	set(cax,'dataaspectratio',[1 1 1]), axis(cax,'off'); set(cax,'NextPlot',next);
	set(get(cax,'xlabel'),'visible','on')
	set(get(cax,'ylabel'),'visible','on')
	%    makemcode('RegisterHandle',cax,'IgnoreHandle',q,'FunctionName','polar');
	% set view to 2-D
	view(cax,2);
	% set axis limits
	if Npan == 2
	    axis(cax,rmax*[-1.15 1.15 0 1.15]);
		line([-rmax rmax],[0 0],'color',tc,'linewidth',1.2);
		line([0 0],[0 rmax],'color',tc);
	else
	    axis(cax,rmax*[0 1.15 0 1.15]);
%	    axis(cax,rmax*[-1 1 -1.15 1.15]);
		line([0 rmax],[0 0],'color',tc,'linewidth',1.2);
		line([0 0],[0 rmax],'color',tc,'linewidth',2);
	end


% FINALY PLOT THE POINTS:
	hold on
%    ALPHABET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
%    ALPHABET=sprintf('%s%s',ALPHABET,lower(ALPHABET));
%    ALPHABET = (1:length(STDs));
%    MarkerSet = 'so^d*+hvp><xvso^d*+hvp><xvso^d*+hvp><xv';
%    ct = repmat(get(gca,'ColorOrder'),ceil(length(STDs)/8),1);
    MarkerSet = repmat(markerstyle,[1 length(STDs)]);
    ct = repmat(markercolor,[length(STDs) 1]);
    for ii = 1 : length(STDs)
%         		pp(ii)=polar(theta(ii),rho(ii));
        pp(ii)=plot(rho(ii)*cos(theta(ii)),rho(ii)*sin(theta(ii)),...
                      'marker',MarkerSet(ii),'markersize',markersize,...
                      'markeredgecolor',ct(ii,:),'markerfacecolor',ct(ii,:),...
                      'LineStyle','none',...
                      'linewidth',2);
        % 		if length(STDs)<=26
        % 			tt(ii)=text(rho(ii)*cos(theta(ii)),rho(ii)*sin(theta(ii)),...
        %                num2str(ii),...
        %                'color','b',...
        %                'fontsize',14);
        % 		elseif length(STDs)<=26*2
        % 			tt(ii)=text(rho(ii)*cos(theta(ii)),rho(ii)*sin(theta(ii)),lower(ALPHABET(ii)),'color','r');
        % 		else
        % 			error('sorry I don''t how to handle more than 52 points labels !');
        % 		end
    end
    
% 	set(tt,'verticalalignment','top','horizontalalignment','left')
%	set(tt,'fontsize',14)
	

%%% OUTPUT
switch nargout
	case 1
		varargout(1) = {pp};
	case 2
		varargout(1) = {pp};
%		varargout(2) = {tt};
		varargout(2) = {ax};		
% 	case 3
% 		varargout(1) = {pp};
% 		varargout(2) = {tt};
% 		varargout(3) = {ax};		
end


end%function


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = disp_optionslist(varargin)

disp('General options:')
dispopt('''Npan''',sprintf('1 or 2: Panels to display (1 for positive correlations, 2 for positive and negative correlations).\n\t\tDefault value depends on CORs'));

disp('RMS axis options:')
dispopt('''tickRMS''','RMS values to plot gridding circles from observation point');
dispopt('''colRMS''','RMS grid and tick labels color. Default: green');
dispopt('''showlabelsRMS''','0 / 1 (default): Show or not the RMS tick labels');
dispopt('''tickRMSangle''','Angle for RMS tick lables with the observation point. Default: 135 deg.');
dispopt('''styleRMS''','Linestyle of the RMS grid');
dispopt('''widthRMS''','Line width of the RMS grid');
dispopt('''titleRMS''','0 / 1 (default): Show RMSD axis title');

disp('STD axis options:')
dispopt('''tickSTD''','STD values to plot gridding circles from origin');
dispopt('''colSTD''','STD grid and tick labels color. Default: black');
dispopt('''showlabelsSTD''','0 / 1 (default): Show or not the STD tick labels');
dispopt('''styleSTD''','Linestyle of the STD grid');
dispopt('''widthSTD''','Line width of the STD grid');
dispopt('''titleSTD''','0 / 1 (default): Show STD axis title');
dispopt('''limSTD''','Max of the STD axis (radius of the largest circle)');

disp('CORRELATION axis options:')
dispopt('''tickCOR''','CORRELATON grid values');
dispopt('''colCOR''','CORRELATION grid color. Default: blue');
dispopt('''showlabelsCOR''','0 / 1 (default): Show or not the CORRELATION tick labels');
dispopt('''styleCOR''','Linestyle of the COR grid');
dispopt('''widthCOR''','Line width of the COR grid');
dispopt('''titleCOR''','0 / 1 (default): Show CORRELATION axis title');

disp('Marker display options:')
dispopt('''markersize''','Size of marker');
dispopt('''markerstyle''','Style of marker (indicated by character)');
dispopt('''markercolor''','Color of marker (indicated by letter or 3-element RGB vector');

end%function

function [] = dispopt(optname,optval)
	disp(sprintf('\t%s',optname));
	disp(sprintf('\t\t%s',optval));
end
./base_functions/taylordiag_change.m0000700442630740010010000004263614203004460020516 0ustar  adelmanDomain Users% TAYLORDIAG Plot a Taylor Diagram with arrows showing change
%
% [hp ht axl] = taylordiag_change(STDs,RMSs,CORs,['option',value])
%
% Plot a Taylor diagram from statistics of different series.
%
% INPUTS:
%	STDs: Standard deviations
%	RMSs: Centered Root Mean Square Difference 
%	CORs: Correlation
%
%	Each of these inputs have dimensions Nx2; first column corresponds to the origin of the arrows and second column to where they point. First
%	row indicates the reference series for the diagram. For exemple
%	STDs(1,:) is the standard deviation of the reference serie and STDs(2:N,:) 
%	are the standard deviations of the two other series.
% 
% OUTPUTS:
% 	hp: returns handles of plotted points
%	ht: returns handles of the text legend of points
%	axl: returns a structure of handles of axis labels
%
% LIST OF OPTIONS:
%	For an exhaustive list of options to customize your diagram, please call the function
%	without arguments:
%		>> taylordiag
%
%	
% REF: 	K. Taylor 
%		Summarizing multiple aspects of model performance in a single diagram
%		Journal of Geophysical Research-Atmospheres, 2001, V106, D7.
%
% Adapted from rev. by Guillaume Maze on 2010-02-10: Help more helpful ! Options now displayed by call.
% Copyright (c) 2008 Guillaume Maze. 
% http://codes.guillaumemaze.org
% All rights reserved.

% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% 	* Redistributions of source code must retain the above copyright notice, this list of 
% 	conditions and the following disclaimer.
% 	* Redistributions in binary form must reproduce the above copyright notice, this list 
% 	of conditions and the following disclaimer in the documentation and/or other materials 
% 	provided with the distribution.
% 	* Neither the name of the Laboratoire de Physique des Oceans nor the names of its contributors may be used 
%	to endorse or promote products derived from this software without specific prior 
%	written permission.
%
% THIS SOFTWARE IS PROVIDED BY Guillaume Maze ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Guillaume Maze BE LIABLE FOR ANY 
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
% BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 
% STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%




function varargout = taylordiag_change(varargin)

%%
if nargin == 0
	disp_optionslist;
	return
else
	narg = nargin - 3;
	if mod(narg,2) ~=0 
		error('taylordiag.m : Wrong number of arguments')
	end
end

STDs = varargin{1};
RMSs = varargin{2};
CORs = varargin{3};

%% CHECK THE INPUT FIELDS:
apro = 100;
di   = fix(RMSs*apro)/apro - fix(sqrt(STDs.^2 + repmat(STDs(1,:),[size(STDs,1) 1]).^2 - 2*STDs.*repmat(STDs(1,:),[size(STDs,1) 1]).*CORs)*apro)/apro;
if find(di~=0)
%	help taylordiag.m
	ii = find(di~=0);
	if length(ii) == length(di)
		error(sprintf('taylordiag.m : Something''s wrong with ALL the datas\nYou must have:\nRMSs - sqrt(STDs.^2 + STDs(1)^2 - 2*STDs*STDs(1).*CORs) = 0 !'))
	else
		error(sprintf('taylordiag.m : Something''s wrong with data indice(s): [%i]\nYou must have:\nRMSs - sqrt(STDs.^2 + STDs(1)^2 - 2*STDs*STDs(1).*CORs) = 0 !',ii))		
	end
end
		
%% IN POLAR COORDINATES:
rho   = STDs;
theta = real(acos(CORs));
dx    = max(rho(1,:));	% Observed STD

%%



%% BEGIN THE PLOT HERE TO GET AXIS VALUES:
hold off
cax = gca;
tc = get(cax,'xcolor');
%ls = get(cax,'gridlinestyle');
ls = '-'; % DEFINE HERE THE GRID STYLE
next = lower(get(cax,'NextPlot'));

%% LOAD CUSTOM OPTION OF AXE LIMIT:
nopt = narg/2; foundrmax = 0;
for iopt = 4 : 2 : narg+3
	optvalue = varargin{iopt+1};
	switch lower(varargin{iopt}), case 'limstd', rmax = optvalue; foundrmax=1; end
end

% make a radial grid
hold(cax,'on');
if foundrmax==0
	maxrho = max(max(abs(rho)));
else
	maxrho = rmax;
end
hhh = line([-maxrho -maxrho maxrho maxrho],[-maxrho maxrho maxrho -maxrho],'parent',cax);
set(cax,'dataaspectratio',[1 1 1],'plotboxaspectratiomode','auto')
v = [get(cax,'xlim') get(cax,'ylim')];
ticks = sum(get(cax,'ytick')>=0);
delete(hhh);


% check radial limits and ticks
rmin = 0; 
if foundrmax == 0;
	rmax = v(4); 
end
rticks = max(ticks-1,2);
if rticks > 5   % see if we can reduce the number
    if rem(rticks,2) == 0
        rticks = rticks/2;
    elseif rem(rticks,3) == 0
        rticks = rticks/3;
    end
end
rinc  = (rmax-rmin)/rticks;
tick  = (rmin+rinc):rinc:rmax;

%% LOAD DEFAULT PARAMETERS:
if find(CORs<0)
	Npan = 2; % double panel
else
	Npan = 1;
end
tickRMSangle  = 135;	
showlabelsRMS = 1;
showlabelsSTD = 1;
showlabelsCOR = 1;
colSTD = [0 0 0];
colRMS = [0 .6 0];
colCOR = [0 0 1];
tickCOR(1).val = [1 .99 .95 .9:-.1:0];
tickCOR(2).val = [1 .99 .95 .9:-.1:0 -.1:-.1:-.9 -.95 -.99 -1];
widthCOR = .8;
widthRMS = .8;
widthSTD = .8;
styleCOR = '-.';
styleRMS = '--';
styleSTD = ':';
titleRMS = 1;
titleCOR = 1;
titleSTD = 1;
tickRMS = tick; rincRMS = rinc;
tickSTD = tick; rincSTD = rinc;
linewidth = 2;
linecolor = [0 0 1];


%% LOAD CUSTOM OPTIONS:
nopt = narg/2;
for iopt = 4 : 2 : narg+3
	optname  = varargin{iopt};
	optvalue = varargin{iopt+1};
	switch lower(optname)
		
		case 'tickrms'
			tickRMS = sort(optvalue);
			rincRMS = (max(tickRMS)-min(tickRMS))/length(tickRMS);				
		case 'showlabelsrms'
			showlabelsRMS = optvalue;
		case 'tickrmsangle'
			tickRMSangle = optvalue;
		case 'colrms'
			colRMS = optvalue;
		case 'widthrms'
			widthRMS = optvalue;
		case 'stylerms'
			styleRMS = optvalue;
		case 'titlerms'
			titleRMS = optvalue;

		case 'tickstd'
			tickSTD = sort(optvalue);
			rincSTD = (max(tickSTD)-min(tickSTD))/length(tickSTD);	
		case 'showlabelsstd'
			showlabelsSTD = optvalue;
		case 'colstd'
			colstd = optvalue;
		case 'widthstd'
			widthSTD = optvalue;
		case 'stylestd'
			styleSTD = optvalue;
		case 'titlestd'
			titleSTD = optvalue;
		case 'npan'
			Npan = optvalue;	
			
		case 'tickcor'
			tickCOR(Npan).val = optvalue;
		case 'colcor'
			colCOR = optvalue;
		case 'widthcor'
			widthCOR = optvalue;
		case 'stylecor'
			styleCOR = optvalue;
		case 'titlecor'
			titleCOR = optvalue;
		case 'showlabelscor'
			showlabelsCOR = optvalue;
                case 'linewidth'
                        linewidth = optvalue;
                case 'linecolor'
                        linecolor = optvalue;
	end
end
		
	
%% CONTINUE THE PLOT WITH UPDATED OPTIONS:
		
% define a circle
    th = 0:pi/150:2*pi;
    xunit = cos(th);
    yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
    inds = 1:(length(th)-1)/4:length(th);
    xunit(inds(2:2:4)) = zeros(2,1);
    yunit(inds(1:2:5)) = zeros(3,1);
% plot background if necessary
   if ~ischar(get(cax,'color')),
%		ig = find(th>=0 & th<=pi);
		ig = 1:length(th);
      patch('xdata',xunit(ig)*rmax,'ydata',yunit(ig)*rmax, ...
            'edgecolor',tc,'facecolor',get(cax,'color'),...
            'handlevisibility','off','parent',cax);
   end

% DRAW RMS CIRCLES:
	% ANGLE OF THE TICK LABELS
    c82 = cos(tickRMSangle*pi/180); 
	s82 = sin(tickRMSangle*pi/180); 
    for ic = 1 : length(tickRMS)
		i = tickRMS(ic);
		iphic = find( sqrt(dx^2+rmax^2-2*dx*rmax*xunit) >= i ,1);
		ig = find(i*cos(th)+dx <= rmax*cos(th(iphic)));
		hhh = line(xunit(ig)*i+dx,yunit(ig)*i,'linestyle',styleRMS,'color',colRMS,'linewidth',widthRMS,...
                   'handlevisibility','off','parent',cax);	
		if showlabelsRMS
	        text((i+rincRMS/20)*c82+dx,(i+rincRMS/20)*s82, ...
	            ['  ' num2str(i)],'verticalalignment','bottom',...
	            'handlevisibility','off','parent',cax,'color',colRMS,'rotation',tickRMSangle-90)
		end
    end
	
% DRAW DIFFERENTLY THE CIRCLE CORRESPONDING TO THE OBSERVED VALUE
%      hhh = line((cos(th)*dx),sin(th)*dx,'linestyle','--','color',colSTD,'linewidth',1,...
%                   'handlevisibility','off','parent',cax);


% DRAW STD CIRCLES:
	% draw radial circles
    for ic = 1 : length(tickSTD)
		i = tickSTD(ic);
        hhh = line(xunit*i,yunit*i,'linestyle',styleSTD,'color',colSTD,'linewidth',widthSTD,...
                   'handlevisibility','off','parent',cax);
		if showlabelsSTD
			if Npan == 2
				if length(find(tickSTD==0)) == 0
					text(0,-rinc/20,'0','verticalalignment','top','horizontalAlignment','center',...
		            'handlevisibility','off','parent',cax,'color',colSTD);
				end
		        text(i,-rinc/20, ...
		             num2str(i),'verticalalignment','top','horizontalAlignment','center',...
		            'handlevisibility','off','parent',cax,'color',colSTD)
			else
				if length(find(tickSTD==0)) == 0
					text(-rinc/20,rinc/20,'0','verticalalignment','middle','horizontalAlignment','right',...
		            'handlevisibility','off','parent',cax,'color',colSTD);
				end
		        text(-rinc/20,i, ...
		             num2str(i),'verticalalignment','middle','horizontalAlignment','right',...
		            'handlevisibility','off','parent',cax,'color',colSTD)
			end
		end
    end
    set(hhh,'linestyle','-') % Make outer circle solid

% DRAW CORRELATIONS LINES EMANATING FROM THE ORIGIN:
	corr = tickCOR(Npan).val;
	th  = acos(corr);
    cst = cos(th); snt = sin(th);
    cs = [-cst; cst];
    sn = [-snt; snt];
    line(rmax*cs,rmax*sn,'linestyle',styleCOR,'color',colCOR,'linewidth',widthCOR,...
         'handlevisibility','off','parent',cax)

	% annotate them in correlation coef
	if showlabelsCOR
		rt = 1.05*rmax;		
	    for i = 1:length(corr)
	        text(rt*cst(i),rt*snt(i),num2str(corr(i)),...
	             'horizontalalignment','center',...
	             'handlevisibility','off','parent',cax,'color',colCOR);
	        if i == length(corr)
	            loc = int2str(0);
				loc = '1';
	        else
	            loc = int2str(180+i*30);
				loc = '-1';
	        end
	    end
	end

% AXIS TITLES
	axlabweight = 'bold';
	ix = 0;
	if Npan == 1
		if titleSTD
			ix = ix + 1;
			ax(ix).handle = ylabel('Standard deviation','color',colSTD,'fontweight',axlabweight);
		end
		
		if titleCOR
			ix = ix + 1;
			clear ttt
			pos1 = 45;	DA = 6;
			lab = 'Correlation';
% 			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.15*rmax;	ii = 0;
% 			for ic = 1 : length(c)
% 				ith = c(ic);
% 				ii = ii + 1;
% 				ttt(ii)=text(dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
% 				set(ttt(ii),'rotation',ith-90,'color',colCOR,'horizontalalignment','center',...
% 					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
% 			end
                        ttt = text(dd*cos((pos1+DA)*pi/180),dd*sin((pos1+DA)*pi/180),lab);
                        set(ttt,'rotation',pos1 - 90,'color',colCOR,'fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight)
			ax(ix).handle = ttt;
		end
		
		if titleRMS
			ix = ix + 1;
			clear ttt
			pos1 = tickRMSangle+(180-tickRMSangle)/2; DA = 5; pos1 = 85;
			lab = 'RMSD';
% 			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.05*tickRMS(1); 
			dd = .95*tickRMS(2);
% 			ii = 0;
% 			for ic = 1 : length(c)
% 				ith = c(ic);
% 				ii = ii + 1;
% 				ttt(ii)=text(dx+dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
% 				set(ttt(ii),'rotation',ith-90,'color',colRMS,'horizontalalignment','center',...
% 					'verticalalignment','top','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
% 			end
			ttt=text(dx+dd*cos((pos1+DA)*pi/180),dd*sin((pos1+DA)*pi/180),lab);
                        set(ttt,'rotation',pos1 - 90,'color',colRMS,'fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight)
			ax(ix).handle = ttt;
		end
		
		
	else
		if titleSTD
			ix = ix + 1;
			ax(ix).handle =xlabel('Standard deviation','fontweight',axlabweight,'color',colSTD);
		end
		
		if titleCOR
			ix = ix + 1;
			clear ttt
			pos1 = 90;	DA = 5;
			lab = 'Correlation';
			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.1*rmax;	ii = 0;
			for ic = 1 : length(c)
				ith = c(ic);
				ii = ii + 1;
				ttt(ii)=text(dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
				set(ttt(ii),'rotation',ith-90,'color',colCOR,'horizontalalignment','center',...
					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
			end
			ax(ix).handle = ttt;
		end
		
		if titleRMS
			ix = ix + 1;
			clear ttt
			pos1 = 160; DA = 10;
			lab = 'RMSD';
			c = fliplr(linspace(pos1-DA,pos1+DA,length(lab)));
			dd = 1.05*tickRMS(1); ii = 0;
			for ic = 1 : length(c)
				ith = c(ic);
				ii = ii + 1;
				ttt(ii)=text(dx+dd*cos(ith*pi/180),dd*sin(ith*pi/180),lab(ii));
				set(ttt(ii),'rotation',ith-90,'color',colRMS,'horizontalalignment','center',...
					'verticalalignment','bottom','fontsize',get(ax(1).handle,'fontsize'),'fontweight',axlabweight);
			end
			ax(ix).handle = ttt;
		end
	end
	

% VARIOUS ADJUSTMENTS TO THE PLOT:
	set(cax,'dataaspectratio',[1 1 1]), axis(cax,'off'); set(cax,'NextPlot',next);
	set(get(cax,'xlabel'),'visible','on')
	set(get(cax,'ylabel'),'visible','on')
	%    makemcode('RegisterHandle',cax,'IgnoreHandle',q,'FunctionName','polar');
	% set view to 2-D
	view(cax,2);
	% set axis limits
	if Npan == 2
	    axis(cax,rmax*[-1.15 1.15 0 1.15]);
		line([-rmax rmax],[0 0],'color',tc,'linewidth',1.2);
		line([0 0],[0 rmax],'color',tc);
	else
	    axis(cax,rmax*[0 1.15 0 1.15]);
%	    axis(cax,rmax*[-1 1 -1.15 1.15]);
		line([0 rmax],[0 0],'color',tc,'linewidth',1.2);
		line([0 0],[0 rmax],'color',tc,'linewidth',2);
	end


% FINALY PLOT THE POINTS:
	hold on
%    ALPHABET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
%    ALPHABET=sprintf('%s%s',ALPHABET,lower(ALPHABET));
%    ALPHABET = (1:length(STDs));
%    MarkerSet = 'so^d*+hvp><xvso^d*+hvp><xvso^d*+hvp><xv';
%    ct = repmat(get(gca,'ColorOrder'),ceil(length(STDs)/8),1);
    ct = repmat(linecolor,[size(STDs,1) 1]);
    for ii = 1:size(STDs,1)
%         		pp(ii)=polar(theta(ii),rho(ii));
        pp((3*ii) - 2) = line(rho(ii,:).*cos(theta(ii,:)),rho(ii,:).*sin(theta(ii,:)),'LineStyle','-','linewidth',linewidth,'Color',linecolor);
        diff_pts = diff(rho(ii,:).*exp(1i*theta(ii,:)),1,2);
        currvec_length = abs(diff_pts);
        currvec_angle = imag(log(diff_pts));
        curr_endpts = (rho(ii,2)*exp(1i*theta(ii,2))) + [((currvec_length/3)*exp(1i*(currvec_angle + (5*pi/6)))) 0];
        pp((3*ii) - 1) = line(real(curr_endpts),imag(curr_endpts),'LineStyle','-','linewidth',linewidth,'Color',linecolor);
        curr_endpts = (rho(ii,2)*exp(1i*theta(ii,2))) + [((currvec_length/3)*exp(1i*(currvec_angle - (5*pi/6)))) 0];
        pp(3*ii) = line(real(curr_endpts),imag(curr_endpts),'LineStyle','-','linewidth',linewidth,'Color',linecolor);
        % 		if length(STDs)<=26
        % 			tt(ii)=text(rho(ii)*cos(theta(ii)),rho(ii)*sin(theta(ii)),...
        %                num2str(ii),...
        %                'color','b',...
        %                'fontsize',14);
        % 		elseif length(STDs)<=26*2
        % 			tt(ii)=text(rho(ii)*cos(theta(ii)),rho(ii)*sin(theta(ii)),lower(ALPHABET(ii)),'color','r');
        % 		else
        % 			error('sorry I don''t how to handle more than 52 points labels !');
        % 		end
    end
    
% 	set(tt,'verticalalignment','top','horizontalalignment','left')
%	set(tt,'fontsize',14)
	

%%% OUTPUT
switch nargout
	case 1
		varargout(1) = {pp};
	case 2
		varargout(1) = {pp};
%		varargout(2) = {tt};
		varargout(2) = {ax};		
% 	case 3
% 		varargout(1) = {pp};
% 		varargout(2) = {tt};
% 		varargout(3) = {ax};		
end


end%function


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = disp_optionslist(varargin)

disp('General options:')
dispopt('''Npan''',sprintf('1 or 2: Panels to display (1 for positive correlations, 2 for positive and negative correlations).\n\t\tDefault value depends on CORs'));

disp('RMS axis options:')
dispopt('''tickRMS''','RMS values to plot gridding circles from observation point');
dispopt('''colRMS''','RMS grid and tick labels color. Default: green');
dispopt('''showlabelsRMS''','0 / 1 (default): Show or not the RMS tick labels');
dispopt('''tickRMSangle''','Angle for RMS tick lables with the observation point. Default: 135 deg.');
dispopt('''styleRMS''','Linestyle of the RMS grid');
dispopt('''widthRMS''','Line width of the RMS grid');
dispopt('''titleRMS''','0 / 1 (default): Show RMSD axis title');

disp('STD axis options:')
dispopt('''tickSTD''','STD values to plot gridding circles from origin');
dispopt('''colSTD''','STD grid and tick labels color. Default: black');
dispopt('''showlabelsSTD''','0 / 1 (default): Show or not the STD tick labels');
dispopt('''styleSTD''','Linestyle of the STD grid');
dispopt('''widthSTD''','Line width of the STD grid');
dispopt('''titleSTD''','0 / 1 (default): Show STD axis title');
dispopt('''limSTD''','Max of the STD axis (radius of the largest circle)');

disp('CORRELATION axis options:')
dispopt('''tickCOR''','CORRELATON grid values');
dispopt('''colCOR''','CORRELATION grid color. Default: blue');
dispopt('''showlabelsCOR''','0 / 1 (default): Show or not the CORRELATION tick labels');
dispopt('''styleCOR''','Linestyle of the COR grid');
dispopt('''widthCOR''','Line width of the COR grid');
dispopt('''titleCOR''','0 / 1 (default): Show CORRELATION axis title');

disp('Line (arrow) display options:')
dispopt('''linewidth''','Width of line');
dispopt('''linecolor''','Color of line (indicated by letter or 3-element RGB vector');

end%function


function [] = dispopt(optname,optval)
	disp(sprintf('\t%s',optname));
	disp(sprintf('\t\t%s',optval));
end
./base_functions/taylor_stats_txy.m0000700442630740010010000000420614172647112020511 0ustar  adelmanDomain Usersfunction [STD,RMS,CORR] = taylor_stats_txy(ref,tf,w,norma)

% [STD,RMS,CORR] = taylor_stats_txy(reference_field,test_field,weights,norma)
% REF and TF    data(time,X)
% W             data weights W(X)
% NORMA:        1 = normalized to REF; 0 (default) = absolute
%
% weights: fraction of total (non-zero) area, such that sum(weights) = 1!
% assumes constant spacing in time! (i.e.: time weights are all the same!)
% centered RMS: (1) time-means of 'TF' and 'REF' will be subtracted (at each grid point)!
%               (2) spatial-mean of time-mean will be subtracted!    
%
% felix.w.landerer@jpl.nasa.gov [last modified: Nov-17, 2014]

defval('norma',0);
defval('w',ones(size(ref,2),size(ref,3)));

if ndims(ref) ~= ndims(tf)
    disp('Input fields REF and TF must have same dimensions!');
    return
elseif ndims(ref) == 3
    disp('collapsing dimensions of input fields to (time,space)');
    ref = ref(:,:);
    tf = tf(:,:);
elseif ndims(ref) <=1 || ndims(ref) > 3
    disp('check dimensions of input fields!')
end

w = w(:);
I = find(w); %only use non-zero-weight pixels

reft = zeros(size(ref,1),size(I,1));
tft = reft;
wt = reft;

% 
for t=1:size(ref,1);
    reft(t,:) = ref(t,I);
    tft(t,:) = tf(t,I);
    wt(t,:) = w(I);
end

% subtract time-mean at each data grid-point:
% NOTE:  - only do this if # of time samples > 1!
%        - if nt=1, proceed with removing spatial mean only
%          (i.e., compute stats for time-invariant field)
if size(tft,1) > 1
tft = tft - repmat(mean(tft,1),[size(tft,1),1]);
reft = reft - repmat(mean(reft,1),[size(reft,1),1]);
end

% subtract weighted spatial-mean from all data grid-points:
% NOTE:  - this step can be redundant in some configs, but do nonetheless
%          to make script work for all cases
tft = tft - mean(tft,1) * wt(1,:)';
reft = reft - mean(reft,1) * wt(1,:)';

% make one long data vector:
tft = tft(:);
reft = reft(:);
wt = wt(:);

ref_std = sqrt(1/sum(wt) * sum(reft.^2 .* wt));

STD = sqrt(1/sum(wt) * sum(tft.^2 .* wt));
RMS = sqrt(1/sum(wt) * sum((tft-reft).^2 .* wt));
CORR = 1/sum(wt) * sum(reft .* tft .* wt) / (ref_std * STD);

if norma
    STD = STD / ref_std;
    RMS = RMS / ref_std;
end

return