function [alpha, beta, beta2, gamma] = process_single_charge_meas(filename,QUAD, drive_freq, plot_flag)
%PROCESS_SINGLE_CHAREG_MEAS import log file recording times of charge
%measurements, return values of coefficients from EQ 2 in G1600699

%% get the data

%import times of excitations
[gps_start, gps_stop, description] = importfile(filename);

addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/
gwd = GWData();
gwd.make_kerberos_ready;
gwd.site_info(2).server = 'nds2.ligo.caltech.edu';
gwd.site_info(2).port = 31200;
chans = {'H1:GDS-CALIB_STRAIN',strcat('H1:SUS-', QUAD , '_L3_MASTER_OUT_DC_DQ'),strcat('H1:SUS-', QUAD, '_L3_DRIVEALIGN_L_OUT_DQ')};
%sus_chans={strcat('H1:SUS-', QUAD , '_L3_MASTER_OUT_DC_DQ'),strcat('H1:SUS-', QUAD, '_L3_DRIVEALIGN_L_OUT_DQ')};
dt = 50;
for ii=1:length(gps_start)
    [data{ii}, t, info] = gwd.fetch(gps_start(ii)+12, dt, chans);%the +12 is
    %because there is delay between when the exitation starts and when I thought
    %it would start when I wrote the time logging.  Should be fixed in excitation script form now on.
end

%% process data

Fs=16384;  %sampling frequency
stepSize = 10; %1 second
overlap = 0.5;
NFFT = stepSize*Fs;
w = hanning(NFFT);

%calculate coherences and TFs
for ii=1:length(gps_start)
    this_data=data{ii};
    [bias_coh_all(:,ii),f] = mscohere(this_data(:,2), this_data(:,1),w,[],[],Fs);
    [length_coh_all(:,ii),f] = mscohere(this_data(:,3), this_data(:,1),w,[],[],Fs);
    [bias_tf_all(:,ii),f] = tfestimate(this_data(:,2), this_data(:,1),w,[],[],Fs);
    [length_tf_all(:,ii),f] = tfestimate(this_data(:,3), this_data(:,1),w,[],[],Fs);
end

%pick out the coherences and TF at the drive freq
drive_f_index=find(f>drive_freq,1);

bias_coh=bias_coh_all(drive_f_index,:);
length_coh=length_coh_all(drive_f_index,:);
bias_tf=bias_tf_all(drive_f_index,:);
length_tf=length_tf_all(drive_f_index,:);

%get these strain/V into meters/volt, correct for sign on x arm
if strcmp(QUAD, 'ETMX')
    sign_flip = -1;
elseif strcmp(QUAD, 'ITMX')
    sign_flip = -1;
else
    sign_flip = 1;
end

TF_bias_drive=bias_tf(1:2)*sign_flip*4000;
TF_signal_drive=length_tf(3:5)*sign_flip*4000;


%% calibrate 

%clibrate bias drive into volts applied to ESD
DAC_volts_per_cnts=20/2^18; %volts per count out of DAC
HV_gain=40; %volts out of HV driver per Volts out of DAC  %the bias sees this gain even if we are in LVLN mode
LV_gain=2;
%bias low pass from page 10 of D1500016
RC1=30e3*3e-6;
RC2=10e3*10e-9;
bias_filter=zpk([],[1/RC1,1/RC1, 1/RC2],(RC1*RC1*RC2)^-1);
bias_filt_frsp=evalfr(bias_filter, 1j*2*pi*drive_freq);  %get the frequency response of the filter at the excitation frequency
TF_bias_drive_m_per_V=TF_bias_drive./(DAC_volts_per_cnts*HV_gain*bias_filt_frsp);  %meters of DARM per volt applied to bias electrode

%calibrate signal drive into volts applied to ESD
EUL2ESD=0.25; %drive was measured at DRIVEALIGN_L_OUT_DQ
% %ESD was in LOW voltage mode with low pass on, DC gain of 2 (page 7 of
% %D1500016)
TF_signal_drive_m_per_V=TF_signal_drive./(EUL2ESD*DAC_volts_per_cnts*LV_gain);
signal_offset=-1*50000*EUL2ESD*DAC_volts_per_cnts*LV_gain;  %in votls

%use sus model to convert DARM meters to Newtons on test mass
addpath('/ligo/svncommon/SusSVN/sus/trunk/QUAD/Common/FilterDesign/MatFiles/')
addpath('/ligo/svncommon/SusSVN/sus/trunk/Common/MatlabTools/')
load('dampingfilters_QUAD_2014-11-10_LLO_model.mat','susModel', 'freq')
tstL2tstL_mpN = squeeze(susModel.dampedf(susModel.dampedout.tst.disp.L, susModel.dampedin.tst.drive.L,:));
%find excitation freq
m_per_Newton=tstL2tstL_mpN(find(freq>drive_freq,1));

signal_N_per_V=abs(TF_signal_drive_m_per_V./m_per_Newton).*sign(cos(angle(TF_signal_drive_m_per_V./m_per_Newton))); %am I screwing up signs? probably
bias_N_per_V=abs(TF_bias_drive_m_per_V./m_per_Newton).*sign(cos(angle(TF_bias_drive_m_per_V./m_per_Newton)));

%% find coefficents
 %find applied bias and signal offsets
 for ii=1:length(gps_start)
    this_data=data{ii};
    bias_DC(ii)= mean(this_data(:,2))*DAC_volts_per_cnts*HV_gain;
    signal_DC(ii)=mean(this_data(:,3))*DAC_volts_per_cnts*LV_gain;
 end
%beta+beta2= drive bias no DC bias or signal  beta-beta2= drive signal no
%DC bias or signal

beta=(bias_N_per_V(2)+signal_N_per_V(1))/2;
beta2=(bias_N_per_V(2)-signal_N_per_V(1))/2;
 
%alpha+gamma= 2*slope of bias_N_per_V with Vb
alpha_plus_gamma1=(bias_N_per_V(1)-bias_N_per_V(2))/(2*bias_DC(1))  %since I forgot to turn of the linearization this meas is not good
alpha_plus_gamma2=-1*(signal_N_per_V(3)-signal_N_per_V(1))/(2*bias_DC(5))

%gamma-alpha=
gamma_minus_alpha=(signal_N_per_V(2)-signal_N_per_V(1))/(2*signal_DC(4));

gamma=(gamma_minus_alpha+alpha_plus_gamma1)/2;
alpha=(alpha_plus_gamma1-gamma_minus_alpha)/2;
%% plots for sanity checks
if plot_flag 
    %make spectra

    Pxx=[];
    for ii=1:length(gps_start)
        this_data=data{ii};
        [Pstrain(:,ii),f]=pwelch(this_data(:,1)-mean(this_data(:,1)),w,[],[],Fs);
        [Pbias(:,ii),f]=pwelch(this_data(:,2)-mean(this_data(:,2)),w,[],[],Fs);
        [PL(:,ii),f]=pwelch(this_data(:,3)-mean(this_data(:,3)),w,[],[],Fs);
    end

    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','none');

    blah=data{1};
    figure(1)
    plot(t, blah(:,2)-mean(blah(:,2)))

    figure(2)
    semilogy(f, sqrt(Pstrain(:,1)),...
        f, sqrt(Pstrain(:,2)),...
        f, sqrt(Pstrain(:,3)),...
        f, sqrt(Pstrain(:,4)),...
        f, sqrt(Pstrain(:,5)))
    xlim([5,30])
    legend(description)
    title('Strain')

    figure(3)
    semilogy(f, sqrt(Pbias(:,1)),...
        f, sqrt(Pbias(:,2)),...
        f, sqrt(Pbias(:,3)),...
        f, sqrt(Pbias(:,4)),...
        f, sqrt(Pbias(:,5)))
    xlim([5,30])
    legend(description)
    title('BIAS drive')

    figure(4)
    semilogy(f, sqrt(PL(:,1)),...
        f, sqrt(PL(:,2)),...
        f, sqrt(PL(:,3)),...
        f, sqrt(PL(:,4)),...
        f, sqrt(PL(:,5)))
    xlim([5,30])
    legend(description)
    title('DRIVEALIGN OUT')

    figure(10)
    semilogy(f, bias_coh_all(:,1),...
        f, bias_coh_all(:,2),...
        f, bias_coh_all(:,3),...
        f, bias_coh_all(:,4),...
        f, bias_coh_all(:,5))
    xlim([5,30])
    ylim([0.01,1])
    grid on
    legend(description)
    title('Coherence of Strain with bias drive')

    figure(11)
    semilogy(f, length_coh_all(:,1),...
        f, length_coh_all(:,2),...
        f, length_coh_all(:,3),...
        f, length_coh_all(:,4),...
        f, length_coh_all(:,5))
    xlim([5,30])
    ylim([0.01,1])
    grid on
    legend(description)
    title('Coherence of Strain with length drive')
end


