files = glob("/home/greg/Documents/Summer Research/rinex files/ma*")
poop=rinexobs(files[6])

plt.figure(figsize=(14,14))
ax1 = plt.subplot(211)
ax1.xaxis.set_major_formatter(fmt)

plt.plot(2.85*(poop[:,23,'P2','data']*1.0E9/3.0E8-poop[:,23,'C1','data']*1.0E9/3.0E8)[10:],
         '.',markersize=3,label='pr tec')
plt.plot(2.85E9*((poop[:,23,'L1','data'])/f1-(poop[:,23,'L2','data'])/f2)[10:],
         '.',markersize=3,label='ph tec')

plt.title('mah13 sv23, biased')
plt.xlabel('time')
plt.ylabel('TECu')
plt.legend()
plt.grid()
plt.show()


"""
Explanation: These graphs have all been produced according to appendix A of "Effect of GPS System Biases on Differential
Group Delay Measurements" available here: http://www.dtic.mil/dtic/tr/fulltext/u2/a220284.pdf. It is referenced in Anthea Coster's 1992 paper "Ionospheric Monitoring System". Equation 2.3 in the paper states that TEC in TECu is equal to the propogation time difference between the two frequencies. In order to find the propogation time, I divided the distance by the speed of light. I did the same thing with the phase but first I had to multiply the cycles (which is available from the rinex file) by the wavelength. 
Basically what I want to explore with this notebook is combining pseudorange and phase data to get better TEC plots. Pseudorange is an absolute measurement but is subject to a lot of noise. Phase is low noise but has an unknown cycle ambiguity. According to the previously mentioned paper, the cycle ambiguity in phase data can be estimated by finding the average difference between the pseudorange and phase. In this script I converted pseudorange into cycles by dividing by wavelength. Then I calculated TEC according to equation 2.3 and plotted it next to the pseudorange calculated TEC and the biased phase calculated TEC.
Shown above is only a short time period free of cycle slips and loss of lock. In order to implement this averaging strategy, I think you have to re-average every time there is a cycle slip or missing data.
End of explanation
"""


sl=200
plt.figure(figsize=(15,15))
ax1=plt.subplot(211)
ax1.xaxis.set_major_formatter(fmt)

plt.plot(2.85E9*(poop[:,23,'P2','data']/3.0E8
                 -poop[:,23,'C1','data']/3.0E8),'b.',label='prtec',markersize=3)

for i in range(int(len(poop[:,23,'L1','data'])/sl)):

    phtec = 2.85E9*(poop[poop.labels[sl*i:sl*(i+1)],23,'L1','data']/f1
                   -poop[poop.labels[sl*i:sl*(i+1)],23,'L2','data']/f2)
          
    prtec = 2.85E9*(poop[poop.labels[sl*i:sl*(i+1)],23,'P2','data']/3.0E8
                   -poop[poop.labels[sl*i:sl*(i+1)],23,'C1','data']/3.0E8)
    
    b = np.average((phtec-prtec)[np.logical_not(np.isnan(phtec-prtec))])
    plt.plot(phtec-b,'r-',linewidth=3,label='')
    

plt.axis([poop.labels[10],poop.labels[10000],-50,50])
plt.title('bias corrected phase data')
plt.xlabel('time')
plt.ylabel('TECu')
plt.grid()
plt.legend()
plt.show()



"""
Explanation: This plot is uncorrected, it is a remake of the plot in Anthea's email on Wed, Jun 15, 2016 at 7:06 AM. The phase calculated TEC appears to match the pseudorange calculated TEC in parts, but other parts are offset by a bias. The next script and plot shows how the phase data can be shifted to line up with the average of the pseudorange data. It does so according to how large the slice value is. Basically it slices up the whole data set into a specifiable range of points and does the averaging on each range of points individually in order to fix phase lock cycle slips. This isn't the best way of doing it, the program should check for loss of lock and shift the phase data in chunks based on that. I am going to read more about cycle slips now.
End of explanation
"""


f1 = 1575.42E6 
f2 = 1227.6E6  

svn = 23
L1 = -1*3.0E8*poop[:,svn,'L1','data']/f1  #(1a)
L2 = -1*3.0E8*poop[:,svn,'L2','data']/f2  #(1b)
P1 = poop[:,svn,'C1','data']              #(1c)
P2 = poop[:,svn,'P2','data']              #(1d)


#wide lane combination
wld = 3.0E8/(f1-f2)
Ld = (f1*L1-f2*L2)/(f1-f2)                #(3)
prd = (f1*P1+f2*P2)/(f1+f2)               #(4)
bd = (Ld-prd)/wld                         #(5)

#wide lane cycle slip detection
bdmean = bd[1]
rms = 0
nans_in_a_row = 0
bdindex=2
ph_con_arc=0
avb=[]
rmsb=[]
lowers=[]
uppers=[]
inarc=False

while(nans_in_a_row<10):
    
    if not np.isnan(bd[bdindex]):
        
        if(bd[i]>4*np.sqrt(rms)):
            avb.append(bdmean)
            rmsb.append(np.sqrt(rms))
            bdmean=bd[i+1]
            rms=0
            ph_con_arc=0
            uppers.append(bdindex)
            inarc=False
            
        else:
            if not inarc:
                inarc=True
                lowers.append(bdindex)
            ph_con_arc+=1
            nans_in_a_row=0
            rms = rms+((bd[bdindex]-bdmean)**2-rms)/(ph_con_arc) #(8b)
            bdmean = bdmean+(bd[bdindex]-bdmean)/(ph_con_arc)    #(8a)
            
    else:
        if inarc:
            avb.append(bdmean)
            rmsb.append(rms)
            ph_con_arc=0
            bdmean=bd[bdindex+1]
            rms=0
            inarc=False
            uppers.append(bdindex)
        nans_in_a_row += 1
    bdindex+=1
            
print('ended at index {}'.format(bdindex-nans_in_a_row))

lengths=[a-b for a,b in zip(uppers,lowers)]
eim = [np.sqrt(a)/np.sqrt(b-1) for a,b in zip(rmsb,lengths)]
        

#ionospheric combination
LI = L1-L2      #(6)
PI = P2-P1      #(7)

#ionospheric cycle slip detection
# get x and y vectors
mask=~np.isnan(PI)
x = np.arange(len(PI[mask]))
y = PI[mask]

# calculate polynomial
z = np.polyfit(x, y, 6)
f = np.poly1d(z)

# calculate new x's and y's
x_new = np.linspace(x[0], x[-1], len(x))
Q = f(x_new)

residual = LI[mask]-Q



print(lowers,uppers)
plt.figure(figsize=(16,16))
plt.subplot(211).xaxis.set_major_formatter(fmt)
plt.plot(2.85E9*(poop[:,23,'P2','data']/3.0E8
                 -poop[:,23,'C1','data']/3.0E8),'b.',label='prtec',markersize=3)
plt.plot(2.85E9*((poop[:,23,'L1','data'])/f1-(poop[:,23,'L2','data'])/f2)[10:],
         'g.',markersize=3,label='ph tec')
for i in range(len(lowers)):

    phtec = 2.85E9*(poop[poop.labels[lowers[i]:uppers[i]],23,'L1','data']/f1
                   -poop[poop.labels[lowers[i]:uppers[i]],23,'L2','data']/f2)
          
    prtec = 2.85E9*(poop[poop.labels[lowers[i]:uppers[i]],23,'P2','data']/3.0E8
                   -poop[poop.labels[lowers[i]:uppers[i]],23,'C1','data']/3.0E8)
    
    b = np.average((phtec-prtec)[np.logical_not(np.isnan(phtec-prtec))])
    plt.plot(phtec-b,'r-',linewidth=2,label='')
    
plt.xlim([poop.labels[1],poop.labels[8000]])
plt.legend()
plt.grid()
plt.show()


"""
Explanation: try some stuff out from "An Automatic Editing Algorithm for GPS Data" by Blewitt
End of explanation
"""


