function out = read_odm_fix_DSD_et(filename)

%%% reads rd raw files from ODM-470 disdrometer


%%%%%%%%%%%%%%%% drops file ... ?
% Date:10.09.2018 23:59:08
% 233 000 000 030 233 024 000 001 
% 233 000 000 039 233 032 000 004 
% 232 240 000 037 233 016 000 006 
% 232 244 000 030 233 024 000 003 
% 232 244 000 030 233 008 000 009 
% 232 252 000 030 233 024 000 000 
% 232 244 000 036 233 008 000 000 
% 232 192 000 037 232 240 000 000 
% 232 224 000 030 232 244 000 009 
% 232 248 005 074 233 012 000 001 
% 232 252 000 030 233 024 000 000 
% 232 248 001 224 233 028 000 000 
% 232 232 000 044 233 000 000 000 
% 232 244 000 032 233 008 000 003 
% 233 020 000 030 233 052 000 000 

%%%%%%%%%%%%%%%%%% rd file: README
%  date     time  volt  urel_anem nsn sbin nrn rbin
% 10092018 235708 4.99 01.78      0034 009 0034 009

%         BS occupied Bins of snow
%         DS number of Drops per bin of snow
%         TS Total snow drops, should equal sum(DR)
% %%% BS 013 014 015 016 017 018 019 020 034 
% %%% DS 008 004 002 006 003 003 006 001 001 
% %%% TS 000289 000062 000042 000134 000034 000044 000167 000018 000018 

%         BR occupied Bins of snow
%         DR number of Drops per bin of snow
%         TR Total snow drops, should equal sum(DR)
% %%% BR 013 014 015 016 017 018 019 020 034 
% %%% DR 008 004 002 006 003 003 006 001 001 
% %%% TR 000289 000062 000042 000134 000034 000044 000167 000018 000018 

%%%%%%%% Snow and rain rows look identical... so we could just ignore snow!


% read_odm_DSD
%
% read datafiles from the ODM470 disdrometer and return drop-size
% distribution info
%
% k.drushka // march 2016
%
% reads in the rd file 

%disp('------------------------------------------');
%disp('  processing DSD data for these times:    ');

% //// conversion between rain bin(bin) and dropsize and rain rate: -- see Klepp 2015
% in matlab, log(x) is ln(x) is the natural log.
binnum=1:128;
D_p=(exp(binnum/94*log(10))-1+exp((binnum+1)/94*log(10))-1)/2; % mm (eq. 1 in Klepp)

%%% matches what is in Christian's/Kyla's fortran code:
%%% notes... log(10) in fortran is the natural log.
% radius = D/2 =      rmm(i)=(exp((i/94.)*log(10.))-1.
%      *          +exp(((i+1)/94.)*log(10.))-1.)/4.
     
     
l=0.120 ; % used to be 120 mm; % length of catchment volume
w=0.022 ; % used to be d = 22 mm;  % diameter of catchment volume
dt=60;  % measurement interval, in s
% % % Mparticle=4/3*pi*1000*(D_p./200).^3;     % mass of drops, in kg (Klepp eq. 3)
Mparticle=4/3*pi*1000*(D_p./2000).^3;     % mass of drops, in kg (Klepp eq. 3 - modified! should be 2000 not 200 in the denominator, to convert from diam[mm] to radius[m])
% vfall= 9.65 - 10.3*exp(-1.2*D_p*10/2);  % terminal fall velocity , for D in cm, vfall in m/s (Klepp eq. 4)
vfall = 9.65-10.3*exp(-.6.*D_p);   %% Takes d in mm and calculates vfall in m/s
% // for snow ... not used //
% Mparticle_snow=0.0000107*(D_p).^(3.1);     % mass of snow drops, in kg (Klepp eq. 5)
% vfall_snow= 7.33.*D_p.^(0.78);  % terminal fall velocity of snow, in m/s (Klepp eq. 4)

% open file  
fid=fopen(filename,'r');

% read line by line
cnt=0;
while 1
    cnt=cnt+1;
    % line 1: header
    hdr=textscan(fid,'%s%s%f%f%u%u%u%u',1);
%     disp(hdr);
    if isempty(hdr{1})
        break
    end
    thist=datenum([cell2mat(hdr{1}) cell2mat(hdr{2})],'ddmmyyyyHHMMSS');
%     disp(datestr(thist));
    Vref=hdr{3};  % reference voltage
    urel_anem=hdr{4};  % wind speed
    Nsnowpart=hdr{5}; % # particles of snow
    Nsnowbins=hdr{6}; % # occupied bins (classes) of snow
    Nrainpart=hdr{7}; % # particles of rain
    Nrainbins=hdr{8}; % # occupied bins (classes) of rain
    % if not enough drops, skip this
%     Nmin=1;
%     if (Nrainpart+Nsnowpart)<Nmin
%         cnt=cnt-1;  % don't count this
%         for l=1:8
%             fgetl(fid);  % read the rest of this entry
%         end
%         continue
%     end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if Nsnowpart > 0
        % line 2: occupied classes of snow
        % --- skip. we only have rain.
        %thisline=fgetl(fid);
%         disp('found snow...');
        BS=textscan(fid,['%s' repmat('%u',1,Nsnowbins)],1);
        BS{1}=[];
        BS=double([BS{:}]);

        
        % line 3: number of drops per bin in snow
        % --- skip. we only have rain.
        %thisline=fgetl(fid);
        DS=textscan(fid,['%s' repmat('%u',1,Nsnowbins)],1);
        DS{1}=[];
        DS=double([DS{:}]);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % bin detected rain and snow onto full diameter bin array
        %[dum,dumi] = ismember(BS,binnum); % classi gives the index of the classes
        nsnow=zeros(size(binnum));
        %nsnow(dumi)=DS;  % this is the number of snow in each bin
        nsnow(BS) = DS;
        
        
        % check: sum of snowbins should equal Nsnowbins
        if sum(DS)~=Nsnowpart
            fprintf(1,'mismatch in rain particles! (cnt=%d)\n',cnt)
        end
        
        TS=textscan(fid,['%s' repmat('%u',1,Nsnowbins)],1);
        
    else
        % skip lines 2 and 3 if rain is not detected
        fgetl(fid);
        fgetl(fid);
        % skip line 4 anyway
        fgetl(fid);
        nsnow = zeros(1,length(binnum));
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if Nrainpart > 0
        % line 5: occupied classes of rain
%         disp('found rain...');
        BR=textscan(fid,['%s' repmat('%u',1,Nrainbins)],1);
        BR{1}=[];
        BR=double([BR{:}]);

        % line 6: number of drops per bin in rain
        DR=textscan(fid,['%s' repmat('%u',1,Nrainbins)],1);
        DR{1}=[];
        DR=double([DR{:}]);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % bin detected rain and snow onto full diameter bin array
        %[dum,dumi] = ismember(BR,binnum) % classi gives the index of the classes
        nrain=zeros(size(binnum));
        %nrain(dumi)=DR  % this is the number of drops in each bin
        nrain(BR) = DR;
    
        % check: sum of rainbins should equal Nrainbins
        if sum(DR)~=Nrainpart
            fprintf(1,'mismatch in rain particles! (cnt=%d)\n',cnt)
        end
        
        TR=textscan(fid,['%s' repmat('%u',1,Nrainbins)],1);

    else
        % skip lines 5 and 6 if rain is not detected
        fgetl(fid);
        fgetl(fid);
        % skip line 7 anyway
        fgetl(fid);
        nrain = zeros(1,length(binnum));
    end

    % line 8: empty
    fgetl(fid);
    
    % skip line 9 because ?
    fgetl(fid);
    
    %% precip rate and concentrations. 
    %%% snow is getting classified for some reason, we may never know why.
    %%% Use "all" drops as snow+rain. Because we know they are rain. 
    % psd = concentrations (count per m^3 volume) (eq. 2 in Klepp)
    % prate = precipitation rate in mm/hr
    psdall = (nrain+nsnow)./(l*w*dt*sqrt(urel_anem^2 + vfall.^2)); % count per m^3
    psdrain = (nrain)./(l*w*dt*sqrt(urel_anem^2 + vfall.^2)); % count per m^3
    psdsnow = (nsnow)./(l*w*dt*sqrt(urel_anem^2 + vfall.^2)); % count per m^3
    praterain = nansum(3600 * psdrain .* vfall .* Mparticle); % rain rate in mm/hr (Klepp eq. 6)
    pratesnow = nansum(3600 * psdsnow .* vfall .* Mparticle); % rain rate in mm/hr (Klepp eq. 6)
    prateall = nansum(3600 * psdall .* vfall .* Mparticle); % rain rate in mm/hr (Klepp eq. 6)
    
    

    
    %%% before, this had nrain... which was just count.
%     prate = nansum(3600 * nrain .* vfall .* Mparticle); 

    %% German fortran code... can't make much sense of this
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           rrzeit=rrzeit+uhist(iibin(k))*Mtr(iibin(k))*Vfall(iibin(k))
%      *             /60./(((0.022/2.))**2.*0.12*3.1415)
%      *             *3.6/100.
%      
%                    rzeit(ibin(k))=uhist(ibin(k))/iuspek(ibin(k))*
%      *             60.*0.0000456/0.022*sqrt(Vfall(ibin(k))**2+ff**2)
%             endif

% %%% in kyla's processing as of May 2020
% rrzeit=rrzeit+uhist(iibin(k))*Mtr(iibin(k))*Vfall(iibin(k))
% *             /60./(((0.022/2.))**2.*0.12*3.1415)
% *             *3.6/100.
            
%             rra=rra+3600.*n_a(nr(i))*Vfall(nr(i))*Mtr(nr(i))
%             nt0(nr(i))=r(nr(i))/(3600.*Vfall(nr(i))*Mtr(nr(i)))
%             lwca=lwca+(nint(rnanz(nr(i)))*Mtr(nr(i)))/
%      *           ((0.022-2.*dcm(nr(i))/100.)*0.12*Vel*mz)*1000.
%                  lwcu=lwcu+nt0(nr(i))*Mtr(nr(i))*1000.
% 
%      
%       dsd_refl(nr(i))= n_a(nr(i))*(dmm_mean(nr(i))**6)
% c     *dmm_interval(nr(i))
% c       dsd_refl(nr(i))=dsd_norm(nr(i))*(dmm_mean(nr(i))**6)
%       dsd_refl_sum=dsd_refl_sum+dsd_refl(nr(i))


    % store data
    out.t(cnt,1)                            =thist;    
    out.t(out.t==0)                         =nan;
    out.voltage(cnt,1)                      =Vref;
    out.urel_anem(cnt,1)                    =urel_anem;
    out.nsumrain(cnt,1)                     =nansum(nrain);
    out.nrain(cnt,1:length(binnum))         =nrain;
    out.nsumsnow(cnt,1)                     =nansum(nsnow);
    out.nsnow(cnt,1:length(binnum))         =nsnow;
    out.nsumall(cnt,1)                      =nansum(nrain+nsnow);
    out.nall(cnt,1:length(binnum))          =nrain+nsnow;
    out.psdrain(cnt,1:length(binnum))       =psdrain;
    out.psdsnow(cnt,1:length(binnum))       =psdsnow;
    out.psdall(cnt,1:length(binnum))        =psdall;
    out.praterain(cnt,1)                    =praterain;
    out.pratesnow(cnt,1)                    =pratesnow;
    out.prateall(cnt,1)                     =prateall;

    1;  %%% I don't know what this line does?
end 

out.Mparticle=Mparticle;
out.vfall=vfall;
out.binnum=binnum;
out.D_p=D_p;
out.nbins = length(binnum);
out.ntimes = length(out.t);
% disp('   done processing DSD data               ');
% disp('------------------------------------------');
