clear all
close all
numfiles=1;
fs=10e6;   % Sampling rate
ts=1/fs;   % Sampling time
ts=4e-3;
%filepath='C:\Locald~1\backup~2\FladeIsblinkData\AccumulationData\jun2_06\';
%gps_filename='acc_jul21_gps'; %don't need above filepath if gps_filename
gps_filename='nmea20060604';   %nmea.... is in current directory.

data_array=13001:50000;
numsam=length(data_array);
fftpow=nextpow2(numsam);
nfft=2^fftpow;
fx=(0:nfft-1)*fs/nfft;
fx=fx(1:nfft/2);
gausshigh=1-exp((-fx.^2)/(2*4.772085352e5^2));    %set up gaussian high pass filter weights for Dick
gausshigh(1,1:88)=1;    %drop filter (i.e. set = 1) for first 88 samples (r = -5 - 2.0176 m)
tm=4e-3;
tau=tm*fx/1.5e9;
r=((3e8/sqrt(3.17))*tau/2);  %sqrt(1.6) for standard density profiles, sqrt(3.17)for solid ice, sqrt(1) for free space
% note: in above calc, use surface correction factor of ~ 7.8503 m for solid
% ice, ~9.9 m for standard density profile and ~14.0 m for free space depth.
r=r-7.8503;     %surface correction factor (for r=0) at Flade Isblink assuming using sqrt(3.17)above
rconstrho=r;    %save r uncorrected for density variation for comparison with r corrected
%the following 5-line, near-surface correction for r using the Flade Isblink index of refraction
%profile works only assuming a surface correction factor of 7.8503 m is used
%so that r(1,230)=0, and that sqrt(3.17) = 1.78 index of refraction for
%solid ice was also used above.
deltar=(r(1,5999)-r(1,1))/(5999-1);
load('C:\Program Files\MATLAB71\work\refraction_index_accum_Flade_Isblink.mat');
for ii=[231:5999]
    r(1,ii)=deltar*1.78/refraction_index_FI(1,ii)+r(1,ii-1);
end

accumdata=[];
accumdata2=[];
time_data=[];
%for z=[0:1 3:11]
for z=[29:38]  %modify this according to the following (line 30) file no.s
    z

%fn = ['C:\Research\Prism\Summit\Accum\accum2005\jul21\acc_jul21_',num2str(z)];
fn = ['G:\FladeIsblinkData\AccumulationData\proc\jun4_',num2str(z)];
load(fn);

maxlen1=size(filtfin,1);
numrec1=size(filtfin,2);
gausshigh=gausshigh(1:maxlen1);
gausshigha=gausshigh'*ones(1,numrec1);
filtfin2=filtfin.*1./gausshigha;    %remove high pass filter built into radar (see Rohit thesis)
datad=diff(filtfin2);    %create 1st diffs; filtfin for regular gaussian high pass, filtfin2 with high pass removed (i.e. using reciprocal filter)

accumdata=[accumdata datad];
accumdata2=[accumdata2 filtfin2];   %filtfin for regular gaussian high pass, filtfin2 with high pass removed (i.e. using reciprocal filter)
time_data=[time_data timearr];

numfiles=numfiles+1;
end

datamed2=medfilt2(accumdata,[1 20]); %nominally set to [1 20], getting good results with [1 200].
windowSize=1;   %this line and next do a windowSize-point moving average on datamed2's columns (3-5 seemed to give best result)
datamed2=filter(ones(1,windowSize)/windowSize,1,datamed2);
datamed22=medfilt2(accumdata2,[1 1]); %nominally set to [1 20], getting good results with [1 200].
datamed22=datamed22.^2/50;  %voltage power in watts
disp('End Processing');
maxlen=size(accumdata,1);

num_rec=size(datamed2,2);
r=r(1:maxlen);

figure(1)
h=gca;
%image([],r,300*datamed2)
%image([],r,8e7*datamed2)   %nominal setting for accum data at Flade was 8e6
imagesc([],r,8e7*datamed2,[-2.0,2.0]) %plots voltage diffs in color with specified color range
%imagesc([],r,8e7*datamed2,[-0.8e4,-0.4e4]) %plots voltage diffs in color with specified color range
%imagesc([],r,8e7*datamed2,[-0.25e5,0.5e5]) %plot in color with specified color range
colorbar    %show color range
grid
figure(2)
imagesc([],r,10*log10(datamed22)+30, [-30,5]) %plots voltage power in dBm, plot in color with specified color range (dBm)
colorbar    %show color range
grid
set(h,'XTick',[0:100:1000])   %when commented out/ommitted, uses default (previously  used [0:200:2400])
%colormap(1-gray)  %when used, makes plot black and white
%axis([0 2500 -0.1 3])
ylabel('Depth (meters)')
%xlabel('Along track')


sizeimage = length(time_data);
% 3 m  max depth of 3 m
%startx=00;
%space=500;   % X axis spacing (in number of records) between GPS labels
%space=100;   % X axis spacing (in number of records) between GPS labels
%dis1 = 3.25; % Location on y-axis (in meters) where longitude GPS label is tagged
%pos = 0.075; % Offset from longitude position where lat and dist labels are placed (+pos for lat, -pos for dist) 
%fnt = 8;
%offset = round(space/5);
%sc=1;
%st=1;
%sc=4;
%sc=5;

% 30m
%startx=0;
%space=100;  %specifies lat/lon/dist label spacing on x axis (in records)
%dis1 = 32.25; %nominal 32.25, formerly 32.0 and 32.5 (originally)
%pos = 1;
%fnt = 8;
%offset = round(space/5);
%sc=.75; %scaling of lat/lon/dist/elev labels relative to lon labels formerly 1 (originally)
%st=1;
%correction = 10;  %= 15 for 2400, 10 for 1000, 2 for 100 and 200
% XL = xlim ;  %axis 'auto x'
%axis ([XL(1) XL(2) -0.1 30])  %allows matlab to autoscale the x axis
%axis(['auto x', 'auto x', -0.1 30])
%axis([0 1000 -5  30])  %specifies x axis min/max in records & y axis min/max in meters
% axis 'auto x'
%title('0.5-2 GHz Radar')
%title('June 4, 2006 (files 244-253) 6 km trek')
%title('June 4, 2006 (files 37-46)')
% % 210m
startx=00;
%space=500;
 space=100;   % X axis spacing (in number of records) between GPS labels
 dis1 = 212.5;
 pos = 1;
 fnt = 8;
 offset = round(space/5);
 sc=1;
 st=1;
 %axis([0 2500 -0.1 30])
 title('0.5-2 GHz Radar')
% axis([0 2500 180 210])
axis([0 1000 -5 210])


[latitude,longitude,elev,time]=getGPS(time_data,gps_filename);
latitude=latitude*180/pi;
longitude=longitude*180/pi;

for num=1:space:sizeimage

    if longitude(num) < 0
           k=text(startx+num-offset-correction,dis1,[num2str(-longitude(num),6),'W']);
    else
           k=text(startx+num-offset-correction,dis1,[num2str(longitude(num),6),'E']);
    end
    set(k,'fontsize',[fnt]); 
    set(k,'fontweight','bold');

    if latitude(num) < 0
        h=text(startx+num-offset-correction,dis1-(sc*pos),[num2str(-latitude(num),6),'S']);
    else
        h=text(startx+num-offset-correction,dis1-(sc*pos),[num2str(latitude(num),6),'N']);
    end
    set(h,'fontsize',[fnt]);
    set(h,'fontweight','bold');
    
    d=gcdist(latitude(num),longitude(num),latitude(st),longitude(st));
    d=d*1000;
    h2=text(startx+num-offset-correction,dis1+(sc*pos),[sprintf('%6.1f m',d)]);
    set(h2,'fontsize',[fnt]);
    set(h2,'fontweight','bold');
    el=elev(num); % also change lines 158-161 below
    h3=text(startx+num-offset-correction,dis1+(2*sc*pos),[sprintf('%6.2f m',el)]);  
    set(h3,'fontsize',[fnt]);
    set(h3,'fontweight','bold');
end

if sizeimage-num >= space*3/4  % if enough space after exit above, put in last x labels
	if longitude(sizeimage) < 0
        k=text(sizeimage-offset-correction+1,dis1,[num2str(-longitude(sizeimage),6),'W']);
    else
        k=text(sizeimage-offset-correction+1,dis1,[num2str(longitude(sizeimage),6),'E']);
    end    
   	set(k,'fontsize',[fnt]);
	set(k,'fontweight','bold');

    if latitude(num) < 0
    	h=text(sizeimage-offset-correction+1,dis1-(sc*pos),[num2str(-latitude(sizeimage),6),'S']);    
    else
    	h=text(sizeimage-offset-correction+1,dis1-(sc*pos),[num2str(latitude(sizeimage),6),'N']);
    end
    
	set(h,'fontsize',[fnt]);
	set(h,'fontweight','bold');
    d=gcdist(latitude(sizeimage),longitude(sizeimage),latitude(st),longitude(st));
    d=d*1000;
    h2=text(sizeimage-offset-correction+1,dis1+(sc*pos),[sprintf('%6.1f m',d)]);
	set(h2,'fontsize',[fnt]);
	set(h2,'fontweight','bold');
    el=elev(num);  % also change lines 130 - 133 abov%e
    h3=text(sizeimage-offset-correction+1,dis1+(2*sc*pos),[sprintf('%6.2f m',el)]);  
    set(h3,'fontsize',[fnt]);
    set(h3,'fontweight','bold');

end

