
% Read engineering Data produced
% by the PRIMA Metrology laser stabilization LCU's (lprmls) 
% Created: 10 July 2011, SLE , after adding 2 columns in the engineering file:
           % pmChannelA: power meter channel A in V
           % pmChannelB: power meter channel B in V

            % 14/7/2011: the phase of the lock-in is directly
            % given in deg using hard coded  conversion 18deg/V


% Structure of the file:
% Date/Start time: 2011-07-06T22.27.16
% TAC filename: pmlss_IL.tac
% Loop status: CLOSED
% Wavelength: 1319.1788
% Laser diagnostic: Laser OK
% DC set. temp: 49.50
% DC act. temp: 49.60
% EOM set. temp: 30.00
% EOM act. temp: 33.40
% I2 set. temp: 70.00
% I2 act. temp: 70.00
% Lock-in ref. detected/frequency: 24999600.0000000000
% Lock-in sensitivity: 100 mVrms / -7 dBm
% Lock-in filter slope: 12 db/oct
% Lock-in output Ch1: 0.027465820
% Lock-in output Ch2: 0.007934570
% Number of Samples: 20000
% Definition of the columns:
% RelativeTime,errorSignalX,phaseSignalY,tempCmd,piezoCmd,pmChannelA,pmChan
% nelB
% 
% 0.000,0.036011,0.001831,1.128141,5.587613

%
%----------------------
% Get the right file
%----------------------


clear all;close all;pack

[FileName, PathName] = uigetfile('MetFreqStab*.txt', '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
fgets(fid);% line 3
[LoopStatus,count]=fscanf(fid,'Loop status: %s',[1,1]);
fgets(fid);%line 4
[wavelength,count]=fscanf(fid,'Wavelength: %f',[1,1]);
fgets(fid);% line 5
%[LaserDiagnostic,count]=fscanf(fid,'Laser diagnostic: %c',[1,10]);  % Laser diagnostic: Laser OK
[LaserDiagnostic,count]=textscan(fid,'Laser diagnostic: %s %s')  % Laser diagnostic: Laser OK
LaserDiagnostic=[ cell2mat(LaserDiagnostic{1}) cell2mat(LaserDiagnostic{2})];

    % Jump line 6 to line 20
      for i=6:20;
          fgets(fid);
      end


fprintf('\n')
fprintf(' **************** \n')
fprintf(' Laser Frequency Stabilization Data \n')
fprintf(' ==> Filename: %s \n ',FileName);
fprintf(' ==> date_time: %s \n ',date_time_txt);
fprintf(' ==> Loop Status: %s \n ',LoopStatus);
fprintf(' ==> Wavelength in nm: %f \n ',wavelength);
fprintf(' ==> Laser Diagnostic: %s \n ',LaserDiagnostic);

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

seconds = data(1,:)'; % in seconds
errorX = data(2,:)'; % in Volts
phase=data(3,:); % in deg
tempCmd=data(4,:); % In V
piezoCmd=data(5,:); % In V
powerMeterA=data(6,:); % In V
powerMeterB=data(7,:); % In V

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sensitivity=25*1.e-9; % in V/Hz=25V/GHz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

errorHz=errorX/Sensitivity; %in Hz
ErrorX_PV=max(errorX)-min(errorX);

c=299792458;%speed of light in m/s
nu=c/(wavelength*1.e-9); % in Hz

error_rel=errorHz./nu;


%----------------------------
%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]    
%-------------------------------        
% Plot Recorded Data
%-------------------------------    
    h=figure(3);set(h,'Position',[2*SizeFig(1) scrsz(4)/2 SizeFig]);
    subplot(2,1,1)
         set(gca,'FontSize',13);
         plot(seconds,errorX, 'b.-');grid
         xlabel('time in sec');ylabel('error signal in V');
        legend('errorX')
        title({[FileName '  ' LoopStatus '   ' LaserDiagnostic],['Sensitivity= ' num2str(Sensitivity) ' V/Hz'],['ErrorX PV= ' num2str(ErrorX_PV) ' V   ' num2str(ErrorX_PV/Sensitivity/1.e3) ' kHz' ],['ErrorX Std= ' num2str(std(errorX)) ' V   ' num2str(std(errorX)/Sensitivity/1.e3) ' kHz' ]})       
   subplot(2,1,2)
        set(gca,'FontSize',13);
        plot(seconds,phase, 'b.-');grid
        xlabel('time in sec');ylabel('phase in deg');
        legend('phase')
        
%     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,tempCmd,'r.-');grid
    xlabel('time in sec');ylabel('Temp Command in V');
    title({[FileName '  ' LoopStatus '   ' LaserDiagnostic] ,['wavelength= ' num2str(wavelength) ' nm ']})

subplot(2,1,2)
    set(gca,'FontSize',13);
    plot(seconds,piezoCmd,'r.-');grid
    xlabel('time in sec');ylabel('Piezo Command in V');

 h=figure(5);set(h,'Position',[1 1 SizeFig]);
    set(gca,'FontSize',13);
    plot(seconds,error_rel, 'g.-');grid
    xlabel('time in sec');ylabel('relative error');
    legend('relative error')
    title({[FileName '  ' LoopStatus '   ' LaserDiagnostic] ,['wavelength= ' num2str(wavelength) ' nm '],['Sensitivity= ' num2str(Sensitivity) ' V/Hz'],['Std= ' num2str(std(error_rel)) ' unitless'] })

h=figure(6);set(h,'Position',[SizeFig(1) 1 SizeFig]);
    set(gca,'FontSize',13);
    plot(seconds,powerMeterA, 'r.-');grid
    xlabel('time in sec');ylabel('Power Meter SHG');
    legend('SHG power')
    title({[FileName '  ' LoopStatus '   ' LaserDiagnostic] ,['Std= ' num2str(std(powerMeterA)) ' V'] })
     

h=figure(6);set(h,'Position',[SizeFig(1) 1 SizeFig]);
   subplot(2,1,1); set(gca,'FontSize',13);
    plot(seconds,powerMeterA, 'r.-');grid
    xlabel('time in sec');ylabel('Relative SHG power in V');
    legend('SHG power')
    title({[FileName '  ' LoopStatus '   ' LaserDiagnostic] ,['SHG:Std/mean= ' num2str(std(powerMeterA)/mean(powerMeterA)*100) ' %'] })
    subplot(2,1,2) 
    plot(seconds,powerMeterB, 'k.-');grid
    xlabel('time in sec');ylabel('Relative IR power in V');
    legend('IR power (if connected...)')


     
    
%-----------------------------------------
% Compute PSD and plot the results
%-----------------------------------------
fs=1/mean(diff(seconds));
nfft=length(seconds);
noverlap = round(0.9*nfft);
[PSD_S,fre] = psdbk(errorX,fs,nfft,noverlap);

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

% STD calculation 
f_cut=1; 
    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_S,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_S,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('    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);
    subplot(2,1,1);
    set(gca,'FontSize',13);
    loglog(fre,sqrt(PSD_S),'r');grid
    xlabel('Frequency in Hz');ylabel('PSD  in V/sqrt(Hz)');
   title({[FileName '  ' LoopStatus '   ' LaserDiagnostic]})
    subplot(2,1,2)
    set(gca,'FontSize',13);
    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) 'V']})
%      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);



