% mvp_mean_sigma.m
%
% first stage in leveling liwi mvp data
% creating a time averaged sigma theta on pressure profile
% of twist/liwi (A B or C) MVP data
% john separated off steps from mags test.m   Nov 2000

%It is an editted version of John's average_profiles_4.m
% routine to average a set of MVP grdfmt data profiles
%this version only sums the pressure bins that have data

% make a mean mvp sigma theta profile
disp(['Now we''re going to make a mean LIWI/TWIST mvp sigma theta profile',sprintf('\n')]);

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

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

% load each mvp profile

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

% plot each individual mvp sigma profile

figure(1)
plot(sigthave,pgrid);

if(h==1)
axis ij
hold on
zoom on
grid on
nprofs=0*pgrid;
sumsgth=nprofs;
sumsgth2=nprofs;
sumdpdt=nprofs;
end

% sum up only the scans of each profile that have data

j=find(~isnan(pave));

nprofs(j)=nprofs(j)+1;

sumsgth(j)=sumsgth(j)+sigthave(j);
sumdpdt(j)=sumdpdt(j)+dpdtave(j);
sumsgth2(j)=sumsgth2(j)+sigthave(j).*sigthave(j); 

end; %if exist
end; %for h=1:length(stas)

disp([sprintf('\n'),'some divide by zeros here where no data in any profile']);
disp(['press return to continue',sprintf('\n')]);

% divide sums by # of obs
pause

sigthave=sumsgth./nprofs;
sgthvar=sumsgth2./nprofs-sigthave.*sigthave;
sgthstd=sqrt(sgthvar);sgtherr=sgthstd./(sqrt(nprofs-1));
dpdtave=sumdpdt./nprofs;

%rename output variables

sgthmean=sigthave;
pmean=pgrid;
dpdtmean=dpdtave;

disp([sprintf('\n'),'Here we will plot mean mvp profile in figure 1']);
disp(['press return to continue',sprintf('\n')]);
pause


% plot mean MVP sigma profile with error bars

figure(1)
plot(sgthmean,pmean,'r');
plot(sgthmean+sgtherr,pmean,'m');plot(sgthmean-sgtherr,pmean,'m');
ylabel('Pressure');
xlabel('Sigma Theta');
title('Input MVP stations = blue Mean = Red Error bars = Magenta  With HRPMean Extension = black *');

outfile=(['/j2/data/mvp/twist_prgs/mvp_mean_sigma_',ABC,'.mat']);


eval(['save ',outfile,' sgthmean pmean']);

