function [moments,AZ,EL] = magicMoments(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 5 km
    endRange=ceil(5000/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
end

% Get other range parameters from meta-data
%rng_cell_zero = 99; 
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 = inf;
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);

make_cf_file = true;

%% for cfradial file
latitude = 37.02045856274912;
longitude = 0;
altitude = 0;

if make_cf_file
    cf_filename = "test_cfr2.nc";
   
    h5create(cf_filename,'/sweep_0001/latitude',1); h5write(cf_filename,'/sweep_0001/latitude',37.02045856274912);
    h5create(cf_filename,'/sweep_0001/longitude',1); h5write(cf_filename,'/sweep_0001/longitude',-76.3366642324842);
    h5create(cf_filename,'/sweep_0001/altitude',1); h5write(cf_filename,'/sweep_0001/altitude',40);
    h5create(cf_filename,'/sweep_0001/time', numNcips); h5write(cf_filename,'/sweep_0001/time',linspace(1,numNcips,numNcips));
    h5create(cf_filename,'/sweep_0001/elevation', numNcips); h5write(cf_filename,'/sweep_0001/elevation',linspace(10,10,numNcips));
    h5create(cf_filename,'/sweep_0001/azimuth',numNcips); h5write(cf_filename,'/sweep_0001/azimuth',linspace(-55,55,numNcips));

    h5create(cf_filename,'/sweep_0001/sweep_mode',1 ,Datatype='string'); h5write(cf_filename,'/sweep_0001/sweep_mode',"azimuth_surveillance");

    h5create(cf_filename,'/sweep_0001/fixed_angle',1); h5write(cf_filename,'/sweep_0001/fixed_angle',10.0'); % just what we have for angle now; multiple scans in file = multiple angles.

    h5create(cf_filename,'/sweep_0001/sweep_start_ray_index',1); h5write(cf_filename,'/sweep_0001/sweep_start_ray_index',0); % just what we have for angle now; multiple scans in file = multiple angles.
    h5create(cf_filename,'/sweep_0001/sweep_end_ray_index',1); h5write(cf_filename,'/sweep_0001/sweep_end_ray_index',numNcips); % just what we have for angle now; multiple scans in file = multiple angles.
end

% range bin distances in meters.
%range = linspace(startrange, stoprange, stoprange)  
%az slices. Either min/max it and create it or just o with it

%% Loop through input file one CIP (coherent integration period) at a time
radialbin=1;
while ~eof && cipNum < max_scan+start_scan
    %figure(100);

    %fprintf('Processing cip %1d\n',cipNum);
    
    % Perform DC bias correction
    myCip.correctDcBias();
    myCip.calibrate();

    [cipDataDc,cipMetaData] = myCip.getCipCube();
    
    % Perform pulse compression
    myCip.compress();

    %myCip.dopplerFilter();   NO!
    myCip.normalize();

    [cipDataPc,cipMetaData] = myCip.getCipCube();
    
    nBurst = nBurst + 1;
    cipStack(:,:,nBurst) = squeeze(cipDataPc(dispChan,:,:));
    % Perform Doppler filtering
  
    % [cipDataDf,cipMetaData] = myCip.getCipCube();
    % zero out invalid range of data.
    %cipStack(1:100,:,nBurst) = NaN;
    % 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
        censor = false;

        feaux = squeeze(getRadarMomentsNCip(OURadarStruct, cipStack, censor));

        moments(:,numFinePsd,:) = feaux;
        nBurst = 0; % this will reinit sumCoarsePsd to 0
        %pause(0.2);

    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

radialbin = 1;
%% Plot additional results
if plot
    minRange = 100;
    maxRange = 5000;
    %choice = 'SR';
    %hds determine where channel data is coming from (radChan.sumxr is an enum)
    if strcmp(lr_or_sr,'LR')
        dispChan=uint32(radChan.sumLr);
        %     endRange=cipMetaData.cells.detEnd(dispChan);
        %hds process 8km
        maxRange = 8000;
        if( cipMetaData.wxMode == true)
            maxRange = 14000;
            minRange = 2000;
        end
        endRange=ceil(maxRange/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
    else
        dispChan=uint32(radChan.sumSr);
        % process only up to 2.2 km
        endRange=ceil(5000/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
        minRange = 100;
        maxRange = 5000;
        %     endRange=cipMetaData.cells.detEnd(dispChan);
    end
    startRange=cipMetaData.cells.rangeZero(dispChan)+2;
    if (cipMetaData.wxMode)

        %scaling=[7 30];
        scaling=[-5 30];
    else
        scaling=[8 50];
    end
 %   range=[startRange endRange];
 

    %% Plot reflectivity
    figure(88);
 
    set(gca,'YDir','normal');
    colormap(cmaps('zmap')); % dmap is weird.
    hcb=colorbar;
    title('Weather Radial Reflectivity Map');
    ylabel(hcb,'dBZ');

    xlabel('Radar Azimuth (deg relative to True North)')
    ylabel('Range (m)')
   
    %xlim(linspace(-55,55,numNcips))
    %xlim([-55 55])
    xlim([1 numNcips])
   
    ylim([minRange maxRange]);

    clim([0 60])

    imagesc('XData',1,'YData',[minRange maxRange] ,'CData',moments(startRange:endRange,:,1));
  
    if make_cf_file
        h5create(cf_filename,'/sweep_0001/dbz',size(moments(startRange:endRange,:,1)));
        h5write(cf_filename,'/sweep_0001/dbz',moments(startRange:endRange,:,1));
        h5create(cf_filename,'/sweep_0001/range',endRange-startRange);
        h5write(cf_filename,'/sweep_0001/range',linspace(minRange, maxRange,endRange-startRange));
    end
    
    
    %% Plot Radial Velocity

    figure(99);


    set(gca,'YDir','normal');
    colormap(cmaps('rgmap',64));
    hcb=colorbar;
    title('Weather Radial Velocity Map (Spectral Estimation)');
    ylabel(hcb,'Radial Velocity (mph)');
 
    xlabel('Radar Azimuth (deg relative to True North)')
    ylabel('Range (m)')

    xlim([1 numNcips])
    ylim([minRange maxRange])

    clim([-60 60]);

    imagesc('XData',1,'YData',[minRange maxRange],'CData',moments(startRange:endRange,:,2));
  
    clim([-60 60])

    %% Plot Spectral Moment (W) 
    figure(77);
   
    set(gca,'YDir','normal');
    %colormap(cmaps('rgmap',64));
    colormap(cmaps('rapmap',64));
    hcb=colorbar;
    title('Weather Spectral Moment Map (Spectral Estimation)');
    ylabel(hcb,'?');
    %xlabel('Azimuth Sector')
    % xlabel('Non-coherent Processing Interval')
    xlabel('Radar Azimuth (deg relative to True North)')
    ylabel('Range (m)')

    xlim([1 numNcips])
    clim([0 40]);
    ylim([minRange maxRange])
   
    radialbin = 1;
  
    imagesc('XData',1,'YData',[minRange maxRange],'CData',moments(startRange:endRange,:,3));
  
    clim([0 40])
  
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;

