function [passstart1] = newCRB( datevec,sat, passstart1 )
%Main Function
%Paths


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/tengdin/reu11/';
options.temp_path='/home/tengdin/reu11/';
options.idm_source_path='/hao/aim2/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
MLTPlot = 0; %=1 if you want the plots showing orbit

%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/tengdin/reu11';
script_path='/home/tengdin/reu11/';
temp_path='/home/tengdin/reu11/';
idm_source_path=('/hao/aim2/dmsp/dmsp_idm_4s/');
hemisphere = 'N';


yrstr= num2str(datevec(1)); %file name ingredients
yearstring = yrstr(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,'/');
idm_source_path = strcat(idm_source_path,idm_folder);

day = datevec(3);
if day < 10;
    dayS = strcat('0',num2str(day));
else
    dayS = num2str(day);
end

year = datevec(1);
yearS = int2str(year);
month = datevec(2);
monthS = sprintf('%2.2d',month);
year2S = yearS(3:4);
daynum = dayno(datevec);
daynumS = sprintf('%3.3d',daynum);
satS = sprintf('%2.2d',sat);
PlotTitle = ([yearS,' ',monthS,' ',dayS] );

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');
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 = exist(idm_path,'file');

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

disp('Reading IDM file...');

[data1_idm,instr_date ] = read_idm(idm_path,idm_file); %read in file
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
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)),...
    '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));

data1 = [(IDMF.UT)./86400 IDMF.MLAT IDMF.MLT IDMF.VtFl];
cvar_range = [-1 1];
max_height = 15;
varargin = {};
initialized = 0;

if MLTPlot ==1
for k = 1:length(passend)
PlotTitle1 = [PlotTitle,' ','sec:',num2str(data(passstart(k),1) ),'-',num2str(data(passend(k)))];
    if ~initialized
        [out,track,lines] = satelliteTrackPlot(data1(passstart(k):passend(k),:), ...
            cvar_range, max_height, PlotTitle1, varargin);
        initialized = 1;
    else
        [out,track,lines] = satelliteTrackPlot(data1(passstart(k):passend(k),:), ...
            cvar_range, max_height, PlotTitle1, 'Update');
    end
    pause
    close all
end
end

plotVy(IDMF,passstart,passend,PlotTitle,year,daynum,idm_plot,idm_Bound) %plot each plot individually

end


function [data,instr_date ] = read_idm(idm_path,idm_file)
filename = idm_path;
fname = idm_file;
fid = fopen(filename,'r');
instr_date = fgets(fid);
data = dlmread(filename,'',2,0);
fclose(fid);
end

function plotVy(IDMF,passstart,passend,PlotTitle,year,daynum,idm_path,idm_Bound)
%%save Vy values before detrending
data = IDMF.data;
VyValues = data(:,2);

o=omniBox(year); %call the omni data object in
o.setfile
o.readyear

for k = 1:length(passend)
    clf
    y_prime = debase(data(passstart(k):passend(k),1),data(passstart(k):passend(k),2),15); %detrend Vy values
    
     
    %plot(data(passstart(k):passend(k),1),data(passstart(k):passend(k),10),'r.-') %%plot MLAT
    
    MLT = data(passstart(k):passend(k),10);
    
%     for i = 1:length(MLT)
%         if MLT(i) < 15 && MLT(i) > 9 %don't plot if the pass doesn't cross the polar regions
%             disp('Satelite Pass not clean')
%         else
            
         
            BoundArr = [];
            BzArray = [];
            ByArray = [];
            j=1;
            for i = passstart(k):70:passend(k)
                [t, BzArray(j)] = o.inquire(daynum,data((i),1),'Bz');
                [t, ByArray(j)] = o.inquire(daynum,data((i),1),'By');
                j=j+1;
            end
            
            Bz = mean(BzArray);
            By = mean(ByArray);
            clear ('t', 'BzArray', 'ByArray','j')
            
            if Bz < 0 && abs(By) < abs(Bz)
                
                yy = smooth(y_prime,10);
                
                subplot(2,1,1)
                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(data(passstart(k):passend(k),1),0,'--')
                xlabel('UT (sec)')
                ylabel('Vy (km/sec)')
                title([PlotTitle,' ','sec:',num2str(data(passstart(k),1) ),'-',num2str(data(passend(k))), '  '...
                    ,'Bz:',num2str(Bz),'   ','By:',num2str(By)]);
                
                hold off
                
                subplot(2,1,2)
                plot(data(passstart(k):passend(k),1),data(passstart(k):passend(k),10),'r.-') %%plot MLT
                xlabel('UT (sec)')
                ylabel('MLT (hrs)')
                pause
                
                
                passMLT = data(passstart(k):passend(k),10);
                passGLAT = data(passstart(k):passend(k),12);
                passUT = data(passstart(k):passend(k),1);
                for i = 1:length(yy) %get first boundary
                    if i>20
                        if yy(i) < 0 && yy(i+3) < 0
                            BoundArr(1,1) = IDMF.year;
                            BoundArr(1,2) = IDMF.doy;
                            BoundArr(1,3) = passUT(i);
                            BoundArr(1,4) = yy(i);
                            BoundArr(1,5) = passMLT(i);
                            BoundArr(1,6) = passGLAT(i);
                            BoundArr(1,7) = Bz;
                            BoundArr(1,8) = By;
                            break;
                        end
                    end
                end
                backData = sortrows(data(passstart(k):passend(k),:),-1);
                yy = sortrows(yy,-1);
                backData(:,2) = yy;
                
                for i = 1:length(backData(:,1)) %get second boundary starting from the other side
                    if 1>20
                        if backData(i,2) < 0 && backData(i+2,2) < 0
                            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,12);
                            BoundArr(2,7) = Bz;
                            BoundArr(2,8) = By;
                            break;
                        end
                    end
                end
                
                
                fplot=strcat(idm_path, num2str(data(passstart(k),1),'%4.0f'), '_', num2str(data(passend(k),1),'%4.0f') );
                 
                dlmwrite(strcat(idm_Bound,'Boundaries','.txt'),BoundArr,'delimiter',' ','-append'); %write boundaries file
                
                set(gcf,'PaperPositionMode','auto');
                print(gcf,'-dtiff',[fplot '_plot']);
                clear fplot
               
            else
                disp('Bz positive or By large, moving to next pass')
            end
%         end
%     end
end


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 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