function varargout = CRBPlot(varargin)
% CRBPLOT M-file for CRBPlot.fig
%      CRBPLOT, by itself, creates a new CRBPLOT or raises the existing
%      singleton*.
%
%      H = CRBPLOT returns the handle to a new CRBPLOT or the handle to
%      the existing singleton*.
%
%      CRBPLOT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CRBPLOT.M with the given input arguments.
%
%      CRBPLOT('Property','Value',...) creates a new CRBPLOT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before CRBPlot_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to CRBPlot_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help CRBPlot

% Last Modified by GUIDE v2.5 04-Aug-2010 21:09:10

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @CRBPlot_OpeningFcn, ...
                   'gui_OutputFcn',  @CRBPlot_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

function get_CRB_peakFlux( datevec, sat, handles )
%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/liamk/data/cenergy/';
options.log_path='/home/liamk/data/log/';
options.script_path='/home/liamk/IDL_Files/';
options.cenergy_path='/home/liamk/IDL_Files/';
options.temp_path='/home/liamk/IDL_Files/';
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/liamk/data/cenergy/';
log_path='/home/liamk/data/log/';
script_path='/home/liamk/IDL_Files/';
cenergy_path='/home/liamk/IDL_Files/';
temp_path='/home/liamk/IDL_Files/';
idm_source_path='/hao/aim3/dmsp/dmsp_idm_4s/';



%file name ingredients
year = datevec(1); yearS = int2str(year);
month = datevec(2); monthS = sprintf('%2.2d',month);
day = datevec(3); dayS = sprintf('%2.2d',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');
ssj_bin_file = strcat('j',satS,year2S,'/j4f',satS,year2S,daynumS);
ssj_txt_file = strcat('cenergy_all_f',satS,'_',yearS,monthS,dayS,'.txt');

%make sure we have data for these days
ssj = exist(strcat(ssj_binary_path, ssj_bin_file),'file');
idm = exist(strcat(idm_source_path, idm_file),'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 = read_idm(strcat(idm_source_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 = createCenergyText([datevec(1),datevec(2),datevec(3),0,0,0], sat);

    disp('Reading DMSP ASCII file...');
    pause(2)
    
    data1_ssj = read_Cenergy(strcat(ssj_daily_path,ssj_txt_file));
else %Assume the ASCII file already exists and read from that
    
    ssj_txt_file = strcat('cenergy_all_f',satS,'_',yearS,monthS,dayS,'.txt');
    disp(ssj_txt_file);
    data1_ssj = read_Cenergy(strcat(ssj_daily_path, ssj_txt_file));
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 = data1_ssj(:,8); MLT = data1_ssj(:,10);
if(NH)
kkN = find(MLAT < MLAT_CUTOFF);
end
if(SH)
kkS = find(MLAT > -MLAT_CUTOFF);
end

% 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); EEFlx_N(kkN) = NaN; 
IEFlx_N = data1_ssj(:,14); 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
pe = 1;
ps = 1;

for i = istart:length(Vy)
    this = eF(i);
    last = eF(i-1);
    next = eF(i+1);
    
    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 (( pe - ps ) > 800)
        reasonable = 1;
    else 
        reasonable = 0;
    end
    
    if (gotstart && gotend && reasonable)
     
        passMLAT = MLAT(ps:pe);
        [value,pp] = max(abs(passMLAT));
        pp = ps+pp;
        %1st side of oval
        
        passVy = Vy(ps:pp);
        passeF = eF(ps:pp);
        passiF = iF(ps:pp);
        passMLAT = MLAT(ps:pp);
        passtime = UTTime(ps:pp);
        
      
        passVyC = correctBaseline(passVy);
        crb = CRB(passVyC,1);
        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
        %}
        set(gcf,'CurrentAxes',handles.axesLT);        
        plot(passtime,passVyC,'k.'); hold on
        plot(passtime(crb),passVyC(crb), '-ro'); hold off
        title('Dawn Side Cross-Track Velocity');
        legend('V_{y} km/s');
        
        set(gcf,'CurrentAxes',handles.axesLB);
        plot(passtime,passeF, 'k-'); hold on
        plot(passtime(eFpeak), passeF(eFpeak) , '-ro'); hold off
        title('Dawn Side Electron Flux');
        legend('Electron flux: mW/m^{2}','1/e points and maximum point'); hold off

        %2nd side of oval
        
        passVy = Vy(pp:pe);
        passeF = eF(pp:pe);
        passiF = iF(pp:pe);
        passMLAT = MLAT(pp:pe);
        passtime = UTTime(pp:pe);
        
        passVyC = correctBaseline(passVy);
        crb = CRB(passVyC,2);
        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;
        
        set(gcf,'CurrentAxes',handles.axesRT)
        plot(passtime, passVyC,'k.'); hold on
        plot(passtime(crb),passVyC(crb), '-ro'); hold off
        title('Dusk Side Cross-Track Velocity');
        legend('V_{y} km/s');
        
        set(gcf,'CurrentAxes',handles.axesRB)
        plot(passtime, passeF, 'k-'); hold on
        plot(passtime(eFpeak), passeF(eFpeak) , '-ro');
        title('Dusk Side Flux: Peak and 1/e points');
        legend('Electron flux: mW/m^{2}','1/e points and maximum point'); hold off
        
        waitfor(handles.nextbutton,'UserData',1);
        set(handles.nextbutton,'UserData',0); 
        clear passMLAT passVy passVyC passeF passtime passiF
    end
       

end


        
function VyCorr = correctBaseline( inVy )
inVy2 = inVy(find(~isnan(inVy)));
VyCorr = inVy - mean(inVy2); 
return


function crbout = CRB(inVy, DuskDawn)
[Max, Maxi] = nanmax(inVy);
[Min, Mini] = nanmin(inVy);
if(DuskDawn == 2) 
    [value, index] = nanmin(abs(inVy(Mini:Maxi)));
    crbout = index+Mini;
    out = sprintf('Min: %f Max: %f index: %f CRB: %f',Mini,Maxi,index,crbout);   
end
if(DuskDawn == 1) 
    [value, index] = nanmin(abs(inVy(Maxi:Mini)));
    crbout = index+Maxi;
    out = sprintf('Max: %f Min: %f index: %f CRB: %f',Maxi,Mini,index,crbout);
end
disp(out);
clear inVy

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 data = read_Cenergy(filename)
%there are 53300 lines in a single day of DMSP SSJ4 data output from
%cenergy.pro
data = NaN(86400,16);
fid = fopen(filename,'rt');
if (fid ~= -1)
        for i = 1:2
            tline = fgets(fid); 
        end          
        i = 1;
        j = 1;
        h = 2;
        while ~feof(fid)
            tline = fgets(fid);
            a = sscanf(tline,'%f');
            if (length(a) == 16) 
                if (i ~= (a(2)+1)) %missing data check; a(2) is time column
                    warn = sprintf('cenergy.pro output file is missing data beginning at line %d', j+h);
                    disp(warn);
                    i = a(2)+1;
                end
            data(i,:) = a; %jth row of data
            j = j+1;
            i = i+1;
            else 
            tline = fgets(fid); %clear header     
            h = h+1;
            end
        end
    fclose(fid);
end


function data = read_idm(filename)
data = NaN(86400,13);
fid = fopen(filename,'rt');
if (fid ~= -1)
        tline = fgets(fid); %clear 1 line header
        i=1;
        last=1;
        while ~feof(fid)
            tline = fgets(fid);
            a = sscanf(tline,'%f'); 
            if (length(a) == 13)
                if (last ~= a(1)-4) %missing data check; a(1) is time column
                    warn = sprintf('idm output file inconsistant between UT secs: %i and %i', last,a(1));
                    disp(warn);
                    last = a(1);
                    i=a(1); 
                    
                else
                    last = a(1);
                    i=a(1);
                    
                end
                data(i,:) = a; %ith row of data
            else
                disp(tline);
            end
                
                
        end
end
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 runs 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;

% --- Executes just before CRBPlot is made visible.
function CRBPlot_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to CRBPlot (see VARARGIN)

% Choose default command line output for CRBPlot
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes CRBPlot wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = CRBPlot_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in nextbutton.
function nextbutton_Callback(hObject, eventdata, handles)
% hObject    handle to nextbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(hObject,'UserData', 1);

% --- Executes on button press in gobutton.
function gobutton_Callback(hObject, eventdata, handles)
% hObject    handle to gobutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
get_CRB_peakFlux([2001 01 29 0 0 0], 13, handles)

