function [wimedian,wiresidual,wimad] = get_mad(obs_val,npmad,n1,minnmad)
%%
%% This program is written for calculating the median, residual, and the
%% median absolute deviation from obs and the TIEGCM values 
%% BA Emery, March 2011

wimedian = NaN;
wiresidual = NaN;
wimad = NaN;

clear wi_npmad;

nmad = 0;
for i=n1:n1+npmad-1 
% if (obs_val(i) ~= NaN)	% This structure fails!
 if ~isnan(obs_val(i))
  nmad = nmad + 1;
  wi_npmad(nmad) = obs_val(i);
 end
end
% Calculate median, residual, and mad
%if (nmad > 0)
% Want half the points before finding the median to avoid edges (esp in obs)
%if (nmad > npmad/2)
 if (nmad > minnmad)
 nmid = fix((nmad+1)/2);
%n1,nmad,wi_npmad
 wimedian = median(wi_npmad);
% y=mad(x,0 or 1) calculates the mean or median absolute deviations;  NaNs are treated as missing
% for a vector y, mad(x,1) = median(abs(x-median(x)))
 wiresidual = wi_npmad(nmid)-wimedian;
 wimad = mad(wi_npmad,1);
% Check assertion (passed!)
% wimad2 = median(abs(wi_npmad-wimedian));
% if (abs(wimad-wimad2) > 0.1)
%wimad,wimad2
% end
end


