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;
plotem = 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)
lrmode = false;
%% 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));
    lrmode = true;
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);

OURadarStruct.lr=lrmode;

%[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,numPulsesPerCip); % dimensions are range x doppler bin
cipStack = zeros(wx_rng_cell_max,myCip.cipMetaData.numPris,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 = false;

%% for cfradial file
% Hampton location 
latitude = 37.02046216662675;
longitude = -76.33666891140884;

altitude = 0;

%fixme use netcdf create groups, etc.

%% build cfradial 1/2 'common' data for a sweep. using netcdf functions.

if make_cf_file
    cf_filename = "test_ref.nc";

%     ncid = netcdf.create(cf_filename,"netcdf4");
%     childGroupId = netcdf.defGrp(ncid,"latitude");
%    

    h5create(cf_filename,'/latitude',1); h5write(cf_filename,'/latitude',37.02045856274912);
    h5create(cf_filename,'/longitude',1); h5write(cf_filename,'/longitude',-76.3366642324842);
    h5create(cf_filename,'/altitude',1); h5write(cf_filename,'/altitude',40);
    h5create(cf_filename,'/time', numNcips); h5write(cf_filename,'/time',linspace(1,numNcips,numNcips));
    h5create(cf_filename,'/elevation', numNcips); h5write(cf_filename,'/elevation',linspace(10,10,numNcips));
    h5create(cf_filename,'/azimuth',numNcips); h5write(cf_filename,'/azimuth',linspace(-55,55,numNcips));

    h5create(cf_filename,'/sweep_mode',1 ,Datatype='string'); h5write(cf_filename,'/sweep_mode',"azimuth_surveillance");
    h5create(cf_filename,'/sweep_number',1); h5write(cf_filename,'/sweep_number',1); % just what we have for angle now; multiple scans in file = multiple angles.

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

    h5create(cf_filename,'/sweep_start_ray_index',1); h5write(cf_filename,'/sweep_start_ray_index',0); % just what we have for angle now; multiple scans in file = multiple angles.
    h5create(cf_filename,'/sweep_end_ray_index',1); h5write(cf_filename,'/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;
one_start_notification = true;
currentel = 0.0;
target = 10.0;
%PPI -- derive from json 
targets = [2.0 4.0 6.0 10.0 15.0 ]

%targets = [0 2 4] % clutter capture every other ppi elevation.
%targets = [1 3 5]
lastel = -10; 

% RHI specifically from Hampton University. 
rhi_ele = linspace(0,45,1);
rhi_az = [331 334 338 343];

one_start_notification = true;

radialmetadata = cipMetaData;

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();

    % should be able to do a 'mycip.calibrate' at least for leakage.
    % the next step is determining if the pulse loc is correct.
    % otherwise the whole CIP + 1 needs a barrel/circshift left. MPS.

    [cipDataPc,cipMetaData] = myCip.getCipCube();

    if(nBurst == numCipsPerNcip/2)
        radialmetadata = cipMetaData;
    end
    nBurst = nBurst + 1;
    %cipStack(1:65,:,:) = 0;
    cipStack(:,:,nBurst) = squeeze(cipDataPc(dispChan,:,:));
    %cipStack(1:92,:,:) = 0;

        % 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

        % should pull ALL elevation targets as slices, and skip
        % transitional elevations.
        if ismember(round(radialmetadata.cipCenterElevation_d), targets)

            if one_start_notification
                fprintf("starting to process elevation: %f at %d\n", round(radialmetadata.cipCenterElevation_d), cipNum);
                one_start_notification = false;
            end
            %figure(33); plot(abs(cipDataDc(1,:,1)))
            %findchangepts(abs(cipDataDc(1,:,1)),"MaxNumChanges",2,"Statistic",'mean')
            % 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)

            % isn't this a *sign* thing rather than a matrix flip?  
            if (meteorology_mode)
                finePsd = fliplr(finePsd)
            end
            censor = false;

            feaux = squeeze(getRadarMomentsNCip(OURadarStruct, cipStack, censor));
            %fprintf("radial angle (uncorrected): %f at %d\n", abs(radialmetadata.cipCenterAzimuth_d), cipNum);

            moments(:,numFinePsd,:) = feaux;
        else
            currentel = round(radialmetadata.cipCenterElevation_d);
            if currentel ~= lastel && currentel ~= target
                % record end angle and time, plot current slice/ create / add to file.
                fprintf("skipping: %f at %d\n", currentel,cipNum);
                lastel = currentel;
                one_start_notification = true;

            end



            %pause(0.2);

        end

        nBurst = 0; % this will reinit sumCoarsePsd to 0
    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 plotem
    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 = 19000;
            minRange = 2000;
        end
        endRange=ceil(maxRange/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
        startRange=cipMetaData.cells.rangeZero(dispChan)+OURadarStruct.startOffset_LR; % LR adjustment +36;
        
    else
        dispChan=uint32(radChan.sumSr);
        % process only up to 2.2 km

        minRange = 100;
        maxRange = 12000;
        endRange=ceil(maxRange/cipMetaData.rangeCellSize_m(dispChan)+cipMetaData.cells.rangeZero(dispChan));
      
        %     endRange=cipMetaData.cells.detEnd(dispChan);
       startRange=cipMetaData.cells.rangeZero(dispChan)+OURadarStruct.startOffset_SR; % LR adjustment +36;

    end
    %startRange=cipMetaData.cells.rangeZero(dispChan)+40; % LR adjustment +36;
    if (cipMetaData.wxMode)

        %scaling=[7 30];
        scaling=[-5 30];
    else
        scaling=[8 50];
    end
 %   range=[startRange endRange];
 
 fprintf("start/end: %d/%d\n", startRange, endRange);


  figure(66);
    title('raw plot single radial:');
    yvals = moments(:,1,1);
    testvals = linspace(1,433,433)';
    xlim([1 433])
    plot(yvals);
    hold on;
    devcheck = movstd(yvals,8);
    plot(devcheck);
    hold on;
%     [changes,~] = findchangepts(yvals,'MaxNumChanges',4,'Statistic','std');
%     [pvals,plocs] =findpeaks(yvals,"MinPeakHeight",70); % FIXME. Magic #. 
%     startRange  = changes(1); % should be ~100m ; seems to be more like 396; 
% startRange = 60;
%     fprintf("change point: %d\n", changes(1));
%     plot(changes(1),yvals(changes(1)),'b+');
%     text(changes(1),yvals(changes(1)),'Start of good data');

    %plot(1:433, squeeze(moments(1:433,1,2)'));
    %plot(1:433, squeeze(moments(1:433,1,3)'));
    %grid;


    %startRange = 58;
    %% 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 numFinePsd])
   
    ylim([minRange maxRange]);

    clim([0 60])

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

   % for each datatype write slices in el to file; pull this out into
   % functions; makecffile, adddata, setmetadata

   %% populate horizontal reflectivity sweep
    if make_cf_file
        h5create(cf_filename,'/reflectivity_horizontal',size(moments(startRange:endRange,:,1)));
        h5write(cf_filename,'/reflectivity_horizontal',moments(startRange:endRange,:,1));
        h5create(cf_filename,'/range',endRange-startRange);
        h5write(cf_filename,'/range',linspace(minRange, maxRange,endRange-startRange));
%        netcdf.close(ncid);
    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 numFinePsd])
    ylim([minRange maxRange])

    clim([-60 60]);

    imagesc('XData',1,'YData',[minRange maxRange],'CData',moments(startRange:endRange,1:numFinePsd,2));
  
    %% too much?
    clim([-30 30])

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

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

