
% Read engineering Data produced
% by the PRIMA Metrology Phase-Meter LCU's (lprmac and lprma2) 
% Created: April 06, SLE 
%          23/8/06: SLE, new header definition, 
%                   add print out of mean DeltaL    
%          30/8/06: SLE, modify computation relative Std and plots
%          31/8/06: SLE, avoid pb of rmspsd failing depending on nfft
%          4/12/06: SLE add warning for NaN in case of block lost
%          21/3/07: SLE update to read 3 new columns in data file (raw data
%          from the phase meter
%          11/12/08: SLE, add  ('Interpreter','none') in the figures 
%                         Organise all figures on the screen and update
%                         units.  Status analysis done at the end 
%          19/11/10: SLE, improve computation time and display of the
%                         status errors. Loop through bit numbers (1-10) instead
%                         of through the whole data length.
%                         Status analysis done at the begining
%          30/11/10  SLE  Plot Acceleration; 

% Structure of the file:
% Date/Start time: 2006-04-06T08.56.11
% acquisition/loop rate (Hz): 100.00
% Configuration:
% CHA: Instrument
%   lp3 (+38.65 MHz)
%   lp1 (+38.00 MHz)
% CHB: FSU B
%   lp4 (-39.55 MHz)
%   lp2 (-40 MHz)
% Laser Frequency (nm): 1319.00
% Frequency Shift (MHz): 78.00
% Index of refraction: n=1
% Number of Samples: 3950
% Definition of the data columns:
% Relativetime,status,DeltaL_PM,Epsilon_FC,DeltaL_REF, Phi_sum(digit), Phi_FC (digit),Number of samples
%
% 0.003087,0x013F,-0.016268763,-7092406.000000000,0.000000000, 0.000000000,0.000000000,25
% 0.033000, NaN,NaN,NaN,NaN   

% NaN-> if block lost occured ; Attention 
%                             In that case fscanf will stop !!!
%                             only the data before the NaN will be loaded
%                             feature of matlab on MS-windows that can be
%                             overtaken only by looping on line... (too
%                             slow. If this happens, a warning occurs
%
%----------------------
% Get the right file
%----------------------


clear all;close all;pack

[FileName, PathName] = uigetfile('lprma*PhaseMeter*', 'Click on the data File to process'); 
eval([ 'cd ''',PathName,'''']);
fid = fopen(FileName,'r');


%-----------------------
% Read Header
%----------------------
    %first line
[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 test of line 2
%[fs,count]=fscanf(fid,'acquisition/loop rate (Hz): %f',[1,1]);
[fs,count]=fscanf(fid,'%% Sampling frequency in Hz: %f',[1,1]);
%Sampling frequency in Hz:
    % Jump line 3 to line 10
      for i=3:10;
          fgets(fid);
      end
[lambda,count]=fscanf(fid,'%% Laser Frequency (nm): %f',[1,1]);
lambda=lambda*1.e-9;% wavelength in meter

fgets(fid);%line11
[delta_nu,count]=fscanf(fid,'%% Frequency Shift (MHz): %f',[1,1]);
delta_nu=delta_nu*1.e6;% frequency shift in Hz
fgets(fid);%line12
[n,count]=fscanf(fid,'%% Index of refraction: n=%f',[1,1]);
fgets(fid);%line13
[N_samples,count]=fscanf(fid,'%% Number of Samples: %d',[1,1]);
fgets(fid);%line14->description of columns
fgets(fid);%line 16 (empty)
fgets(fid);%line 17  start of data

fprintf('\n')
fprintf(' **************** \n')
fprintf(' Phase Meter Data \n')
fprintf(' ==> Filename: %s \n    Sampling frequency= %.0f Hz \n    Number of samples= %.0f \n',FileName,fs,N_samples);


%---------------------------
% Read data vs time 
%---------------------------
[data,count]=fscanf(fid,'%f,%i,%f,%f,%f,%f,%f,%f',[8,inf]);
fclose(fid);

seconds = data(1,:)'; % in sec 
status_dec = data(2,:)'; % status delivered by the phase meter

                    %Raw data delivered by the phase meter
Phi_sum=data(6,:)'; % summed phase delivered by the phasemeter expressed in digits (1 digit = 2?/1024 rad)
Phi_FC=data(7,:)';  % error compensation fringe counter delivered by the phasemeter expressed in digit (1 digit =2? rad)
NSamples_PM=data(8,:)'; % numbser of samples delivered by the phase meter

                       % Differential OPD Pre-computed by the LCU 
DeltaL_PM= data(3,:)'; % Differential OPD computed by the LCU in m -> should wrapp after a "PM overflow"
epsilon_FC = data(4,:)';% Error correction computed by LCU in m --> should NOT wrapp after a "FC overflow"
DeltaL_REF= data(5,:)'; % Reference OPD value as set on the pmacqGui in m

DeltaL=DeltaL_PM+epsilon_FC-DeltaL_REF; % corrected Differential OPD in m


c=299792458;%speed of light in m/s
nu=n*c/lambda;

                       % Verification that DeltaL_PM and epsilon_FC
                       % are computed according to the formulas specified
                       % as below
                       %
DeltaL_PMRaw=c/2/nu*Phi_sum/1024./(NSamples_PM+1); % in meters
epsilon_FC_Raw=c/2/nu*delta_nu/nu*Phi_FC; % in meters
DeltaL_Raw=DeltaL_PMRaw+epsilon_FC_Raw-DeltaL_REF; % in meters

fprintf('    Sampling= %e Hz  P_V Sampling= %e sec \n',1/mean(diff(seconds)),max(diff(seconds))-min(diff(seconds)));
fprintf('    DeltaL Mean= %e m\n    DeltaL PV= %e m\n    DeltaL rms= %e m\n',mean(DeltaL),max(DeltaL)-min(DeltaL),std(DeltaL));


% check for block lost: NaN in data data or bad deltaL ref
% and generate warning window
if length(DeltaL_REF)~=N_samples
    fprintf(' WARNING for NaN: Check for block lost or bad DeltaL_ref \n')
    warndlg('WARNING for NaN: Check for block lost or bad DeltaL_ref','!! Warning !!')
end


%----------------------------
%prepare display of figures:
%----------------------------
 Nfig=8;
 scrsz = get(0,'ScreenSize'); %[1 1 width height]
 SizeFig=[scrsz(3)/(Nfig/2) scrsz(4)/2.5];%[width height]


%-------------------------
% Check for Status Error
%-------------------------
    %Status describtion         Default value when OK 
    %450kHz Ref detected        1 (LSB)
    %650kHz Ref detected        1
    %450k Probe detected        1
    %650k Probe detected        1
    %PLL Locked                 1
    %200kHz Probe detected      1
    %PM Overflow                0
    %FC Overflow                0
    %Trigger detected           1
    %Reset detected             0 (MSB)

    status_OK_dec=bin2dec('0100111111'); % Status OK corresponds to dec number 319 or hex 13F
    index_Pb=find(status_dec~=status_OK_dec);% find indexes where status has a problem
    Nb_Pb=length(index_Pb);
    fprintf('    Status Not OK:  %.0f \n',Nb_Pb);

    

%  Analyse in details the origin of the status error
%-------------------------------------------------------
    
    if (Nb_Pb~=0)
        
      fprintf(' Analyse status...\n',Nb_Pb);
      Pb=[];
          %generate "Pb" matrix [length(status_dec) lines x 10 columns]: each line = status
          %in binary;  First  column=LSB, last column=MSB (reverse writing
          % order)

     for i=1:1:10   %scan from LSB(bit=1) to MSB(bit=10)
           Pb(:,i)=bitget(uint16(status_dec),repmat(i,[length(status_dec) 1])); 
                     % repmap is a vector of length equal to status_dec
                     %containing only the bit number to be checked
     end
      
 
    % Check the nature of the status problem
    i_450kHz_Ref_not_detected=find(Pb(:,1)~=1); % LSB
    i_650kHz_Ref_not_detected=find(Pb(:,2)~=1);
    i_450kHz_Probe_not_detected=find(Pb(:,3)~=1);
    i_650kHz_Probe_not_detected=find(Pb(:,4)~=1);
    i_PLL_not_locked=find(Pb(:,5)~=1);
    i_200kHz_not_detected=find(Pb(:,6)~=1);
    i_PM_overflow=find(Pb(:,7)~=0);
    i_FC_overflow=find(Pb(:,8)~=0);
    i_trigger_not_detected=find(Pb(:,9)~=1);
    i_reset_detected=find(Pb(:,10)~=0);        % MSB

    % Display Cumulative number of errors 
    h=figure(1);set(h,'Position',[1 scrsz(4)/2 SizeFig]); %[left bottom width height]
    subplot(2,1,1)
    set(gca,'FontSize',13);
    bar([10:-1:1],[length(i_reset_detected),length(i_trigger_not_detected),length(i_FC_overflow),length(i_PM_overflow),length(i_200kHz_not_detected),length(i_PLL_not_locked),length(i_650kHz_Probe_not_detected),length(i_450kHz_Probe_not_detected),length(i_650kHz_Ref_not_detected),length(i_450kHz_Ref_not_detected)],'r');
    title({[FileName],['Number of samples= ',num2str(N_samples)],['Bad Status = ' num2str(Nb_Pb)]})
    xlabel('Bit number: MSB=10 LSB=1');
    ylabel('Cumulative number of errors');
    text(10,0,'Reset detected','Rotation',90,'Fontsize',16);
    text(9,0,'trigger not detected','Rotation',90,'Fontsize',16);
    text(8,0,'FC overflow','Rotation',90,'Fontsize',16);
    text(7,0,'PM overflow','Rotation',90,'Fontsize',16);
    text(6,0,'200kHz not detected','Rotation',90,'Fontsize',16);
    text(5,0,'PLL not Locked','Rotation',90,'Fontsize',16);
    text(4,0,'650kHz Probe not detected','Rotation',90,'Fontsize',16);
    text(3,0,'450kHz Probe not detected','Rotation',90,'Fontsize',16);
    text(2,0,'650kHz Ref not detected','Rotation',90,'Fontsize',16);
    text(1,0,'450kHz Ref not detected','Rotation',90,'Fontsize',16);
    
    %Display time of occurance of each errors
    subplot(2,1,2)
    set(gca,'FontSize',13);
    plot(seconds(i_reset_detected),10*ones(1,length(i_reset_detected)),'ro')
    hold on;
    plot(seconds(i_trigger_not_detected),9*ones(1,length(i_trigger_not_detected)),'rx')
    plot(seconds(i_FC_overflow),8*ones(1,length(i_FC_overflow)),'r+')
    plot(seconds(i_PM_overflow),7*ones(1,length(i_PM_overflow)),'r*')
    plot(seconds(i_200kHz_not_detected),6*ones(1,length(i_200kHz_not_detected)),'rs')
    plot(seconds(i_PLL_not_locked),5*ones(1,length(i_PLL_not_locked)),'rd')
    plot(seconds(i_650kHz_Probe_not_detected),4*ones(1,length(i_650kHz_Probe_not_detected)),'rv')
    plot(seconds(i_450kHz_Probe_not_detected),3*ones(1,length(i_450kHz_Probe_not_detected)),'r>')
    plot(seconds(i_650kHz_Ref_not_detected),2*ones(1,length(i_650kHz_Ref_not_detected)),'rp')
    plot(seconds(i_450kHz_Ref_not_detected),1*ones(1,length(i_450kHz_Ref_not_detected)),'rh')
    hold off
    ylim([ 1 10])
    text(seconds(index_Pb(1)),10.2,'Reset detected','Fontsize',11);
    text(seconds(index_Pb(1)),9.2,'trigger not detected','Fontsize',11);
    text(seconds(index_Pb(1)),8.2,'FC overflow','Fontsize',11);
    text(seconds(index_Pb(1)),7.2,'PM overflow','Fontsize',11);
    text(seconds(index_Pb(1)),6.2,'200kHz not detected','Fontsize',11);
    text(seconds(index_Pb(1)),5.2,'PLL not Locked','Fontsize',11);
    text(seconds(index_Pb(1)),4.2,'650kHz Probe not detected','Fontsize',11);
    text(seconds(index_Pb(1)),3.2,'450kHz Probe not detected','Fontsize',11);
    text(seconds(index_Pb(1)),2.2,'650kHz Ref not detected','Fontsize',11);
    text(seconds(index_Pb(1)),1.2,'450kHz Ref not detected','Fontsize',11);
    xlabel('Time in sec')
    ylabel('Bit number MSB=10 LSB=1')
    grid
   % title({[FileName],['Number of samples= ',num2str(N_samples)],['Bad Status = ' num2str(Nb_Pb)]})
    
    end
 
    
    
    
%-------------------------------        
% Plot Recorded Data
%-------------------------------    
    
    
    
    h=figure(3);set(h,'Position',[2*SizeFig(1) scrsz(4)/2 SizeFig]);
    subplot(3,1,1)
         set(gca,'FontSize',13);
         plot(seconds,DeltaL, 'b.-',seconds(index_Pb),DeltaL(index_Pb),'ro');grid
         xlabel('time in sec');ylabel('DeltaL in m');
        legend('DeltaL','Status not OK')
        title({[FileName],['Fs=',num2str(fs), ' Hz'],['Bad Status = ' num2str(Nb_Pb)],['PV= ' num2str(max(DeltaL)-min(DeltaL)) ' m'],['Std= ' num2str(std(DeltaL)) ' m']})
    subplot(3,1,2)
        set(gca,'FontSize',13);
        speed=diff(DeltaL)./diff(seconds);
        %********
       % plot(seconds(1:end-1),speed, 'b.-',seconds(index_Pb),speed(index_Pb),'ro');grid
        plot(seconds(1:end-1),speed, 'b.-');grid
         xlabel('time in sec');ylabel('speed DeltaL in m/s');
        legend('speed DeltaL','Status not OK')
     subplot(3,1,3)
        set(gca,'FontSize',13);
        acceleration=diff(speed)./diff(seconds(1:end-1));
        %plot(seconds(1:end-2),acceleration, 'b.-',seconds(index_Pb),acceleration(index_Pb),'ro');grid
         plot(seconds(1:end-2),acceleration, 'b.-');grid
         xlabel('time in sec');ylabel('Acceleration DeltaL in m/s2');
        legend('Acceleration DeltaL','Status not OK')
    
    %subplot(2,1,1)
% 
%     set(gca,'FontSize',13);
%     plot(seconds,DeltaL, 'b.',seconds,DeltaL_Raw,'go', seconds(index_Pb),DeltaL(index_Pb),'r.');grid
%     xlabel('time in sec');ylabel('DeltaL in m');
%     legend('DeltaL','DeltaL Raw','Status not OK')
%     title({['Phasemeter   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Bad Status = ' num2str(Nb_Pb)],['Std= ' num2str(std(DeltaL)) ' m']})
% 

% subplot(2,1,2)
%     set(gca,'FontSize',13);
%     plot(seconds,DeltaL-DeltaL_Raw, 'b.-');grid
%     xlabel('time in sec');ylabel('DeltaL in m');
%     legend('DeltaL-DeltaL Raw')
%     figName = strrep(FileName, '.txt','.jpg');
%     print('-djpeg90','-r300',figName);




h=figure(4);set(h,'Position',[3*SizeFig(1) scrsz(4)/2 SizeFig]);
subplot(2,1,1)
    set(gca,'FontSize',13);
    plot(seconds,DeltaL_PM,'b.',seconds,DeltaL_PMRaw,'g.');grid
    xlabel('time in sec');ylabel('OPD in m');
    h1=legend('DeltaL_PM','DeltaL_PMRaw');set(h1,'Interpreter','none')
    title({[FileName],['Fs=',num2str(fs), ' Hz'],['Bad Status = ' num2str(Nb_Pb)],['Std= ' num2str(std(DeltaL)) ' m']})
subplot(2,1,2)
    set(gca,'FontSize',13);
    plot(seconds,DeltaL_PM-DeltaL_PMRaw,'b');grid
    xlabel('time in sec');ylabel('OPD in m');
    h1=legend('DeltaLPM-DeltaLPMRaw');set(h1,'Interpreter','none')
%     figName = strrep(FileName, '.txt','.jpg');
%     print('-djpeg90','-r300',figName);

h=figure(5);set(h,'Position',[1 1 SizeFig]);
subplot(2,1,1)
    set(gca,'FontSize',13);
    plot(seconds,epsilon_FC,'b.',seconds,epsilon_FC_Raw,'g.');grid
    xlabel('time in sec');ylabel('OPD in m');
    h1=legend('epsilon_FC','epsilon_FC_Raw');set(h1,'Interpreter','none')
    title({[FileName],['Fs=',num2str(fs), ' Hz'],['Bad Status = ' num2str(Nb_Pb)],['Std= ' num2str(std(DeltaL)) ' m']})
subplot(2,1,2)
    set(gca,'FontSize',13);
       plot(seconds,epsilon_FC-epsilon_FC_Raw,'b.-');grid
    xlabel('time in sec');ylabel('OPD in m');
    h1=legend('epsilon_FC-epsilon_FC_Raw');set(h1,'Interpreter','none')


h=figure(6);set(h,'Position',[SizeFig(1) 1 SizeFig]);
    set(gca,'FontSize',13);
    subplot(2,1,1)
    plot(seconds(1:end-1),diff(seconds)*1.e6,'g.-');grid
    xlabel('time in sec');ylabel('Sampling time in microsec');
    legend('Stability of Sampling')
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['mean sampling time' num2str(mean(diff(seconds)*1000)) 'msec'],['PV sampling time' num2str(max(diff(seconds)*1000)-min(diff(seconds)*1000)) 'msec']})
    subplot(2,1,2)
    plot(seconds(1:end-1),diff(seconds)/5e-6,'kd-',seconds(1:end-1),floor(diff(seconds)/5e-6),'m.-',seconds,NSamples_PM,'r.-');grid
    xlabel('time in sec');ylabel('Number of samples');
    legend('Expected Nb of samples','Expected Integer Nb of samples','Nb of samples Provided by PM data ')
    
h=figure(7);set(h,'Position',[2*SizeFig(1) 1 SizeFig]);
subplot(2,1,1)
    set(gca,'FontSize',13);
    plot(seconds,Phi_sum,'b.');grid
    xlabel('time in sec');ylabel('Digits (1 digit=2pi/1024 rad=0.6nm)');
    h1=legend('Phi_sum in digit');set(h1,'Interpreter','none')
    title({[FileName],['Fs=',num2str(fs), ' Hz'],['Bad Status = ' num2str(Nb_Pb)],['RAW PHASE METER DATA']})
subplot(2,1,2)
    set(gca,'FontSize',13);
    plot(seconds,Phi_FC,'b.-');grid
    xlabel('time in sec');ylabel('Digits (1 digit=2pi rad=659nm)');
    h1=legend('Phi_FC in digit');set(h1,'Interpreter','none')


%   [C,F] = mscohere(NSamples_PM(1:end-1),(diff(seconds)/5e-6),[],[],[],fs);
%   figure(5); semilogx(F,C,'k.-')
%   ylabel('Coherence');
%   legend('Expected Nb of samples')

%     figName = strrep(FileName, '.txt','.jpg');



     
    
%-----------------------------------------
% Compute PSD and plot the results
%-----------------------------------------
%[PSD_DeltaL,fre,std_f,f_plot]=PSDcalc_sam(fs,DeltaL);
nfft=length(DeltaL);
noverlap = round(0.9*nfft);
[PSD_DeltaL,fre] = psdbk(DeltaL,fs,nfft,noverlap);

fprintf(' Compute PSDs... \n')

% STD calculation 
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,PSD_DeltaL,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,PSD_DeltaL,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

fprintf('    DeltaL rms above f_cut= %.0f Hz= %e m\n',f_cut,rms_above_f_cut');

h=figure(8);set(h,'Position',[3*SizeFig(1) 1 SizeFig]); set(gca,'FontSize',13);
set(gca,'FontSize',13);
    subplot(2,1,1);
    loglog(fre,sqrt(PSD_DeltaL),'r');grid
    xlabel('Frequency in Hz');ylabel('PSD DeltaL in m/sqrt(Hz)');
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['Std over (' num2str(fre(1)) '-' num2str(fre(end)) ') Hz= ' num2str(std(DeltaL)) 'm']})

    subplot(2,1,2)
    plot(f_plot,std_f_relativeToFcut*100,'m.-');grid
    xlabel('Frequency in Hz');
    ylabel('% to Std above fcut');
    ii=find(std_f_relativeToFcut<0.99);
    axis([0 f_plot(ii(end)) 0 100])
     %Title of separated figure: title({['Phasemeter   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cut) 'm']})
    title({['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cut) 'm']})
%      figName = strrep(['PSD' FileName], '.txt','.jpg');
%      print('-djpeg90','-r300',figName);   


    
    
% save results in a File
%  fid2 = fopen('./resultsPhaseMeter.txt','at+');
%  fprintf(fid2,'%33s %0.f %.0f %.0f %e %e \n',FileName,Nb_Pb,fs,N_samples,max(DeltaL)-min(DeltaL),std(DeltaL));
%  fclose(fid2);



