function get_CRB_peakFlux( datevec, sat )
%Main Function
%Paths

%parameters/constants
MLAT_CUTOFF=45;
NH=1; %calculate for northern hemisphere
SH=1; %ditto southern hemisphere


global options
options.omni_source_path='/hao/polar1/emery/aurora_i/xuw/combined_output/eln_tot_72_10_combnx7_hrly'; 
options.ssj_binary_path='/hao/aim3/dmsp/dmsp_ssj/';
options.omni_daily_path='/home/liamk/data/omni/';
options.ssj_daily_path='/home/tengdin/reu11/';
options.log_path='/home/liamk/data/log/';
options.script_path='/home/tengdin/reu11/';
options.cenergy_path='/home/tengdin/reu11/';
options.temp_path='/home/tengdin/reu11/';
options.idm_source_path='/hao/aim3/dmsp/dmsp_idm_4s/';

%paths
omni_source_path='/hao/polar1/emery/aurora_i/xuw/combined_output/eln_tot_72_10_combnx7_hrly'; 
ssj_binary_path='/hao/aim3/dmsp/dmsp_ssj/';
omni_daily_path='/home/liamk/data/omni/';
ssj_daily_path='/home/tengdin/reu11/';
log_path='/home/tengdin/reu11';
script_path='/home/tengdin/reu11/';
cenergy_path='/home/tengdin/reu11/';
temp_path='/home/tengdin/reu11/';
idm_source_path=('/hao/aim3/dmsp/dmsp_idm_4s/');


yrstr= num2str(datevec(1)); %%code that we may not need
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 = strcat(num2str(day));
end

%file name ingredients
year = datevec(1); 
yearS = int2str(year);
month = datevec(2); 
monthS = sprintf('%2.2d',month);
%dayS = sprintf('%2.Valid options:-r [--revision] a2d',day);
year2S = yearS(3:4);
daynum = dayno(datevec); 
daynumS = sprintf('%3.3d',daynum);
satS = sprintf('%2.2d',sat);

%file names
%idm_file = strcat(yearS,'/',year2S,monthS,'/l',satS,year2S,daynumS,'a.txt')
idm_file =strcat('l',satS,year2S,daynumS,'a.txt');
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 );

Cenergy_file_path = strcat(ssj_daily_path, ssj_txt_file);

%make sure we have data for these days
ssj = exist(strcat(ssj_binary_path, ssj_bin_file),'file');
idm = exist(idm_path,'file');
%if (ssj == 0)
%    disp('no ssj data available! Exiting...');
%    return
%end
if (idm == 0)
    disp('no idm data available! Exiting...');
    return
end    

%see if we have text ssj data so we don't run cenergy
ssj = exist(strcat(ssj_daily_path,ssj_txt_file));
if (ssj == 2)
    disp('text SSJ data exists');
    quickLoad = 1;
else
    disp('text SSJ data does not exist');
    quickLoad = 0;
end

%read the IDM 4s data from text files
%-------------------------------------



disp('Reading IDM file...');
    
[data1_idm,IDMFile] = read_idm(idm_path,idm_file);


% read the DMSP SSJ and Ephemeris data (cenergy.pro output)
%----------------------------------------------------------


if(~quickLoad) %Run cenergy.pro and generate the ASCII file
    disp('Running cenergy.pro...');
    pause(2)
    
    ssj_txt_file = createCenergyText2([datevec(1),datevec(2),datevec(3),0,0,0], sat);

    disp('Reading DMSP ASCII file...');
    pause(2)
    
    [data1_ssj, DayFile] = read_Cenergy(Cenergy_file_path);
else %Assume the ASCII file already exists and read from that
    
    ssj_txt_file = strcat('cenergy_all_f',satS,'_',yearS,monthS,dayS,'.txt');
    disp(Cenergy_file_path);
    [data1_ssj,DayFile] = read_Cenergy(Cenergy_file_path);
end

% make all missing values NaN
%----------------------------

k = find(abs(data1_ssj) < 10^(-30)); 
data1_ssj(k) = NaN;
k = find(data1_idm == 999.0000); 
data1_idm(k) = NaN;

% make all points that correspond to a maglat below cutoff point NaN
%-------------------------------------------------------------------
kkN=[];
MLAT = IDMFile.MLAT; 
MLT = IDMFile.MLT;
if(NH)
kkN = find(MLAT < MLAT_CUTOFF);
end
if(SH)
kkS = find(MLAT > -MLAT_CUTOFF);
end
%MLAT(kkN) = NaN; %delete MLAT values outside cutoff
%MLAT(kkS) = NaN;

% extract the relavent variables from the data files
%---------------------------------------------------

% IDM data
%all quantities are taken for the midpoint of the bin


idm_Vy_N = data1_idm(:,2); 
idm_Vy_N(kkN) = NaN;  %km/s 
%idm_sVy = data1_idm(:,4); 
%idm_sVy(kk1) = NaN; %stdev Vy
%idm_Vz = data1_idm(:,3); %km/s
%idm_sVz = data1_idm(:,5); %stdev Vz
%idm_pot = data1_idm(:,7); %average electrostatic potential in kV
%idm_mlat110 = data1_idm(:,8); %corrected magnetic latitude at 110 km trace
%along field line
%idm_mlatssp = data1_idm(:,9); %corrected mag lat at sub satellite point 
%idm_mlt = data1_idm(:,10); %magentic local time
%idm_invlat = data1_idm(:,11); %invariant lat
%idm_lat = data1_idm(:,12); %geographic lat
%idm_lon = data1_idm(:,13); %geographic lon

EEFlx_N = data1_ssj(:,11); %In_ele_eflx(mW/m2)
EEFlx_N(kkN) = NaN; 
IEFlx_N = data1_ssj(:,14); %Int_ion_eflx
IEFlx_N(kkN) = NaN; 

UTTime = 0:length(EEFlx_N);

clear('data1_ssj','data1_idm');
%correct the Vy baseline, find CRB and flux peaks for each pass
%-------------------------------------

%set pass start, pole and end markers
if ( ~isnan(idm_Vy_N(1)) ) %in case we start in the middle of a pass 
    k = find(isnan(idm_Vy_N),1); %find first point beyond cutoff and start from there
    istart = k;    
else
    istart = 2;
end


Vy = idm_Vy_N;
eF = EEFlx_N;
iF = IEFlx_N;
%clear idm_Vy_N EEFlx_N IEFlx_N

gotstart=0; gotend=0;
j = 1; %pass index
for i = istart:length(Vy)
    this = eF(i);
    last = eF(i-1);
    try
    next = eF(i+1);
    catch
        if ( (isnan(this)) && (~isnan(last)) )
        pe = i; 
        passend( j ) = i; 
        gotend = 1;
        else
            return;
        end
    end

    if ( (~isnan(this)) && (isnan(last)) )
        if (~isnan(next))
        ps = i;
        passstart( j ) = i; 
        gotstart = 1;
        end
    end
    if ( (isnan(this)) && (~isnan(last)) )
        if (isnan(next))
        pe = i; 
        passend( j ) = i; 
        gotend = 1;
        end
    end
    
    if (gotstart && gotend)
     
        %passMLAT = MLAT(ps:pe);
        %[value,pp] = max(abs(MLAT(ps:pe))); %fixed indexing here
        %pp = ps + pp - 1;
        %1st side of oval
        
        passVy = Vy(ps:pe);
        passeF = eF(ps:pe);
        passiF = iF(ps:pe);
        passMLAT = MLAT(ps:pe);
        passtime = UTTime(ps:pe);
        
      
        passVyC = correctBaseline(passVy);
        k = find(~isnan(passVyC));
        VyCt = [transpose(passtime),passVyC];
        VyCt = VyCt(k,:);
        crb = CRB(passVyC);
        eFpeak = findFluxPeak(passeF);
        %iFpeak = findFluxPeak(passiF);
         
        %{
        pass.Vy = passVy; pass.eF = passeF; pass.iF = passiF;
        pass.MLAT = passMLAT; pass.time = passtime; pass.VyC = passVyC;
        pass.crb = crb; pass.eFpeaks = eFpeaks; pass.iFpeaks = iFpeaks;
        

        filename = sprintf('passDawn%2d.mat',j);
        save(filename, 'pass');
        clear pass
        %}
        
        axes();
        
        plot(VyCt(:,1), VyCt(:,2), 'k-'); hold on
        plot(passVyC(crb),passtime(crb), '-ro'); 
        
        pause
        
        axes();
        
        plot(passtime,passeF, 'k-'); hold on
        plot(passtime(eFpeak), passeF(eFpeak) , '-ro');

        
        pause
        
        %2nd side of oval
        
                          
        passVy = Vy(ps:pe);
        passeF = eF(ps:pe);
        passiF = iF(ps:pe);
        passMLAT = MLAT(ps:pe);
        passtime = UTTime(ps:pe);
        
        passVyC = correctBaseline(passVy);
        crb = CRB(passVyC);
        eFpeak = findFluxPeak(passeF);
        iFpeak = findFluxPeak(iF);
        
        %{ 
        pass.Vy = passVy; pass.eF = passeF; pass.iF = passiF;
        pass.MLAT = passMLAT; pass.time = passtime; pass.VyC = passVyC;
        pass.crb = crb; pass.eFpeaks = eFpeaks; pass.iFpeaks = iFpeaks;
        
        filename = sprintf('passDusk%2d.mat',j);
        save(filename, 'pass');
        clear pass
        %}
        
        gotstart = 0; gotend = 0; j=j+1;
        
%         axes();
%         
%         plot(passVyC, passtime, 'k-'); hold on
%         plot(passVyC(crb),passtime(crb), '-ro'); 
%         
%         pause
%         
%         axes();
%         
%         plot(passeF, passtime, 'k-'); hold on
%         plot(passeF(eFpeak), passtime(eFpeak), '-ro');
%         title('Electron flux: Peak and 1/e points');
%         legend('Electron flux: mW/m^{2}','1/e points and maximum point');
      
        pause
     
    end

end


        
function VyCorr = correctBaseline( Vy )
VyCorr = Vy - mean(Vy);
return


function crb = CRB(Vy)
[value, index] = min(abs(Vy));
crb = index;
return

function peak = findFluxPeak(Flx)
array = [transpose(1:length(Flx)),Flx];
k = find(~isnan(array(:,2)));
array=array(k,:);
%clear Flx k
array = sortrows(array,-2);
max = mean(array(1:10,2));
level = max*exp(-1);
array(:,2) = array(:,2)-level;
array(:,2) = abs(array(:,2));
array = sortrows(array,2);
array = array(1:10,:);
array = sortrows(array,1);
dawn = array(1,1);
dusk = array(10,1);
max = floor((dawn+dusk)/2);
peak = [dawn, max, dusk];
return
    

function [data1_ssj,DayFile] = read_Cenergy(Cenergy_file_path)
fname = Cenergy_file_path;
fid = fopen(fname,'r');
instr_date = fgets(fid);
defs = fgets(fid);
header = textscan(defs,'%s', 'delimiter', ',');
data = dlmread(fname,'',3,0);
DayFile = struct('data',data,'dataNames',header,'instrument', str2double(instr_date(2:3)),...
    'year',str2double(instr_date(5:8)), 'month',str2double(instr_date(9:10)),'day',str2double(instr_date(11:12)),...
    'number', data(:,1),'utsec',data(:,2),'glat',data(:,3),'glon',data(:,4),'alt',data(5,:),'mlat110',data(:,8),...
    'mlt',data(:,10),'Int_ele_eflx',data(:,11),'Int_ion_eflx',data(:,14));
data1_ssj = data;
    fclose(fid);

function [data,IDMFile] = 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);
IDMFile = 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));
fclose(fid);

function yesno = lastDayCheck(year, daynum)
if ((mod(year,4) == 0)) %leap year
    switch daynum 
        case {31,60,91,121,152,182,213,244,274,305,335,366} 
            yesno = 1;
        otherwise
            yesno = 0;
    end
else %non-leap year
    switch daynum 
        case{31,59,90,120,151,181,212,243,273,304,334,365}
            yesno = 1;
        otherwise
            yesno = 0;
    end
end

%Function which run IDL program cenergy.pro and outputs the name of the
%ascii file of all of the SSJ data for a particular day
function fileStr = createCenergyText(dateVec,sat)
global options
%options members: 
%omni_daily_path ssj_daily_path omni_source_path cenergy_path script_path 
%temp_path ssj_binary_path

    year = dateVec(1); yearS = int2str(year);
    month = dateVec(2); monthS = sprintf('%2.2d',month);
    day = dateVec(3); dayS = sprintf('%2.2d',day);
    daynumber = dayno([year,month,day,0,0,0]);
    satS = int2str(sat);
    
    %Directories
    temp_path = options.temp_path;
    ssj_daily_path = options.ssj_daily_path;
    script_path = options.script_path;
    log_path = options.log_path;
    
    %Files
    
    fileStr = strcat('cenergy_all_f',satS,'_',yearS,monthS,dayS,'.txt');
    logStr = strcat(log_path,'cenergy_f',satS,'_',yearS,monthS,dayS,'_full.log');
    
    lastdayswitch=lastDayCheck(year, daynumber);
    
    cmd = sprintf('%sdocenergy.sh %i %i %i %i %i %i config.txt',script_path, ...
        year, month, day, daynumber, sat, lastdayswitch);
    
    disp(cmd);
    exitstat = system(cmd)
    if (exitstat ~= 0)
        outstr = strcat('shell script exited with non-zero status check',logStr);
        disp(outstr)
    end
       
% --- Simple functions used to convert between date and array index
function out = calculateIndex(startDate, endDate, firstDate, stepPerD)
%lmk- startDate and endDate must have format [YYYY,MM,DD, 0, 0, 0] (standard matlab
%vector date). This only works with entire days. Returns
%[Startindex,endindex]. firstDate is earliest date in database.
%All use vector date format except stepPerD (steps per day) which must be integer.

sY = startDate(1);
eY = endDate(1);
fY = firstDate(1);

startD = dayno(startDate); endD = dayno(endDate); firstD = dayno(firstDate);

startI = 1; %start from firstDate at index 1
endI = stepPerD;

corr = 0; %leap year correction


    for y = fY:sY  
        if ((mod(y,4) == 0)), corr = 1; else corr = 0; end
        
        if (fY == sY)
           startI = startI + (startD - firstD)*stepPerD;
        elseif (y == sY)
           startI = startI + startD*stepPerD;          
        elseif (y == fY)
           startI = startI + ((365+corr)-firstD)*stepPerD;
        else
           startI = startI + (365+corr)*stepPerD;
        end
    end
    
    for y = fY:eY  
        if ((mod(y,4) == 0)), corr = 1; else corr = 0; end
        
        if (eY == fY)
            endI = endI + (endD - firstD)*stepPerD;
        elseif (y == eY)
            endI = endI + endD*stepPerD;          
        elseif (y == fY)
            endI = endI + ((365+corr)-firstD)*stepPerD;
        else
           endI = endI + (365+corr)*stepPerD;
        end
    end
    
    out = [startI, endI];
    
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;
