plothlts_spectra.m

%% Script overview - plothlts_spectra
% Function is to take spectra for a large triple suspension for each OSEM location & DOF (both with and without damping) 
%
%% Version History:-
% S. Aston,  Version 3,  28/03/2012, Initial release
% J. Kissel, Version 4,  01/04/2012, Added UTC time, Updated channel list for L1, Displayed progress timing/comments  
% S. Aston,  Version 5,  04/04/2012, Added option to save Matlab data file and tidy-up
% S. Aston,  Version 6,  06/04/2012, Updated channel list for L1 SUS PR2 DQ
% S. Aston,  Version 7,  09/04/2012, Added start 'now' option, bonus plot sanity checks, changed save data directory
% S. Aston,  Version 8,  27/04/2012, Tidy-up and removed bonus plots from script (added to plotallhlts_spectra.m)
% S. Aston,  Version 9,  29/05/2012, Added auto-generated channel list and generating plot filenames, plus a tidy-up
% S. Aston,  Version 10, 17/09/2012, Fixed site switching for X1 & X2, & plot fillpage
% S. Aston,  Version 11, 18/09/2012, Tidy-up and moved measurement GPS start times to plotallhlts_spectra.m
% S. Aston,  Version 12, 19/09/2012, Tidy-up configuration parameters and used tconvert to calculate UTC from GPS time
% S. Aston,  Version 13, 29/10/2012, Added FillPage and IDfig (including measurement author info)
% S. Aston,  Version 14, 03/12/2012, Remove individual plot's, added H2 case, fixed RCG2.6 lower data rate, new sensor calibration  
% S. Aston,  Version 15, 20/12/2012, Reverted Doff to Ref and Don to Res for compatibilty with previous results

%% Initialise and tidy-up
close all
clear all
clc

%% Configuration parameters
printFigs = 1; % Save PDF plots?
saveData = 1;  % Save data file (*.mat)?

ifo            = 'L1';          % Interferometer name ('X1','X2','L1','H1','H2')
opticID        = 'PR3';         % Optic name or "BSFM" if on a test stand ['BSFM','ETMX','ITMY','MC2']
susType        = 'HLTS';        % Type of suspension ('HLTS')
meas.sensCalib = true;          % true or false -- are sensors calibrated?
meas.author    = 'S.ASton';     % Let us know who did all the work!
automatic_gps_saving = 1;       % If true, let the script switching the damping on and off

svnDir = '/ligo/svncommon/SusSVN/sus/trunk/'; % LLO and everything except LHO triple test stand
% svnDir = '/ligo3/svncommon/SusSVN/sus/trunk/'; % ONLY LHO triple test stand

resultsDirName = 'Results';     % Results directory name
figFileTag = '2020-04-21_1300'; % Output filename, use to differentiate plots of different collections of data (date_time)

susID = [ifo 'SUS' opticID];
figDir = [svnDir susType '/' ifo '/' opticID '/'];

%% Measurement parameters

    numAvgs = 5;
    freqRes = 0.01;
    measDuration = 500;
    freqRange = [0.01 900];
    damping_on_gps=1085533588;         % GPS Start Time Damping ON
    damping_off_gps=1085534588;         % GPS Start Time Damping OFF

%% New Feature added Feb 13th 2013 
% If automatic_gps_saving is true, it turns on the damping, saves the gps time, wait for measDuration, turns
% off damping and wait for measDuration

if automatic_gps_saving == false
    startTimeRef.gps = damping_off_gps;  
    startTimeRes.gps = damping_on_gps;   
else
    [damping_on_gps damping_off_gps]=switch_damp_spectra(ifo,opticID,measDuration) %switches damping ON then OFF and saves the gps time associated
    startTimeRef.gps = damping_off_gps;  % GPS Start Time Damping OFF
    startTimeRes.gps = damping_on_gps ;  % GPS Start Time Damping ON
end

%% OSEM sensor calibration (ref https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=4756)
% aosemCalib = 3.2*10^-8; % AOSEM (ref https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=3526)
% bosemCalib = 4.1*10^-8; % BOSEM (ref https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=3526)

if meas.sensCalib
                  % <-Sensor-> ;
                  %    OSEM    ;
                  % Sensitivity; 
                  %   [m/um]   ; 
    calibration = ( (1/1e6) )  ;
else
                  % <----------------- Sensor Chain ---------------->; 
                  %     OSEM                     SatAmp              ;
                  %  Sensitivity                TransImp.    ADC     ; 
                  %   [mm/uA]    [m/mm]  [uA/A]   [A/V]   [V/sens ct]; 
    calibration = ( (0.7/76.29)*(1/1e3)*(1e6/1)*(1/240e3)*(40/2^16) );
end

susID = [ifo 'SUS' opticID];
figDir = [svnDir susType '/' ifo '/' opticID '/'];

%% GPS to UTC time conversion (tconvert)
tctimeRef = num2str(startTimeRef.gps);
tctimeRes = num2str(startTimeRes.gps);

tcstringRef = ['tconvert -f %Y-%m-%d" "%H:%M:%S" "UTC ' tctimeRef];
tcstringRes = ['tconvert -f %Y-%m-%d" "%H:%M:%S" "UTC ' tctimeRes];

[statusRef,utctimeRef] = system(tcstringRef);
[statusRes,utctimeRes] = system(tcstringRes);

startTimeRef.utc = strcat(utctimeRef);
startTimeRes.utc = strcat(utctimeRes);

%% Determine which test-stand or ifo and auto generate channel list 
if (ifo == 'X1') | (ifo == 'X2') 
    
    channelList = {[ifo,':SUS-HXTS_M1_OSEMINF_T1_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M1_OSEMINF_T2_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M1_OSEMINF_T3_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M1_OSEMINF_LF_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M1_OSEMINF_RT_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M1_OSEMINF_SD_OUT_DQ'];...
                   
                   [ifo,':SUS-HXTS_M1_DAMP_L_IN1_DQ'];...
                   [ifo,':SUS-HXTS_M1_DAMP_T_IN1_DQ'];...
                   [ifo,':SUS-HXTS_M1_DAMP_V_IN1_DQ'];...
                   [ifo,':SUS-HXTS_M1_DAMP_R_IN1_DQ'];...
                   [ifo,':SUS-HXTS_M1_DAMP_P_IN1_DQ'];...
                   [ifo,':SUS-HXTS_M1_DAMP_Y_IN1_DQ'];...
               
                   [ifo,':SUS-HXTS_M2_OSEMINF_LL_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M2_OSEMINF_LR_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M2_OSEMINF_UL_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M2_OSEMINF_UR_OUT_DQ'];...
               
                   [ifo,':SUS-HXTS_M2_WIT_L_DQ'];...
                   [ifo,':SUS-HXTS_M2_WIT_P_DQ'];...
                   [ifo,':SUS-HXTS_M2_WIT_Y_DQ'];...
               
                   [ifo,':SUS-HXTS_M3_OSEMINF_LL_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M3_OSEMINF_LR_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M3_OSEMINF_UL_OUT_DQ'];...
                   [ifo,':SUS-HXTS_M3_OSEMINF_UR_OUT_DQ'];...
                   
                   [ifo,':SUS-HXTS_M3_OSEMWIT_L_DQ'];...
                   [ifo,':SUS-HXTS_M3_OSEMWIT_P_DQ'];...
                   [ifo,':SUS-HXTS_M3_OSEMWIT_Y_DQ']};
               
end

if (ifo == 'L1') | (ifo == 'H1') | (ifo == 'H2')
    
    channelList = {[ifo,':SUS-',opticID,'_M1_OSEMINF_T1_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_OSEMINF_T2_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_OSEMINF_T3_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_OSEMINF_LF_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_OSEMINF_RT_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_OSEMINF_SD_OUT_DQ'];...
                   
                   [ifo,':SUS-',opticID,'_M1_DAMP_L_IN1_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_DAMP_T_IN1_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_DAMP_V_IN1_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_DAMP_R_IN1_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_DAMP_P_IN1_DQ'];...
                   [ifo,':SUS-',opticID,'_M1_DAMP_Y_IN1_DQ'];...
               
                   [ifo,':SUS-',opticID,'_M2_OSEMINF_LL_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M2_OSEMINF_LR_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M2_OSEMINF_UL_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M2_OSEMINF_UR_OUT_DQ'];...
               
                   [ifo,':SUS-',opticID,'_M2_WIT_L_DQ'];...
                   [ifo,':SUS-',opticID,'_M2_WIT_P_DQ'];...
                   [ifo,':SUS-',opticID,'_M2_WIT_Y_DQ'];...
               
                   [ifo,':SUS-',opticID,'_M3_OSEMINF_LL_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M3_OSEMINF_LR_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M3_OSEMINF_UL_OUT_DQ'];...
                   [ifo,':SUS-',opticID,'_M3_OSEMINF_UR_OUT_DQ'];...
                   
                   [ifo,':SUS-',opticID,'_M3_WIT_L_DQ'];...
                   [ifo,':SUS-',opticID,'_M3_WIT_P_DQ'];...
                   [ifo,':SUS-',opticID,'_M3_WIT_Y_DQ']};
               
end

%% Get Data from DAQ channels
disp(['Getting ' num2str(measDuration) ' secs of data'])
disp(['Getting Damping OFF data, ' num2str(measDuration) ' secs of data starting at ' num2str(startTimeRef.gps) ' ...'])
tic
rawDataRef = get_data(channelList,'raw',startTimeRef.gps,measDuration);
disp(['   Finished in ',num2str(toc),' seconds.'])

disp(['Getting Damping ON data, ' num2str(measDuration) ' secs of data starting at ' num2str(startTimeRes.gps) ' ...'])
tic
rawDataRes = get_data(channelList,'raw',startTimeRes.gps,measDuration);
disp(['   Finished in ',num2str(toc),' seconds.'])

%% Calculate additional measurement parameters
dt = 1/rawDataRef(1).rate;
time = dt:dt:measDuration;
time = time(:);

nyquistFreq = rawDataRef(1).rate / 2;
nFFTs = rawDataRef(1).rate / freqRes;

%% Analysis and plotting
disp('Beginning data analysis and plotting ...')

% Set Default Plot Parameters
set(0,'DefaultLineLineWidth',2) % Thick Lines
set(0,'DefaultAxesFontSize',16) % Nice big font
set(0,'DefaultAxesXGrid','on')  % Grid on by default
set(0,'DefaultAxesYGrid','on')
set(0,'DefaultAxesZGrid','on')

for iChan = 1:length(channelList) 
    % Detrend data (simply remove mean)
    detrendDataRef(iChan,:) = detrend(rawDataRef(iChan).data);
    detrendDataRes(iChan,:) = detrend(rawDataRes(iChan).data);
    
    % Calculate PSD via P.Welch method
    [psdDataRef(iChan,:), w] = pwelch(detrendDataRef(iChan,:),hann(nFFTs),[],nFFTs,rawDataRef(1).rate,'onesided');
    [psdDataRes(iChan,:), w] = pwelch(detrendDataRes(iChan,:),hann(nFFTs),[],nFFTs,rawDataRef(1).rate,'onesided');

    % Scale data
    scaleDataRef(iChan,:) = sqrt(psdDataRef(iChan,:));
    scaleDataRes(iChan,:) = sqrt(psdDataRes(iChan,:));

    % Calibrate Data for M1 (BOSEMs)
    if isempty(cell2mat(regexpi(channelList(iChan), '_M1_'))) == 0
        calibrateDataRef(iChan,:) = scaleDataRef(iChan,:) * calibration;
        calibrateDataRes(iChan,:) = scaleDataRes(iChan,:) * calibration;
    
        % Set data directory and filename
        sagID = 'M1';
        dataDir = [svnDir susType '/' ifo '/' opticID '/SAG' sagID '/' resultsDirName '/'];
        figFileDesc = char(channelList(iChan));
        fileName = [figFileTag '_' susID '_Spectra_' sagID '_0p1to900Hz_' figFileDesc];
    end
    
    % Calibrate Data for M2 (AOSEMs)
    if isempty(cell2mat(regexpi(channelList(iChan), '_M2_'))) == 0
        calibrateDataRef(iChan,:) = scaleDataRef(iChan,:) * calibration;
        calibrateDataRes(iChan,:) = scaleDataRes(iChan,:) * calibration;
    
        % Set data directory and filename
        sagID = 'M2';
        dataDir = [svnDir susType '/' ifo '/' opticID '/SAG' sagID '/' resultsDirName '/'];
        figFileDesc = char(channelList(iChan));
        fileName = [figFileTag '_' susID '_Spectra_' sagID '_0p1to900Hz_' figFileDesc];
    
    end
    
    % Calibrate Data for M3 (AOSEMs)
    if isempty(cell2mat(regexpi(channelList(iChan), '_M3_'))) == 0
        calibrateDataRef(iChan,:) = scaleDataRef(iChan,:) * calibration;
        calibrateDataRes(iChan,:) = scaleDataRes(iChan,:) * calibration;
    
        % Set data directory and filename
        sagID = 'M3';
        dataDir = [svnDir susType '/' ifo '/' opticID '/SAG' sagID '/' resultsDirName '/'];
        figFileDesc = char(channelList(iChan));
        fileName = [figFileTag '_' susID '_Spectra_' sagID '_0p1to900Hz_' figFileDesc];
    
    end
    
    % Plot power spectra
    figure(iChan)
    loglog(w, calibrateDataRef(iChan,:), w, calibrateDataRes(iChan,:),'-r')
    axis([1e-2,1e+3,1e-13,1e-5])
    set(gca,'XTick',10.^(-2:3),'YTick',10.^(-13:-5)) 
    set(gcf,'PaperOrientation','landscape')
    orient rotated;
    legend(['Damping OFF (' startTimeRef.utc ')'],['Damping ON (' startTimeRes.utc ')'],'Location','SouthWest') 
    xlabel('Frequency (Hz)')
    ylabel('Displacememt (m / Hz^1^/^2)','interpreter', 'tex')
    title ({[susID ' (' susType ') Amplitude Spectral Density'];figFileDesc},'interpreter', 'none')
        
    % Save pdf figures
    if printFigs
        figure(iChan)
        FillPage('w')
        IDfig(meas.author)
        
        if (ifo == 'X1')
            plotName(iChan,:) = {[figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID strrep(figFileDesc,'X1:SUS-HXTS','') '_Spectra.pdf']};
        end
        
        if (ifo == 'X2')
            plotName(iChan,:) = {[figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID strrep(figFileDesc,'X2:SUS-HXTS','') '_Spectra.pdf']};
        end
        
        if (ifo == 'L1') | (ifo == 'H1') | (ifo == 'H2')
            plotName(iChan,:) = {[figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID strrep(figFileDesc,[ifo ':SUS-',opticID],'') '_Spectra.pdf']};
        end
               
        saveas(gcf,plotName{iChan,:})
    end
end

disp('    Done plotting.')

disp('Merging .pdfs ... ')

%% PDF figure merge (M1)
mergeOrderM1 = [1 2 3 4 5 6 7 8 9 10 11 12]; % M1 spectra merge order (T1, T2, T3, LF, RT, SD, L, T, V, R, P, Y)
sagID = 'M1';

if printFigs
    mergeCommand = ['pdfmerge ' figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID '_' sagID '_ALL_Spectra.pdf'];
    for iFile = mergeOrderM1
        mergeCommand = [mergeCommand ' ' plotName{iFile}];
    end
        
    % SMA - 01/27/2017 - Added pdfmerge workaround for SL7 (https://services.ligo-la.caltech.edu/FRS/show_bug.cgi?id=7240)
    savedPath = getenv('LD_LIBRARY_PATH');
    setenv('LD_LIBRARY_PATH', '');
    system(mergeCommand);
    setenv('LD_LIBRARY_PATH', savedPath);
    
    for iFile = mergeOrderM1
        system(['rm ' plotName{iFile}]);
    end
    
end

%% PDF figure merge (M2)
mergeOrderM2 = [13 14 15 16 17 18 19]; % M2 spectra merge order (LL, LR, UL, UR, L, P, Y)
sagID = 'M2';

if printFigs
    mergeCommand = ['pdfmerge ' figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID '_' sagID '_ALL_Spectra.pdf'];
    for iFile = mergeOrderM2
        mergeCommand = [mergeCommand ' ' plotName{iFile}];
    end
        
    % SMA - 01/27/2017 - Added pdfmerge workaround for SL7 (https://services.ligo-la.caltech.edu/FRS/show_bug.cgi?id=7240)
    savedPath = getenv('LD_LIBRARY_PATH');
    setenv('LD_LIBRARY_PATH', '');
    system(mergeCommand);
    setenv('LD_LIBRARY_PATH', savedPath);
    
    for iFile = mergeOrderM2
        system(['rm ' plotName{iFile}]);
    end
        
end

%% PDF figure merge (M3)
mergeOrderM3 = [20 21 22 23 24 25 26]; % M3 spectra merge order (LL, LR, UL, UR, L, P, Y)
sagID = 'M3';

if printFigs
    mergeCommand = ['pdfmerge ' figDir 'SAG' sagID '/' resultsDirName '/' figFileTag '_' susID '_' sagID '_ALL_Spectra.pdf'];
    for iFile = mergeOrderM3
        mergeCommand = [mergeCommand ' ' plotName{iFile}];
    end
        
    % SMA - 01/27/2017 - Added pdfmerge workaround for SL7 (https://services.ligo-la.caltech.edu/FRS/show_bug.cgi?id=7240)
    savedPath = getenv('LD_LIBRARY_PATH');
    setenv('LD_LIBRARY_PATH', '');
    system(mergeCommand);
    setenv('LD_LIBRARY_PATH', savedPath);
    
    for iFile = mergeOrderM3
        system(['rm ' plotName{iFile}]);
    end
        
end

%% Save the Data
disp('Saving Data .mat... ')

if saveData
    saveFileName = [svnDir susType '/Common/Data/' figFileTag '_' susID '_M1M2M3_Spectra'];
    save([saveFileName '.mat'],'channelList','susID','susType','figFileDesc','calibrateDataRef','calibrateDataRes','w')
end

disp('   Totally done.')