clear
DATEYR = '2021-';
DAYNUM = '250' ;

%set the data file path name
%data_dir =
%strcat('K:\ify\sol\soon_to_be_deleted_(spring-2013)\Instruments\AndorSpectrometer-old\Data\')
%data_dir = strcat('S:\ify\sol\instruments\AndorSpectrometer\');
data_dir = strcat('M:\Andor_data_from_Autumn_2020\');

load(strcat(data_dir,DATEYR,DAYNUM,'.mat'))

[nrofmeas filelength]=size(daymat);

qqqcnt = 0;
qqqbkg = 0;
ibkg1 = 0;
for i=1: nrofmeas
    n=i;
    tag = daymat(n,10);
    %tag=0->bkg, tag=1->data
    if(i==1)
        DATE1=['Date: ' int2str(daymat(n,1)) '-' int2str(daymat(n,2)) '-' int2str(daymat(n,3)) ' UTC Time: ' int2str(daymat(n,4)) ':' int2str(daymat(n,5))]
    end
    DATE2=['Date: ' int2str(daymat(n,1)) '-' int2str(daymat(n,2)) '-' int2str(daymat(n,3)) ' UTC Time: ' int2str(daymat(n,4)) ':' int2str(daymat(n,5))]
    % if background, set as such
    if(tag==0)
        ibkg1 = n;
        continue
    else
        qqqcnt = qqqcnt + 1;
        dat = daymat(i,12:end);
    end
    %find the next background. Or if the first doesn't exist, start over
    if(ibkg1 == 0)
        continue
    else
        nend = n + 20
        if (nend >= nrofmeas)
            nend = nrofmeas;
        end
        for kkk= n+1:nend
            tag = daymat(kkk,10);
            ibkg2 = kkk;
            if(tag==0)
                continue
            end
        end
        if(tag~=0)
            ibkg2=ibkg1;
        end
    end
    half = (ibkg2-ibkg1)/2.0 ;
    ibkg = ibkg1;
    if((n-ibkg1) > half)
        ibkg = ibkg2;
    end
    %subtract off nearest background, all corrected to counts per sec(niumber of intervals * sec/interval)
    dat = dat/(daymat(i,7)*daymat(i,8)) - daymat(ibkg,12:end)/(daymat(ibkg,7)*daymat(ibkg,8));
    %median filter data
    dat=medfilt1(dat,3);
    
    a31(qqqcnt) = sum(dat(200:223))/(223-200+1);
    a42(qqqcnt) = sum(dat(690:705))/(705-690+1);
    aNa(qqqcnt) = sum(dat(993:1024))/(1024-993+1);
    dat31 = dat-a31(qqqcnt);
    dat42 = dat-a42(qqqcnt);
    datNa = dat-aNa(qqqcnt);
       
    a31(qqqcnt) = sum(dat31(231:289));
    a42(qqqcnt) = sum(dat42(631:670));
    aNa(qqqcnt) = sum(datNa(918:987));
        
    t(qqqcnt) = (daymat(n,4))+(daymat(n,5))/60 ;
    if(t(qqqcnt) < 10)
        t(qqqcnt) = t(qqqcnt) + 24;
    end
end
plot(t,a31,'b')
hold
plot(t,a42,'r')
plot(t,aNa,'k', 'Linewidth', 1.5)
hold
title(strcat(DATE1,' to : ',DATE2));
xlabel('UT hours')
ylabel('OH M(3,1) & (4,2) Q-branch Intensities [C/s]')
legend('OH M(3,1)', 'OH M(4,2)', 'Na 1638nm(cloud)')
ylim([0 inf]);
mat = [ a31', a42', aNa'];
cmat = corrcoef(mat);
qqq = [' (3,1)-(4,2) correlation coeff = ',num2str(cmat(1,2))]

qqq = [' (3,1)-Na correlation coeff = ',num2str(cmat(1,3))]

qqq = [' (4,2)-Na correlation coeff = ',num2str(cmat(2,3))]
