function [moments,AZ,EL] = MpdrWxMoments(Filename,Pathname,lr_or_sr,numCipsPerNcip)
% parse a data recording file for MPDR radar and process the data

tic

radarFile = 'mpdr2PlusRadarConfig.yaml';
UDP=false;
plot = true;
myRadar = mpdrRadar(radarFile);
if exist('Filename',"var")
    if exist('Pathname',"var")
%         myRadar.openDataRecordFile(Filename,Pathname);
        myCip=cipCube.factory(myRadar, UDP, Filename, Pathname);
    else
%         myRadar.openDataRecordFile(Filename);
        myCip=cipCube.factory(myRadar, UDP, Filename);
    end
else
    myCip=cipCube.factory(myRadar, UDP);
end

% Set flag to process long range or short range waveform data
if ~exist('lr_or_sr',"var")
    lr_or_sr = 'LR';
%     dlgTitle    = 'Long Range / Short Range Selection';
%     dlgQuestion = 'Display LR or SR returns?';
%     lr_or_sr = questdlg(dlgQuestion,dlgTitle,'LR','SR', 'LR');
end

% Set number of coherent integrations to be averaged for each non-coherent integration
if ~exist('numCipsPerNcip',"var")
    numCipsPerNcip = 1;
end

% Specify maximum number of coherent integration periods to process (inf to
% process all of them)
% maxCips = numCipsPerNcip*2;
maxCips = inf; %numCipsPerNcip*2000;

% Initialize CIP object
%myCip=cipCube.factory(myRadar, UDP);

% Load the first CIP and read CIP data and meta-data
%[cipMetaData,eof] = myCip.loadCip(870); % set starting CIP if desired
%[cipMetaData,eof] = myCip.loadCip([],13.86432); % check starting time if desired
[cipMetaData,eof] = myCip.loadCip();
% [cipData,cipMetaData] = myCip.getCipCube();

fileTimeSec = cipMetaData.priTime_s*cipMetaData.fileSizePris;
numPulsesPerCip = cipMetaData.numPris;
numCips = cipMetaData.fileSizeCips;
numNcips = floor(numCips/numCipsPerNcip);
ncipTimeSec = cipMetaData.priTime_s*cipMetaData.numPris*numCipsPerNcip;
fps = 30; % MP4 video frames per second
%framesPerNcip = floor(ncipTimeSec*fps)
framesPerNcip_realtime = floor(ncipTimeSec*fps);
framesPerNcip_user = 2*fps; % alternative delay between Ncip updates during MP4 video playback
framesPerNcip = framesPerNcip_realtime;

fprintf('File contents: %1.1f sec, %1d cips (%1d pulses per cip), %1d ncips (%1d cips per ncip)\n',fileTimeSec,numCips,numPulsesPerCip,numNcips,numCipsPerNcip);
% pause(3)

%% Determine whether to process long range or short range waveform and select range cells
if strcmp(lr_or_sr,'LR')
    dispChan=uint32(radChan.sumLr);
    % In Wx mode, the LR Sum channel data has 4 duplicate values per valid
    % sample.  The valid samples are moved to the 1st quarter of the array
    % when read and the rest of the samples are not valid.
    wx_rng_cell_max = floor(size(myCip.dataCube,2));%/4
    %     endRange=cipMetaData.cells.detEnd(dispChan);
    %     endRange=ceil(8000/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
else
    dispChan=uint32(radChan.sumSr);
    wx_rng_cell_max = size(myCip.dataCube,2);
    % process only up to 2.2 km
    %    endRange=ceil(2200/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
end

% Get other range parameters from meta-data
rng_cell_zero = cipMetaData.cells.rangeZero(dispChan);
m_per_rng_cell = cipMetaData.rangeCellSize_m(dispChan);
OURadarStruct = radarStructFromCipMetaData(myRadar, dispChan);

%[rng_m_start,rng_m_end] = rng_cell_to_m(rng_cell_zero,m_per_rng_cell,rng_cell_start,rng_cell_end)
% [rng_m_start,rng_m_end] = rng_cell_to_m(rng_cell_zero,m_per_rng_cell,rng_cell_zero,wx_rng_cell_max);
% [rng_cell_start,rng_cell_end] = rng_m_to_cell(rng_cell_zero,m_per_rng_cell,0,17000);

mph_per_m_per_sec = 2.23694;
% pri_us = 4249/32e6*1e6; % 132.78125 us (not used, get from cipMetaData if needed)
% cip_sec = 0.01593375 sec
thresh_sidelobe_dB = -15; % Censor bins with power more than 15 dB below peak bin power
thresh_noise_dB = 4;
pwr_offset_dB = -37;
meteorology_mode = 0; % set to 1 to change sign of velocity data to follow meteorology convention (i.e. pos Doppler == neg velocity)
save_mode = 0; % set to one to save DC bias corrected data and pulsed compressed data in MAT files

cipNum = 1;
nBurst = 0;
start_scan = 1; %6; %30;
max_scan = 100;
numFinePsd = 0; % accumulated number of power spectral density estimates in NCIP
% sumCoarsePsd = zeros(size(myCip.dataCube,[2 3])); % dimensions are range x doppler bin
sumCoarsePsd = zeros(wx_rng_cell_max,120); % dimensions are range x doppler bin
cipStack = zeros(wx_rng_cell_max,myCip.cipMetaData.numPrisWx,numCipsPerNcip);
store_enable = 0; % set to 0 until want to save partial data set when save_mode = 1
moments = zeros(wx_rng_cell_max, numNcips, 3);
% fprintf('File contents: %1.1f sec, %1d cips (%1d pulses per cip), %1d ncips (%1d cips per ncip)\n',fileTimeSec,numCips,numPulsesPerCip,numNcips,numCipsPerNcip);

%% Loop through input file one CIP (coherent integration period) at a time
while ~eof && cipNum < max_scan+start_scan
    fprintf('Processing cip %1d\n',cipNum);
    
    % Perform DC bias correction
    myCip.correctDcBias();
    [cipDataDc,cipMetaData] = myCip.getCipCube();
    
    % Perform pulse compression
    myCip.compress();
    [cipDataPc,cipMetaData] = myCip.getCipCube();
    %limitrate;
    nBurst = nBurst + 1;
    cipStack(:,:,nBurst) = squeeze(cipDataPc(dispChan,:,:));
    % Perform Doppler filtering
%     myCip.dopplerFilter();
%     [cipDataDf,cipMetaData] = myCip.getCipCube();
    
    
    % Perform per NCIP (non-coherent CIP processing)
    if nBurst == numCipsPerNcip %87 for nominal Wx mode processing
        
        % Increment number of fine PSD estimates (Pxx) that are calculated
        % after each non-coherent integration period (NCIP) is complete
        numFinePsd = numFinePsd+1;
        
        % Flip to treat motion toward radar as negative radial velocity
        % (this is convention used in meteorology, not by Helios)
        if (meteorology_mode)
            finePsd = fliplr(finePsd)
        end
        feaux = squeeze(getRadarMomentsNCip(OURadarStruct, cipStack, true));
        
        moments(:,numFinePsd,:) = feaux;
        nBurst = 0; % this will reinit sumCoarsePsd to 0
        %pause(0.5)
        
    end % if cipNum >= start_scan, process cip data
    
    cipNum = cipNum + 1;
    % load the next CIP...
    [~,eof] = myCip.loadCip();
    
end % while ~eof && cipNum < max_scan+start_scan
toc

%% Plot additional results
if plot
    figure
    rng_axis=rng_cell_zero:wx_rng_cell_max;
   % Vr = squeeze(moments(:,:,2));
   %Vr = squeeze(moments(:,1,:)); % just one az slize, though?
   Vr = moments(:,1,1);
   
        xlim([290 326]);
    set(gca,'YDir','normal');
    colormap(cmaps('rgmap',32));
    hcb=colorbar;
    title('Weather Radial Velocity Map (Spectral Estimation)');
    ylabel(hcb,'Radial Velocity (mph)');
    %xlabel('Azimuth Sector')
    % xlabel('Non-coherent Processing Interval')
    xlabel('Radar Azimuth (deg relative to True North)')
    ylabel('Range (km)')
    %caxis([-15 15])
    caxis([-30 30]);
    h1 = imagesc([0:(numFinePsd-1)] - floor(numFinePsd/2),rng_axis*m_per_rng_cell,Vr);

    
end

%keyboard

%%------------------------------------------------------------------------------------------------------------
function [rng_m_start,rng_m_end] = rng_cell_to_m(rng_cell_zero,m_per_rng_cell,rng_cell_start,rng_cell_end)
rng_m_start = (rng_cell_start-rng_cell_zero)*m_per_rng_cell;
rng_m_end = (rng_cell_end-rng_cell_zero)*m_per_rng_cell;

function [rng_cell_start,rng_cell_end] = rng_m_to_cell(rng_cell_zero,m_per_rng_cell,rng_m_start,rng_m_end)
rng_cell_start = round(rng_m_start/m_per_rng_cell)+rng_cell_zero;
rng_cell_end = round(rng_m_end/m_per_rng_cell)+rng_cell_zero;

