%% First stage is to match the calibration networks match (for cal-verify this will be disabled).
% Open Channel 1 data -- data collected with signal divider channel one connected to sum channel.  This is the master data set all other data will be corrected to this one. (other choice is to force all to zero phase)

% function MpdrChanMatchCal(Filename,Pathname)
% MpdrQuickLook(Filenane,Pathname) -- parse a data recording file for MPDR radar
%   filename -- filename of the file to parse
CollectFftWisdom = false;


[aFilename, aPathname]= uigetfile('*.dat',...
    'Select Data-Recording file(s) to process',...
    'multiselect','on');
if (isfloat(aFilename))
    error('No files selected');
end

if (ischar(aFilename))
    aFilename={aFilename}; % make it always a cell array...
end

% loop through files selected...
NumFiles = length(aFilename);
for fileIndex=1:NumFiles
    
    myRadar = mpdrRadar('.\Config\mpdr2RadarConfig.yaml');
    
    myRadar.openDataRecordFile(aFilename{fileIndex},aPathname);

%% Set up pulse compressor...

    
    
    % create pulse compressors (only first time through loop)...
    if fileIndex==1
        cipMetaData=myRadar.cipMetaData;
        
        LrPc = PulseCompressor(waveformTypes.LfmUpChirp,...
            myRadar.cipMetaData.bandWidth(uint32(radChan.sumLr)),...
            myRadar.cipMetaData.pulseLength_s(uint32(radChan.sumLr)),...
            myRadar.cipMetaData.sampleRate_sps,...
            myRadar.cipMetaData.priSize);
        
        SrPc = PulseCompressor(waveformTypes.LfmUpChirp,...
            myRadar.cipMetaData.bandWidth(uint32(radChan.sumSr)),...
            myRadar.cipMetaData.pulseLength_s(uint32(radChan.sumSr)),...
            myRadar.cipMetaData.sampleRate_sps,...
            myRadar.cipMetaData.priSize);
        
        % Pulse compress the CIPs
        LrChans=radChan.getLrChans();
        SrChans=radChan.getSrChans();
        FigOrder=radChan.getFigOrder();
        LrIndex=uint32(radChan.sumLr);
        SrIndex=uint32(radChan.sumSr);
        
        % load calibration coefficients...A
        load ("CalCorr.mat");
        
        %     preallocate statistics...
        MagAvg = zeros(cipMetaData.numChannels,NumFiles,"double");
        MagStd = zeros(cipMetaData.numChannels,NumFiles,"double");
        PhaseAvg = zeros(cipMetaData.numChannels,NumFiles,"double");
        PhaseStd = zeros(cipMetaData.numChannels,NumFiles,"double");
        MeanMonRatio=zeros(4,NumFiles,"double"); % magic #4 = (Daz + Del) * (Sr + Lr)
        MaxPcTime = 0;
    end
    

%

    %   read the whole file at once...
    [inputData,myReadSize]=fread(myRadar.fd,"int16");
    
    % truncate to an even PRI length...
    priLength = 2 * cipMetaData.numChannels * cipMetaData.priSize;
    numPris= floor(myReadSize/priLength);
    % truncate down to even number of PRIs...
    inputData = inputData(1:(priLength*numPris));
    inputData = reshape(inputData,2,[]);
    cipData=complex(inputData(1,:),inputData(2,:));
    cipData = reshape(cipData,...
        cipMetaData.numChannels,...
        cipMetaData.priSize,[]);

%

    
    if CollectFftWisdom && fileIndex==1
        fftw("planner","exhaustive");
    end
    tic
    % pulse compress long/short range channels...
    
    cipData(LrChans,:,:)=LrPc.compress(cipData(LrChans,:,:));
    cipData(SrChans,:,:)=SrPc.compress(cipData(SrChans,:,:));
    
    PcTime=toc;
    MaxPcTime=max(PcTime,MaxPcTime);
    fprintf("%d pris PC time: %f sec MaxTime: %f\n",numPris,PcTime,MaxPcTime);
    
    if CollectFftWisdom && fileIndex==1
        FftwWisdom = fftw("dwisdom");
        save FftwWisdom FftwWisdom -mat;
    end
    %% Plot PC results in dB
    %     plotCube(Data,cipMetaData,['Pulse compressed channels: ' cipMetaData.Filename],3);

%
%% get access to the peaks and examine the differences

    
    cipData=reshape(cipData,cipMetaData.numChannels,[]);
    
    pcDataSum=cipData([LrIndex SrIndex],:);

% normilize to max value...

    maxValue = max(abs(pcDataSum),[],2);

% put it into log form:

    pcDataSum_dB=(20*log10(abs(pcDataSum)./maxValue));
    
    pcDataBlob_dB=zeros(size(cipData));
    pcDataBlob_dB(LrChans,:)=(20*log10(abs(cipData(LrChans,:))./maxValue(1,1)));
    pcDataBlob_dB(SrChans,:)=(20*log10(abs(cipData(SrChans,:))./maxValue(2,1)));
    

% find the peaks, Each peak  has to be > - 5 dB of the largest value,
% plus be a local maximum (greater than the cell before it and greater than the cell after it)

    % P1 Signal is > - 5dB/peak value
    p1=(pcDataSum_dB(:,2:end-1) > -5.0);
    % p2 Signal is > prior cell
    p2=(pcDataSum_dB(:,1:end-2)<pcDataSum_dB(:,2:end-1));
    % p3 Signal is > next cell
    p3=(pcDataSum_dB(:,2:end-1)<pcDataSum_dB(:,3:end));
    % p4 is the logical and of P1-3 (aligning for index differences)
    p4=p1&p2&(p3-1);
    % find the non-zero cells
    LrPeaks = find(p4(1,:))+1;
    SrPeaks = find(p4(2,:))+1;
    fprintf("First Detection: %d  %d \n",LrPeaks(1,1),SrPeaks(1,1));

%
%% Calibrate the data, convert to db...

    cipData=cipData.*calCorr;
    
    maxValue= zeros(cipMetaData.numChannels,1);
    maxValue(LrChans,1)=max(abs(cipData(LrIndex,:)),[],2);
    maxValue(SrChans,1)=max(abs(cipData(SrIndex,:)),[],2);
    
    calPcDataBlob_db = zeros(size(cipData));
    calPcDataBlob_dB=(20*log10(abs(cipData./maxValue)));

%

    peakSamples = zeros(8,max(length(LrPeaks),length(SrPeaks)));
    peakSamples(LrChans,:)=cipData(LrChans,LrPeaks(1,:));
    peakSamples(SrChans,:)=cipData(SrChans,SrPeaks(1,:));
    
    peakMonRatio= [imag(peakSamples(uint32(radChan.dAzLr),:)./peakSamples(uint32(radChan.sumLr),:));...
        imag(peakSamples(uint32(radChan.dElLr),:)./peakSamples(uint32(radChan.sumLr),:));...
        imag(peakSamples(uint32(radChan.dAzSr),:)./peakSamples(uint32(radChan.sumSr),:));...
        imag(peakSamples(uint32(radChan.dElSr),:)./peakSamples(uint32(radChan.sumSr),:))];
    
    
    [peakPhase_r,unCalPeakMag]=cart2pol(real(peakSamples),imag(peakSamples));
    unCalPeakPhase_d = rad2deg(peakPhase_r);
    peakSamples(LrChans,:)=cipData(LrChans,LrPeaks(1,:));
    peakSamples(SrChans,:)=cipData(SrChans,SrPeaks(1,:));
    [peakPhase_r,calPeakMag]=cart2pol(real(peakSamples),imag(peakSamples));
    
    calPeakPhase_d = rad2deg(peakPhase_r);

%

    MagAvg(:,fileIndex) = mean(calPeakMag,2);
    MagStd(:,fileIndex) = std(calPeakMag(:,:),0,2);
    PhaseAvg(:,fileIndex) = mean(calPeakPhase_d,2);
    PhaseStd(:,fileIndex) = std(calPeakPhase_d,0,2);
    MeanMonRatio(:,fileIndex) = mean(peakMonRatio,2);
    for ii=FigOrder
        fprintf("%-40s Chan:%-12s Mag:%-10.2g Sig:%-10.2g Var:%-7.4f P(d):%+-10.2f PSig:%-10.2f",...
            cipMetaData.Filename,...
            radChan.getString(ii), ...
            MagAvg(ii,fileIndex), ...
            MagStd(ii,fileIndex), ...
            MagStd(ii,fileIndex)/MagAvg(ii,fileIndex), ...
            PhaseAvg(ii,fileIndex), ...
            PhaseStd(ii,fileIndex));
        if (ii == uint32(radChan.sumLr))
            fprintf(" MAz:%-7.4f MEl:%-7.4f\n", MeanMonRatio(1,fileIndex), MeanMonRatio(2,fileIndex));
        elseif (ii == uint32(radChan.sumSr))
            fprintf(" MAz:%-7.4f MEl:%-7.4f \n", MeanMonRatio(3,fileIndex), MeanMonRatio(4,fileIndex));
        else
            fprintf("\n");
        end
    end
    
    clear myRadar;
end

%
%% Plot magnitude between files

if(NumFiles>1)
    lFig=figure;
    set(lFig,'position',[40 50 1500 600]);
    set(lFig,'Name','Magnitude between files');
    tiledlayout(2,4);
    
    maxValue=max(MagAvg(:,:),[],2);
    rOrder=10.^(floor(log10(maxValue)));
    plotLimit=ceil(maxValue./rOrder).*rOrder;
    
    for ii=FigOrder
        nexttile
        plot((MagAvg(ii,:)));
        title(radChan.getString(ii));
        xlabel('Record Files');
        ylabel('counts');
        ylim([0 plotLimit(ii)]);
        xlim([1 NumFiles]);
    end

%
%% Plot Phase between files

    lFig=figure;
    set(lFig,'position',[40 50 1500 600]);
    set(lFig,'Name','Phase between files');
    tiledlayout(2,4);
    
    for ii=FigOrder
        nexttile
        plot((PhaseAvg(ii,:)));
        title(radChan.getString(ii));
        xlabel('Record Files');
        ylabel('degrees');
        ylim([-180 180]);
        xlim([1 NumFiles]);
    end

%
%% Plot Monopulse Ratio between files

    lFig1=figure;
    set(lFig1,'position',[40 50 1500 600]);
    set(lFig1,'Name','Monopulse Ratio');
    tiledlayout(2,2);
    
    maxValue=max(max(MeanMonRatio(:,:)));
    minValue=min(min(MeanMonRatio(:,:)));
    rOrderHigh=10.^(floor(log10(abs(maxValue))));
    rOrderLow=10.^(floor(log10(abs(minValue))));
    plotLimitH=ceil(maxValue./rOrderHigh)*rOrderHigh;
    plotLimitL=floor(minValue./rOrderLow).*rOrderLow;
    
    ratioTitle={'DAz LR','DEl LR','DAz SR', 'DEl SR'};
    
    for ii=1:4
        nexttile
        plot((MeanMonRatio(ii,:)));
        title(ratioTitle{ii});
        xlabel('Record Files');
        ylabel('Imag(Delta/sum)');
        ylim([plotLimitL plotLimitH]);
        xlim([1 NumFiles]);
    end
end