%% script to compute pressures along a constant heigh surface in the SLV. 

close all;
clear all;

R=287;
g=9.81;
fname='ISFS1.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));
latb=str2double(ST(3));
lonb=str2double(ST(4));
base=elev;

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15; basetemp=temp;
pres=data(:,13);
pref=pres;


figure(1); clf
subplot(2,1,1); hold on;
hf=fill([datenum(2010,12,2,14,33,0),datenum(2010,12,3),datenum(2010,12,3),datenum(2010,12,2,14,33,0)],[880 880 860 860],[1 1 200/255],'linestyle','none');
hf=fill([datenum(2010,12,3,14,33,0),datenum(2010,12,4),datenum(2010,12,4),datenum(2010,12,3,14,33,0)],[880 880 860 860],[1 1 200/255],'linestyle','none');

%set(hf,'facealpha',.5);
plot(times,pres,'k','linewidth',2);
grid on; 
datetick('x');
subplot(2,1,2)
plot(times,temp,'k','linewidth',2);
grid on; hold on;
datetick('x');
ISFS1pres=pres; ISFS1times=times;



%

fname='ISFS7.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));
dz=base-elev;
lat=str2double(ST(3));
lon=str2double(ST(4));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1


data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=temp;%(temp+basetemp)/2;
pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));

figure(1)
subplot(2,1,1)
%plot(times,pres);hold on;
plot(times,PR,'-b','linewidth',2);
subplot(2,1,2)
plot(times,temp,'linewidth',1); hold on;

figure(2)
hold on;
%plot(times,(pref-PR)./dx);
datetick('x')

ISFS7pres=PR;
ISFS7times=times;

%%
fname='ISFS2.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(6));
dz=base-elev;
lat=str2double(ST(4));
lon=str2double(ST(5));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=(temp+basetemp)/2;

pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1)
%plot(times,pres,'r');hold on;
plot(times,PR,'--r');
subplot(2,1,2)
plot(times,temp,'r'); hold on;



ISFS2pres=PR;
ISFS2times=times;

figure(2)
hold on;
%plot(times,(PR-pref)./dx,'r');
datetick('x')
%% 
fname='ISFS4.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(6));
dz=base-elev;
lat=str2double(ST(4));
lon=str2double(ST(5));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=(temp+basetemp)./2;
pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1)
%plot(times,pres,'c');hold on;
plot(times,PR,'--c');
subplot(2,1,2)
plot(times,temp,'c'); hold on;


ISFS4pres=PR;
ISFS4times=times;

figure(2)
hold on;
%plot(times,(PR-pref)./dx,'c');
datetick('x')

 
%% 
fname='ISFS5.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(6));
dz=base-elev;
lat=str2double(ST(4));
lon=str2double(ST(5));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=(temp+basetemp)./2;
pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1)
%plot(times,pres,'c');hold on;
plot(times,PR,'--c');
subplot(2,1,2)
plot(times,temp,'c'); hold on;


ISFS5pres=PR;
ISFS5times=times;

figure(2)
hold on;
%plot(times,(PR-pref)./dx,'c');
datetick('x')


%% 
fname='ISFS6.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(6));
dz=base-elev;
lat=str2double(ST(4));
lon=str2double(ST(5));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=(temp+basetemp)./2;
pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1)
%plot(times,pres,'c');hold on;
plot(times,PR,'--c');
subplot(2,1,2)
plot(times,temp,'c'); hold on;


ISFS6pres=PR;
ISFS6times=times;

figure(2)
hold on;
%plot(times,(PR-pref)./dx,'c');
datetick('x')

%% 
fname='ISFS3.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};

elev=str2double(ST(5));
dz=base-elev;
lat=str2double(ST(3));
lon=str2double(ST(4));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

temp=data(:,7)+273.15;
tmean=(temp+basetemp)./2;
pres=data(:,13);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1)
%plot(times,pres,'c');hold on;
plot(times,PR,'--c');
subplot(2,1,2)
plot(times,temp,'c'); hold on;


ISFS3pres=PR;
ISFS3times=times;

figure(2)
hold on;
%plot(times,(PR-pref)./dx,'c');
datetick('x')

%% 
tref=times;
fname='ISSSF.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));
dz=base-elev;
lat=str2double(ST(3));
lon=str2double(ST(4));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1
data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));
datestr(times(1:5))
temp=data(:,7)+273.15;
temp=interp1(times,temp,tref);
tmean=temp;%(temp+basetemp)./2;
pres=data(:,11);
pres=interp1(times,pres,tref);

PR=pres.*(exp((-g./(R*tmean))*dz));
figure(1)
subplot(2,1,1); hold on;
%plot(tref,pres,'m');hold on;
plot(tref,PR,'-m','linewidth',2);
xlim([datenum(2010,12,3) datenum(2010,12,3,15,0,0)])


ISSpres=PR;
ISStimes=tref;

subplot(2,1,2); hold on;
plot(tref,temp,'m','linewidth',2); hold on;
xlim([datenum(2010,12,3) datenum(2010,12,3,15,0,0)])

figure(2)
hold on;
pgrad=(pref-PR)./dx;
[b,c]=butter(5,1./5);
pgradsm=filtfilt(b,c,pgrad);

plot(tref,pgradsm,'m','linewidth',2);
hold on;
plot(tref,pgrad,'-k','linewidth',1)
datetick('x')


%% 
fname='ISS3.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 

ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));
dz=base-elev;
lat=str2double(ST(3));
lon=str2double(ST(4));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));
datestr(times(1:5))
temp=data(:,7)+273.15;
temp=interp1(times,temp,tref);
tmean=(temp+basetemp)./2;
pres=data(:,11);
pres=interp1(times,pres,tref);


PR=pres.*(exp((-g./(R*tmean))*dz));


figure(1)
subplot(2,1,1); hold on;

%plot(tref,pres,'g');hold on;
plot(tref,PR,'-g','linewidth',2);
ylim([868.5 878.5])
%datetick('x')
grid on;
xticks=[datenum(2010,12,2,18,0,0):3/24:datenum(2010,12,3,18,0,0)];
set(gca,'xtick',xticks,'xticklabel',{datestr(xticks,'HH mmm-DD')});
xlim([datenum(2010,12,2,18,0,0) datenum(2010,12,3,18,0,0)])
set(gca,'box','on','layer','top')

ISS3pres=PR;
ISS3times=tref;

subplot(2,1,2)
plot(tref,temp,'g'); hold on;
set(gca,'layer','top','box','on');
grid on;
xticks=[datenum(2010,12,2,18,0,0):3/24:datenum(2010,12,3,18,0,0)];
set(gca,'xtick',xticks,'xticklabel',{datestr(xticks,'HH mmm-DD')});
xlim([datenum(2010,12,2,18,0,0) datenum(2010,12,3,18,0,0)])
set(gca,'box','on','layer','top')

figure(2)
hold on;
%plot(tref,(PR-pref)./dx,'g');
%datetick('x')
plot(times,0*times,'--k');
grid on;
xticks=[datenum(2010,12,2,0,0,0):3/24:datenum(2010,12,3,18,0,0)];
set(gca,'xtick',xticks,'xticklabel',{datestr(xticks,'HH mmm-DD')});
xlim([datenum(2010,12,2,0,0,0) datenum(2010,12,3,18,0,0)])
set(gca,'box','on','layer','top')


%% FPS site
%PARM = MON,DAY,YEAR,HR,MIN,TMZN,TMPF,RELH,SKNT,GUST,DRCT,QFLG,PRES,SOLR,P05I,TSOI,TGND,NETS,NETL,TSNC,WSNC,USTD,VSTD,WSTD,TSTD,HFLX,USTR,WRAT,OTOT,WARN,SPDS,DIRS,VOLT,DWPF

fname='FPS.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 

ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(7));
dz=base-elev;
lat=str2double(ST(5));
lon=str2double(ST(6));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));
datestr(times(1:5))
temp=data(:,7)+273.15;
temp=interp1(times,temp,tref);
tmean=(temp+basetemp)./2;
pres=data(:,13);
pres=interp1(times,pres,tref);


PR=pres.*(exp((-g./(R*tmean))*dz));


figure(1)
subplot(2,1,1); hold on;

%plot(tref,pres,'g');hold on;
plot(tref,PR,'--m','linewidth',4);
ylim([868.5 878.5])
%datetick('x')
grid on;
xticks=[datenum(2010,12,2,18,0,0):3/24:datenum(2010,12,3,18,0,0)];
set(gca,'xtick',xticks,'xticklabel',{datestr(xticks,'HH mmm-DD')});
xlim([datenum(2010,12,2,18,0,0) datenum(2010,12,3,18,0,0)])
set(gca,'box','on','layer','top')

FPSpres=PR; FPSpres(end)=FPSpres(end-1);
FPStime=tref;


%% Filter out 12 hr+ signal

nhr=6;
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

nhr2=.25;
w2=1./(nhr.*12);
wn2=w2./.5;
[b2,c2]=butter(5,wn2,'low');

PRprime=zeros(5,length(tref));
PRprime(1,:)=filtfilt(b,c,ISFS1pres);
%PRprime(1,:)=filtfilt(b,c,PRprime(1,:));
PRprime(2,:)=filtfilt(b,c,ISFS2pres);
PRprime(3,:)=filtfilt(b,c,ISFS3pres);
PRprime(4,:)=filtfilt(b,c,ISFS4pres);
PRprime(5,:)=filtfilt(b,c,ISFS5pres);
PRprime(6,:)=filtfilt(b,c,ISFS6pres);
PRprime(7,:)=filtfilt(b,c,ISFS7pres);
PRprime(8,:)=filtfilt(b,c,ISSpres);
PRprime(9,:)=filtfilt(b,c,ISS3pres);
PRprime(10,:)=filtfilt(b,c,FPSpres);



%now low pass it to get rid of some features
[b2,c2]=butter(5,1/2,'low');
% for i=1:10
%    PRprime(i,:)=filtfilt(b2,c2,PRprime(i,:)); 
% end

%tidx=find(tref>=datenum(2010,12,2) & tref<=datenum(2010,12,2,15,0,0));
%PRprime=PRprime(:,tidx);
%tref=tref(tidx);

f=figure(300);clf;set(f,'color','w')
title('North South Transect');
plot(tref,PRprime(1,:),'b','linewidth',2);
hold on;
%plot(tref,PRprime(2,:),'k','linewidth',2);
%plot(tref,PRprime(3,:),'b');
%%plot(tref,PRprime(4,:)+.5,'g');
%plot(tref,PRprime(5,:)+.5,'k');
%plot(tref,PRprime(6,:),'b');
plot(tref,PRprime(8,:)+.5,'r','linewidth',2);
plot(tref,PRprime(7,:)+1,'c','linewidth',2);
%plot(tref,PRprime(9,:),'b');
%plot(tref,PRprime(10,:)+1.5,'m');
xlim([datenum(2010,12,2,4,0,0) datenum(2010,12,2,18,0,0)]);
datetick('x','keeplimits');
%legend('ISFS1-Playa (North)','ISFS4-West Valley','ISS-Center','ISFS7- Riverton');
legend('ISFS1-Playa (North)','ISS-Center','ISFS7- Riverton');
set(gca,'fontsize',10,'fontweight','bold','linewidth',2)
xlabel('Time [UTC]','fontsize',10,'fontweight','bold');
ylabel('hPa','fontsize',10,'fontweight','bold');


figure(888);clf
title('East-West Transect');
hold on;
plot(tref,PRprime(3,:)-.5,'g');
plot(tref,PRprime(4,:)+.5,'g');

datetick('x')
xlim([datenum(2010,12,2,0,0,0) datenum(2010,12,3,0,0,0)]);
legend('ISFS3','ISFS4');

figure(301);clf;
title('East-West Transect');
%plot(tref,PRprime(1,:),'b','linewidth',2);
hold on;
%plot(tref,PRprime(2,:),'k','linewidth',2);
%plot(tref,PRprime(3,:)-.5,'g');
%plot(tref,PRprime(4,:)+.5,'g');
plot(tref,PRprime(5,:),'k');
plot(tref,PRprime(8,:)+.5,'r','linewidth',2);
plot(tref,PRprime(6,:)+1,'b');
%plot(tref,PRprime(7,:)+1.5,'c','linewidth',2);
%plot(tref,PRprime(9,:),'b');
%plot(tref,PRprime(10,:)+1.5,'m');
datetick('x')
xlim([datenum(2010,12,2,4,0,0) datenum(2010,12,3,0,0,0)]);
legend('ISFS5','ISS','ISFS6');

%clockwise rotation?
f=figure(302); clf; set(f,'color','w');
%sp1=subplot(2,1,1);
title('ISFS Pressure Perturbations','fontsize',10,'fontweight','bold');
%plot(tref,PRprime(1,:),'b','linewidth',2);
hold on;
plot(tref,PRprime(1,:)-.5,'b','linewidth',2);
plot(tref,PRprime(2,:),'r','linewidth',2);
plot(tref,PRprime(3,:)+.5,'k','linewidth',2);
plot(tref,PRprime(4,:)+1,'g','linewidth',2);
plot(tref,PRprime(5,:)+1.5,'m','linewidth',2);
plot(tref,PRprime(6,:)+2,'c','linewidth',2);
plot(tref,PRprime(7,:)+2.5,'y','linewidth',2);
legend('ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7','location','bestoutside');
box on;
grid on;
xlim([datenum(2010,12,2,0,0,0) datenum(2010,12,3,6,0,0)]);

datetick('x','keeplimits')

set(gca,'fontsize',10,'fontweight','bold','linewidth',1)
xlabel('Time [UTC]','fontsize',10,'fontweight','bold');
ylabel('hPa','fontsize',10,'fontweight','bold');
%spos=get(sp1,'pos');

return
%% Spectral analysis
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(PRprime)/L);
disp(nffts)

tidx=find(tref>=datenum(2010,12,2) & tref <= datenum(2010,12,3));

%apply the analysis to the detrended data... which is the PRprime data that
%has variations longer than 12 hours removed
figure(302)
sp2=subplot(2,1,2);
spos2=get(sp2,'pos');
spos2(3)=spos(3);
set(sp2,'pos',spos2);
%f=figure(900);clf; set(f,'color','w');
clrs={'-b','-r','-k','-g','-m','-c','-y','--k'};
for i=1:8;%10
[pxx,ff,pxxc]=pwelch(detrend(PRprime(i,tidx)),hanning(L),round(4*L/5),L,fs,'ConfidenceLevel',.95);
loglog(ff,ff.*pxx,clrs{i},'linewidth',2)
hold on;
xlabel('Hz','fontsize',10,'fontweight','bold');
ylabel('','fontsize',10,'fontweight','bold');
t12=1./(60*60*12); %12 hours
t3=1./(3600*3); %4hours
t1=1./(3600*1.2); %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 .12],'--k','linewidth',3);
plot([t3 t3],[10^-4 .12],'--k','linewidth',3);
text(t1,10^-3,' 1.2 hr ','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','linewidth',2);
grid on;
ticker=get(gca,'xtick');
ylim([10^-4 .12])
xlim([10^-5 .5*10^-2])
ylabel('Power Spectral Density (hpa^2/s)');
end
%legend('ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7','ISS');



return


%% 
fname='KOGD.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(6));
dz=base-elev;
lat=str2double(ST(4));
lon=str2double(ST(5));
dx=deg2km(distance(latb,lonb,lat,lon))*10^3;%distance in meters to ISFS1

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));
temp=data(:,7)+273.15;
tmean=(temp);
pres=data(:,14); pres(pres>1100)=nan;
pres=inpaint_nans(pres);
PR=pres.*(exp((-g./(R*tmean))*dz));
figure(800)
subplot(2,1,1)
%plot(times,pres,'c');hold on;
plot(times,PR,'--c');
subplot(2,1,2)
plot(times,temp,'c'); hold on;

w=1/(12*8);
wn=w./.5;
[b2,c2]=butter(5,wn,'high');

hatprime=filtfilt(b2,c2,pres);
figure; plot(times,hatprime,'linewidth',2); datetick('x');
hold on;
plot(tref,PRprime(8,:),'r');

