%TO DO: MAKE QC INDICATOR AND NOT 30 SUM MATCH getDirunal_nc
%This is because getDirunal is no longer used to produce PBLH NetCDFs,
%so, the QC vector is using our old way of determining the QC indicator
%instead of negative SNR indicator. This should still match for consistency

function[diurnalCycle,smoothDiurnalCycle,polyEvalTimes,diurnalCycleAbsFactors,snrRaster,actuallevelsTimeRange,actualprofilesTimeRange,matrixMapR,matrixMapC,minheight,maxheight,minutesvector,heightvector,allminutes,allheights,emergeTime,emergeHeight,minutestamps,resolution,timeResolution,removal,smoothResolutionCurve,QCvector2]=getDiurnal_csv(m,d,y,levelheights,snr,dates,avgtime)
    %Run doloop on the raw data to separate profiles for a specified day
    [actuallevels,actualprofiles,newHours,newMinutes,~]=makeProfiles(m,d,y,levelheights,snr,dates,avgtime);   
    %convert the times of each profile start into number of minutes since
    %midnight and store them in minutestamps if they are between 9AM and
    %8PM. If the timestamp is out of that range, put NaN in its position.
    

    
    [minutestamps]=NaN.*ones(size(newHours,1),1);
    for i=1:size(newHours,1)
        stamp=newHours(i)*60+newMinutes(i);
        if stamp>=540 && stamp<=1200
            minutestamps(i)=stamp;
        end
    end
        
    
    %Build a vector that has the same size as minutestamps. Fill all values
    %with NaN. Change the value to 1 if minutestamps has a value at that
    %position. 
    profsInRangeIdx=NaN.*ones(size(minutestamps,1),1);
    for i=1:size(minutestamps)
        if ~isnan(minutestamps(i))
            profsInRangeIdx(i)=1;
        end
    end


    %Use the index of timestamps in range to pull the profiles and levels
    %associated with those timestamps.
    [actuallevelsTimeRange]=actuallevels;
    [actualprofilesTimeRange]=actualprofiles;

    for i=1:size(actuallevels,2)
        actuallevelsTimeRange(:,i)=actuallevelsTimeRange(:,i).*profsInRangeIdx;
        actualprofilesTimeRange(:,i)=actualprofilesTimeRange(:,i).*profsInRangeIdx;
    end




    %For levels and profiles, we have to be careful to only remove rows
    %that are entirely NaN because some profiles/levels may have a NaN
    %within the profile, and we do not want to remove that profile.
    actuallevelsTimeRange(~any(~isnan(actuallevelsTimeRange), 2),:)=[];
    actualprofilesTimeRange(~any(~isnan(actualprofilesTimeRange), 2),:)=[];
    
    


   %Now that all of our timestamps, profiles, and levels are in range, we
   %remove the data outside of range.
    minutestamps=minutestamps(~isnan(minutestamps));
   % disp(minutestamps);

    maxheight=max(actuallevelsTimeRange(:));
    minheight=min(actuallevelsTimeRange(:));

    [heightvector]=minheight:1:maxheight;
    try
    [minutesvector]=minutestamps(1,1):1:minutestamps(size(minutestamps,1),1);
    catch
        disp("Not enough data for this date.");
    end

    numheights=size(heightvector,2);
    numminutes=size(minutesvector,2);

    allheights=1:1:4000;
    allminutes=540:1:1200;

    %Build empy matrix for interpolated SNR matrix based on number of
    %meters between lowest and highest recording and number of minutes
    %between first and last recording
    [snrRaster]=NaN.*ones(numheights,numminutes);

    %Building 2 matrices of size actuallevelsTimeRange that will contain the row
    %and column indices of where each element is plotted in the
    %interpolated SNR matrix
    %Also inputs values from actuallevelsTimeRange into SNR matrix
    matrixMapR=NaN.*ones(size(actuallevelsTimeRange,1),size(actuallevelsTimeRange,2));
    matrixMapC=NaN.*ones(size(actuallevelsTimeRange,1),size(actuallevelsTimeRange,2));
    for i =1:size(actuallevelsTimeRange,1)
        colIndex=find(minutesvector==minutestamps(i));  
        for j=1:size(actuallevelsTimeRange,2)
            if ~isnan(actuallevelsTimeRange(i,j))
                rowIndex=find(heightvector==actuallevelsTimeRange(i,j));
                snrRaster(rowIndex,colIndex)=actualprofilesTimeRange(i,j);
                matrixMapR(i,j)=rowIndex;
                matrixMapC(i,j)=colIndex;
            end
        end
    end

    %linear interpolation in height and time
    snrRaster0 = fillmissing(snrRaster,'linear');
    snrRaster0 = fillmissing(snrRaster0','linear')';
    [snrRaster]=snrRaster0;
    

    %Finding profiles that have high starting heights or low ending heights
    %to remove overly extrapolated values
    highStarts=NaN.*ones(size(actuallevelsTimeRange,1),1);
    lowEnds=NaN.*ones(size(actuallevelsTimeRange,1),2);

    for i=1:size(actuallevelsTimeRange,1)
        if actuallevelsTimeRange(i,1)~=minheight
            highStarts(i,1)=i;
        else
            highStarts(i,1)=0;
        end
        [max0,idx]=max(actuallevelsTimeRange(i,:));
        if max0~=maxheight
            lowEnds(i,1)=i;
            lowEnds(i,2)=idx;
        else
            lowEnds(i,1)=0;
            lowEnds(i,2)=9999;
        end
    end

    %remove overly extrapolated values
    for i=1:size(highStarts,1)
        nancolr=0;
        nancoll=0;
        if highStarts(i)>0
            nanrow=matrixMapR(i,1)-50;
            nancol=matrixMapC(i,1);
            if i<size(highStarts,1)
                nancolr=matrixMapC(i+1,1);
            end
            if i>1
                nancoll=matrixMapC(i-1,1);
            end
            snrRaster(1:nanrow,nancol)=NaN;
            if nancolr>0
                snrRaster(1:nanrow,nancol:nancolr-1)=NaN;
            end
            if nancoll>0
                snrRaster(1:nanrow,nancoll+1:nancol)=NaN;
            end
        end
        if lowEnds(i,1)>0
            nanrow=matrixMapR(i,lowEnds(i,2))+50;
            nancol=matrixMapC(i,lowEnds(i,2));
            if i<size(lowEnds,1)
                [~,idx]=max(actuallevelsTimeRange(i+1,:));
                nancolr=matrixMapC(i+1,idx);
            end
            if i>1
                [~,idx]=max(actuallevelsTimeRange(i-1,:));
                nancoll=matrixMapC(i-1,idx);
            end
            snrRaster(nanrow:end,nancol)=NaN;
            if nancolr>0
                snrRaster(nanrow:end,nancol:nancolr-1)=NaN;
            end
            if nancoll>0
                snrRaster(nanrow:end,nancoll+1:nancol)=NaN;
            end
        end
    end
    
    

    %emergence time
    %720=12PM
    if minutesvector(size(minutesvector,2))>=720
        twelvepmInd=find(minutesvector==720);
    else
        twelvepmInd=(size(minutesvector,2));
    end

    emergeIndR=NaN;
    emergeIndC=NaN;

    for i=1:size(snrRaster,1)
        rastRow=snrRaster(i,1:twelvepmInd);
        numNans=sum(isnan(rastRow));
        %make sure the row we are looking at for emergence time is filled
        %at least 35% with values
        if numNans<0.65*(minutesvector(1,twelvepmInd)-minutesvector(1,1))
            emergeIndR=i;
            break
        end
    end
    if ~isnan(emergeIndR)
        [~,emergeIndC]=max(snrRaster(emergeIndR,1:twelvepmInd));
    else
        for i=1:size(minutesvector,2)
            if minutesvector(i)>720 && minutesvector(i)<840
                emergeIndC=i;
                break
            end
        end
    end

    startProf=find(matrixMapC==emergeIndC);
    while isempty(startProf) && emergeIndC<size(minutesvector,2)
        emergeIndC=emergeIndC+1;
        startProf=find(matrixMapC==emergeIndC);
    end


    emergeTime=minutesvector(emergeIndC);
    if ~isnan(emergeIndR)
        emergeHeight=heightvector(emergeIndR);
    else
        emergeHeight=heightvector(1);
    end

    diurnalCycle=NaN.*ones(1,size(actuallevelsTimeRange,1));
    QCvector=NaN.*ones(1,size(actuallevelsTimeRange,1));

    difs2=[];
    for i=1:size(actuallevelsTimeRange,1)
       difs=[];
       for j=2:size(actuallevelsTimeRange,2)
            dif=actuallevelsTimeRange(i,j)-actuallevelsTimeRange(i,j-1);
            difs=[difs dif];     
       end
       difs2=[difs2 nanmean(difs)];
    end
    [resolution]=nanmean(difs2);


    within2HrsEtime=max(find(minutestamps<=(emergeTime+120)));



    %applying algorithm to profiles
    for i=startProf(1):size(actuallevelsTimeRange,1)
        if i<=within2HrsEtime
            hookstat=0;
        else
            hookstat=1;
        end
        disp('profnum');
        disp(i);
        disp('hookstat');
        disp(hookstat);
        sigvec=snrRaster(:,matrixMapC(i,1))';
        signonan=sigvec(~isnan(sigvec));
        heightvector2=heightvector(~isnan(sigvec));
        


        if sum(signonan<0) > 0.1*size(signonan,2) && heightvector2(1)>1000
            QCvector(i)=4;
        elseif sum(signonan<0) > 0.1*size(signonan,2)
            QCvector(i)=2;
        elseif heightvector2(1)>1000
            QCvector(i)=3;
        else
            QCvector(i)=1;
        end

        
        if heightvector2(1)>1000 || heightvector2(size(heightvector2,2))<2000
            disp('Insufficient Height Range');
            disp(heightvector2(1)>1000);
            disp(heightvector2(size(heightvector2,2)));
            signonan=NaN;
        end

        noninterpsig=actualprofilesTimeRange(i,:);
        %noninterpsig=noninterpsig(~isnan(noninterpsig));


        %applying algorithm to interpolated profiles from times of
        %measurement
        %40 and 30 are the percentages applied to the absolute max of a
        %profile to determine the abs factor
        if (heightvector2(1,1)-emergeHeight)<1000
            [PBLHdia,absFactor,negIndicator,not30Ind]=runAlgorithm(signonan,heightvector2,noninterpsig,hookstat,30,40);
        else
            PBLHdia=NaN;
            absFactor=-99;
        end


        if i==startProf(1)
            diurnalCycle(i)=emergeHeight; 
            diurnalCycleAbsFactors(i)=-99;
        else
            diurnalCycle(i)=PBLHdia;
            diurnalCycleAbsFactors(i)=absFactor;   
        end                      
    end


    difs=[];
    for i=1:size(minutestamps,1)-1
       dif=minutestamps(i+1,1)-minutestamps(i,1);
       difs=[difs dif];     
    end
    [timeResolution]=nanmean(difs);

    if timeResolution<=48
        consecNanCheck=round(120/timeResolution);
    else
        consecNanCheck=2;
    end

    nanct=0;
    removal=0;
    for i=startProf(1):size(diurnalCycle,2)
        if nanct>consecNanCheck
            diurnalCycle=NaN;
            removal=1;
            break
        end
        if isnan(diurnalCycle(1,i))
            nanct=nanct+1;
        else
            nanct=0;
        end
    end

    numNansInCycle=sum(isnan(diurnalCycle));
    if numNansInCycle>=.75*size(diurnalCycle,2)
        diurnalCycle=NaN;
        removal=1;
    end






    if size(diurnalCycle,2)>5
        forPolyCurve=diurnalCycle(1,:);
        if rem(emergeTime,60)==0
            nearDiv60start=emergeTime;
        else
            nearDiv60start = emergeTime + (60 - rem(emergeTime,60));
        end
        lastnonan=size(diurnalCycle,2);
        for i = size(diurnalCycle,2):-1:1
            if ~isnan(diurnalCycle(1,i))
                lastnonan=i;
                break
            end
        end

        lasttime=minutestamps(lastnonan,1);

        nearDiv60end = lasttime - rem(lasttime,60) ;
        polyEvalTimes=(nearDiv60start:60:nearDiv60end);

        minForPoly=minutestamps(:,1);
        nanIndexes = isnan(forPolyCurve);

        minForPoly(nanIndexes) = [];
        forPolyCurveNoNan=forPolyCurve;
        forPolyCurveNoNan(nanIndexes) = [];
        smfact=3;

        [q,~,mu] = polyfit(minForPoly', forPolyCurveNoNan, smfact);
        [smoothDiurnalCycle] = polyval(q,polyEvalTimes,[],mu);
    else
        smoothDiurnalCycle=NaN;
        polyEvalTimes=NaN;
    end
    
    
    

    if any(smoothDiurnalCycle<0) || any(smoothDiurnalCycle>4000)
        smoothDiurnalCycle=NaN;
    end
    

    forResPoly=difs2*2;    
    minForPoly=minutestamps(:,1);
    nanIndexes = isnan(forResPoly);
    minForPoly(nanIndexes) = [];
    forResPolyNoNan=forResPoly;
    forResPolyNoNan(nanIndexes) = [];
    smfact=3;

    [q,~,mu] = polyfit(minForPoly, forResPolyNoNan, smfact);
    [smoothResolutionCurve] = polyval(q,polyEvalTimes,[],mu);

    QCvector2=NaN.*ones(1,size(polyEvalTimes,2));
    for i = 1:size(polyEvalTimes,2)   
        [~,index] = min(abs(minutestamps-polyEvalTimes(i)));
        QCvector2(i) = QCvector(index);
    end

        
        
    
end
