%Shane Crockett - last updated 17 JUN 98

c = zeros(1,65);
c(1:20) = (.5*10^17:.5*10^17:10*10^17)./10^17;
b = -1.34*(100000000^2)./(f.^2);
[a] = b*c;
a = a(:,1:20)./(.1599);																	% Delta t
w = abs(round(a));																%integer values of delta t

%Determining the instantaneous power of the dechirped spectrogram

for x = 1:20
   for cc = 1:65;
      B1 = shift(supcar(cc,:),w(cc,x));
      DE(cc,:) = B1;
   end
   merit(x) = sum(sum(DE).^4);
end


%Refining the figure of merit by varying the trial TEC by a step ten fold finer around the highest FOM

optimum = find(merit == max(merit));
refine = zeros(1,65);
interval = (1/30);%*10^17;
if optimum == 1
   refine(1:31) = c(optimum):interval:c(optimum + 2);
elseif optimum == 20
   refine(1:31) = c(optimum - 2):interval:c(optimum);      
else
   refine(1:31) = c(optimum - 1):interval:c(optimum + 1);
end
[anew] = b*refine;
anew = anew(:,2:30)./(.1599);																	% Delta t
w1 = abs(round(anew));

%Determining the instantaneous power of the dechirped spectrogram for the refined values of TEC

for x = 1:29
   for cc = 1:65;
      RB1 = shift(supcar(cc,:),w1(cc,x));
      RDE(cc,:) = RB1;
   end
   rmerit(x) = sum(sum(RDE).^4);
end
TEC = [c(1:optimum-1) refine(2:30) c(optimum+1:20)];
FOM = [merit(1:optimum-1) rmerit(1:29) merit(optimum + 1:20)];
plot(TEC,FOM/max(FOM), 'bd')
newoptimum = find(rmerit == max(rmerit));

for cc = 1:65;
      B1 = shift(supcar(cc,:),w1(cc,newoptimum));
      DChirp(cc,:) = B1;															%Optimally dechirped spectrum is called DChirp   
end


%Shifting back the entire spectrogram
%Plotting the optimally dechirped, time-shifted spectrogram

reshift = w1(65,newoptimum);
DChirp = [DChirp(:,length(DChirp) - reshift:length(DChirp)) DChirp(:,1:length(DChirp)-reshift-1)];
sr = image(t, f, log(DChirp));

titl = [filename ' Shifted with Optimum TEC'];
title(titl)

xlabel('Time (sec)')
ylabel('Frequency (MHz)')

set(sr, 'CDataMapping', 'scaled')
set(gca, 'YDir', 'normal', 'YTicklabel', [30 35 40 45 50 55])%, 'XTickLabel', [0 50 100 150 200 250 300 350 400])  
set(gcf, 'Position', [50 25 400 350])

colorbar

% Administrative to free space
clear a b c w x cc B1 DE anew w1 merit optimum refine interval RB1 RDE rmerit opt reshift sr supcar