function bins = daily_rain_hist(myDate);
%
% calculates histograms of rain rates along each azimuth, but only
% when RATE_KDP exists.  Limits range from 17 km to 70 km.  Looks
% only at the first three elevations (SUR, only)
%

direc = '/export/d1/dynamo/data/qc_v1/cfradial/moments/merged/sur';
set(0,'DefaultTextInterpreter','none');

global rt_direc;
rt_direc = '/export/d1/dynamo/data/qc_v1/cfradial';

rates = {'RATE_HYBRID' 'RATE_KDP' 'RATE_KDP_ZDR' 'RATE_ZH' 'RATE_Z_ZDR'};
elevs = {'e1' 'e2' 'e3' 'e4'};

% limit range of interest (17km to 70km)
ming = 113;
maxg = 466;

bins = struct();
bins.e1 = struct();
bins.e2 = struct();
bins.e3 = struct();
bins.e4 = struct();

for nn=1:numel(rates);
   bins.e1.(rates{nn}) = zeros([200,360]);
   bins.e2.(rates{nn}) = zeros([200,360]);
   bins.e3.(rates{nn}) = zeros([200,360]);
   bins.e4.(rates{nn}) = zeros([200,360]);
end;

f_name = get_sorted_file_list(myDate,'moments','merged',{'sur'});

[ ATT ] = get_hdf5_param_atts( char(f_name(1,:)) );

klim = size(f_name,1);
for nf=1:klim;
   [ ATT ] = get_hdf5_param_atts( char(f_name(nf,:)) );
   SWP = get_and_scale_hdf5_data( char(deblank(f_name(nf,:))), ATT, ...
         [rates { 'time_coverage_end' 'time_coverage_start' 'time' 'range' ...
         'sweep_start_ray_index' 'sweep_end_ray_index' 'sweep_number' ...
         'fixed_angle' 'azimuth' 'elevation' 'antenna_transition'}]);
   fprintf('\n%5d  ', nf);

   % replace zero-degree az with 360
   A = find(SWP.azimuth == 0 );
   SWP.azimuth(A) = 360;

   % consider rain rates only where there is a value for RATE_KDP

   mask = SWP.RATE_KDP ./ SWP.RATE_KDP;  % uses properties of NaN
                                         % to create identity/NaN
                                         % mask

   el_lim = min([numel(elevs)  numel(SWP.sweep_end_ray_index)]);

   for nn=1:numel(rates);
       b4 = SWP.sweep_end_ray_index(el_lim);
       SWP.(rates{nn}) = mask .* SWP.(rates{nn});
       fprintf('%10.2f  ', max(max(SWP.(rates{nn})(ming:maxg,1:b4))));
   end;

   for nelev=1:el_lim;
       b1 = SWP.sweep_start_ray_index(nelev) + 1;
       b2 = SWP.sweep_end_ray_index(nelev) + 1;
       for nbeam=b1:b2;
          if( SWP.antenna_transition(nbeam) == 1 || SWP.elevation(nbeam) < 0.0)
          %        skip the transition beam
             continue; 
          end;
          az  = round(SWP.azimuth(nbeam));
          for nrate=1:numel(rates);
              for ngate=ming:maxg;  % limit range
                  if( isnan(SWP.(rates{nrate})(ngate,nbeam)))
                      continue;
                  end;
                  thisbin = min([ceil(SWP.(rates{nrate})(ngate,nbeam)) 200]);
                  %                  fprintf('%10.4f  %10d  %10d\n', SWP.(rates{nrate})(ngate,nbeam),thisbin,az);
                  bins.(elevs{nelev}).(rates{nrate})(thisbin,az) = ...
                      1 + bins.(elevs{nelev}).(rates{nrate})(thisbin,az);
              end;
          end;
       end;
   end;
end;
fprintf('\n');
