%function [passstart1] = newCRB( o,datevec,sat, passstart1 )
 function [passstart1] = newCRB( o,year,daynum,sat, passstart1,hemisphere,idm_hem_year_bndy)
% datevec = yy mm dd hh mm ss where yymm is used for IDM directory, 
%   and daynum is calculated via the function here called dayno which uses datenum
%  yy mm dd then in plot title, which could be year dayno instead
%  month is needed in idm_source_path and idm_Bound name
%Main Function
%Paths

% SH 98199 problems (18 July)

% PlotBoundaries showed 96017 SH (18274-18958) had bad sunward convection (peak .05km/s)
%  Therefore, deleted this from lS139601Boundaries_Bz.txt

global options
options.omni_source_path='/hao/polar1/emery/aurora_i/xuw/combined_output/eln_tot_72_10_combnx7_hrly';
options.omni_daily_path='/home/liamk/data/omni/';
options.log_path='/home/liamk/data/log/';
options.script_path='/home/holden/reu12/';
options.temp_path='/home/holden/reu12/';
options.idm_source_path='/hao/aim3/dmsp/dmsp_idm_4s/';
%options.hemisphere = 'N';
fileMod = 0; %in case we have cleaned up the files and put a new in front of the name
%  Cleaning up is deleting the header line of initial UT, Kp, etc (can have *****)
%    Code works with and without the 1st of 811 lines (so are #orbits/day less lines)
year = 
if (year==2005)
 fileMod = 1;
end

Boundaries = 1; % =1 if you want to make the boundaries file - set =0 if want all orbits when MLTPlot=1
MLTPlot = 0; % =1 if you want the plots showing orbit  
% makePlots should be 0 if MLTPlot=1 since MLTPlot puts all orbits on single polar plot 
%  each day (most good if Boundaries=1, or all available if Boundaries=0)
makePlots = 1; % =1 if you want plots made
if (MLTPlot==1 & makePlots==1)
MLTPlot,makePlots
makePlots=0
end
badstart = 1;
badend = 1;

%paths
omni_source_path='/hao/polar1/emery/aurora_i/xuw/combined_output/eln_tot_72_10_combnx7_hrly';
omni_daily_path='/home/liamk/data/omni/';
log_path='/home/holden/reu12/';
script_path='/home/holden/reu12/';
temp_path='/home/holden/reu12/';
idm_source_path=('/hao/aim3/dmsp/dmsp_idm_4s/');

satS = sprintf('%2.2d',sat);

%year = datevec(1);
%month = datevec(2);
%daynum = dayno(datevec);

% Find month from daynum and year
month = findmonth(year,daynum);	

% Instead of instr_date from read_idm, use daynum for IDMF.doy and year for IDMF.year
%  NOTE: IDMF.year is the last 2 digits of a year, while year is 4 digits

yearS = int2str(year);
year2S = yearS(3:4);
monthS = sprintf('%2.2d',month);	% Need monthS for idm_folder and idm_Bound 
daynumS = sprintf('%3.3d',daynum);

%yearstring = yearS(3:4);
%if datevec(2) < 10;
%    idm_folderMo = strcat('0',num2str(datevec(2))));
%else
%    idm_folderMo = strcat(num2str(datevec(2)));
%end
%idm_folder = strcat(yearstring,idm_folderMo,'/');
if (year==2005)
 idm_folder = strcat(yearS,'/');
else
 idm_folder = strcat(year2S,monthS,'/');
end
idm_source_path = strcat(idm_source_path,idm_folder);

% Do not use day (of month) or dayS
%day = datevec(3);
%if day < 10;
%    dayS = strcat('0',num2str(day));
%else
%    dayS = num2str(day);
%end
%PlotTitle = ([yearS,' ',monthS,' ',dayS] );

PlotTitle = ([yearS,' ',daynumS,' ',hemisphere,'H'] );

% Temporary read of Weimer Vy ONLY if 96078, 05105, 05289
ifwei = 0;
%if (year==1996 && daynum==78)
if (year==2005 && daynum==105 && hemisphere=='N')
%if (year==2005 && daynum==289 && hemisphere=='S')
ifwei = 1;
% Read Weimer Vy information for a day of good IDM orbits
if (hemisphere=='N')
% NH plot UTs #1-end except: 37600-39200(#1-347-333)
%filname = '/home/emery/tieconvec/Nw05_Ealong_B_Vicross_1996_078.txt';
npwei = 6;
nskpwei=0;
filname = '/home/emery/tieconvec/Nw05_Ealong_B_Vicross_2005_105.txt';
npwei = 4;
nskpwei=0;
else
% SH plot UTs: 22800-24200 (#32-end), 29000-30400 (#33-end), 41200-42600 (#37-end)
%filname = '/home/emery/tieconvec/Sw05_Ealong_B_Vicross_1996_078.txt';
nskpwei=36;
npwei = 3;
filname = '/home/emery/tieconvec/Sw05_Ealong_B_Vicross_2005_289.txt';
npwei = 8;
nskpwei=0;
end
WEI=fopen(filname,'r');
for nn=1:npwei
 topwei = fgetl(WEI);
 weitop = sscanf(topwei,'%d %d %d %f %f %f %f');
 nwei(nn) = weitop(1);
 wyear = weitop(2);
 wday = weitop(3);
 wbz(nn) = weitop(4);
 wby(nn) = weitop(5);
 for j=1:nwei(nn)
  linwei = fgetl(WEI);
  if (j>nskpwei)
   i = j-nskpwei;
  weidat = sscanf(linwei,'%f');
  weiuts(i,nn) = weidat(2);
  weimlat(i,nn) = weidat(3);
  weimlt(i,nn) = weidat(4);
  weivy(i,nn) = weidat(11);
  end
 end
end
if (hemisphere=='N')
 nwei(2) = nwei(2)-14;
end
nw = 1;
utmidwei = 0.5*(weiuts(1,nw)+weiuts(nwei(nw)-nskpwei,nw));
end	% of read Weimer 2005 Vy for 96078

if fileMod == 1
  idm_file = strcat('new','l',satS,year2S,daynumS,'a.txt');
  idm_name = strcat('l',satS,year2S,daynumS);
else
  idm_file = strcat('l',satS,year2S,daynumS,'a.txt');
  idm_name = strcat('l',satS,year2S,daynumS);
end

%ssj_bin_file = strcat('j',satS,year2S,'/j4f',satS,year2S,daynumS);
%ssj_txt_file = strcat(satS,'_',yearS,monthS,dayS,'.txt');	% bad idea to use month/day
idm_path = strcat(idm_source_path,idm_file )
idm_plot = strcat(idm_source_path,'l',satS,year2S,daynumS,'a');
%idm_Bound = strcat(idm_source_path,'l',satS,year2S,monthS);
idm_Bound = strcat(idm_source_path,'l',hemisphere,satS,year2S,monthS,'Boundaries_Bz.txt');

idm = exist(idm_path,'file');

if (idm == 0)
    disp('no idm data available! Exiting...');
    return
end

%disp('Reading IDM file...');
% NOTE:  If reading new files w/o the first line for each block of 810 data,
%   CANNOT use instr_date since it does not exist!  
%  Therefore, change read_idm and/or elim use of istr_date!
%[data1_idm,instr_date ] = read_idm(idm_path,idm_file); %read in file
 [data1_idm] = read_idm(idm_path,idm_file); %read in file
%  for the first line in each block of 810, columns 2-9 are:
%  00999.00999.00999.00999.00999.00999.00999.00999.00, while the first col is date
%  so k=find(..) sets the header to missing NaN to skip with isnan (no problem if header missing in the 'new' files)
k = find(data1_idm == 999.0000); %label bad data
data1_idm(k) = NaN;
k = find(~isnan(data1_idm(:,2))); %%get rid of empty values
data1_idm = data1_idm(k,:);

[passstart,passend,data,passstart1] = parsepasses(data1_idm,hemisphere,passstart1); %create indices for passes

%%define structure with handles for each entity
%  Delete use of instr_date
%IDMF = struct('data',data,'instrument',str2double(instr_date(2:3)),'year',str2double(instr_date(4:5)),...
%    'doy',str2double(instr_date(6:8)),'minute', str2double(instr_date(9:10)),'sec', str2double(instr_date(11:12)),...

%%%%%%%%%%%%%%%%%%%%%%%%%%Define Params in IDM Read (HzFl,VtFl in km/s)%%%%%%%%%%%%%%%%%%%%%%%%%%
IDMF = struct('data',data,...
    'UT',data(:,1),'HzFl',data(:,2),'VtFl',data(:,3),'StdDevHz',data(:,4),'StDevVt',data(:,5),'ElctroPot',data(:,7),...
    'MLAT',data(:,9),'MLT',data(:,10),'GLAT',data(:,12),'GLON',data(:,13));

%%%if MLTPlot ==1	% Moved plotting to 2 other plot areas

if MLTPlot ==1
% VtFl = vertical flow but want HzFl and preferably smoothed 20 sec horizontal flow!
%data1 = [(IDMF.UT)./86400 IDMF.MLAT IDMF.MLT IDMF.VtFl];
%cvar_range = [-1 1];
cvar_range = [-3 3];
max_height = 15;
varargin = {};
initialize = 1;
%  8/11 Barbara Emery:  Set normalization so is normalized to the magnitude and not to the colorbar if cvar_range(1)<0. and cvar_range(2)>0.
  if cvar_range(1)<0 && cvar_range(2)>0
    cbarmax = max(abs(cvar_range(1)),cvar_range(2));
  else
    cbarmax = 0.
  end
end



%save Vy values before detrending
data = IDMF.data;
VyValues = data(:,2);

% Read OMNI file of 5-min data once only for a year or month of data, not all the time!
%  But, if this is not done always, then o is undefined - unless it is passed in!
%o=omniBox(year); %call the omni data object in
%o.setfile
%o.readyear

% Loop for day IDM file
% 1996 smin passes for NH winter? especially cut off before completing the 
%  second cut through the aurora because of light ion contamination(?)
%  Therefore, can try to find the second CRB and equatorward peak if ends between 60-70 mlat
%    BUT DO NOT DETREND - Skip entirely if cut off >|70| mlat
detrendmlat = 60.;
badmlat = 70.;
% Do not keep Vy sunward maxima if Vyx<50m/s (0.05km/s)
vymin = 0.05;

dark_green = [0 0.5 0];
fsize = 12;

for k = 1:length(passend)
% nend is the same as length of passUT, data(passtart(k):passend(k),1)), etc
nend = passend(k) - passstart(k) + 1; 
  
                        passUT = data(passstart(k):passend(k),1);
                        passMLAT = data(passstart(k):passend(k),9);
                        passMLT = data(passstart(k):passend(k),10);
                        passGLAT = data(passstart(k):passend(k),12);
                        passGLON = data(passstart(k):passend(k),13); 

%  Find 20 sec smoothed Vy values and detrend here before any checking in order to plot all available orbits on 1 day if MLTPLOT=1 and 
ifdetrend = 1;
% 1996 smin NH pass especially cuts off before completing the 
%  second cut through the aurora because of light ion contamination(?)
%  Therefore, can try to find the second CRB and equatorward peak if ends between 60-70 mlat
%    BUT DO NOT DETREND - Skip entirely if cut off >|70| mlat 
     if abs(passMLAT(1)) > detrendmlat | abs(passMLAT(nend)) > detrendmlat  
	 ifdetrend = 0;
     end
% 1998 SH passes often large negative at -45mlat ~6MLT, but are ~0 ~-55mlat and ~0 at -45mlat ~18 MLT
%   DO NOT DETREND SH passes under these circumstances 
vy1 = VyValues(passstart(k));
vy2 = VyValues(passend(k)); 
     if passMLAT(1)<0. && vy1<-0.25 && abs(vy2)<0.20
	 ifdetrend = 0;
     end
        if (ifdetrend==1)          
                    y_prime = debase(data(passstart(k):passend(k),1),data(passstart(k):passend(k),2),15); %detrend Vy values
	else
	 	y_prime = VyValues(passstart(k):passend(k));
	end
% smooth(y) does a moving smoothing over default of 5 values (20 sec for 4 sec data); 
% smooth(y,span) where span MUST BE ODD - hence smooth(y,10) is converted to the default 5
%                    yy = smooth(y_prime,10);
                    yy = smooth(y_prime,5);

if MLTPlot==1 && Boundaries==0
data1 = [passUT/3600. passMLAT passMLT yy];  % will label start of pass in hours=passUT/3600.
    if initialize==1
   [out,track,lines] = satelliteTrackPlot(data1,cvar_range, max_height, PlotTitle, initialize,cbarmax,{});
        initialize = 0;
    else
% 'Update' does not seem to work, so passed in initialize and cbarmax
   [out,track,lines] = satelliteTrackPlot(data1,cvar_range, max_height, PlotTitle, initialize,cbarmax,'Update');
    end
     hold on;
end				% if MLTPlot==1 to plot all available orbits

% Do not find boundaries if so many gaps have < 200 points in the polar pass (200*4s=800s or 47.04 arc length)
%  Check to see if have at least 150 4-s points in the pass (20-s 1.176deg*30=35.28 deg if no gaps)                
%                if length(data(passstart(k):passend(k),1)) > 20 
               if nend > 150 
   
% Check to see that begin or end of pass is >|70|mlat
      if abs(passMLAT(1)) < badmlat && abs(passMLAT(nend)) < badmlat
                 
            BoundArr = [];
            VswArray = [];
            DswArray = [];
            BzArray = [];
            ByArray = [];
            j=1; %counter for Bz
% Do not set m since this is used to see if the first boundary has been found or not
%            m=1; %counter for checking area under boundaries
            mark = 0;
            for i = passstart(k):70:passend(k)
                if  data((i),1) > 600
                    [t, BzArray(j)] = o.inquire(daynum,(data((i),1) - 600),'Bz'); %find Bz values for the time that the spacecraft
                    [t, ByArray(j)] = o.inquire(daynum,(data((i),1) - 600),'By'); %made orbit
                    [t, VswArray(j)] = o.inquire(daynum,(data((i),1) - 600),'vsw');
                    [t, DswArray(j)] = o.inquire(daynum,(data((i),1) - 600),'dsw');
                    j=j+1;
%  Assume although j increases, can have missing values of Vsw or Dsw especially, but mean correct
                end
            end
            
            Bz = mean(BzArray);
            By = mean(ByArray);
	    Vsw = mean(VswArray);
	    Dsw = mean(DswArray);
            clear ('t', 'BzArray', 'ByArray','j', 'VswArray', 'DswArray')
            
            if Bz < 0 && (1*(abs(By))) < abs(Bz)
%            if Bz < 0 && (3*(abs(By))) < abs(Bz)


	        
% Find boundaries if correct IMF conditions (and >600s or 150 pts) but only save to a file if Boundaries==1                   
%    if Boundaries ==1 
badstart = 0;
badend = 0;
% Give MLT range of hours ~6 dawn and ~18 dusk for NH and SH for beginning  and end of pass ~+/-45 mlat
if(strcmp(hemisphere,'N'))
 begmlt1 = 15;
 begmlt2 = 22;
 endmlt1 = 2;
 endmlt2 = 9;
elseif(strcmp(hemisphere,'S'))
 begmlt1 = 2;
 begmlt2 = 9;
 endmlt1 = 15;
 endmlt2 = 22;
end
      if data(passstart(k),10) < begmlt1 | data(passstart(k),10) > begmlt2 
	badstart=1;
      end
      if data(passend(k),10) < endmlt1 | data(passend(k),10) > endmlt2 
	badend=1; 
      end
% Print to screen if pass is a possibility - skip print since UT is in title
%check_bnds_utbegin = passUT(1) 
%end_mlat = passMLAT(nend)

% Find the index of the maximum mlat in order to find the maximum sunward velocities on either side of this index
xmlat = 0.5*(badmlat+detrendmlat);
ixmlat = 0;
for i=1:nend
  if (abs(passMLAT(i)) > xmlat)
    xmlat = abs(passMLAT(i));
    ixmlat = i;
  end
end

% Start of loop for boundaries
if (ixmlat>0)
% NEED to find the max Vy here and then slide polewards to find FIRST neg
        yy1x = -100.;
	for i=1:ixmlat
	  if yy(i) > yy1x & abs(passMLAT(i)) > 50.
	     yy1x = yy(i);
	     i1x = i;
          else
	  end
	end
%yy1x
%ut1x = passUT(i1x)
%mlat1x = passMLAT(i1x)
% Now find first negative polewards
	i = i1x+1;
        while yy(i) > 0 && i < ixmlat
	  i = i + 1;
        end
%      if (length(yy)-i)>20
%yy(i-1),yy(i),yy(i+1)
%  yy(i-1) is positive and yy(i) is the first negative	
        i10 = i; 
	if abs(yy(i-1)) < abs(yy(i)) 
		i10 = i-1; 
        end     
% Toss pass if yy1x<vymin
        if yy1x<vymin 
	 badstart=1;
        end
% Toss pass if i1x close to a gap
	if (i1x > 3)
	if passUT(i1x)-passUT(i1x-3)>20. | passUT(i1x+3)-passUT(i1x)>20. 
	 badstart=1;
        end
	else
	 badstart=1;
        end
% Toss pass if i10 close to a gap
	if (i10 > 3)
	if passUT(i10)-passUT(i10-3)>20. | passUT(i10+3)-passUT(i10)>20. 
	 badstart=1;
        end
	else
	 badstart=1;
        end
%                        chkar = [];
%                        for i = 1:(length(yy)-35) %get first boundary and write to Boundary array 
%%subtract 20 so that algorithm won't choke
%                            if i>20
%                                if yy(i) < 0  %this loop checks to make sure that all the next 30 points are also negative
%                                    for p = i:(i+35)
%                                        chkar(p) = yy(p);
%                                        if yy(p)>0
%                                            mark = 1;
%                                        else
%                                        end
%                                    end
%                                    chknum = min(chkar);
%                                    if mark == 0 && chknum < -.10
ib = 1;
%                                       BoundArr(ib,1) = IDMF.year;
%                                       BoundArr(ib,2) = IDMF.doy;
                                        BoundArr(ib,1) = year;
                                        BoundArr(ib,2) = daynum;
                                        BoundArr(ib,3) = Bz;
                                        BoundArr(ib,4) = By;
                                        BoundArr(ib,5) = Vsw;
                                        BoundArr(ib,6) = Dsw;
i = i1x;
                                        BoundArr(ib,7) = passUT(i);
                                        BoundArr(ib,8) = yy1x;
                                        BoundArr(ib,9) = passMLAT(i);
                                        BoundArr(ib,10) = passMLT(i);
                                        BoundArr(ib,11) = passGLAT(i);
                                        BoundArr(ib,12) = passGLON(i);
i = i10;
                                        BoundArr(ib,13) = passUT(i);
                                        BoundArr(ib,14) = yy(i);
                                        BoundArr(ib,15) = passMLAT(i);
                                        BoundArr(ib,16) = passMLT(i);
                                        BoundArr(ib,17) = passGLAT(i);
                                        BoundArr(ib,18) = passGLON(i);

%                                       s = i;
                                        m = passstart(k) + i;
%                                       break;
%                                    end
%                                end
%                            end
%                            mark = 0;
%                        end		% for i=1:(length(yy)-35)
%	end		% for (length(yy)-i)>20
                     
                        backData = sortrows(data(passstart(k):passend(k),:),-1);
                        y = flipud(yy); %flip the data
                        backData(:,2) = y; %use detrended and averaged data

% NEED to find the max Vy here in ~80% of length/2 and then slide polewards to find FIRST neg
        yy1x = -100.;
	for i=1:nend-ixmlat
	  if y(i) > yy1x & abs(backData(i)) > 50.
	     yy1x = y(i);
	     i1x = i;
          else
	  end
	end
%yy1x
%ut1x = backData(i1x,1)
%mlat1x = backData(i1x,9)
% Now find first negative polewards
	i = i1x+1;
% for 96009 NH, never found negative after positive peak on dayside
        while y(i) > 0 && i < nend-ixmlat
	  i = i + 1;
        end
%      if (length(y)-i)>20
%y(i-1),y(i),y(i+1),y(i+2)
%  yy(i-1) is positive and yy(i) is the first negative	
        i10 = i; 
	if abs(y(i-1)) < abs(y(i)) 
		i10 = i-1; 
	end                             
% Toss pass if yy1x<vymin
        if yy1x<vymin 
	 badend=1;
        end
% Toss pass if i1x close to a gap
	if (i1x > 3)
	if backData(i1x-3,1)-backData(i1x,1)>20. | backData(i1x,1)-backData(i1x+3,1)>20. 
	 badend=1;
	end
	else
	 badend=1;
        end
% Toss pass if i10 close to a gap
	if (i10 > 3)
	if backData(i10-3,1)-backData(i10,1)>20. | backData(i10,1)-backData(i10-3,1)>20. 
	 badend=1;
	end
	else
	 badend=1;
        end
%                        for i = 1:(length(backData(:,1))-35) %get second boundary starting from the other side
%                            if i>20 %skip the first 20 points
%                                if backData(i,2) < 0
%                                    for p = i:(i+35)
%                                        chkar(p) = backData(p,2);
%                                        if backData(p,2) > 0
%                                            mark = 1; %don't use this crossing
%                                        else
%                                        end
%                                    end
%                                    chknum = min(chkar);
%                                    if mark == 0 && chknum <-.10
%                                        BoundArr(2,1) = IDMF.year;
%                                        BoundArr(2,2) = IDMF.doy;
%                                        BoundArr(2,3) = backData(i,1);
%                                        BoundArr(2,4) = backData(i,2);
%                                        BoundArr(2,5) = backData(i,10);
%                                        BoundArr(2,6) = backData(i,9);
%                                        BoundArr(2,7) = backData(i,12);
%                                        BoundArr(2,8) = backData(i,13);
%                                        BoundArr(2,9) = Bz;
%                                        BoundArr(2,10) = By;
%                                        BoundArr(2,11) = Vsw;
%                                        BoundArr(2,12) = Dsw;
ib = 2;
%                                       BoundArr(ib,1) = IDMF.year;
%                                       BoundArr(ib,2) = IDMF.doy;
                                        BoundArr(ib,1) = year;
                                        BoundArr(ib,2) = daynum;
                                        BoundArr(ib,3) = Bz;
                                        BoundArr(ib,4) = By;
                                        BoundArr(ib,5) = Vsw;
                                        BoundArr(ib,6) = Dsw;
i = i1x;
                                        BoundArr(ib,7) = backData(i,1);
                                        BoundArr(ib,8) = yy1x;
                                        BoundArr(ib,9) = backData(i,9);
                                        BoundArr(ib,10) = backData(i,10);
                                        BoundArr(ib,11) = backData(i,12);
                                        BoundArr(ib,12) = backData(i,13);
i = i10;
                                        BoundArr(ib,13) = backData(i,1);
                                        BoundArr(ib,14) = backData(i,2);
                                        BoundArr(ib,15) = backData(i,9);
                                        BoundArr(ib,16) = backData(i,10);
                                        BoundArr(ib,17) = backData(i,12);
                                        BoundArr(ib,18) = backData(i,13);
                                        n=i;
%                                        break;
%                                    end
%                                end
%                            end			% for i=1:length
%                            mark = 0;
%                        end               % for i=1:(length(backData(:,1))-35)
%	end		% for (length(y)-i)>20
 else
   badstart=1;
   badend=1;    
 end		% for ixmlat>0 
          
% Make plots here in order to put CRB on plot and peaks if found

% satelliteTrackPlots have smallest length for min cvar_range and NOT for 0 which is wrong!
%   Color bar is not very useful anyway, since length should correspond to magnitude
if MLTPlot ==1 && badstart==0 && badend==0
data1 = [passUT/3600. passMLAT passMLT yy];  % will label start of pass in hours=passUT/3600.
%PlotTitle1 = [PlotTitle,' ','sec:',num2str(data(passstart(k),1) ),'-',num2str(data(passend(k)))];
%   [out,track,lines] = satelliteTrackPlot(data1,cvar_range, max_height, PlotTitle1, varargin);
    if initialize==1
%        [out,track,lines] = satelliteTrackPlot(data1(passstart(k):passend(k),:), ...
%            cvar_range, max_height, PlotTitle1, varargin);
   [out,track,lines] = satelliteTrackPlot(data1,cvar_range, max_height, PlotTitle, initialize,cbarmax,{});
        initialize = 0;
    else
   [out,track,lines] = satelliteTrackPlot(data1,cvar_range, max_height, PlotTitle, initialize,cbarmax,'Update');
%        [out,track,lines] = satelliteTrackPlot(data1(passstart(k):passend(k),:), ...
%            cvar_range, max_height, PlotTitle1, 'Update');
    end
%    pause
%    close all
     hold on;
end				% if MLTPlot==1
        
    if makePlots ==1 && badstart==0 && badend==0
%    if makePlots ==1
%   title_label = sprintf('%4.1f MLT1 %s sec:%5.0f-%5.0f Bz=%5.1f By=%5.1f %4.1f MLT2',data(passstart(k),10),PlotTitle,data(passstart(k),1),data(passend(k),1),Bz,By,data(passend(k),10))
%  do not need MLT1,2 in plot title if plot MLT also
    title_label = sprintf('%s sec:%5.0f-%5.0f Bz=%5.1f By=%5.1f',PlotTitle,data(passstart(k),1),data(passend(k),1),Bz,By)
% top plot
	H1=figure;
%                       subplot(2,1,1)
	if ifwei==1
		subplot ('Position',[0.1,0.10,0.83,0.80])
	else
		subplot ('Position',[0.1,0.50,0.83,0.40])
	end
                        plot(data(passstart(k):passend(k),1),VyValues(passstart(k):passend(k)),'r.-'); %plot original values (not detrended)
                        hold on
                        plot(data(passstart(k):passend(k),1),yy,'g.-');
                        plot(data(passstart(k):passend(k),1),y_prime,'b.-'); %%plotting UT vs. Vy corrected

% Plot Weimer 2005 Vy for 96078 (is NOT at the same UTs as passstart(k):passend(k)!)
if ifwei==1 && utmidwei>passUT(1) && utmidwei<passUT(nend) 
 clear utwei;
 clear vywei;
  for i=1:nwei(nw)-nskpwei
    utwei(i) = weiuts(i,nw);
    vywei(i) = weivy(i,nw);
  end
	plot(utwei,vywei,'.-','Color',dark_green);
	text (passUT(nend),BoundArr(1,8),'Weimer05','Color',dark_green,'FontSize',fsize,'HorizontalAlignment','Right','fontweight','bold');
        nw = nw + 1;
	if (nw<=npwei)
           utmidwei = 0.5*(weiuts(1,nw)+weiuts(nwei(nw)-nskpwei,nw));
	else
	   utmidwei = -99999.
        end
end
grid;
if badstart==0
	plot(BoundArr(1,7),BoundArr(1,8),'om','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,14),'oc','MarkerSize',10);
   else
	plot(BoundArr(1,7),BoundArr(1,8),'xm','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,14),'xc','MarkerSize',10);
end
if badend==0
	plot(BoundArr(2,7),BoundArr(2,8),'om','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,14),'oc','MarkerSize',10);
   else
	plot(BoundArr(2,7),BoundArr(2,8),'xm','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,14),'xc','MarkerSize',10);
end
                        plot(data(passstart(k):passend(k),1),0,'--')
                        ylabel('Vy (km/sec)','fontsize',fsize,'fontweight','bold')
	                title (title_label,'fontsize',fsize,'fontweight','bold');
	if (ifwei==1)
                       xlabel('UT (sec)','fontsize',15,'fontweight','bold')
			hold off
	else
	set (gca, 'XTickLabel',{' ';});
%                       hold off
% bottom plots MLAT and MLT
%                        subplot(2,1,2)
% Plot MLAT in the middle
		subplot ('Position',[0.1,0.30,0.83,0.20])
	hold on;
                        plot(data(passstart(k):passend(k),1),data(passstart(k):passend(k),9),'r.-') %%plot MLAT
grid;
if badstart==0 
	plot(BoundArr(1,7),BoundArr(1,9),'om','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,15),'oc','MarkerSize',10);
   else
	plot(BoundArr(1,7),BoundArr(1,9),'xm','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,15),'xc','MarkerSize',10);
end
if badend==0
	plot(BoundArr(2,7),BoundArr(2,9),'om','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,15),'oc','MarkerSize',10);
   else
	plot(BoundArr(2,7),BoundArr(2,9),'xm','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,15),'xc','MarkerSize',10);
end

%                       xlabel('UT (sec)','fontsize',fsize,'fontweight','bold')
                        ylabel('MLAT (deg)','fontsize',fsize,'fontweight','bold')
	set (gca, 'XTickLabel',{' ';});
% Plot MLT axis
		subplot ('Position',[0.1,0.10,0.83,0.20])
	hold on;
                        plot(data(passstart(k):passend(k),1),data(passstart(k):passend(k),10),'b.-'); %%plot MLT
grid;
if badstart==0 
	plot(BoundArr(1,7),BoundArr(1,10),'om','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,16),'oc','MarkerSize',10);
   else
	plot(BoundArr(1,7),BoundArr(1,10),'xm','MarkerSize',10);
	plot(BoundArr(1,13),BoundArr(1,16),'xc','MarkerSize',10);
end
if badend==0
	plot(BoundArr(2,7),BoundArr(2,10),'om','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,16),'oc','MarkerSize',10);
   else
	plot(BoundArr(2,7),BoundArr(2,10),'xm','MarkerSize',10);
	plot(BoundArr(2,13),BoundArr(2,16),'xc','MarkerSize',10);
end
                       xlabel('UT (sec)','fontsize',fsize,'fontweight','bold')
                       ylabel('MLT (h)','fontsize',fsize,'fontweight','bold')

%                        pause                         
%                        fplot=strcat(idm_plot, num2str(data(passstart(k),1),'%4.0f'), '_', num2str(data(passend(k),1),'%4.0f') );
%                        set(gcf,'PaperPositionMode','auto');
%                        print(gcf,'-dtiff',[fplot '_plot']);
%                        clear fplot

		hold off;
    end		% if ifwei==1 (end of else part)
    end		% if makePlots==1

% NOTE:  With suppression of condition that i be 20 away from ends, m and n should always exist unless ixmlat=0 or ilx=0 (ie, highest abs(mlat)<65deg)
%criteria for saving the boundaries
%                            if exist('m','var') == 1 && exist('n','var') == 1 


%      if data(passstart(k),10) > begmlt1 && data(passstart(k),10) < begmlt2 && backData(1,10) > endmlt1 && backData(1,10) < endmlt2 

%  Check to see if have skimming pass (>600 sec  or 38.22 deg arc from start or end to CRB)
%  Changed from 600 sec or 35.28 deg arc 
%                               if data(m,1) - data(passstart(k),1) <650 && backData(1,1) - backData(n,1) <650%don't save a pass that is irregular
                                    %area1 = trapz(data(passstart(k):m-1,1),yy(1:s)); %take area before the first crossing
                                    %area2 = trapz(data(m:(passstart(k) + length(yy) - n),1),yy(s:(length(yy) - n))); %take area between crossings
                                    %ind = area1/area2;
% write boundaries file
    if Boundaries ==1 && badstart==0 && badend==0
%  dlmwrite(idm_Bound,BoundArr,'delimiter',' ','-append');	% Skip writing in yrmo
			dlmwrite(idm_hem_year_bndy,BoundArr,'delimiter',' ','-append'); 
    end
%                                    dlmwrite(strcat(idm_Bound,'Boundaries_Bz','.txt'),BoundArr,'delimiter',' ','-append'); %write boundaries file
%                                else
%	skim_label = sprintf('%d %f skip following skimming(?) pass where time>650s from CRB to start or end',daynum,data(passstart(k),1)); 
%                                   dlmwrite(idm_Bound,BoundArr,'delimiter',' ','-append'); 
%                                   disp('skimming or irregular pass where time>650s from CRB to start or end, not saved')
%                               end
%                            else
%                                disp('begin/end of pass outside of 4-h limits')
%                            end
                             
%                         else
%                             disp('boundary not found, moving to next pass')
%                         end                           

%                    else
%                    end		% if Boundaries==1 

                else
%                    disp('Bz positive or By large, moving to next pass')
                end             % if Bz>0 or By large
                         else
%                             disp('begin/end of pass poleward of |70|mlat limits')
                         end
                else
%                    disp('pass is not correct size (<150 pts or 600s), skipping')
                end		% if length > 20 (changed to 150)
end		% for k=1:length(passend)

end		% function


%function [data,instr_date ] = read_idm(idm_path,idm_file)
% Do not get istr_date from first line anymore (new... files don't have it!)
function [data] = read_idm(idm_path,idm_file)
filename = idm_path;
fname = idm_file;
fid = fopen(filename,'r');
%instr_date = fgets(fid);
% This starts at the 3rd line of the record (skips first 2 data lines for new.. files)
data = dlmread(filename,'',2,0);
fclose(fid);
end

function [passstart,passend,data,passstart1] = parsepasses(data1_idm,hemisphere,passstart1)

%parse out missing values
%requires 4 column input data, expects time in days in first
%column, mlat in second column

%REQUIRES UNITS OF DAYS FOR TIME COLUMN, TO CONVERT FROM UT
%SECONDS DO XUT./86400;

%NS is string with values 'N' or 'S'

data = data1_idm;
if passstart1 ~=0 %add in leftover values from previous day
    data = passstart1 + data;
end

g =  ~isnan(data(:,1)) & ~isnan(data(:,9)); %check to make sure all values are good
data = data(g,:);


%parse out N polar parts of passes
gn = data(:,9) >= 45;
gs = data(:,9) <= -45;

if(strcmp(hemisphere,'N'))
    hemidata = data(gn,:);
    hemi = gn;
elseif(strcmp(hemisphere,'S'))
    hemidata = data(gs, :);
    hemi = gs;
end
hemidata(:,1) = hemidata(:,1)./86400;
n = 1; passstart(1) = 1;
for k = 1:length(hemidata(:,1))-1
    deltat = hemidata(k+1,1) - hemidata(k,1);
    if deltat >= 1/24;
        passend(n) = k;
        n=n+1;
        passstart(n) = k+1;
    end
end
data = data(hemi,:);

%save leftover data
if length(passstart) > length(passend)
    passstart1 = data(passstart(end):end,:);
end
end

function out = dayno(datevec)
%lmk - Expects date vector format [YYYY MM DD 0 0 0]
dateNum = datenum(datevec)-datenum([datevec(1) 1 1 0 0 0]) + 1;
out = dateNum;
end

function out = findmonth(year,daynum)
ndmon = [31,28,31,30,31,30,31,31,30,31,30,31];
iy4 = int32(year/4);
if (iy4*4==year) 
 ndmon(2) = 29;
else
end
i = 0;
nm = 0;
while i<daynum
  nm = nm + 1;
  i = i + ndmon(nm);
end
out = nm; 
end

function y_prime = debase(x1,y1,N)
%DEBASE - removes a linear baseline from an array based on
%a linear fit with preferential weight given to the first and last N
%points

x2 = [x1(1:N);x1(end-N:end)];
y2 = [y1(1:N);y1(end-N:end)];

coeff = polyfit(x2,y2,1);
m = coeff(1); b = coeff(2);

fun = @(x)(m.*x + b);

line = fun(x1);

y_prime = y1-line;

end
