%script to read and process Sebastian's HOBO data
close all;
clear all;
scrz=get(0,'screensize');

fname='hobo_grandeur_v2.csv';
hobo_data=importdata(fname,',',8);
GDdata=hobo_data.data;
clear hobo_data
GDhgts=[1300,1305,1364,1398,1454,1513,1545,1603,1649,1696,1749,1801,1846,1899,1950,1997,2049,2104,2151,2209,2251,2306,2354,2405,2455,2498];

GDdata(GDdata==-99.99)=nan;
%data=inpaint_nans(data);
%MONTH,DAY,YEAR,HOUR (UTC),MINUTE
month=GDdata(:,1);
day=GDdata(:,2);
year=GDdata(:,3);
hour=GDdata(:,4);
mins=GDdata(:,5);
GDdata(:,1:5)=[];
GDtemp=GDdata;
GDtemp=inpaint_nans(GDtemp);
clear GDdata;
GDhobotime=datenum(year,month,day,hour,mins,zeros(length(mins),1));

%%
fname='hobo_harker_v2.csv';
hobo_data=importdata(fname,',',8);
HKdata=hobo_data.data;
clear hobo_data
HKhgts=[1312,1312,1352,1394,1474,1519,1563,1594,1657,1707,1757,1813,1857,1905,1958,2005,2056,2111,2156,2209,2260,2307,2374,2410,2448,2503];
HKdata(HKdata==-99.99)=nan;
%data=inpaint_nans(data);
%MONTH,DAY,YEAR,HOUR (UTC),MINUTE
clear month day year hour min
month=HKdata(:,1);
day=HKdata(:,2);
year=HKdata(:,3);
hour=HKdata(:,4);
mins=HKdata(:,5);
HKdata(:,1:5)=[];
HKtemp=HKdata;
HKtemp=inpaint_nans(HKtemp);
clear HKdata;
HKhobotime=datenum(year,month,day,hour,mins,zeros(length(mins),1));


%%
% fname='hobo_kennecott_t_v2.csv';
% hobo_data=importdata(fname,',',8);
% KTdata=hobo_data.data;
% clear hobo_data
% KThgts=[2159,2060,1986,1945,1896,1848,1807,1765,1712,1673,1608,1551,1503,1461,1395,1354,1328];
% KTdata(KTdata==-99.99)=nan;
% KTdata=inpaint_nans(KTdata);
% %MONTH,DAY,YEAR,HOUR (UTC),MINUTE
% clear month day year hour min
% month=KTdata(:,1);
% day=KTdata(:,2);
% year=KTdata(:,3);
% hour=KTdata(:,4);
% mins=KTdata(:,5);
% KTdata(:,1:5)=[];
% KTdata(:,1:16)=[];
% KTtemp=KTdata;
% KTtemp(KTtemp>30)=nan;
% KTtemp=inpaint_nans(KTtemp,5);
% clear KTdata;
% KThobotime=datenum(year,month,day,hour,mins,zeros(length(mins),1));



%% Now load the RASS data to see how these fit together

%first load the 1 hourly data (sounding grid)
load('double_blend.mat','blended','BlendTC','DPth','RHth','hgts','master');
%blended = Potential Temperature Time Height on Sounding Grid (uses RASS,
%surface obs, and soundings).
% DPth = Dew points time height
% RHth= Relative humidity time heights
% master = matlab serial time vector
% hgts= heights;


%next load the RASS and RWP gridded data
load('RASS_RWP.mat','RASSpres','THblend','TCblend','THblendRWP','times2','vert2','tv','u','v','shear','RI','dthdz','vert','times');
% times2 = times for RASS data grids
% vert2= heights for RASS grid
% times= times for RWP grids
% vert= range gates for RWP
% THblend = Potential temperature on RASS grid
% tv= RASS virtual temperature on RASS grid
% THblendRWP= Potential temperature on RWP grid
% dthdz= vertical gradient in potential temperature (on RWP grid);
% RI= richardson number (RWP grid)
% u= uwind on RWP grid
% v= vwind on RWP grid
% shear= shear on RWP grid

%convert hobo data to potential temperatuer....

GDpres=interp2(times2,vert2,RASSpres',GDhobotime,GDhgts);
GDpres=permute(GDpres,[2 1]);
GDtheta=(GDtemp+273.15).*(1000./GDpres).^(2/7);
TC_GD=interp2(master,hgts,permute(BlendTC,[2 1])',GDhobotime,GDhgts);
TC_HK=interp2(master,hgts,permute(BlendTC,[2 1])',HKhobotime,HKhgts);

%TC_GD=interp2(times2,vert2,TCblend',GDhobotime,GDhgts);
%TC_HK=interp2(times2,vert2,TCblend',HKhobotime,HKhgts);
TC_GD=permute(TC_GD,[2 1]);
TC_HK=permute(TC_HK,[2 1]);

HK_GD=interp2(HKhobotime,HKhgts(2:end),HKtemp(:,2:end)',GDhobotime,GDhgts);

%now interpolate the T


%%
close all;

for IOP=[1]
    if IOP ==1
        starter=datenum(2010,12,2,0,00,00);
        ender=datenum(2010,12,04,0,00,00);
        %starter=datenum(2010,12,2,18,00,00);
        %ender=datenum(2010,12,04,00,00,00);
    elseif IOP==2
        starter=datenum(2010,12,7,00,00,00);
        ender=datenum(2010,12,11,0,00,00);
    elseif IOP==3
        starter=datenum(2010,12,11,00,00,00);
        ender=datenum(2010,12,15,0,00,00);
    elseif IOP==4
        starter=datenum(2010,12,23,18,00,00);
        ender=datenum(2010,12,26,19,00,00);
    elseif IOP==5
        starter=datenum(2010,12,31,18,00,00);
        ender=datenum(2011,1,10,0,00,00);
    elseif IOP==6
        starter=datenum(2011,1,10,00,00,00);
        ender=datenum(2011,1,18,0,00,00);
    elseif IOP==7
        starter=datenum(2011,1,19,00,00,00);
        ender=datenum(2011,1,23,0,00,00);
    elseif IOP==8
        starter=datenum(2011,1,22,00,00,00);
        ender=datenum(2011,1,27,0,00,00);
    elseif IOP==9
        %starter=datenum(2011,1,25,18,00,00);
        %starter=datenum(2011,1,28,12,00,00);
        %ender=datenum(2011,2,1,00,00,00);
        starter=datenum(2011,1,29,16,0,0);
        ender=datenum(2011,1,30,16,00,00);
    elseif IOP==10
        starter=datenum(2011,2,1,00,00,00);
        ender=datenum(2011,2,6,0,00,00);
    elseif IOP==11 %not really an IOP just the mix out during IOP1
        starter=datenum(2010,12,02,15,00,00);
        ender=datenum(2010,12,03,18,00,00);
    elseif IOP==12
        %cold fropa
        starter=datenum(2010,12,20,00,00,00);
        ender=datenum(2010,12,21,00,00,00);
    end
    midx=find(master>= starter & master<=ender);
    GDidx=find(GDhobotime>= starter & GDhobotime<=ender);
    HKidx=find(HKhobotime>= starter & HKhobotime<=ender);
    
    tidx1=find(times>= starter & times<=ender);
    tidx2=find(times2>= starter & times2<=ender);
    iopnum=strcat('0',num2str(IOP));
    
    cmap=[gray(10); flipud(gray(10))];
    cmap(1:10,3)=1;
    cmap(11:end,1)=1;
    
    
    f=figure(IOP);
    clf
    colormap(cmap);
    %colormap(jet);
    set(f,'position',[scrz(1) scrz(2) scrz(3) scrz(4)],'color','white')
    subplot(2,1,1)
    %pcolor(GDhobotime(GDidx),GDhgts,GDtemp(GDidx,:)');shading interp;
    
    pcolor(GDhobotime(GDidx),GDhgts,GDtemp(GDidx,:)'-TC_GD(GDidx,:)');shading interp;
    %contourf(GDhobotime(GDidx),GDhgts,GDtemp(GDidx,:)'-TC_GD(GDidx,:)',-5:.5:5,'linestyle','none');
    caxis([-5 5]);%10]);
    hold on;
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:1:400,'-','linewidth',1,'linecolor','k')
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:5:400,'-','linewidth',2,'linecolor','k')
    contour(master(midx),hgts,blended(:,midx),200:1:400,'-','linewidth',1,'linecolor','k');
    contour(master(midx),hgts,blended(:,midx),200:5:400,'-','linewidth',2,'linecolor','k');
    
    %caxis([min(min(GDtheta(GDidx,:))) max(max(GDtheta(GDidx,:)))]);
    
    hold on;
    datetick('x','mmm-dd')
    axis tight
    hold on;
    title('Gradeur Temp. Anom.','fontsize',15,'fontweight','bold');
    set(gca,'fontsize',15,'fontweight','bold','layer','top','linewidth',2);
    
    ylim([1290 2500])
    colorbar;
    
    subplot(2,1,2)
    %    pcolor(HKhobotime(HKidx),HKhgts,HKtemp(HKidx,:)');shading interp;
    
    pcolor(HKhobotime(HKidx),HKhgts,HKtemp(HKidx,:)'-TC_HK(HKidx,:)');shading interp;
    caxis([-5 5]);%10]);
    hold on;
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:1:400,'-','linewidth',1,'linecolor','k')
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:5:400,'-','linewidth',2,'linecolor','k')
    contour(master(midx),hgts,blended(:,midx),200:1:400,'-','linewidth',1,'linecolor','k');
    contour(master(midx),hgts,blended(:,midx),200:5:400,'-','linewidth',2,'linecolor','k');
    
    %caxis([min(min(GDtheta(GDidx,:))) max(max(GDtheta(GDidx,:)))]);
    
    hold on;
    datetick('x','mmm-dd')
    axis tight
    hold on;
    title('Harkers Temp. Anom.','fontsize',15,'fontweight','bold');
    set(gca,'fontsize',15,'fontweight','bold');
    ylim([1290 2500])
    colorbar;
    xlabel('Month-Day [UTC]','fontsize',15,'fontweight','bold');
    %export_fig(sprintf('IOP%d_HOBO_COMP.png', IOP),'-m2.5')
    
    
    %%
    figure(500)
    colormap(cmap);
    pcolor(GDhobotime(GDidx),GDhgts,GDtemp(GDidx,:)'-HK_GD(:,GDidx));shading interp;
    
    %pcolor(HKhobotime(HKidx),HKhgts,HKtemp(HKidx,:)'-TC_HK(HKidx,:)');shading interp;
    caxis([-5 5]);
    hold on;
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:1:400,'-','linewidth',1,'linecolor','k')
    %contour(times2(tidx2),vert2,THblend(tidx2,:)',200:5:400,'-','linewidth',2,'linecolor','k')
    contour(master(midx),hgts,blended(:,midx),200:1:400,'-','linewidth',1,'linecolor','k');
    contour(master(midx),hgts,blended(:,midx),200:5:400,'-','linewidth',2,'linecolor','k');
    
    %caxis([min(min(GDtheta(GDidx,:))) max(max(GDtheta(GDidx,:)))]);
    
    hold on;
    datetick('x','mmm-dd')
    axis tight
    hold on;
    title('Grandeur - Harkers','fontsize',15,'fontweight','bold');
    set(gca,'fontsize',15,'fontweight','bold');
    ylim([1290 2500])
    colorbar;
    xlabel('Month-Day [UTC]','fontsize',15,'fontweight','bold');
    
    
    
    
end

%% parse out data near 1600 m to compare fluctuations
close all;

gdhidx=find(GDhgts<=1900 & GDhgts>=1500);
hkhidx=find(HKhgts<=1900 & HKhgts>=1500);
GDmean=mean(GDtemp(GDidx,gdhidx),2);
HKmean=mean(HKtemp(HKidx,hkhidx),2);
figure(50);clf
plot(GDhobotime(GDidx),GDmean);hold on;
%plot(HKhobotime(HKidx),HKtemp(HKidx,hkhidx),'r');
plot(HKhobotime(HKidx),HKmean);
datetick('x')

%now filter these mean profiles
nhr=3;
w1=1/(nhr*12); % cutoff in 1/minutes -5 minute obs-> (60/5)=12
wn=w1/.5; %%normalized by nyquist frq.
[b,c]=butter(5,wn,'high'); %construct the high pass filter


HKprime=filtfilt(b,c,HKmean);
GDprime=filtfilt(b,c,GDmean);


figure(51); clf
plot(GDhobotime(GDidx),GDprime);
hold on;
plot(HKhobotime(HKidx),HKprime,'r');
datetick('x');
xlim([datenum(2010,12,2) datenum(2010,12,3)]);

%dancing hobos

B=ones(1,2)*(1/2);
HKtemp=filter2(B,HKtemp);
GDtemp=filter2(B,GDtemp);

HKGD=interp2(HKhobotime,HKhgts(2:end),HKtemp(:,2:end)',GDhobotime,GDhgts);
HKGD=permute(HKGD,[2 1]);
HKGD=inpaint_nans(HKGD);

%%
sidx=find(GDhobotime(GDidx)==datenum(2010,12,2,5,15,0));
eidx=find(GDhobotime(GDidx)==datenum(2010,12,2,8,35,0));

f=figure(22);clf; set(f,'color','w');
pbar=850*10^2;
tbar1=squeeze(nanmean(273.15+(HKGD(GDidx(sidx:eidx),:,1))));
tbar=nanmean(tbar1);
rhobar=pbar./(287*tbar);
pprime=zeros(size(GDidx(sidx:eidx)),1);
dz=diff(GDhgts);
HKGD(:,1)=HKGD(:,2);
g=9.81;
aa=0;
   subplot(3,1,1:2)
   hold on;
hnew=[1000:10:3000];
adiabat=(-9.8 *10^-3).*hnew;
plot(adiabat +10,hnew,'--g');
plot(adiabat +15,hnew,'--g');
plot(adiabat +20,hnew,'--g');
plot(adiabat +25,hnew,'--g');
plot(adiabat +30,hnew,'--g');
for tt=sidx:1:eidx%1:length(GDidx);
   aa=aa+1;
   subplot(3,1,1:2)
   hold on;
   plot(HKGD(GDidx(tt),:),GDhgts,'k')
       xlim([-6 5]);
       ylim([GDhgts(1) 2500])
%    subplot(1,2,2)
%    hold on;
%    plot(GDtemp(GDidx(tt),:),GDhgts,'r');
%     xlim([-8 8]);
    tsum=0;
    for hh=1:length(GDhgts)-3;
       tnow=((273.15 + HKGD(GDidx(tt),hh)) + (273.15 + HKGD(GDidx(tt),hh+1)))./2;
       tprime=tnow-tbar;
       tsum=(tsum+(tprime.*(dz(hh))));
    end
    pprime(aa)=(-(rhobar.*g)./tbar).*tsum;
    %title(datestr(GDhobotime(GDidx(tt))));
   pause(.5)
end
%%


%%
hold on;
tidx2=find(GDhobotime(GDidx)==datenum(2010,12,2,6,15,0));
tidx3=find(GDhobotime(GDidx)==datenum(2010,12,2,8,10,0));
plot(HKGD(GDidx(tidx2),:),GDhgts,'r','linewidth',2)
plot(tbar1-273.15,GDhgts,'--r','linewidth',3);
plot(HKGD(GDidx(tidx3),:),GDhgts,'r','linewidth',2)
box on;
grid on;
xlabel('Temp {\circ}C','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold')
text(-5.8,2390,'a)','fontsize',12,'fontweight','bold','backgroundcolor',[.7 .7 .7],'edgecolor','k');

%%
subplot(3,1,3);
movegui('northwest')
plot(GDhobotime(GDidx(sidx:eidx)),detrend(pprime),'-k','linewidth',2);datetick('x');
ylabel('Pressure [Pa]','fontsize',10,'fontweight','bold');
axis tight;
ylim([-45 45]);
box on;
grid on;
set(gca,'fontsize',10,'fontweight','bold')
text(datenum(2010,12,2,5,18,0),28,'b)','fontsize',12,'fontweight','bold','backgroundcolor',[.7 .7 .7],'edgecolor','k');

xlabel('Time [UTC]','fontsize',10,'fontweight','bold');
return
%%
figure(505);clf
mh=mean(HKGD(GDidx(sidx:eidx),:),1);
md=median(HKGD(GDidx(sidx:eidx),:),1);
mn=min(HKGD(GDidx(sidx:eidx),:),[],1);
mx=max(HKGD(GDidx(sidx:eidx),:),[],1);

stdh=std(HKGD(GDidx(sidx:eidx),:),1);
plot(mh,GDhgts,'-k','linewidth',2); hold on;
plot(mh+stdh,GDhgts,'--k');
plot(mh-stdh,GDhgts,'--k');
plot(mn,GDhgts,':k');
plot(mx,GDhgts,':k');
return
%% difference plot
f=figure(23);clf;set(f,'color','w');
slosh=(HKGD(GDidx(sidx:eidx),:)-GDtemp(GDidx(sidx:eidx),:));
%now filter slosh
nhr=6;
nsamp=12;
w=1./(nhr.*nsamp);
wn=w./.5;
[b,c]=butter(5,wn,'high')

for hh=1:length(GDhgts);
   sprime(:,hh)=filtfilt(b,c,slosh(:,hh));
end


cmap=[gray(10);flipud(gray(10))];
cmap(1:10,3)=1;
cmap(11:end,1)=1;
colormap(cmap)
pcolor(sidx:eidx,GDhgts,sprime');shading interp;colorbar;caxis([-1.5 1.5]);
hold on;


%%
mslosh=mean(slosh(GDidx,3:18),2);
%low pass it
%[b3,c3]=butter(5,1/6);
%mslosh=filtfilt(b3,c3,mslosh);


f=figure(24); clf; set(f,'Units','Pixels','pos',[100 100 900 800],'color','w');
subplot(2,3,1:3);
plot(GDhobotime(GDidx),(mslosh),'linewidth',2);hold on;
%plot(datenum(2010,12,2,5,0,0),3,'ok');
datetick('x')
xlim([datenum(2010,12,2) datenum(2010,12,3,18,0,0)]);
grid on;
set(gca,'box','on','layer','top','linewidth',2,'fontsize',10,'fontweight','bold');
title('Across Basin Temperature Difference','fontsize',10,'fontweight','bold');
ylabel('Temp [ ^{\circ}C ]','fontsize',10,'fontweight','bold')


%%
nhr=6;
fs=1/(5*60); % 5 minute data -> sampling frequency

% the length of each fft must be long enough to resolve the 6-8 hour
% fluctuations of interest...
L=nhr*12; %number of samples in the FFT (1 minute data)

%number of FFTs
nffts=round(length(mslosh)/L);
disp(nffts)

[pxx,ff,pxxc]=pwelch(detrend(mslosh),hanning(L),round(4*L/5),L,fs,'ConfidenceLevel',.95);
subplot(2,3,5);
loglog(ff,ff.*pxx,'-k','linewidth',2)
hold on;
xlabel('Hz','fontsize',10,'fontweight','bold');
ylabel('Power Spectral Density (K^2/s)','fontsize',10,'fontweight','bold');
t12=1./(60*60*12); %12 hours
t4=1./(3600*4); %4hours
t1=1./(3600*1.25); %1 hour
tkh=1./(60*30); % 30 minutes
N=sqrt((9.81./285).*(15/600)); %brunt-vaisala frequency
solarday=1./(60*60*(9) + 25*60);% 9 hours 25 minutes
%vline(t1,'--k','1.2 Hr');
plot([t1 t1],[10^-4 1],'--k','linewidth',3);
text(t1,10^-2,' 1.2 hr ','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','linewidth',2);
grid on;
ticker=get(gca,'xtick');
ylim([1.5*10^-3 1])
xlim([5*10^-5 2*10^-3])

%%
%figure;
subplot(2,3,4)
%[~,tidx]=min(abs(GDhobotime(GDidx)-datenum(2010,12,2,9,0,0)));

[~,tidx]=min(abs(GDhobotime(GDidx)-datenum(2010,12,2,5,0,0)));
plot(HKGD(GDidx(tidx),:),GDhgts,'linewidth',2);
hold on;
plot(GDtemp(GDidx(tidx),:),GDhgts,'r','linewidth',2);
legend('Harkers','Grandeur')
xlim([-6 6]);
title(datestr(GDhobotime(GDidx(tidx))),'fontsize',10,'fontweight','bold')
grid on;
ylim([min(GDhgts) 2500]);
set(gca,'fontsize',10,'fontweight','bold','box','on','linewidth',2);
xlabel('Temp [ ^{\circ}C ]','fontsize',10,'fontweight','bold')
ylabel('Height [m]','fontsize',10,'fontweight','bold')

%% 
%figure;
subplot(2,3,6);
%[~,tidx]=min(abs(GDhobotime(GDidx)-datenum(2010,12,2,8,00,0)));

[~,tidx]=min(abs(GDhobotime(GDidx)-datenum(2010,12,2,5,40,0)));
plot(HKGD(GDidx(tidx),:),GDhgts,'-b','linewidth',2);
hold on;
plot(GDtemp(GDidx(tidx),:),GDhgts,'-r','linewidth',2);
legend('Harkers','Grandeur')

xlim([-6 6]);
title(datestr(GDhobotime(GDidx(tidx))),'fontsize',10,'fontweight','bold')
grid on;
ylim([min(GDhgts) 2500]);
set(gca,'fontsize',10,'fontweight','bold','box','on','linewidth',2,'yaxislocation','right');
xlabel('Temp [ ^{\circ}C ]','fontsize',10,'fontweight','bold')
ylabel('Height [m]','fontsize',10,'fontweight','bold')

return
%% now load RWP data to look at across basin flows

f2='PCAPS_RWP.nc.nc';
nc_dump(f2);
thetaRWP=nc_varget(f2,'ThetaRWP');
uRWP=nc_varget(f2,'uRWP');
vRWP=nc_varget(f2,'vRWP');
tRWP=nc_varget(f2,'timeRWP');
hRWP=nc_varget(f2,'heightRWP');
rwptidx=find(tRWP>=starter & tRWP <= ender);

figure(60);clf
cmap=[gray(10);flipud(gray(10))];
cmap(1:10,3)=1;
cmap(11:end,1)=1;
colormap(cmap);
pcolor(uRWP(1:10,rwptidx));shading interp; caxis([-2.5 2.5]);colorbar;
hold on;
contour(thetaRWP(1:10,rwptidx),200:1:500,'-w');

figure(24); hold on;
plot(tRWP(rwptidx),detrend(mean(uRWP(4:8,rwptidx),1)),'r');datetick('x')