%Use Amplitude of Calibration line to get AlphaVsTime
%Produce list in ascii format

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Reference G at cal freq chosen
load L1OLG_753424982
%  G    448001x1                  7168016  double array (complex)
%  f    448001x1                  3584008  double array
%  g0        1x1                        8  double array
%  g1        1x1                        8  double array

%Calibration frequency
fcal=927.7;

%Ref time for alphas, betas
tr=753424982;
intv=[-5,5]; %interval to calculate mean line amplitude
%lineref=0.0202;

%Input data file: from SenseMon log file
[T,range,range2,gain,CalLine,Alpha0,Beta0]=...
    textread('L1_SenseMonitor_CumLog_S3.txt','%f%f%f%f%f%f%f', ...
	     'commentstyle','shell');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Algorithm to generate Alphas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Calibration Line Amplitude
%--------------------------


%Reference Line Amplitude 
ir=find(abs(T-tr)==min(abs(T-tr)));
ix=ir+intv;
CalLine0=mean(CalLine(ix));
CalLine0err=std(CalLine(ix));
disp(['CalLine REF=' num2str(CalLine0,3) ' +/- ' ...
      num2str(CalLine0err,2) ' between GPS ' ...
      num2str(T(ix(1))) ' and ' num2str(T(ix(end)))])

%ratio of calline to reference calibration lines
R=CalLine./CalLine0;

%Betas: SenseMon betas, referenced to ref. time
%-------------------------------------------

%break in SenseMon ref for beta
ib=max(find(T<=752560033));
Beta=[Beta0(1:ib)*0.00264/0.001625;Beta0(ib+1:end)];

%normalize to ref time
Beta=Beta/Beta(ir);


%open loop gain at calibration frequency
%----------------------------------------
f0=find(abs(f-fcal)==min(abs(f-fcal)));
ix=[f0-10:f0+10];
G0 = interp1(f(ix),G(ix),fcal);


%Get alpha*beta from quadratic equation
%--------------------------------------

A=zeros(size(T));
B=zeros(size(T));
C=zeros(size(T));

AB=zeros(size(T));

ix=find(R~=0);

A(ix) = R(ix).^2.*abs(G0)^2 -(abs(1+G0))^2./Beta(ix).^2;
B(ix) = 2* R(ix).^2 * real(G0);
C(ix) = R(ix).^2;

%Solution to alpha*beta
AB(ix) = (-B(ix) - sqrt(B(ix).^2-4*A(ix).*C(ix)))./(2*A(ix));

%set all imaginary alphas to zero
ix=find(imag(AB)~=0);
AB(ix)=0;


%Find unphysical values and set them to zero
ilow=find(AB<1/g0);
ihigh=find(AB>1/g1);

AB(ilow)=0;
AB(ihigh)=0;


%Alpha coefficient
%-----------------


ix=find(AB~=0 & Beta~=0);
Alpha=zeros(size(T));
Alpha(ix)=AB(ix)./Beta(ix);

%Plot
%---

figure
plot(T,AB,'.')
xlabel('gps')
ylabel('\alpha\beta')
title('L1 calibration during S3')
shg

pause(.1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Save data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Data=[T AB Alpha Beta R];

Info={'GPS','Alpha*Beta','Alpha','Beta','CalLine/Ref'};
save L1AlphaBeta.mat Data Info

format short g

%save -ascii L1AlphaBeta.txt Data

%Pretty print  the data in an ascii file (slow part!)

disp('writing data...')
fid=fopen('L1AlphaBeta.txt','w');

fprintf(fid,'%s\t','GPS Time');
fprintf(fid,'%s\t','alpha*beta');
fprintf(fid,'%s\t','alpha');
fprintf(fid,'%s\t','beta');
fprintf(fid,'%s\t','CalLine/Ref');


for j=1:size(Data,1)
  fprintf(fid,'%s\n',num2str(Data(j,:)));
end

fclose(fid);