%%Program to calculate the radius of each crossing and
%%plot the radius vsBz strength
 clf
 filename = '/home/emery/tieconvec/1998Jul_1999Jun_Weimer05_CRBs.txt';
 data = dlmread(filename,'');
 data(1:length(data),12) = NaN;
 
 for i=1:2:length(data)
     if (data(i+1,3) < data(i,3))
         (data(i,:))
         (data(i+1,:))
     end
     D = (data(i+1,3) - data(i,3))*7.407;
     angle = (D*360)/(2*pi*(6378.1+850));
     data(i,12) = angle;
     
 end
 

  scatter((data(:,8).*data(:,10)*10^-3),data(:,12),20);            %If you wantBz*SWV
  %scatter(data(:,8),data(:,12),20);                                 %If you want Bz only
  axis([-11 0 10 50])
  hold on
  
  B = ((data(:,8).*data(:,10)*10^-3));
  A = (data(:,12));
  
  %p = [-1.7221, 22.4540];
  p = polyfit(B,A,1)
  y = polyval(p,(0:-1:-20));
  plot((0:-1:-20),y)
  
  pause
  %xlabel('Bz strength(nT)')
  xlabel('Bz*VSW strength(mV/m)')
  ylabel('Crossing angle (degrees)')
  title('Angle vs. Bz*SWV strength for July 1998 through June 1999 Weimer Model 2005')
  %title('Angle vs. Bz*SWV strength for 1996 Weimer Model 2005')
  fplot= ['/hao/aim2/dmsp/dmsp_idm_4s/' '9807_9906Angle_vs_BzSWVWeimer'];
  set(gcf,'PaperPositionMode','auto');
  print(gcf,'-dtiff',fplot);

   