function [quad]=readQuadCentroidBTK(varargin) 

% Reads the engineering files produced by pmpsdGUI
% 4 quadcells record the beam displacement on 
% ip3: FSUA/AMBER/MIDI on Telescope 2
% ip1: FSUA/AMBER/MIDI on Telescope 1
% ip4: FSUB on Telescope 2
% ip2: FSUB on Telescope 1

% Display results and compute statistics of a selected input channel.

% Input parameters
    % filename:  example MetQuadCentroids2007-08-06T13.13.06.000000.txt
    % input_channel: 1,2 3 or 4
    % draw: 'yes' or 'no'
% Output:
    %quad structure containing all info about the selected channel
    
% example:
%[quad]=readQuadCentroidBTK('filename','MetQuadCentroids2007-08-06T13.13.06.txt','input_channel',1,'draw','yes')
%[quad]=readQuadCentroidBTK('input_channel',1,'draw','yes'): a window will
%                                                         put up to
%                                                         selected the file
% SLE August 2007
% december 2007: adapt to new file format (user guiding offsets removed
%                from centroid file.
%                Background removed for data analysis
% March 2009:  engineering file now contains demodX and demodY computed
%              by PRIMET BTK algorithm. Add Plot for radial displacements
%
%
% July 2009:  SLE, the header contains 5 additional lines to show the guiding
%              % Parameters of the controllers
%              % IP3 IM Gain A1 A2 B1 B2 -1 0 0 -1 0 -1.89039 0.89215 -1.99906 0.999058
%              % IP1 IM Gain A1 A2 B1 B2 1 0 0 -1 0 -1.89039 0.89215 -1.99906 0.999058
%              % IP4 IM Gain A1 A2 B1 B2 1 0 0 -1 0 -1.89039 0.89215 -1.99906 0.999058
%              % IP2 IM Gain A1 A2 B1 B2 1 0 0 -1 0.01 0 0 -1 0
%
% November 2010: SLE 
%                     [NoiseAfterQuad,count] removed from file Header
%                     (as of 18/9/2010 UT time~1h00)
%                      readQuadcentroidBTK/correction of september'10
%                     (COM12) renamed->OLD
%                     readQuadcentroidBTK/correction updated
%Feb 2011: SLE
%                       add plot std=f(bandwidth) for the radial
%                       displacement

 
%Default if no argument is given
filename='user_selected';
input_channel=3;
draw=1;

for k=1:2:length(varargin)
    switch lower(varargin{k})
        case 'filename'
            filename=varargin{k+1};
        case 'input_channel'
            input_channel=varargin{k+1};
        case 'draw'
            draw=strcmpi(varargin{k+1},'yes');
        otherwise
            error(sprintf('property ''%s'' not supported',varargin{k}));
    end
end

% Open the right file
if strcmp('user_selected',filename)
    [filename, PathName] = uigetfile('MetQuadCentroids*', 'Click on the data File to process'); 
    eval([ 'cd ''',PathName,'''']);
end;
fid = fopen(filename,'r');


% Read Header
%----------------------
    %first line:     '%% Date/Start time: %s' if string starts with a %
    %                ' Date/Start time: %s' if string starts without a %
    %                '%*c Date/Start time: %s'  to ignore first character

   
    
 %extract data common to all selected quadcell
[date_time_txt,count]=fscanf(fid,'%% Date/Start time: %s',[1,1]);
fgets(fid);% line 2  (alternative line2=fgets(fid), in that case line2 contains the text of line 2
[fs,count]=fscanf(fid,'%% Sampling frequency in Hz: %f',[1,1]);
 %[fs,count]=fscanf(fid,'%% Sampling Rate : %f',[1,1]);
fgets(fid); %line 3
[fc,count]=fscanf(fid,'%% Correction frequency in Hz: %f',[1,1]);
fgets(fid); %line 4
[N_samples,count]=fscanf(fid,'%% Number of Samples: %d',[1,1]);
fgets(fid);%line 6 (Configuration)

IP3.filename=filename;IP1.filename=filename;IP4.filename=filename;IP2.filename=filename;
IP3.SampFreq=fs;IP1.SampFreq=fs;IP4.SampFreq=fs;IP2.SampFreq=fs;
IP3.CorrFreq=fc;IP1.CorrFreq=fc;IP4.CorrFreq=fc;IP2.CorrFreq=fc;


%extract quadcell dependant information
fgets(fid);
[CHA,count]=fscanf(fid,'%% CHA: Instrument= %s',[1,1]);
IP3.IP='IP3';IP1.IP='IP1';
IP3.Instrument=CHA;IP1.Instrument=CHA;
fgets(fid);
[IP3.Active,count]=fscanf(fid,'%% IP3 Active %s',[1,1]);
fgets(fid);
[IP3.Serial,count]=fscanf(fid,'%% IP3 Serial %s',[1,1]);
fgets(fid);
[IP3.Tel,count]=fscanf(fid,'%% IP3 Tel: %s',[1,1]);
fgets(fid);
[IP3.Trans,count]=fscanf(fid,'%% IP3 Trans Gain %f',[1,1]);
fgets(fid);
[IP3.PGA,count]=fscanf(fid,'%% IP3 PGA Gain %f',[1,1]);
fgets(fid);
[IP3.BKG,count]=fscanf(fid,'%% IP3 Background %f %f %f %f',[4,1]);

fgets(fid);
[IP1.Active,count]=fscanf(fid,'%% IP1 Active %s',[1,1]);
fgets(fid);
[IP1.Serial,count]=fscanf(fid,'%% IP1 Serial %s',[1,1]);
fgets(fid);
[IP1.Tel,count]=fscanf(fid,'%% IP1 Tel: %s',[1,1]);
fgets(fid);
[IP1.Trans,count]=fscanf(fid,'%% IP1 Trans Gain %f',[1,1]);
fgets(fid);
[IP1.PGA,count]=fscanf(fid,'%% IP1 PGA Gain %f',[1,1]);
fgets(fid);
[IP1.BKG,count]=fscanf(fid,'%% IP1 Background %f %f %f %f',[4,1]);

fgets(fid);
[CHB,count]=fscanf(fid,'%% CHB: Instrument= %s',[1,1]);
IP4.IP='IP4';IP2.IP='IP2';
IP4.Instrument=CHB;IP2.Instrument=CHB;

fgets(fid);
[IP4.Active,count]=fscanf(fid,'%% IP4 Active %s',[1,1]);
fgets(fid);
[IP4.Serial,count]=fscanf(fid,'%% IP4 Serial %s',[1,1]);
fgets(fid);
[IP4.Tel,count]=fscanf(fid,'%% IP4 Tel: %s',[1,1]);
fgets(fid);
[IP4.Trans,count]=fscanf(fid,'%% IP4 Trans Gain %f',[1,1]);
fgets(fid);
[IP4.PGA,count]=fscanf(fid,'%% IP4 PGA Gain %f',[1,1]);
fgets(fid);
[IP4.BKG,count]=fscanf(fid,'%% IP4 Background %f %f %f %f',[4,1]);

fgets(fid);
[IP2.Active,count]=fscanf(fid,'%% IP2 Active %s',[1,1]);
fgets(fid);
[IP2.Serial,count]=fscanf(fid,'%% IP2 Serial %s',[1,1]);
fgets(fid);
[IP2.Tel,count]=fscanf(fid,'%% IP2 Tel: %s',[1,1]);
fgets(fid);
[IP2.Trans,count]=fscanf(fid,'%% IP2 Trans Gain %f',[1,1]);
fgets(fid);
[IP2.PGA,count]=fscanf(fid,'%% IP2 PGA Gain %f',[1,1]);
fgets(fid);
[IP2.BKG,count]=fscanf(fid,'%% IP2 Background %f %f %f %f',[4,1]);

% added July 2009
fgets(fid);%Line parameter of the controllers
fgets(fid);
[IP3.CtrlParam,count]=fscanf(fid,'%% IP3 IM Gain A1 A2 B1 B2 %f %f %f %f %f %f %f %f %f',[9,1]);
fgets(fid);
[IP1.CtrlParam,count]=fscanf(fid,'%% IP1 IM Gain A1 A2 B1 B2 %f %f %f %f %f %f %f %f %f',[9,1]);
fgets(fid);
[IP4.CtrlParam,count]=fscanf(fid,'%% IP4 IM Gain A1 A2 B1 B2 %f %f %f %f %f %f %f %f %f',[9,1]);
fgets(fid);
[IP2.CtrlParam,count]=fscanf(fid,'%% IP2 IM Gain A1 A2 B1 B2 %f %f %f %f %f %f %f %f %f',[9,1]);
fgets(fid);%Line definition of the columns
fgets(fid);
fgets(fid);

% Read data  
%---------------------------


%[data,count]=fscanf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f',[17,inf]);
[data,count]=fscanf(fid,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f',[25,inf]);

fclose(fid);

data=data';% put back data in columns as in the eng. file

% IP3.data=data(:,2:5);
% IP1.data=data(:,6:9);
% IP4.data=data(:,10:13);
% IP2.data=data(:,14:17);

IP3.data=data(:,2:7);
IP1.data=data(:,8:13);
IP4.data=data(:,14:19);
IP2.data=data(:,20:25);



% Select the input channel and analyse the data with background removal
eval(['quad=', ['IP' num2str(input_channel(1))],[';']]);
eval(['BKG=', ['IP' num2str(input_channel(1))],['.BKG;']]);

quad.seconds=data(:,1);

% WITH BACKGRUOND REMOVAL
quad.Q1=quad.data(:,1)-BKG(1);
quad.Q2=quad.data(:,2)-BKG(2);
quad.Q3=quad.data(:,3)-BKG(3);
quad.Q4=quad.data(:,4)-BKG(4);

%WITHOUT BACKGROUND REMOVAL
% quad.Q1=quad.data(:,1);
% quad.Q2=quad.data(:,2);
% quad.Q3=quad.data(:,3);
% quad.Q4=quad.data(:,4);


quad.demodX=quad.data(:,5);
quad.demodY=quad.data(:,6);

quad.SUM=quad.data(:,1)+quad.data(:,2)+quad.data(:,3)+quad.data(:,4);
quad.X=2/pi*(quad.data(:,2)+quad.data(:,3)-quad.data(:,1)-quad.data(:,4))./quad.SUM;
quad.Y=2/pi*(quad.data(:,1)+quad.data(:,2)-quad.data(:,3)-quad.data(:,4))./quad.SUM;
quad.R=sqrt((quad.X-quad.X(1)).^2+(quad.Y-quad.Y(1)).^2); % radial displacement w.r.t first data point to get size of the cloud (max value)
quad.Pos=sqrt(quad.X.^2+quad.Y.^2); % to get mean position of cloud w.r.t center of quadcell


% Compute PSD's 
nfft=round(length(quad.R));
noverlap = round(0.9*nfft);
[quad.PSDX,fre]=psdbk(quad.X,fs,nfft,noverlap);
[quad.PSDY,fre]=psdbk(quad.Y,fs,nfft,noverlap);
[quad.PSDR,fre]=psdbk(quad.R,fs,nfft,noverlap);
quad.fre=fre;

quad= rmfield(quad, 'data')


% STD calculation vs bandwidth for Radial motion
f_cut=5; 
    index=find(fre>=f_cut);
    f_cut=fre(index(1));    clear index
    % f_cut is replaced by nearest frequency in vector fre
    % otherwise for k=1, 
    % std_f_relativeToFcut(k)=rmspsd(fre,PSD_Data,f_cut,f_cut)/rms_abov
    % e_f_cut will fail
    
% cut frequency for relative STD computation
if fre(end)>=f_cut
rms_above_f_cut=rmspsd(fre,quad.PSDR,f_cut,fre(end)) ; %fre(end)=fs/2

% STD vs BW calculation
f_plot=f_cut:fre(end) ; % if fre(end)=fs/2 is not an integer, f_plot stops at floor(fre(end))!
                        % alternative: loop on fre, but fre can have a lot
                        % of points !!
std_f_relativeToFcut=[];
for k=1:length(f_plot)
 std_f_relativeToFcut(k)=rmspsd(fre,quad.PSDR,f_cut,f_plot(k))/rms_above_f_cut;
end

else fprintf('STD vs BW not possible: largest freq<fcut !  fcut must be changed')
end





% Display the results
fprintf('    Mean Q1 = %e V   PV Q1= %e V   rms Q1= %e V \n',mean(quad.Q1),max(quad.Q1)-min(quad.Q1),std(quad.Q1));
fprintf('    Mean Q2 = %e V   PV Q2= %e V   rms Q2= %e V \n',mean(quad.Q2),max(quad.Q2)-min(quad.Q2),std(quad.Q2));
fprintf('    Mean Q3 = %e V   PV Q3= %e V   rms Q3= %e V \n',mean(quad.Q3),max(quad.Q3)-min(quad.Q3),std(quad.Q3));
fprintf('    Mean Q4 = %e V   PV Q4= %e V   rms Q4= %e V \n',mean(quad.Q4),max(quad.Q4)-min(quad.Q4),std(quad.Q4));
fprintf('    Mean Sum= %e V \n',mean(quad.SUM));
fprintf('\n')
fprintf('    Mean X = %e  p   PV= %e p   rms= %e p \n',mean(quad.X),max(quad.X)-min(quad.X),std(quad.X));
fprintf('    Mean Y = %e  p   PV= %e p   rms= %e p \n',mean(quad.Y),max(quad.Y)-min(quad.Y),std(quad.Y));
fprintf('    Mean R = %e  p   R_Pk= %e p   R_rms= %e p \n',mean(quad.R),max(quad.R),std(quad.R));
fprintf('                                                           R_rms above f_cut= %.0f Hz= %e \n',f_cut,rms_above_f_cut');
fprintf('    Mean position= %e p \n',mean(quad.Pos));

fprintf('    Total recording time= %e sec \n',quad.seconds(end)-quad.seconds(1));
fprintf('    Mean Sampling= %e sec \n',mean(diff(quad.seconds)));
fprintf('    P_V Sampling= %e sec \n',max(diff(quad.seconds))-min(diff(quad.seconds)));


if draw
    
   Nfig=8;
   scrsz = get(0,'ScreenSize'); %[1 1 width height]
   SizeFig=[scrsz(3)/(Nfig/2) scrsz(4)/2.5];%[width height]
   
   %built title for all figures
   line1=quad.filename;
   line2=['fs: ',num2str(quad.SampFreq),' fc: ',num2str(quad.CorrFreq)];
   line3=[quad.IP,' ', quad.Instrument,' ',quad.Serial,' ', quad.Tel];
   line4=['Trans: ',num2str(quad.Trans),' PGA: ',num2str(quad.PGA)];
   line5=['IM Gain A1 A2 B1 B2=',num2str(quad.CtrlParam')];
   
   %Plot figures
   h=figure(1);set(h,'Position',[1 scrsz(4)/2 SizeFig]); %[left bottom width height]
   set(gca,'FontSize',13);
   plot(quad.seconds,[quad.Q1,quad.Q2,quad.Q3,quad.Q4]);grid
   legend('Q1','Q2','Q3','Q4')
   xlabel('Time in sec');ylabel('Voltage in V');
   title({[line1],[line2],[line3],[line4],[line5]});
 
   h=figure(2);set(h,'Position',[SizeFig(1) scrsz(4)/2 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.seconds,quad.SUM);grid
   legend('Sum in V')
   xlabel('Time in sec');ylabel('Voltage in V');
   title({[line1],[line2],[line3],[line4]}); 
   
   h=figure(3);set(h,'Position',[2*SizeFig(1) scrsz(4)/2 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.seconds,[quad.X,quad.Y]);grid
   legend('X','Y')
   xlabel('Time in sec');ylabel('Displacement unitless');
   title({[line1],[line2],[line3],[line4]});
   
    h=figure(4);set(h,'Position',[2*SizeFig(1) scrsz(4)/2-20 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.seconds,[quad.Pos,quad.R]);grid
   legend('Radial Pos. wrt quad center','Radial Pos. wrt 1st pt')
   xlabel('Time in sec');ylabel('Displacement unitless');
   title({[line1],[line2],[line3],[line4]});
   
   
   theta = 0:0.01:(2*pi);
   x = 0.1*cos(theta);
   y = 0.1*sin(theta);
   h=figure(5);set(h,'Position',[3*SizeFig(1) scrsz(4)/2 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.X,quad.Y,'r+',x,y,'k');grid
   axis([-1 1 -1 1]);axis('square')
   xlabel('X');ylabel('Y');
   title({[line1],[line2],[line3],[line4]});
   figName=[ quad.filename(1:35),'-',quad.IP,'.jpg']
   print('-djpeg90','-r300',figName);

   
   
   h=figure(6);set(h,'Position',[1 1 SizeFig]);
   subplot(2,1,1)
        set(gca,'FontSize',13);
        semilogy(quad.fre,[sqrt(quad.PSDX), sqrt(quad.PSDY), sqrt(quad.PSDR)]);grid
        legend('X','Y','R')
        xlabel('freq in Hz');ylabel('PSD Displacement unitless/sqrt(Hz)');
        title({[line1],[line2],[line3],[line4]});
   subplot(2,1,2)
        set(gca,'FontSize',13);
        plot(f_plot,std_f_relativeToFcut*100,'m.-');grid
        xlabel('Frequency in Hz');
        ylabel('% of Std above fcut');
        ii=find(std_f_relativeToFcut<0.99);
        axis([0 f_plot(ii(end)) 0 100])
        legend({['Std above fcut=' num2str(f_cut) 'Hz = ' num2str(rms_above_f_cut) '']})
   
   

   h=figure(7);set(h,'Position',[SizeFig(1) 1 SizeFig]);
   set(gca,'FontSize',13);
   plot(diff(quad.seconds));grid
   xlabel('Index');ylabel('Sampling time in sec');
   title({[line1],[line2],[line3],[line4]});
   
   h=figure(8);set(h,'Position',[2*SizeFig(1) 1 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.seconds,quad.demodX,'r',quad.seconds,quad.demodY,'b',quad.seconds,quad.demodX-quad.demodY,'g');grid
   legend('DemodX','DemodY','DemodX-DemodY=ellipticity')
   xlabel('Time in sec');ylabel('Demod');
   title({[line1],[line2],[line3],[line4],['DemodX&DemodY >threshold1?'],['DemodX-DemodY <threshold2(ellipticity)?']});
 

%    h=figure(10);set(h,'Position',[3*SizeFig(1) 1 SizeFig]);
%    set(gca,'FontSize',13);
%    plot(quad.X,quad.demodX,'r');grid
%    legend('DemodX')
%    xlabel('X');ylabel('DemodX unitless');
%    title({[line1],[line2],[line3],[line4]});
%    
%    h=figure(11);set(h,'Position',[3*SizeFig(1) 1 SizeFig]);
%    set(gca,'FontSize',13);
%    plot(quad.Y,quad.demodY,'b');grid
%    legend('DemodY')
%    xlabel('Y');ylabel('DemodY unitless');
%    title({[line1],[line2],[line3],[line4]});
%    
%    h=figure(12);set(h,'Position',[3*SizeFig(1) 1 SizeFig]);
%    set(gca,'FontSize',13);
%    plot(quad.Pos,quad.demodX-quad.demodY,'m');grid
%    legend('DemodX-DemodY')
%    xlabel('Radial Pos from center');ylabel('DemodX-DemodY unitless');
%    title({[line1],[line2],[line3],[line4]});
   
end





