function [quad]=readQuadCorrection(varargin) 

% Reads the engineering files produced by pmpsdGUI
% MetQuadCorrestions*.txt
% 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
% The corrections sent the the STS are logged in the eng file
% 2 types of corrections: User offsets and Controller offsets

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

% Input parameters
    % filename:  example MetQuadCorrections2007-08-06T13.13.06.000000.txt
    % input_channel: 1,2 3 or 4
    % draw: 'yes' or 'no'
% Output:
    %quadCorr structure containing all info about the selected channel
    
% example:
%[quadCorr]=readQuadCorrection('filename','MetQuadCorrections2007-08-06T13.13.06.txt','input_channel',1,'draw','yes')
%[quadCorr]=readQuadCorrection('input_channel',1,'draw','yes') a window will
%                                                         pop up to
%                                                         select the file
% SLE August 2007

% July 20009:  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



%Default if no argument is given
filename='user_selected';
input_channel=4;
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('MetQuadCorrections*', '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]);
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]);
fgets(fid);
fgets(fid);
fgets(fid);


% new lines added July 2009
fgets(fid);
fgets(fid);
fgets(fid);
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,%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:7);
IP1.data=data(:,8:13);
IP4.data=data(:,14:19);
IP2.data=data(:,20:25);

% Select the input channel and analyse the data
eval(['quad=', ['IP' num2str(input_channel(1))],[';']]);
quad.seconds=data(:,1);

quad.UsrU=quad.data(:,1);
quad.UsrV=quad.data(:,2);
quad.UsrW=quad.data(:,3);
quad.CtrlU=quad.data(:,4);
quad.CtrlV=quad.data(:,5);
quad.CtrlW=quad.data(:,6);
radial_correction=sqrt((quad.UsrU+quad.CtrlU).^2+(quad.UsrW+quad.CtrlW).^2)

% test if corrections along V are really 0:
if (nnz(quad.UsrV)~=0 | nnz(quad.CtrlV)~=0)
    fprintf(' WARNING: Non-zero corrections along V !!  \n')
    warndlg('WARNING: Non-zero corrections along V !!','!! Warning !!')
end;


nfft=round(length(quad.seconds));
[quad.PSDCorrU,fre]=psdbk(quad.UsrU+quad.CtrlU,fc,nfft);
[quad.PSDCorrW,fre]=psdbk(quad.UsrW+quad.CtrlW,fc,nfft);
quad.fre=fre;

quad= rmfield(quad, 'data');


% Display the results
fprintf('    Mean UsrU = %e    PV UsrU= %e    rms UsrU= %e  \n',mean(quad.UsrU),max(quad.UsrU)-min(quad.UsrU),std(quad.UsrU));
fprintf('    Mean UsrW = %e    PV UsrW= %e    rms UsrW= %e  \n',mean(quad.UsrW),max(quad.UsrW)-min(quad.UsrW),std(quad.UsrW));
fprintf('\n')
fprintf('    Mean CtrlU = %e    PV CtrlU= %e    rms CtrlU= %e  \n',mean(quad.CtrlU),max(quad.CtrlU)-min(quad.CtrlU),std(quad.CtrlU));
fprintf('    Mean CtrlW = %e    PV CtrlW= %e    rms CtrlW= %e  \n',mean(quad.CtrlW),max(quad.CtrlW)-min(quad.CtrlW),std(quad.CtrlW));
fprintf('\n')
fprintf('    Mean U = %e    PV U= %e    rms U= %e  \n',mean(quad.UsrU+quad.CtrlU),max(quad.UsrU+quad.CtrlU)-min(quad.UsrU+quad.CtrlU),std(quad.UsrU+quad.CtrlU));
fprintf('    Mean W = %e    PV W= %e    rms W= %e  \n',mean(quad.UsrW+quad.CtrlW),max(quad.UsrW+quad.CtrlW)-min(quad.UsrW+quad.CtrlW),std(quad.UsrW+quad.CtrlW));
fprintf('    Mean Rcorr = %e    PV Rcorr= %e    rms Rcorr= %e  \n',mean(radial_correction),max(radial_correction)-min(radial_correction),std(radial_correction));

fprintf('\n')
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)];
   
   %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.UsrU,quad.UsrW,quad.CtrlU,quad.CtrlW,quad.UsrU+quad.CtrlU,quad.UsrW+quad.CtrlW]);grid
   plot(quad.seconds,quad.CtrlU,'b.',quad.seconds,quad.CtrlW,'bo',quad.seconds,quad.UsrU,'r+',quad.seconds,quad.UsrW,'rd');grid
   legend('CtrlU','CtrlW','UsrU','UsrW')
   xlabel('Time in sec');ylabel('Correction to STS (unitless)');
   title({[line1],[line2],[line3],[line4]});
 
   h=figure(2);
   set(h,'Position',[1*SizeFig(1) 1 SizeFig]);
   set(gca,'FontSize',13);
   plot(quad.UsrU,quad.UsrW,'r+');grid;axis('equal');hold on;plot(quad.UsrU(1),quad.UsrW(1),'kd');hold off;
   xlabel('UsrU');ylabel('UsrW');
   title({[line1],[line2],[line3],[line4]});
   
   h=figure(3);
   set(h,'Position',[SizeFig(1) scrsz(4)/2 SizeFig])
   set(gca,'FontSize',13);
   plot(quad.UsrU+quad.CtrlU,quad.UsrW+quad.CtrlW,'m+');grid;axis('equal')
   xlabel('UsrU+CtrlU');ylabel('UsrW+CtrlW');
   title({[line1],[line2],[line3],[line4]});
   
   h=figure(4);
   set(h,'Position',[2*SizeFig(1) scrsz(4)/2 SizeFig])
   set(gca,'FontSize',13);
   loglog(quad.fre,[sqrt(quad.PSDCorrU), sqrt(quad.PSDCorrW)]);grid
   legend('U','W')
   xlabel('freq in Hz');ylabel('PSD Corrections unitless/sqrt(Hz)');
   title({[line1],[line2],[line3],[line4]});
      
   h=figure(5);
   set(h,'Position',[3*SizeFig(1) scrsz(4)/2 SizeFig])
   set(gca,'FontSize',13);
   plot(quad.seconds(2:end),diff(quad.seconds));grid
   xlabel('Time in sec');ylabel('Sampling time in sec');
   title({[line1],[line2],[line3],[line4]});
   







% h=figure(6);
% set(gca,'FontSize',13);
% axis([-1 1 -1 1])
% axis('equal')
% mov = avifile('samuel5.avi')
% for k=1:100
% 	plot(quad.UsrU(k),quad.UsrW(k),'r+');
%     hold on;
% 	F = getframe(gca);
% 	mov = addframe(mov,F);
% end
% mov = close(mov);
% hold off;

end









