%Read and Plot Boundary Files
close all

% Boundary files have format: 
%year, doy, UT(sec),Vy(km/s),MLT,MLAT,GLAT,Bz,By

% ifwei=1,0 if do,not read the Weimer 2005 files which contain 3 extra params
ifwei = 1;
% 8/11: Problems because crad too small and shift from NH By<0 (duskbanana) to >0
% (dawnbanana) is to dusk and NOT to dawn as previous studies of NH W05 showed for By+/-15

date_title = sprintf('Weimer05 3y');
%date_title = sprintf('SH Weimer05 3y');
%date_title = sprintf('DMSP-F13 3y');
%date_title = sprintf('1996');
%date_title = sprintf('Jul 1998 - Jun 1999');
%date_title = sprintf('2005');

% if fixvsw=1, then fix Vsw and Dsw (and sometimes year) and rewrite revname
fixvsw = 0;

% bmid is Bz value above and below which red or blue is used in polar plots
% 50% of Bz values > and < Bz=-3.25nT
bmid = -3.25;

% ifplotm = 1,0 if do,not plot MLT vs MLAT plots
% ifplotp = 1,0 if do,not plot polar plots
% ifplotb = 1,0 if do,not plot polar plots with fitted circles for 5 neg Bz bins
ifplotm = 0;
ifplotp = 0;
ifplotb = 1;

for ihem=1:2

%fpath = '/hao/aim2/dmsp/dmsp_idm_4s/'; %directory path
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/9807_9906Boundaries_sw'};
% new write for |Bz/By|>1 (old format is ltempS.txt)
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/9601/lS139601Boundaries_Bz'};
if (ihem==1)
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yN'};
fname = {'/home/emery/tieconvec/Boundaries_W05_3yN'};
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Bz_1996N'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_1996Nb.txt';
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_98182_99181_N'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_9899Nb.txt';
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_2005_N'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_2005N.txt';
else
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_3yS'};
fname = {'/home/emery/tieconvec/Boundaries_W05_3yS'};
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Bz_1996S'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_1996Sb.txt';
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_98182_99181_S'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_9899Sb.txt';
%fname = {'/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_2005_S'};
%revname = '/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_Vsw_2005S.txt';
end
% {[fpath '9607/l139607Boundaries_Bz'],[fpath '9608/l139608Boundaries_Bz'],... %here put all of the names&paths
%     [fpath '9609/l139609Boundaries_Bz'],[fpath '9610/l139610Boundaries_Bz'],...      %of the files you want plotted
%     [fpath '9611/l139611Boundaries_Bz'],[fpath '9612/l139612Boundaries_Bz'],...
%     [fpath '9601/l139601Boundaries_Bz'],[fpath '9602/l139602Boundaries_Bz'],...
%     [fpath '9603/l139603Boundaries_Bz'],[fpath '9604/l139604Boundaries_Bz'],...
%     [fpath '9605/l139605Boundaries_Bz'],[fpath '9606/l139606Boundaries_Bz']};

if ifplotm==1
 figure;
 colormap(jet)
end
num = 0;
numarr = [];
Bvec = [];
%subplot(2,1,1)

	h2deg=15.;

for i = 1:length(fname)
    str = char(fname(i));
    filename = strcat( str, '.txt');
    data = dlmread(filename,'');
%   MLT = data(:,5);
%   MLAT = data(:,6);
%   B = sqrt((data(:,8).^2) + (data(:,9).^2) );
% for new format
    B = sqrt((data(:,3).^2) + (data(:,4).^2) );
    Bz = data(:,3);
    By = data(:,4);
    Vsw = data(:,5);
    Dsw = data(:,6);
% for new CRB MLT/MLAT
    MLT = data(:,16);
    MLAT = data(:,15);
% for new peak sunward convection MLT/MLAT
    MLTx = data(:,10);
    MLATx = data(:,9);
    Vyx = data(:,8);
    UTx = data(:,7);
    daynx = data(:,2);
if ifplotm==1
    scatter(MLAT,MLT,40,B,'filled')
    hold on
    scatter(MLATx,MLTx,'+k')
    hold on 
end
    Bvec = [Bvec; B];
    n = length(data);
    numarr(i) = n;
    num = num + n;
end
% Convert mlat and mlt to rho (colat) and theta (angle)
nidm(ihem) = num;
for n=1:nidm(ihem)
%if abs(MLAT(n))>80. 
%n,daynx(n),MLAT(n)
% if (Vyx(n) < 0.1)
%n,Vyx(n)
%end
 byidm(n,ihem) = By(n);
 bzidm(n,ihem) = Bz(n);
 badvsw(n) = 0;
%  Check for bad Vsw or Dsw (fix in newCRB.m instead?)
  if Vsw(n)<200 | Vsw(n)>1000 | Dsw(n)<0 | Dsw(n)>50
   badvsw(n) = 1;
  end
 mlatidm(n,ihem) = 90. - abs(MLAT(n));
 mltidm(n,ihem) = (MLT(n)*h2deg-90.)*pi/180.;
 mlatidmx(n,ihem) = 90. - abs(MLATx(n));
 mltidmx(n,ihem) = (MLTx(n)*h2deg-90.)*pi/180.;
%  Read 4 extra parameters ifwei=1: 
%   #19ctpoten(kV), and in deg: radius from bndyfit/2,offc(+0MLT),dskofc(+6MLT) 
  if ifwei==1
    ctpoten(n,ihem) = data(n,19);
    crad(n,ihem) = data(n,20);
    offc(n,ihem) = data(n,21);
    dskofc(n,ihem) = data(n,22);
  end
end

% Linearly interpolate badvsw values and rewrite into revname
if (fixvsw==1)
for n=1:2:num-1
% change y2 to year for pair of Boundaries
  if (data(n,1)>50 && data(n,1)<100)
   data(n,1) = 1900 + data(n,1);
   data(n+1,1) = 1900 + data(n+1,1);
  end
  if (badvsw(n)==1) 
    nn=n+1;
    while badvsw(nn)==1
     nn = nn + 1;
    end
% badvsw(nn-1)=1, badvsw(nn)=0 and badvsw(n-1)=0
    dayfr1 = data(n-1,13)/86400.;
    dayfr2 = data(nn,2)-data(n-1,2) + data(nn,13)/86400.;
    dayfr = data(n,2)-data(n-1,2) + data(n,13)/86400.;
    fac1 = (dayfr2-dayfr)/(dayfr2-dayfr1);
    fac2 = (dayfr-dayfr1)/(dayfr2-dayfr1);
    data(n,5) = data(n-1,5)*fac1 + data(nn,5)*fac2;
    data(n,6) = data(n-1,6)*fac1 + data(nn,6)*fac2; 
% Set next pair to same value
    data(n+1,5) = data(n,5);
    data(n+1,6) = data(n,6);
    badvsw(n+1) = 0;
  end
end
dlmwrite(revname,data,'delimiter',' ','-append');
end

if ifplotm==1
caxis([0 10])
 %axis([62 81 5 18.5])
xlabel('MLAT','fontsize',15,'fontweight','bold')
ylabel('MLT','fontsize',15,'fontweight','bold')
mtitle = sprintf('%s Boundaries (o=CRB, +=sunward peak)',date_title);
% title('1996 Boundaries','fontsize',15,'fontweight','bold')
% title('1996 Boundaries (o=CRB, +=sunward peak)','fontsize',15,'fontweight','bold')
%title('1998 Jul-Dec Boundaries (o=CRB, +=sunward peak)','fontsize',15,'fontweight','bold')
title(mtitle,'fontsize',15,'fontweight','bold')
colorbar
grid
set(gca,'fontsize',15,'fontweight','bold');
hold off
end	% ifplotm==1

% subplot(2,1,2)
% a = 1:12;
% bar(numarr)
% set(gca,'XTickLabel',{'Jul','Aug','Sep','Oct','Nov','Dec','Jan','Feb','Mar','Apr','May','Jun'})
% ylabel('number of crossings')
% title(strcat('total crossings:',num2str(num)))

% This is the histogram of the B values!
%subplot(2,1,2)
% [n,xout] = hist(Bvec,20);
% bar(xout,n)
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w')
% xlabel('B Value')
% ylabel('number of crossings with each value')
% title('B values of crossings from July 1999 - June 1999')
%colorbar

%pause
%fplot= ['/hao/aim2/dmsp/dmsp_idm_4s/0698_0799Boundaries'];
%fplot= ['/hao/aim2/dmsp/dmsp_idm_4s/9601/S9601Boundaries'];
%fplot= ['/hao/aim2/dmsp/dmsp_idm_4s/Boundaries_1996N'];
%set(gcf,'PaperPositionMode','auto');
%print(gcf,'-dtiff',fplot);

end		% for ihem=1:2 reading Boundaries for both hemispheres


% 8/11 Calculate circle fits for CRB and sunward max Vy values for IDM (and Weimer05)
%	Use 5 (~20%) neg Bz bins: 0to-2,-2to-3,-3to-4,-4to-6,<-6(~15%)
if ifplotb==1
Bzmin = [0 -2 -3 -4 -6];
Bzmax = [-2 -3 -4 -6 -23];
for nbzbin=1:5
 k=nbzbin;
 title_label = sprintf('%s %d<Bz<%d nT',date_title,Bzmin(nbzbin),Bzmax(nbzbin));
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     titr=pi/4;
     xr=65*cos(titr);
     yr=65*sin(titr);
     titb=-3*pi/4;
     xb=65*cos(titb);
     yb=65*sin(titb);
 hf = figure;
 set(hf,'position',[162    41   700   906]); 
%subplot('position',[.1,.1,.80,.80]);
 for ncase=1:4
 clear thetamlt;
 clear rhomlat;
 clear xs;
 clear ys;
 clear bymag;
 clear cp;
 clear thetacrad;
 clear xdksofc;
 clear yoffc;
 n = 0;
 if (ncase==1)
  subplot('position',[.05,.50,.40,.40]);
% Print out title etc
    hold on;
     text(1.15*xr,1.15*yr,title_label,'FontName','Times','FontSize',16,'color','k', 'HorizontalAlignment','Center');
     text(0.85*x,0.85*y,'CRBs','FontName','Times','FontSize',16,'color','b');
%     text(xr,yr,'By<0N By>0S','FontName','Times','FontSize',16,'color','k','HorizontalAlignment','Center');
 end
 if ncase==2
  subplot('position',[.50,.50,.40,.40]);
    hold on;
     text(0.85*x,0.85*y,'maxVys','FontName','Times','FontSize',16,'color','r');
 end
 if ncase==3
  subplot('position',[.05,.05,.40,.40]);
     text(0.85*x,0.85*y,'CRBs','FontName','Times','FontSize',16,'color','b');
%     text(xr,yr,'By>0N By<0S','FontName','Times','FontSize',16,'color','k','HorizontalAlignment','Center');
 end
 if ncase==4
  subplot('position',[.50,.05,.40,.40]);
    hold on;
     text(0.85*x,0.85*y,'maxVys','FontName','Times','FontSize',16,'color','r');
 end

      PolarBackgroundk;
    for ihem=1:2	% ihem=1N,2S
      for i=1:nidm(ihem)
	if (bzidm(i,ihem)<Bzmin(nbzbin) && bzidm(i,ihem)>=Bzmax(nbzbin))
	ifgo = 0;
	if (byidm(i,ihem)<0 && ihem==1) 
	 ifgo = 1;
        end
	if (byidm(i,ihem)>=0 && ihem==2) 
	 ifgo = 1;
        end
	 if (ifgo==1)
	  if (ncase==1)
	   polar(mltidm(i,ihem),mlatidm(i,ihem),'b.');		% CRB
   	   n = n + 1;
	   thetamlt(n) = mltidm(i,ihem);
	   rhomlat(n) = mlatidm(i,ihem);
	   bymag(n) = abs(byidm(i,ihem));
	   if ifwei==1
	     cp(n) = ctpoten(i,ihem);
	     thetacrad(n) = crad(i,ihem);
	     xdskofc(n) = dskofc(i,ihem);
	     yoffc(n) = -offc(i,ihem);
	   end
	  elseif ncase==2
	   polar(mltidmx(i,ihem),mlatidmx(i,ihem),'rx');	% maxVy
   	   n = n + 1;
	   thetamlt(n) = mltidmx(i,ihem);
	   rhomlat(n) = mlatidmx(i,ihem);
	   bymag(n) = abs(byidm(i,ihem));
	  end
         end	% ifgo==1
	ifgo = 0;
	if (byidm(i,ihem)>=0 && ihem==1) 
	 ifgo = 1;
        end
	if (byidm(i,ihem)<0 && ihem==2) 
	 ifgo = 1;
        end
	 if (ifgo==1)
	  if (ncase==3)
	   polar(mltidm(i,ihem),mlatidm(i,ihem),'b.');		% CRB
   	   n = n + 1;
	   thetamlt(n) = mltidm(i,ihem);
	   rhomlat(n) = mlatidm(i,ihem);
	   bymag(n) = abs(byidm(i,ihem));
	   if ifwei==1
	     cp(n) = ctpoten(i,ihem);
	     thetacrad(n) = crad(i,ihem);
	     xdskofc(n) = dskofc(i,ihem);
	     yoffc(n) = -offc(i,ihem);
	   end
	  elseif ncase==4
	   polar(mltidmx(i,ihem),mlatidmx(i,ihem),'rx');	% maxVy
   	   n = n + 1;
	   thetamlt(n) = mltidmx(i,ihem);
	   rhomlat(n) = mlatidmx(i,ihem);
	   bymag(n) = abs(byidm(i,ihem));
	  end
         end	% ifgo==1
        end	% if Bz conds
      end	% i=1:nidm(ihem)
    end		% ihem=1(N):2(S)
  % Convert positions of CRB or maxVy from polar to cartesian
  [xs,ys]=pol2cart(thetamlt,rhomlat);
  % Find best fit circle for positions
  [xfitM(nbzbin),yfitM(nbzbin),RfitM(nbzbin)] = circfit(xs,ys);
  % Convert back to polar
  [MLTcentM(nbzbin),LATcentM(nbzbin)]=cart2pol(xfitM(nbzbin),yfitM(nbzbin));
  % Plot best fit circle
   if ncase==1 | ncase== 3
      rectangle('position',[xfitM(k)-RfitM(k),yfitM(k)-RfitM(k),RfitM(k)*2,RfitM(k)*2],...
          'curvature',[1,1],'linestyle','-','edgecolor','r');
      polar(MLTcentM(k),LATcentM(k),'r.');
     if ifwei==1
      medcp = median(cp);
      medcrad = median(thetacrad);
      meddskofc = median(xdskofc);
      medoffc = median(yoffc);
  % Convert back to polar
  [MLTcentofc,LATcentofc]=cart2pol(meddskofc,medoffc);
  % Plot circle
      rectangle('position',[meddskofc-medcrad,medoffc-medcrad,medcrad*2,medcrad*2],...
          'curvature',[1,1],'linestyle','-','edgecolor','g');
      polar(MLTcentofc,LATcentofc,'g.');
     end
   else
      rectangle('position',[xfitM(k)-RfitM(k),yfitM(k)-RfitM(k),RfitM(k)*2,RfitM(k)*2],...
          'curvature',[1,1],'linestyle','-','edgecolor','b');
      polar(MLTcentM(k),LATcentM(k),'b.');
   end
rfit=RfitM(k);
mltcent = MLTcentM(k);
mlatcent = LATcentM(k);
xcent = xfitM(k);
ycent = yfitM(k);
bymed = median(bymag);
   fit_title = sprintf('%5.2fR %5.2fMLT0 %5.2fMLT6 deg',rfit,-ycent,xcent);
%   fit_title = sprintf('%5.2fR %5.2fMLT0 %5.2fMLT6 deg',rfit,mlatcent,mltcent);
   text (xb,yb,fit_title,'FontName','Times','FontSize',12,'color','k');
   if (ncase==1)
     by_title = sprintf('By<0N By>0S median |By| %5.2fnT',bymed); 
     if ifwei==1
	by_title = sprintf('%s median CP %5.1fkV',by_title,medcp);
     end
     text(xr,yr,by_title,'FontName','Times','FontSize',16,'color','k', 'HorizontalAlignment','Center');
   end
   if ncase==3
     by_title = sprintf('By>0N By<0S median |By| %5.2fnT',bymed); 
     if ifwei==1
	by_title = sprintf('%s median CP %5.1fkV',by_title,medcp);
     end
     text(xr,yr,by_title,'FontName','Times','FontSize',16,'color','k', 'HorizontalAlignment','Center');
   end
	axis equal;
      axis off;
    if (ncase==4)
      hold off;
    end
end	% ncase=1:2 for CRBs and sunward max Vys
end	% nbzbin=1:5
end	% ifplotb==1

% 8/11  Make polar plot of values
if ifplotp==1
for iby=1:2
if (iby==1)
 spec_title = sprintf('%s By<0N By>0S',date_title);
else
 spec_title = sprintf('%s By>0N By<0S',date_title);
end

hf=figure;
 set(hf,'position',[162    41   700   906]); 
subplot('position',[.1,.5,.45,.45]);
title_label = sprintf('%s CRBs (blue,red Bz>,<%f5.2)',spec_title,bmid);
% title_label = sprintf('1996 CRBs (blueN,cyanS Bz>-4, redN,magentaS Bz<-4)');
%title_label = sprintf('1998 Jul-Dec CRBs (blueN,cyanS Bz>-4, redN,magentaS Bz<-4)');
%     XM=1.35;
%      h = polarcontf(trho(:,k),ttheta(:,k),(Epot(:,:,k,is)'),-90:10:70); caxis([-90 70]);
%      axis equal; 
%      axis off;
% Print out title etc
    hold on;
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
%         lab=pi/4;
%         x=65*cos(lab);
%         y=-65*sin(lab);
%         text(x,y,' kV','FontName','Times','FontSize',14,'color','k');
%      set(gca,'layer','bottom');
%      v={'-90';'-80';'-70';'-60';'-50';'-40';'-30';'-20';'-10';'0  ';'10 ';'20 ';'30 ';'40 ';'50 ';'60 ';'70 '};
%      hc=colorbar;      
%      pos=get(hc,'position'); 
%      pos(1)=1.19*pos(1); pos(3)=0.9*pos(3); pos(4)=0.95*pos(4);
%      set(hc,'position',pos,'ylim',[-90 70],'ytick',([-90:10:70]),'yticklabel',v);
      %%%%%%%%%%%%%%%%%%%  to set background (MLAT circle and MLT text)
      PolarBackgroundk;
% Plot IDM CRBs and max sunward drift locations
ihem = 1;
      for i=1:nidm(ihem)
	ifgo = 0;
	if (byidm(i,ihem)<0 && iby==1)
	 ifgo = 1;
	end
	if (byidm(i,ihem)>0 && iby==2)
	 ifgo = 1;
	end
	if (ifgo==1)
	if (bzidm(i,ihem) < bmid)
	 polar(mltidm(i,ihem),mlatidm(i,ihem),'r.');
        else
	 polar(mltidm(i,ihem),mlatidm(i,ihem),'b.');
        end
        end
      end
ihem = 2;	% Remember to flip sign of By
      for i=1:nidm(ihem)
	ifgo = 0;
	if (byidm(i,ihem)<0 && iby==2)
	 ifgo = 1;
	end
	if (byidm(i,ihem)>0 && iby==1)
	 ifgo = 1;
	end
	if (ifgo==1)
	if (bzidm(i,ihem) < bmid)
	 polar(mltidm(i,ihem),mlatidm(i,ihem),'m.');
        else
	 polar(mltidm(i,ihem),mlatidm(i,ihem),'c.');
        end
        end
      end
	axis equal;
% Convert back to polar, where Latcent is like radius >0, MLTcent is neg
%  [x,y]=pol2cart(mltcrb1(k),mlatcrb1(k));
%     plot (x,y,'ok','LineWidth',3,'MarkerSize',14);
%  [x,y]=pol2cart(mltcrb2(k),mlatcrb2(k));
%     plot (x,y,'ok','LineWidth',3,'MarkerSize',14);

      axis off;
subplot('position',[.1,.05,.45,.45]);
%title_label = sprintf('1996 max sunward Vys');
%title_label = sprintf('1998 Jul-Dec max sunward Vys');
title_label = sprintf('%s max sunward drifts',spec_title);
% Print out title etc
    hold on;
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
      PolarBackgroundk;
% Plot IDM CRBs and max sunward drift locations
ihem = 1;
      for i=1:nidm(ihem)
	ifgo = 0;
	if (byidm(i,ihem)<0 && iby==1)
	 ifgo = 1;
	end
	if (byidm(i,ihem)>0 && iby==2)
	 ifgo = 1;
	end
	if (ifgo==1)
	if (bzidm(i,ihem) < bmid)
	 polar(mltidmx(i,ihem),mlatidmx(i,ihem),'rx');
        else
	 polar(mltidmx(i,ihem),mlatidmx(i,ihem),'bx');
        end
        end
      end
ihem = 2;
      for i=1:nidm(ihem)
	ifgo = 0;
	if (byidm(i,ihem)<0 && iby==2)
	 ifgo = 1;
	end
	if (byidm(i,ihem)>0 && iby==1)
	 ifgo = 1;
	end
	if (ifgo==1)
	if (bzidm(i,ihem) < bmid)
	 polar(mltidmx(i,ihem),mlatidmx(i,ihem),'mx');
        else
	 polar(mltidmx(i,ihem),mlatidmx(i,ihem),'cx');
        end
        end
      end
	axis equal;
      axis off;
hold off;
end	% for iby=1:2

% 8/11  Make polar plot of lines from CRB to maxVy
hf=figure;
 set(hf,'position',[162    41   700   906]); 
subplot('position',[.1,.5,.45,.45]);
 spec_title = sprintf('%s By<0N By>0S',date_title);
title_label = sprintf('%s CRB-sunwardmax (blue,red Bz>,<%f5.2)',spec_title,bmid);
% Print out title etc
    hold on;
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
      PolarBackgroundk;
% Plot IDM CRBs and max sunward drift locations
ihem = 1;
      for i=1:nidm(ihem)
       if (byidm(i,ihem) < 0.)
	 theta = [mltidm(i,ihem) mltidmx(i,ihem)];
	 rho = [mlatidm(i,ihem) mlatidmx(i,ihem)];
	if (bzidm(i,ihem) < bmid)
	 polar(theta,rho,'r');
        else
	 polar(theta,rho,'b');
        end
       end
      end
ihem = 2;	% Remember to flip sign of By
      for i=1:nidm(ihem)
       if (byidm(i,ihem) > 0.)
	 theta = [mltidm(i,ihem) mltidmx(i,ihem)];
	 rho = [mlatidm(i,ihem) mlatidmx(i,ihem)];
	if (bzidm(i,ihem) < bmid)
	 polar(theta,rho,'m');
        else
	 polar(theta,rho,'c');
        end
        end
      end
	axis equal;
      axis off;
subplot('position',[.1,.05,.45,.45]);
 spec_title = sprintf('%s By>0N By<0S',date_title);
title_label = sprintf('%s CRB-sunwardmax (blue,red Bz>,<%f5.2)',spec_title,bmid);
% Print out title etc
    hold on;
     tit=3*pi/4;
     x=65*cos(tit);
     y=65*sin(tit);
     text(x,y,title_label,'FontName','Times','FontSize',16,'color','k');
      PolarBackgroundk;
% Plot IDM CRBs and max sunward drift locations
ihem = 1;
      for i=1:nidm(ihem)
       if (byidm(i,ihem) > 0.)
	 theta = [mltidm(i,ihem) mltidmx(i,ihem)];
	 rho = [mlatidm(i,ihem) mlatidmx(i,ihem)];
	if (bzidm(i,ihem) < bmid)
	 polar(theta,rho,'r');
        else
	 polar(theta,rho,'b');
        end
       end
      end
ihem = 2;	% Remember to flip sign of By
      for i=1:nidm(ihem)
       if (byidm(i,ihem) < 0.)
	 theta = [mltidm(i,ihem) mltidmx(i,ihem)];
	 rho = [mlatidm(i,ihem) mlatidmx(i,ihem)];
	if (bzidm(i,ihem) < bmid)
	 polar(theta,rho,'m');
        else
	 polar(theta,rho,'c');
        end
        end
      end
	axis equal;
      axis off;
   hold off;
end	% ifplotp==1
