def check_stinger_position(t04, ref_llh, intervals, out_plot_path, time_offset):
    """
    Plot Stinger position 3D error, 3D sigma, number of SV vs. tow
    Input:    t04 = .T04 files
              ref_llh = a list of reference latitude[degree], longtitude[degree] and height[meter], eg. [37 + 23/60., -(122), 0.] for Auger test
              intervals = list of pairs of on and off time [[],[]]
              out_plot_path = output_dir/plots/x/ = directory where plots go
              time_offset = offset time in plot, unit [sec]
    """
    ref_pos = r_[ref_llh[0], ref_llh[1], ref_llh[2]] 

    # Calculate 3D error, 3D sigma
    #ticol = d[:, k['TIME']]
    d,k = vd2arr(t04, '-d35:16')
    t       = d[:, k['TIME']]
    fixtype = d[:, k['FIXTYPE']]
    svtrk   = d[:, k['NTRK']]
    svused  = d[:, k['NUSED']]
    gpsused = d[:, k['GPSUsd']]
    gloused = d[:, k['GLOUsd']]
    galused = d[:, k['GLOUsd']]

    err_3d = comp_3d_dist( d[:, k['LAT']:k['LAT']+3], ref_pos )
    sigma_3d = sqrt( d[:,k['SigU']]**2 + d[:,k['SigE']]**2 + d[:,k['SigN']]**2 )
    
    # Calculate mask for the period when jammer is on
    interval_included_by_union = numpy.zeros((err_3d.size,), dtype=bool)
    for j,interval in enumerate(intervals):
        interval_included_by_union |= (t>=interval[0])&(t<=interval[1])
    
    # Calculate 68%, 95%, 99% and 100% error when jammer on and off respectively
    jam_on_err = dotdict({})
    jam_off_err = dotdict({})
    idx_on  = find( interval_included_by_union )
    idx_off = find ( interval_included_by_union == False )
    cx3, cy3 = docdf( err_3d[idx_on] )
    jam_on_err.cdf68, jam_on_err.cdf95, jam_on_err.cdf99, jam_on_err.cdf100 = get_cdf_percents( cx3, cy3, [68,95,99,100] )
    cx3, cy3 = docdf( err_3d[idx_off] )
    jam_off_err.cdf68, jam_off_err.cdf95, jam_off_err.cdf99, jam_off_err.cdf100 = get_cdf_percents( cx3, cy3, [68,95,99,100] )

    data_mat = numpy.concatenate([[t], [err_3d], [sigma_3d], [fixtype]], axis=0)
    fname = out_plot_path+'/'+'stinger_pos_err.txt'
    hdr = "Time[sec] Error_3D[m] Sigma_3D[m] FIXTYPE[-]"
    numpy.savetxt(fname, data_mat.T, delimiter=' ', header=hdr)

    def plot_jammer_on_off(data, intervals, y_min=None, y_max=None):
        x = []
        y = []
        plot_y_min = 0
        plot_y_max = max(data) if y_max==None else y_max
        for pair in intervals:
            for z in [pair[0],pair[0],pair[1],pair[1]]:
                x.append(z)
            for z in [plot_y_min/1.05, plot_y_max*1.05, plot_y_max*1.05, plot_y_min/1.05]:
                y.append(z)
        matplotlib.pyplot.fill( [ xi-time_offset for xi in x], y, 'c', alpha=0.3, label='Jammer on')
        plt.legend(loc=1)
        
    # Plot Stinger Position metrics
    msg_3d_error = "Error during jamming:%.1fcm@68%%, %.1fcm@95%%, %.1fcm@99%%"%(100*jam_on_err.cdf68, 100*jam_on_err.cdf95, 100*jam_on_err.cdf99)

    plt.subplots(3,1, layout="constrained")
    ind = numpy.where(t[:]-time_offset > 120)
    # stinger 3D error and sigma
    plt.subplot(3,1,1)
    plt.plot(t-time_offset, 100*err_3d, '.', label='3D error')
    plt.plot(t-time_offset, 100*sigma_3d, '.', label='3D sigma')
    if max(err_3d[ind]) >= max(sigma_3d[ind]):
        y_max_plot = 1.1*max(100*err_3d[ind])
        plot_jammer_on_off(100*err_3d, intervals, y_min=0.0, y_max=y_max_plot)
    else:
        y_max_plot = 1.1*max(100*sigma_3d[ind])
        plot_jammer_on_off(100*sigma_3d, intervals, y_min=0.0, y_max=y_max_plot)   
    plt.ylim([0, y_max_plot])
    plt.title(msg_3d_error)
    plt.ylabel('error[cm]')

    # Number of tracking and used SVs
    plt.subplot(3,1,2)
    plt.plot(t-time_offset, svtrk, '.', label='#SV Tracked')
    plt.plot(t-time_offset, svused, '.', label='#SV Used')
    plt.plot(t-time_offset, gpsused, '.', label='GPS Used')
    plt.plot(t-time_offset, galused, '.', label='GAL Used')
    plt.plot(t-time_offset, gloused, '.', label='GLO Used')
    plot_jammer_on_off(svtrk, intervals)
    plt.ylabel('#SV')
    
    # Fix type
    plt.subplot(3,1,3)
    plt.plot(t-time_offset, fixtype, '.', label='FIXTYPE')
    plot_jammer_on_off(fixtype, intervals)
    plt.ylabel('FIXTYPE')
    plt.ylim([0, 10])

    plt.tight_layout()

    fname = out_plot_path+'/'+'stinger_pos_err.png'
    print('Save', fname)
    save_compressed_png(fname)
    
def cal_double_diff(d, k, rec, sy, freq, track, intervals=[[]]):
    """
    Estimate double differences for GPS L1-C/A and GLO G1C code and carrier-phse
    over the jamming periods.
    eg. dd_code, dd_carr, dd_code_statistics, dd_carr_statistics = cal_double_diff(d, k, '-d35:19', 0, 0, [0], intervals=[[100,200]])

    Input: 
        d = 2D array loaded from .T04 files
        k = key dictionary from .T04 files  
        rec = string '-d35:19'
        sy = satellite type
        freq = block type (frequency), eg. 0
        track = list of track/signal type, eg. [0] 
        intervals = eg. [[]] or [[t0, t1],[t2,t3], ...]
    Returns:
        dd_code, dd_carr, dd_code_statistics, dd_carr_statistics
        dd_code/dd_carr array [tow, ref_sv, target_sv, ref_sv_elevation, target_sv_elevation, double_difference]
        dd_statistics array [Sys, sv, epoches, min_dd, max_dd, mean_dd, std_dd, mad_dd]

    """
    ticol = d[:, k['TIME']]
    svcol = d[:, k['SV']]
    sycol = d[:, k['SAT_TYPE']]
    elcol = d[:, k['EL']] # elevation angle
    rtcol = d[:, k['FREQ']]
    tkcol = d[:, k['TRACK']]
    ancol = d[:, k['ANTENNA']] # antenna id
    mscol = d[:, k['MEAS']] # meas 
    sccol = d[:, k['SLIP_CNT']]
    prcol = d[:, k['RANGE']] # code phase
    cacol = d[:, k['PHASE']] # carrier phase  

    # For 35:19 only currently
    if rec != '-d35:19':
        raise RuntimeError('Not support or unknown record', rec)

    ref_sd_code = []
    ref_sd_carr = []
    dd_code_list = []
    dd_carr_list = []
    dd_code_stat_rows = []
    dd_carr_stat_rows = []

    systr, sv_range = get_sys_name_and_sv_range(sy)

    for ctr,interval in enumerate( intervals ):
        if interval == []:
            interval = [min(ticol), max(ticol)]

        time_interval_str = '_'+str(int(interval[0]))+'_'+str(int(interval[1]))
            
        # Calculate unioned track mask
        if track == []:
            track_included = numpy.ones((ticol.size,), dtype=bool)
        else:
            track_included = numpy.zeros((ticol.size,), dtype=bool)
            for j,track_j in enumerate(track):
                track_included |= (tkcol==track[j])
            
        # Calculate intervals mask
        interval_included = numpy.zeros((ticol.size,), dtype=bool)
        interval_included |= (ticol>=interval[0])&(ticol<=interval[1])
            
        # get all the time stamps for the interval
        ant0_idxes = numpy.where( (sycol==sy)&(ancol==0)&(rtcol==freq)&track_included&interval_included )
        ant1_idxes = numpy.where( (sycol==sy)&(ancol==1)&(rtcol==freq)&track_included&interval_included )
        t_array = []

        # make sure the elevation angle column is valid
        if max(elcol)==0:
            raise RuntimeError("SV elevation angles are unavailable")
            
        # get timestamps for antenna 0
        for t_idx in arange(len(ant0_idxes[0])):
            if t_idx == 0:
                t_array.append(ticol[ant0_idxes[0][t_idx]])
            else:
                if ticol[ant0_idxes[0][t_idx]] != t_array[-1]:
                    t_array.append(ticol[ant0_idxes[0][t_idx]])

        # find the SV having the highest elevation angle
        tmp_subidx = argmax(elcol[ant0_idxes])
        tmp_sv = svcol[ant0_idxes[0][tmp_subidx]]
        # need more than 1 epoch having the sv
        sv_epoches = numpy.where( (svcol[ant0_idxes]==tmp_sv) )

        if len(sv_epoches[0])<(0.5*len(t_array)):
            tmp_idxes = numpy.where( svcol[ant0_idxes]!=tmp_sv )
            tmp_subidx = argmax(elcol[ant0_idxes[0][tmp_idxes]])
            tmp_sv = svcol[ant0_idxes[0][tmp_idxes[0][tmp_subidx]]]
            sv_epoches = numpy.where( (svcol[ant0_idxes]==tmp_sv) )
            if len(sv_epoches[0])<(0.5*len(t_array)):
                raise RuntimeError("The highest elevation SV doesn't have enough epoches")

        ref_sv = tmp_sv
        dt = np.asarray(t_array[1:]) - np.asarray(t_array[:-1])
        data_rate = median(dt[numpy.where(dt[:]!=0)])

        # loop through the time stamps
        for ind, t in enumerate( t_array ):
            # get all the indexes having the same time, system, freq and track
            ant0_subid = numpy.where( (ticol[ant0_idxes]==t) )
            ant1_subid = numpy.where( (ticol[ant1_idxes]==t) )
            ant0_idx = ant0_idxes[0][ant0_subid[0]]
            ant1_idx = ant1_idxes[0][ant1_subid[0]]
            
            if len(ant0_idx)>1 and len(ant1_idx)>1:
                # find the highest elevation sv index
                ref0_subid = numpy.where( (svcol[ant0_idx]==ref_sv) )
                ref1_subid = numpy.where( (svcol[ant1_idx]==ref_sv) )
                if len(ref0_subid[0])==0 or len(ref1_subid[0])==0:
                    continue
                
                ref0_idx = ant0_idx[ref0_subid[0]]
                ref1_idx = ant1_idx[ref1_subid[0]]
                
                msflag = mscol[ref0_idx[0]].astype('int')
                flag1 = mask1 = k['dMEAS_RANGE_VALID']
                flag2 = mask2 = k['dMEAS_SLIP']
                flag3 = mask3 = k['dMEAS_PHASE_VALID']|k['dMEAS_LOCK_POINT']

                if msflag&flag1 != mask1 or msflag&flag2 == mask2 or msflag&flag3 != mask3:
                    continue
                if ind > 0:
                    ant0_subid_prev = numpy.where( (ticol[ant0_idxes]==t_array[ind-1])&(svcol[ant0_idxes]==ref_sv) )
                    ant1_subid_prev = numpy.where( (ticol[ant1_idxes]==t_array[ind-1])&(svcol[ant1_idxes]==ref_sv) )
                    ref0_idx_prev = ant0_idxes[0][ant0_subid_prev[0]]
                    ref1_idx_prev = ant1_idxes[0][ant1_subid_prev[0]]
                    if ((ticol[ref0_idx_prev[0]]-ticol[ref0_idx[0]])-data_rate) > 1e-6:
                        continue
                    if sccol[ref0_idx_prev[0]]!=sccol[ref0_idx[0]]:
                        continue
                    if ((ticol[ref1_idx_prev[0]]-ticol[ref1_idx[0]])-data_rate) > 1e-6:
                        continue
                    if sccol[ref1_idx_prev[0]]!=sccol[ref1_idx[0]]:
                        continue
                ref_elev = elcol[ref0_idx[0]]
                ref_sdcode = prcol[ref1_idx[0]] - prcol[ref0_idx[0]]
                ref_sdcarr = cacol[ref1_idx[0]] - cacol[ref0_idx[0]]

                for sv0_ind, sv0 in enumerate(svcol[ant0_idx]):
                    # skip the reference SV
                    if sv0 == ref_sv:
                        continue
                    for sv1_ind, sv1 in enumerate(svcol[ant1_idx]):
                        if sv0 == sv1:
                            target_sv = sv0
                            target_elev = elcol[ant0_idx][sv0_ind]
                            msflag = mscol[ant0_idx][sv0_ind].astype('int')

                            # code valid
                            if msflag&flag1 == mask1:
                                target_sdcode = prcol[ant1_idx][sv1_ind] - prcol[ant0_idx][sv0_ind]
                                # form double difference
                                dd_code = ref_sdcode - target_sdcode
                                dd_code_list.append([t, ref_sv, target_sv, ref_elev, target_elev, dd_code])
                                ref_sd_code.append(ref_sdcode)

                            # carrier phase valid
                            if (ind > 0):
                                ant0_subid_prev = numpy.where( (ticol[ant0_idxes]==t_array[ind-1])&(svcol[ant0_idxes]==target_sv) )
                                ant1_subid_prev = numpy.where( (ticol[ant1_idxes]==t_array[ind-1])&(svcol[ant1_idxes]==target_sv) )
                                tar0_idx_prev = ant0_idxes[0][ant0_subid_prev[0]]
                                tar1_idx_prev = ant1_idxes[0][ant1_subid_prev[0]]
                    
                                if len(ant0_subid_prev[0]) < 1 or len(ant1_subid_prev[0]) < 1:
                                    continue
                                
                            if ((msflag&flag2) != mask2)and((msflag&flag3) == mask3):
                                if (ind > 0):
                                    if sccol[tar0_idx_prev[0]]!=sccol[ant0_idx][sv0_ind]:
                                        continue
                                    if sccol[tar1_idx_prev[0]]!=sccol[ant1_idx][sv1_ind]:
                                        continue
                                target_sdcarr = cacol[ant1_idx][sv1_ind] - cacol[ant0_idx][sv0_ind]
                                dd_carr = ref_sdcarr - target_sdcarr
                                dd_carr_list.append([t, ref_sv, target_sv, ref_elev, target_elev, dd_carr])
                                ref_sd_carr.append(ref_sdcarr)

    ref_sd_code_array = numpy.asarray(ref_sd_code)
    ref_sd_carr_array = numpy.asarray(ref_sd_carr)
    ref_sd_code_mean = numpy.mean(ref_sd_code_array)
    ref_sd_code_std = numpy.std(ref_sd_code_array)
    ref_sd_carr_mean = numpy.mean(ref_sd_carr_array)
    ref_sd_carr_std = numpy.std(ref_sd_carr_array)
    
    output_dd_code_array = numpy.asarray(dd_code_list)
    output_dd_carr_array = numpy.asarray(dd_carr_list)

    # generate D.D. plots
    # time vs D.D.
    fig1 = plt.figure()
    ax1 = fig1.add_subplot(1, 1, 1)
    fig2 = plt.figure()
    ax2 = fig2.add_subplot(1, 1, 1)
    
    # table for code phase
    filename = 'dd_code_table.csv'
    fname = './'+filename
    f_code = open(fname, 'w')
    writer_code = csv.writer(f_code)
    # plot table for code phase
    fig3, ax3 = plt.subplots(1, 1, dpi=300)
    
    table_code_head_list = [
        'Satellite System','SV', 'Epochs',
        'Min [m]',
        'Max [m]',
        'Mean [m]',
        'Std [m]',
        'Mad [m]'
        ]
    writer_code.writerow(table_code_head_list)
    
    # table file for carrier phase
    filename = 'dd_carr_table.csv'
    fname = './'+filename
    f_carr = open(fname, 'w')
    writer_carr = csv.writer(f_carr)
    # plot table for carrier phase
    fig4, ax4 = plt.subplots(1, 1, dpi=300)

    table_carr_head_list = [
        'Satellite System','SV', 'Epochs',
        'Min [cycles]',
        'Max [cycles]',
        'Mean [cycles]',
        'Std [cycles]',
        'Mad [cycles]'
        ]
    writer_carr.writerow(table_carr_head_list)
        
    tt = str(systr[:-1])+' '+get_band_name(freq)+' D.D. '
    for sv in sv_range:
        if sv == ref_sv:
            continue
            
        plt_sv_id_code = numpy.where( output_dd_code_array[:,2] == sv )
        plt_sv_id_carr = numpy.where( output_dd_carr_array[:,2] == sv )
        dd_code_row = [systr, int(sv)]
        dd_carr_row = [systr, sv]
        if len(plt_sv_id_code[0])>0:
            ax1.plot(output_dd_code_array[plt_sv_id_code[0],0], output_dd_code_array[plt_sv_id_code[0],5], '.')
            code_mean = numpy.mean(output_dd_code_array[plt_sv_id_code[0],5])
            dd_code_row.append( int(len(plt_sv_id_code[0])) )
            dd_code_row.append( "{:.3f}".format(float(min(output_dd_code_array[plt_sv_id_code[0],5]))) )
            dd_code_row.append( "{:.3f}".format(float(max(output_dd_code_array[plt_sv_id_code[0],5]))) )
            dd_code_row.append( "{:.3f}".format(float(code_mean)) )
            dd_code_row.append( "{:.3f}".format(float(numpy.std(output_dd_code_array[plt_sv_id_code[0],5]))) )
            dd_code_row.append( "{:.3f}".format(float(numpy.mean(numpy.absolute( output_dd_code_array[plt_sv_id_code[0],5] - code_mean )))) )
            dd_code_stat_rows.append(dd_code_row)
            writer_code.writerow(dd_code_row)                
            
        if len(plt_sv_id_carr[0])>0:
            ax2.plot(output_dd_carr_array[plt_sv_id_carr[0],0], output_dd_carr_array[plt_sv_id_carr[0],5], '.')
            carr_mean = numpy.mean(output_dd_carr_array[plt_sv_id_carr[0],5])
            dd_carr_row.append( int(len(plt_sv_id_carr[0])) )
            dd_carr_row.append( "{:.3f}".format(float(min(output_dd_carr_array[plt_sv_id_carr[0],5]))) )
            dd_carr_row.append( "{:.3f}".format(float(max(output_dd_carr_array[plt_sv_id_carr[0],5]))) )
            dd_carr_row.append( "{:.3f}".format(float(carr_mean)) )
            dd_carr_row.append( "{:.3f}".format(float(numpy.std(output_dd_carr_array[plt_sv_id_carr[0],5]))) )
            dd_carr_row.append( "{:.3f}".format(float(numpy.mean(numpy.absolute( output_dd_carr_array[plt_sv_id_carr[0],5] - carr_mean )))) )
            dd_carr_stat_rows.append(dd_carr_row)
            writer_carr.writerow(dd_carr_row)

    output_dd_code_stat_array = numpy.asarray(dd_code_stat_rows)
    output_dd_carr_stat_array = numpy.asarray(dd_carr_stat_rows)
    
    # close csv files
    f_code.close()
    f_carr.close()
    
    # D.D. code phase plots
    ax1.set_xlabel('Time [sec]')
    ax1.set_ylabel('m')
    ax1.set_title(tt+'Code Phase {Ref SV#'+str(int(ref_sv))+' S.D. $\mu=$'+'%.2f'%ref_sd_code_mean+',$\sigma=$'+'%.2f'%ref_sd_code_std+'}')
    fname = './'+str(systr[:-1])+'_'+get_band_name(freq)+'_D.D_code.png'
    fig1.savefig(fname, dpi=300)
    plt.close(fig1)
    
    # D.D. carrier phase plots    
    ax2.set_xlabel('Time [sec]')
    ax2.set_ylabel('cycle')
    ax2.set_title(tt+'Carrier Phase {Ref SV#'+str(int(ref_sv))+' S.D. $\mu=$'+'%.2f'%ref_sd_carr_mean+',$\sigma=$'+'%.2f'%ref_sd_carr_std+'}')
    fname = './'+str(systr[:-1])+'_'+get_band_name(freq)+'_D.D_carrier.png'
    fig2.savefig(fname, dpi=300)
    plt.close(fig2)

    # D.D. code and carrier pase tables
    ax3.table(cellText=dd_code_stat_rows, colLabels=table_code_head_list, loc="center")
    ax3.set_title(tt+'Code Phase Statistics, Ref SV#'+str(int(ref_sv)), y=0.85, pad=-14)
    ax3.axis('tight')
    ax3.axis('off')
    fname = './'+str(systr[:-1])+'_'+get_band_name(freq)+'_D.D_code_table.png'
    fig3.savefig(fname, dpi=300)
    plt.close(fig3)
    
    ax4.table(cellText=dd_carr_stat_rows, colLabels=table_carr_head_list, loc="center")
    ax4.set_title(tt+'Carrier Phase Statistics, Ref SV#'+str(int(ref_sv)), y=0.85, pad=-14)
    ax4.axis('tight')
    ax4.axis('off')
    fname = './'+str(systr[:-1])+'_'+get_band_name(freq)+'_D.D_carr_table.png'
    fig4.savefig(fname, dpi=300)
    plt.close(fig4)
    
    return output_dd_code_array, output_dd_carr_array, output_dd_code_stat_array, output_dd_carr_stat_array