% hrp_mean_sigma_2.m
% to create a mean LIWI hrp sigma theta on pressure profile
%
% this step separated off of mags test.m
% by John  Nov 2000

directory1 = (['/n6/hrpdata/oc324/reprocessed/matfmt/']);
% for matching twistb the following stas1 is appropriate (all)
% stas1=[7:13 40:74 84:90 92:144 146:153 164:170 172:218 220:224 241:246 ];
% for matching twista stas1= [49:54 55:79] works well (55 has spike near surface)
% for matching twistc stas1= [52:54 56:100] works well 


[stas1]=input('Enter HRP stations you want to include in the mean  :');
[ABC]=input('Enter the LIWI sensor (a,b,c) you are matching to :','s');

disp('Now we will load and plot the HRP stations');
disp(['press return to continue',sprintf('\n')]);
pause

for h=1:length(stas1)
matfmtfile1=[directory1,'cmb',sprintf('%3.3i',stas1(h)),'.mat'];
disp(['   working on file ',matfmtfile1]);

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

% plot each individual hrp sigma profile

figure(2)
plot(SGTH,PRES);

if(h==1)
  axis ij;
  hold on;
  zoom on;
  grid on;
DUMMY=[10:.5:5000];
NPROFS=0*DUMMY;
SUMSGTH=NPROFS;
SUMSGTH2=NPROFS;
SUMDPDT=NPROFS;
end

% sum up only the scans of each profile that have good SGTH data
% John changed this ... mag was doing this on good PRES data

j=find(~isnan(SGTH));

for k=1:length(j)
l=1+2*(PRES(j(k))-10);
NPROFS(l)=NPROFS(l)+1;
SUMSGTH(l)=SUMSGTH(l)+SGTH(j(k));
SUMSGTH2(l)=SUMSGTH(l)+SGTH(j(k))*SGTH(j(k));
end; %k loop

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


% divide sums by # of obs
kk=find(NPROFS~=0);
SIGTHAVE(kk)=SUMSGTH(kk)./NPROFS(kk);
SGTHVAR(kk)=SUMSGTH2(kk)./NPROFS(kk)-SIGTHAVE(kk).*SIGTHAVE(kk);
SGTHSTD(kk)=sqrt(SGTHVAR(kk));SGTHERR(kk)=SGTHSTD(kk)./(sqrt(NPROFS(kk)-1));

%rename output variables
SGTHMEAN=SIGTHAVE(kk);
PMEAN=DUMMY(kk);
SGTHERR=SGTHERR(kk);
SGTHVAR=SGTHVAR(kk);


disp([sprintf('\n'),'Now we will plot the mean HRP profile']);
disp(['press return to continue',sprintf('\n')]);
pause



% plot mean hrp sigma profile with error bars
figure(2)
plot(SGTHMEAN,PMEAN,'c');axis ij;hold on; zoom on;
plot(SGTHMEAN+SGTHERR,PMEAN,'g');plot(SGTHMEAN-SGTHERR,PMEAN,'g');
ylabel('Pressure');
xlabel('Sigma Theta');
title('Input HRP stations = blue Mean = Cyan Error bars = Green');

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

eval(['save ',outfile,' SGTHMEAN PMEAN']);
