% Read engineering Data produced
% by the PRIMA Metrology Phase-Meter 
% Created: April 06, SLE
%          23/8/06: new header definition , SLE

% 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 columns: Relativetime,status,DC_REF_650kHz,DC_REF_450kHz,DC_PROBE_650kHz,DC_PROBE_450kHz,
%
% 0.003087,0x013F,0.18,0.94,0.18,0.39
% 0.013082,0x013F,0.13,0.92,0.13,0.35
% 0.023081,0x013F,0.29,1.05,0.29,0.48



% Get the right file
%----------------------
clear all
close all
pack

[FileName, PathName] = uigetfile('lprma*Photodiodes*', '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]);
    % 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
fgets(fid);%line15-> after modif from Robert to include description of columns

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
%----------------
[data,count]=fscanf(fid,'%f,%i,%f,%f,%f,%f',[6,inf]);
fclose(fid);

seconds = data(1,:)';
status_dec = data(2,:)';
REF650= data(3,:)';
REF450 = data(4,:)';
PROBE650= data(5,:)';
PROBE450=data(6,:)';

% check for block lost: NaN in data data or fs mismatch at 5% level
% and generate warning window
    if length(seconds)~=N_samples
       fprintf(' WARNING for NaN: Check  block lost  \n')
       warndlg('WARNING for NaN: Check  block lost ','!! Warning !!')
    end

    if (1/mean(diff(seconds))>=fs*(1+0.05) | 1/mean(diff(seconds))<=fs*(1-0.05))
     fprintf(' WARNING for fs: fs written in header and computed from data set do not match  \n')
     warndlg('WARNING for fs: fs written in header and computed from data set do not match ','!! Warning !!')
     end


% Display Statistics in V
fprintf(' Statistics in Volts \n')
fprintf('    REF650 CH#A Mean= %1.3f V    PV= %1.3f V    rms= %1.3f V\n',mean(REF650),max(REF650)-min(REF650),std(REF650));
fprintf('    REF450 FSUB Mean= %1.3f V    PV= %1.3f V    rms= %1.3f V\n',mean(REF450),max(REF450)-min(REF450),std(REF450));
fprintf('    PROBE650 CH#A Mean= %1.3f V  PV= %1.3f V    rms= %1.3f V\n',mean(PROBE650),max(PROBE650)-min(PROBE650),std(PROBE650));
fprintf('    PROBE450 FSUB Mean= %1.3f V  PV= %1.3f V    rms= %1.3f V\n',mean(PROBE450),max(PROBE450)-min(PROBE450),std(PROBE450));
fprintf('\n')

% Display Statistics in W
Volt_per_microwatt=1.2;
fprintf(' Statistics in microWatts \n')
fprintf('    REF650 CH#A Mean= %3.3f micW    PV= %3.3f micW    rms= %3.3f micW\n',mean(REF650)/Volt_per_microwatt,(max(REF650)-min(REF650))/Volt_per_microwatt,std(REF650)/Volt_per_microwatt);
fprintf('    REF450 FSUB Mean= %3.3f micW    PV= %3.3f micW    rms= %3.3f micW\n',mean(REF450)/Volt_per_microwatt,(max(REF450)-min(REF450))/Volt_per_microwatt,std(REF450)/Volt_per_microwatt);
fprintf('    PROBE650 CH#A Mean= %3.3f micW  PV= %3.3f micW    rms= %3.3f micW\n',mean(PROBE650)/Volt_per_microwatt,(max(PROBE650)-min(PROBE650))/Volt_per_microwatt,std(PROBE650)/Volt_per_microwatt);
fprintf('    PROBE450 FSUB Mean= %3.3f micW  PV= %3.3f micW    rms= %3.3f micW\n',mean(PROBE450)/Volt_per_microwatt,(max(PROBE450)-min(PROBE450))/Volt_per_microwatt,std(PROBE450)/Volt_per_microwatt);
fprintf('\n')


% % % Compute PSD's and cumulative std
% fprintf('compute PSDs...\n')
% nfft=length(seconds);
% noverlap = round(0.9*nfft);
% f_cut=5; %frequency above which cumulative std is computed and plotted
% 
% [PSD_REF650,freREF650] = psdbk(REF650,fs,nfft,noverlap);
% [std_f_relativeToFcutREF650,rms_above_f_cutREF650,f_plotREF650]=cumStdPlot(PSD_REF650,freREF650,f_cut);
% 
% [PSD_REF450,freREF450] = psdbk(REF450,fs,nfft,noverlap);
% [std_f_relativeToFcutREF450,rms_above_f_cutREF450,f_plotREF450]=cumStdPlot(PSD_REF450,freREF450,f_cut);
% 
% [PSD_PROBE650,frePROBE650] = psdbk(PROBE650,fs,nfft,noverlap);
% [std_f_relativeToFcutPROBE650,rms_above_f_cutPROBE650,f_plotPROBE650]=cumStdPlot(PSD_PROBE650,frePROBE650,f_cut);
% 
% [PSD_PROBE450,frePROBE450] = psdbk(PROBE450,fs,nfft,noverlap);
% [std_f_relativeToFcutPROBE450,rms_above_f_cutPROBE450,f_plotPROBE450]=cumStdPlot(PSD_PROBE450,frePROBE450,f_cut);




% Analysis of the Status
%-------------------------
    %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);
    
    if (Nb_Pb~=0)
        fprintf('    analyse status... \n');
      %generate "Pb" matrix (Nb_Pb lines  x 10 columns): each line = status in binary
      Pb=[];
      for i=1:Nb_Pb % loop only on indexes with problem detected
          Pb(i,:)=bitget(status_dec(index_Pb(i)),10:-1:1);
      end 
   
 
    % Check the nature of the status problem
    i_reset_detected=find(Pb(:,1)~=0);
    i_trigger_not_detected=find(Pb(:,2)~=1);
    i_FC_overflow=find(Pb(:,3)~=0);
    i_PM_overflow=find(Pb(:,4)~=0);
    i_200kHz_not_detected=find(Pb(:,5)~=1);
    i_PLL_not_locked=find(Pb(:,6)~=1);
    i_650kHz_Probe_not_detected=find(Pb(:,7)~=1);
    i_450kHz_Probe_not_detected=find(Pb(:,8)~=1);
    i_650kHz_Ref_not_detected=find(Pb(:,9)~=1);
    i_450kHz_Ref_not_detected=find(Pb(:,10)~=1);
    
    figure(1)
    bar([1:10],[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({['Photodiode   ', date_time_txt],['Number of samples= ',num2str(N_samples)],['Bad Status = ' num2str(Nb_Pb)]})
    xlabel('Bit number MSB=1')
    text(1,0,'Reset detected','Rotation',90,'Fontsize',16);
    text(2,0,'trigger not detected','Rotation',90,'Fontsize',16);
    text(3,0,'FC overflow','Rotation',90,'Fontsize',16);
    text(4,0,'PM overflow','Rotation',90,'Fontsize',16);
    text(5,0,'200kHz not detected','Rotation',90,'Fontsize',16);
    text(6,0,'PLL not Locked','Rotation',90,'Fontsize',16);
    text(7,0,'650kHz Probe not detected','Rotation',90,'Fontsize',16);
    text(8,0,'450kHz Probe not detected','Rotation',90,'Fontsize',16);
    text(9,0,'650kHz Ref not detected','Rotation',90,'Fontsize',16);
    text(10,0,'450kHz Ref not detected','Rotation',90,'Fontsize',16);
     
    figure(2)
    plot(seconds(index_Pb(i_reset_detected)),1*ones(1,length(i_reset_detected)),'ro')
    hold on;
    plot(seconds(index_Pb(i_trigger_not_detected)),2*ones(1,length(i_trigger_not_detected)),'rx')
    plot(seconds(index_Pb(i_FC_overflow)),3*ones(1,length(i_FC_overflow)),'r+')
    plot(seconds(index_Pb(i_PM_overflow)),4*ones(1,length(i_PM_overflow)),'r*')
    plot(seconds(index_Pb(i_200kHz_not_detected)),5*ones(1,length(i_200kHz_not_detected)),'rs')
    plot(seconds(index_Pb(i_PLL_not_locked)),6*ones(1,length(i_PLL_not_locked)),'rd')
    plot(seconds(index_Pb(i_650kHz_Probe_not_detected)),7*ones(1,length(i_650kHz_Probe_not_detected)),'rv')
    plot(seconds(index_Pb(i_450kHz_Probe_not_detected)),8*ones(1,length(i_450kHz_Probe_not_detected)),'r>')
    plot(seconds(index_Pb(i_650kHz_Ref_not_detected)),9*ones(1,length(i_650kHz_Ref_not_detected)),'rp')
    plot(seconds(index_Pb(i_450kHz_Ref_not_detected)),10*ones(1,length(i_450kHz_Ref_not_detected)),'rh')
    hold off
    ylim([ 1 10])
    text(seconds(index_Pb(1)),1.2,'Reset detected','Fontsize',11);
    text(seconds(index_Pb(1)),2.2,'trigger not detected','Fontsize',11);
    text(seconds(index_Pb(1)),3.2,'FC overflow','Fontsize',11);
    text(seconds(index_Pb(1)),4.2,'PM overflow','Fontsize',11);
    text(seconds(index_Pb(1)),5.2,'200kHz not detected','Fontsize',11);
    text(seconds(index_Pb(1)),6.2,'PLL not Locked','Fontsize',11);
    text(seconds(index_Pb(1)),7.2,'650kHz Probe not detected','Fontsize',11);
    text(seconds(index_Pb(1)),8.2,'450kHz Probe not detected','Fontsize',11);
    text(seconds(index_Pb(1)),9.2,'650kHz Ref not detected','Fontsize',11);
    text(seconds(index_Pb(1)),10.2,'450kHz Ref not detected','Fontsize',11);
    xlabel('Time in sec')
    ylabel('Bit number MSB=1')
    grid
    title({[FileName],['Number of samples= ',num2str(N_samples)],['Bad Status = ' num2str(Nb_Pb)]})
   
    
    
     end


% Plot Data
%-------------------
figure(3);
plot(seconds,REF650,'g',seconds,REF450,'c',seconds,PROBE650,'r',seconds,PROBE450,'m');
%add levels correspondong to warning for saturation  (1.6 V) and warning for low flux (60mV)
hold on;plot(seconds,1.6*ones(1,length(seconds)),'r',seconds,0.06*ones(1,length(seconds)),'r');hold off;grid
xlabel('time in sec');ylabel('Voltage in V');
legend('REF 650','REF 450','PROBE 650','PROBE 450');
title({['Photodiode   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Status not OK= ' num2str(Nb_Pb)]})

figure(4)
    plot(seconds(1:end-1),diff(seconds)*1.e3,'k.');grid
    xlabel('time in sec');ylabel('Sampling time in msec');
    legend('Stability of Sampling')
    title({[FileName],['Fs=',num2str(fs), ' Hz']})

    
figure(5); set(gca,'FontSize',12);
    subplot(2,1,1);
    loglog(freREF650,sqrt(PSD_REF650),'g');grid
    xlabel('Frequency in Hz');ylabel('PSD DeltaL in m/sqrt(Hz)');
    legend('REF650 CHA')
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['Std over (' num2str(freREF650(1)) '-' num2str(freREF650(end)) ') Hz= ' num2str(std(REF650)) 'Volts']})

    subplot(2,1,2)
    plot(f_plotREF650,std_f_relativeToFcutREF650*100,'g.-');grid
    xlabel('Frequency in Hz');
    ylabel('% to Std above fcut');
    legend('REF650 CHA')
    ii=find(std_f_relativeToFcutREF650<0.99);
    axis([0 f_plotREF650(ii(end)) 0 100])
     %Title of separated figure: title({['Photodiode   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutREF650) 'Volts']})
    title({['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutREF650) 'Volts']})
    
    
figure(6); set(gca,'FontSize',12);
    subplot(2,1,1);
    loglog(freREF450,sqrt(PSD_REF450),'c');grid
    xlabel('Frequency in Hz');ylabel('PSD DeltaL in m/sqrt(Hz)');
    legend('REF450 FSUB')
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['Std over (' num2str(freREF450(1)) '-' num2str(freREF450(end)) ') Hz= ' num2str(std(REF450)) 'Volts']})

    subplot(2,1,2)
    plot(f_plotREF450,std_f_relativeToFcutREF450*100,'c.-');grid
    xlabel('Frequency in Hz');
    ylabel('% to Std above fcut');
    legend('REF450 FSUB')
    ii=find(std_f_relativeToFcutREF450<0.99);
    axis([0 f_plotREF450(ii(end)) 0 100])
     %Title of separated figure: title({['Photodiode   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutREF450) 'Volts']})
    title({['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutREF450) 'Volts']})
    
figure(7); set(gca,'FontSize',12);
    subplot(2,1,1);
    loglog(frePROBE650,sqrt(PSD_PROBE650),'r');grid
    xlabel('Frequency in Hz');ylabel('PSD DeltaL in m/sqrt(Hz)');
    legend('PROBE650 CHA')
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['Std over (' num2str(frePROBE650(1)) '-' num2str(frePROBE650(end)) ') Hz= ' num2str(std(PROBE650)) 'Volts']})

    subplot(2,1,2)
    plot(f_plotPROBE650,std_f_relativeToFcutPROBE650*100,'r.-');grid
    xlabel('Frequency in Hz');
    ylabel('% to Std above fcut');
    legend('PROBE650 CHA')
    ii=find(std_f_relativeToFcutPROBE650<0.99);
    axis([0 f_plotPROBE650(ii(end)) 0 100])
     %Title of separated figure: title({['Photodiode   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutPROBE650) 'Volts']})
    title({['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutPROBE650) 'Volts']})
    
figure(7); set(gca,'FontSize',12);
    subplot(2,1,1);
    loglog(frePROBE450,sqrt(PSD_PROBE450),'m');grid
    xlabel('Frequency in Hz');ylabel('PSD DeltaL in m/sqrt(Hz)');
    legend('PROBE450 FSUB')
    title({[FileName],['Fs=',num2str(fs), ' Hz; Bad Status = ' num2str(Nb_Pb)],['Std over (' num2str(frePROBE450(1)) '-' num2str(frePROBE450(end)) ') Hz= ' num2str(std(PROBE450)) 'Volts']})

    subplot(2,1,2)
    plot(f_plotPROBE450,std_f_relativeToFcutPROBE450*100,'m.-');grid
    xlabel('Frequency in Hz');
    ylabel('% to Std above fcut');
    legend('PROBE450 FSUB')
    ii=find(std_f_relativeToFcutPROBE450<0.99);
    axis([0 f_plotPROBE450(ii(end)) 0 100])
     %Title of separated figure: title({['Photodiode   ', date_time_txt],['Fs=',num2str(fs), ' Hz'],['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutPROBE450) 'Volts']})
    title({['Std above fcut=' num2str(f_cut) 'Hz : ' num2str(rms_above_f_cutPROBE450) 'Volts']})
    
    
