%%Program to calculate the radius of each crossing and
%%plot the radius vsBz strength
 clf
 filename = '/home/emery/tieconvec/1998Jul_1999Jun_Weimer05_CRBs.txt';
%filename = '/home/emery/tieconvec/1996Jan_Dec_Weimer05_CRBs.txt';
%filename = '/home/emery/tieconvec/9807_9906Boundaries_sw.txt';
%filename = '/hao/aim2/dmsp/dmsp_idm_4s/9807_9906Boundaries_sw.txt';
%filename = '/home/emery/tieconvec/1996Boundaries_sw.txt';
%filename = '/hao/aim2/dmsp/dmsp_idm_4s/1996Boundaries_sw_lat.txt';
 data = dlmread(filename,'');
 data(1:length(data),12) = NaN;
 
 n = 0;
 for i=1:2:length(data)

     D = (data(i+1,3) - data(i,3))*7.407;
     angle = (D*360)/(2*pi*(6378.1+850));
     data(i,12) = angle;
% Weimer problems (reverse convection in polar cap near noon for 96252 both orbits so max-min time is reversed
    if (data(i+1,3) < data(i,3))
       (data(i,:))
       (data(i+1,:))
% 1996 low lat problems
%      if (angle>60.)
%       angle
%       data(i+1,8)
%       data(i,8)     
     else
      n = n + 1;
% Bz
      Bzdata(n) = data(i,8);
% BzVsw
%      Bzdata(n) = data(i,8)*data(i,10)*0.001; 
      Diameterdata(n) = angle;
%     end
     end
        
 end
 
min_diam=min(Diameterdata)
max_diam=max(Diameterdata)
min_bz=min(Bzdata)
max_bz=max(Bzdata)
%  scatter((data(:,8).*data(:,10)*10^-3),data(:,12),20);
  scatter(Bzdata,Diameterdata,40,'ko','filled','LineWidth',3);

% axis for Bz
  xmin = -20;
% axis for BzVsw
% xmin = -12;
  axis([xmin 0 10 60])
  grid
  hold on

% Find the first order (1) polynomial fit y=p(2)+p(1)*x (x=Bz, y=Diameter)  
  p = polyfit(Bzdata,Diameterdata,1)
  y=polyval(p,(0:-0.5:-20));
  plot((0:-0.5:-20),y,'Color','k','LineWidth',3);

% Find the cross-correlation of x and y (use transpose with a quote to get in cols, not rows)
  cc_n = xcorr(Bzdata,Diameterdata,'coeff');
n
% cc is in middle of n*2-1, so for n=186 is 186 of 371 size 
% 08/11:  These values changed slightly - why??  - get corrcoeff working instead!
  cross_correl = cc_n(n)

% corrcoeff fails here but works in function coeff_ls ???
% [r,pr,rlo,rhi] = corrcoeff(Bzdata',Diameterdata');
% cross_correl = r(1,2)
%  probability = pr(1,2);
%  rlo05percent = rlo(1,2);
%  rhi95percent = rhi(1,2);
[coeff,prob,rlo05,rhi95,a,b,c,d,e,f,g,h,i,j,avy1,avy2,avr,mdy1,mdy2,mdr] = coeff_ls(Bzdata',Diameterdata');
  cross_correl = coeff
  cctitle = sprintf('r=%5.2f',cross_correl);
% NOTE, cross_correl from [r,p,rlo,rhi]=corrcoeff(x,y), where cc=r(1,2), is DIFFERENT from
%   midpoint of xcorr (and smaller - like -0.6000 instead of -0.8861 for 9807-9906 Bz IDM)
 
fsiz = 15;
bzp9 = xmin*0.9;

% text for cross-correlation for Bz 
 xlabel('Bz (nT)','fontsize',fsiz,'fontweight','bold')
% text for cross-correlation for BzVsw 
%xlabel('BzVsw (mV/m)','fontsize',fsiz,'fontweight','bold')

text (bzp9,15.,cctitle,'fontsize',fsiz,'fontweight','bold')

 ylabel('Crossing angle (degrees)','fontsize',fsiz,'fontweight','bold')

%title('Weimer05 1996','fontsize',fsiz,'fontweight','bold')
%fplot= ['1996_Angle_vs_Bz_W05'];
%title('F13 IDM 1996','fontsize',fsiz,'fontweight','bold')
%fplot= ['1996_Angle_vs_Bz_IDM'];
%fplot= ['1996_Angle_vs_BzVsw_IDM'];
% title('F13 IDM July 1998 to June 1999','fontsize',fsiz,'fontweight','bold')
% fplot= ['9807_9906_Angle_vs_Bz_IDM'];
% fplot= ['9807_9906_Angle_vs_BzVsw_IDM'];
  title('Weimer05 July 1998 to June 1999','fontsize',fsiz,'fontweight','bold')  
  fplot= ['9807_9906_Angle_vs_Bz_W05'];

% set(gca,'fontsize',fsiz,'fontweight','bold');
  set(gca,'xlim',[xmin 0],'ylim',[10 60],'fontsize',fsiz,'fontweight','bold');
  set(gcf,'PaperPositionMode','auto');
  print(gcf,'-dtiff',fplot);

   
