% This Plots the processed Deep Ice Radar Data 
% Change following parameters to match the data location time:
% inputfolder, gpsfile, figname, titletext, starttime, endtime - other
% parameters may be changed as needed
% A png file with 600 dpi resolution is saved with 'figname' 
% !! Change 'figname' to prevent overwriting existing file!!
% 09/25/2014
% Mike Wolovick, Tej Dhakal, 
% LDEO, Columbia University

clc
clear all
close all
tic

TxInd1us = 321;        
TxInd10us = 321; 

% Parameters: 
inputfolder='.../Flight5/Pt20-22\'; %Set this to data loacation path
gpsfile='.../glnd3_ft05_200hz_gps_v1.mat';  % Set this to gps mat file for the flight

filekey='FL05_1us_Ch2*.mat';  % search key.  Should contain "*"

figname=[inputfolder,'G3FL04_20140718_1DSAR_1usCh2_Pt5toPt6_test.png']; %% change this for every plot to avoid overwriting old one
titletext='GR3, FL05, ch2, 1\mus,Eqip Line, Pt20 to Pt22';

if(filekey(6:7) == '10')
    TxInd = TxInd10us;
else
    TxInd = TxInd1us;
end


clims=[-60,30];                    % Change limits to change image contrast
ylims=[-200, 5000];                   % m
cmap='gray';                       % valid colormap string
invertcaxis=1;                     % 0 or 1
elevcorrection=1;                  % 0 or 1
timeoffset= 0;%-14;                    % s added to gps timestamp
c_air=3e8;                         % m/s
c_ice=1.68e8;                      % m/s
cutterH=1;                         % horizontal decimation
starttime=  '19-07-2014 14:44:00';   % UTC time of data to start plotting at.  Format: 'dd-mm-yyyy HH:MM:SS'
endtime=    '19-07-2014 14:45:00';     % UTC time of data to end plotting at.  Format: 'dd-mm-yyyy HH:MM:SS'
pagesize=[15,8];                   % inches
resolution=600;                    % dpi
% Find input files:
cd(inputfolder) 
inputfiles=dir(filekey);
numfiles=length(inputfiles);
if numfiles==0 
    error('Could not locate input files.')
end
 
%Check if gps data exist:
if exist('gpsfile','var')
     
    % Load gps data:
    load(gpsfile,'Lat','Lon','Alt','dateUTC','sodUTC');
    GPSUTC = 86400*((dateUTC + sodUTC/86400) - datenum('01-01-1970 00:00:00','dd-mm-yyyy HH:MM:SS'));
end

% Compute start and end echoTimeUTC:
startutc=86400*(datenum(starttime,'dd-mm-yyyy HH:MM:SS')-datenum('01-01-1970 00:00:00','dd-mm-yyyy HH:MM:SS'));
endutc=86400*(datenum(endtime,'dd-mm-yyyy HH:MM:SS')-datenum('01-01-1970 00:00:00','dd-mm-yyyy HH:MM:SS'));

% Loop through and plot input files:
lastind=0;
isfirstfile=1;
for ii=1:numfiles-1
    
    % Load echoTimeUTC time:
    load(inputfiles(ii).name,'echoTimeUTC')
   
    % load data:
    if isempty(echoTimeUTC)==0
        if echoTimeUTC(end)>startutc && echoTimeUTC(1)<endutc
            load(inputfiles(ii).name,'CG','TWT')
               disp(inputfiles(ii).name);
        else
            continue
        end
    else
        continue
    end
    % Interpolate gps:
    if exist('gpsfile','var')
        thislat=interp1(GPSUTC,Lat,echoTimeUTC,'linear');
        thislon=interp1(GPSUTC,Lon,echoTimeUTC,'linear');
        thiselev=interp1(GPSUTC,Alt,echoTimeUTC,'linear');
    end
    % Define first file:
    if isfirstfile==1
        if exist('gpsfile','var')
            refelev=thiselev(1);
        end
        mintime=86400*floor(echoTimeUTC(1)/86400);
        firsttime=(echoTimeUTC(1)-mintime)/86400;
        isfirstfile=0;
    end
    % Perform elevation correction:
    if elevcorrection==1 && exist('gpsfile','var')
        DepthScale=-(TWT - TWT(TxInd))*c_ice/2;
        for jj=1:length(echoTimeUTC)
            elevoffset=thiselev(jj)-refelev;
            indoffset=elevoffset/((TWT(2)-TWT(1))*c_air/2);
            CG(:,jj)=interp1([1:size(CG,1)]',CG(:,jj),[1:size(CG,1)]'+indoffset,'linear');
        end
    else
        DepthScale=(TWT - TWT(TxInd))*c_ice/2;
    end
    
    % Determine data within elevation range:
    InRange= abs(DepthScale)>ylims(1)& abs(DepthScale)<ylims(2);
    
    % plot data:
    figure(1)
    if invertcaxis==1
        imagesc((echoTimeUTC-mintime)/86400,DepthScale(InRange),-20*log10(abs(CG(InRange,1:cutterH:end))),-clims(2:-1:1))
    else
        imagesc((echoTimeUTC-mintime)/86400,DepthScale(InRange),20*log10(abs(CG(InRange,1:cutterH:end))), clims); 
    end
    hold on
    eval(['colormap(',cmap,')'])
    if isempty(echoTimeUTC)==0
        xlim([firsttime,(echoTimeUTC(end)-mintime)/86400])
    end
    xlabel('Flight Time, UTC')
    datetick2('x','HH:MM','keeplimits')
    if elevcorrection
        set(gca,'ydir','normal')
        ylabel('Relative Elevation, m')
        ylim([-ylims(2),-ylims(1)])
    else
        set(gca,'ydir','reverse')
        ylabel('Range in Ice, m')
        ylim(ylims)
    end
    title(titletext)
    set(gca,'TickDir','out')
    % Draw figure:
    drawnow
%     shading interp
end

% Save figure:
set(gcf,'PaperSize',pagesize)
set(gcf,'PaperPosition',[0,0,pagesize])
print('-dpng',figname,['-r',num2str(resolution)])
toc