% level.m
% mfile to take each liwi grdfmt mvp file
% compare it to a mean sigma theta profile
%    made from a combination of mean mvp and mean hrp data
% and level densities in each mvp profile to match the mean
%  inputting from /j2/data/mvp/twist(abc)/grdfmt/grd*.mat
%  outputting to /j2/data/mvp/twist(abc)/levfmt/lev*.mat


[ABC]=input('   enter twist/liwi MVP directory a b or c :   ','s');
stas=input('Which stations do you want to process? eg. [1:248] for twistb  :   ');

% mean hrp mvp sigma file
% John wants to used one mean profile = leveling_hrp_mvp_mean.mat to level
% all sensors (originally we did each sensor to it's own extended mean)
%meanfile=(['/j2/data/mvp/twist_prgs/hrp_mvp_mean_sort_',ABC,'.mat']);

meanfile=(['/j2/data/mvp/twist_prgs/leveling_hrp_mvp_mean.mat']);
eval(['load ',meanfile]);
pgrid=pmean;

% input directory
directory = (['/j2/data/mvp/twist',ABC,'/grdfmt']);

% output directory
outdir=(['/j2/data/mvp/twist',ABC','/levfmt']);

% load each mvp profile again

for h=1:length(stas)
matfmtfile=[directory,'/grd',sprintf('%4.4i',stas(h)),'.mat'];
disp(['   working on file ',matfmtfile]);

if(exist(matfmtfile)==2);
eval(['load ',matfmtfile]);

%figure(1)
%plot(sigthave,pgrid,'b');
%if (h==1)
%  axis ij; hold on; grid on; zoom on;
%  xlabel('sigma');
%  ylabel('pressure');
%end

% do low pass filter here first !!!!!
% [b,a]=mybutter(order,1/frequency(m),half the sample rate)
% I did all this xx stuff to get sgtest to be the same length
% as all the others ...
%  but it still is NOT monotonic ... which makes the rest bomb.

xx=find(~isnan(sigthave));   % this same as j above
[b,a]=mybutter(4,1/20,.5);
[sgtest(xx)]=myfiltfilt(b,a,sigthave(xx)');
sgtest(1:(xx(1)-1))=sgtest(xx(1));
sgtest(xx(end)+1:length(sigthave))=sgtest(xx(end));
sgtest=sort(sgtest);
sgthmean=sort(sgthmean);

% do this to make them monotonic 
% so that no two consecutive values are the same

xxx=diff(sgtest);
yyy=diff(sgthmean);
fxxx=find(xxx==0);
fyyy=find(yyy==0);
for z=1:length(fxxx)
 sgtest(fxxx(z)+1)=sgtest(fxxx(z))+.00001;
end
for z=1:length(fyyy)
  sgthmean(fyyy(z)+1)=sgthmean(fyyy(z))+.00001;
end
tst_p=tab1(sgthmean,pgrid,sgtest,'nonan');

% but now tst_p is still not monotonic
% so tab1 will not work
% do this to make tst_p monotonic

tst_p=sort(tst_p);
jjj=diff(tst_p);
fjjj=find(jjj==0);
for z=1:length(fjjj)
 tst_p(fjjj(z)+1)=tst_p(fjjj(z))+.00001;
end

displacement=pgrid-tst_p;
plev=tst_p;
ulev=tab1(tst_p,uave,pgrid,'nonan');
vlev=tab1(tst_p,vave,pgrid,'nonan');
wlev=tab1(tst_p,wave,pgrid,'nonan');
thetalev=tab1(tst_p,thetave,pgrid,'nonan');
slev=tab1(tst_p,s_ave,pgrid,'nonan');



figure(1)
plot(sigthave,tst_p,'k');
if h<10
  axis ij; hold on; grid on; zoom on;
end

xlabel('sigma theta');
ylabel('pressure');
title('Leveled MVP sigma theta profiles');

%figure(4)
%plot(new_disp);
%grid on; hold on; zoom on;


% Now I need to OUTPUT these new files
% One per station per instrument (abc)

outfile=[outdir,'/lev',sprintf('%4.4i',stas(h)),'.mat'];
disp(['    outputting to file : ',outfile,sprintf('\n')]);

eval(['save ',outfile,' pgrid displacement plev ulev vlev wlev thetalev slev startdaytime stopdaytime']);



end
end
