diff --git a/CEP/Calibration/ExpIon/src/fitClockTEC.py b/CEP/Calibration/ExpIon/src/fitClockTEC.py
index 463645ce0d8e01b2d75d8e812ea4559cbb3299b0..c61a29f8257b35c929f147d710dea7c745ca9cbd 100644
--- a/CEP/Calibration/ExpIon/src/fitClockTEC.py
+++ b/CEP/Calibration/ExpIon/src/fitClockTEC.py
@@ -390,7 +390,7 @@ def getClockTECBaselineFit(ph,amp,freqs,SBselect,polIdx,stIdx,useOffset=False,st
stepdTEC=0.005
else:
stepdTEC=0.001 # faster wrapping for LBA
-
+ stepDelay=3
succes=False
initprevsol=False
@@ -421,9 +421,9 @@ def getClockTECBaselineFit(ph,amp,freqs,SBselect,polIdx,stIdx,useOffset=False,st
iD1=-4
iD2=4
else:
- iTEC1=-2
- iTEC2=2
- iD1=-300
+ iTEC1=-1.5
+ iTEC2=1.5
+ iD1=-50
iD2=300
print "First",iTEC1,iTEC2,iD1,iD2
@@ -485,7 +485,7 @@ def add_to_h5_func(h5file,data,name='test'):
-def getAll(ionmodel,refstIdx=0,doClockTEC=True,doRM=False,add_to_h5=True,stationSelect='',label='fit',SBselect='all',allBaselines=True,useOffset=False,initFromPrevious=False,flagBadChannels=False,flagcut=1.5,chi2cut=300.,removePhaseWraps=False,combine_pol=False,fixedClockforCS=False,timerange='all',CStec0=False):
+def getAll(ionmodel,refstIdx=0,doClockTEC=True,doRM=False,add_to_h5=True,stationSelect='',label='fit',SBselect='all',allBaselines=True,useOffset=False,initFromPrevious=False,flagBadChannels=False,flagcut=1.5,chi2cut=30000.,removePhaseWraps=False,combine_pol=False,fixedClockforCS=False,timerange='all',CStec0=False):
global tecarray
global clockarray
global offsetarray
@@ -578,27 +578,19 @@ def getAll(ionmodel,refstIdx=0,doClockTEC=True,doRM=False,add_to_h5=True,station
#halfwraps=np.remainder(np.round(np.absolute(wraps[stationIndices]*2)),2)==1
#print "found halfwraps for",np.array(stations)[halfwraps]
if CStec0:
-
- #get slope of average TECvs lattitude
- lats=ionmodel.piercepoints.cols.positions[:][timerange[0]:timerange[1],0,:,1]
- stat_select=ionmodel.stat_select[:]
- lats=lats[:,stationIndices[stat_select]]
- lats-=lats[:,[0]] #relative -> no offset
- lats=np.average(lats,axis=0)
- avgtec=np.average(tecarray[timerange[0]:timerange[1],stationIndices,pol]-tecarray[timerange[0]:timerange[1],[0],pol],axis=0)+steps[0]*(np.round(wraps[stationIndices])-np.round(wraps[0]))
-
-
- slope= np.dot(1./(np.dot(lats.T,lats)), np.dot(lats.T,avgtec.T)).flatten()
- print "slope",slope
-
-
-
- offsets=-1* (avgtec-slope*lats)*2.*np.pi/steps[0]
- print "average TEC",avgtec
+ pos=myion.station_positions[:]
+ lats=np.degrees(np.arctan2(pos[:,2],np.sqrt(pos[:,0]*pos[:,0]+pos[:,1]*pos[:,1])))
+ lats-=lats[0]
+ lats=lats[stationIndices]
+ TEC=tecarray[timerange[0]:timerange[1],stationIndices,pol]-tecarray[timerange[0]:timerange[1],[0],pol]+steps[0]*(np.round(wraps[stationIndices])-np.round(wraps[0]))
+ slope=np.dot(1./(np.dot(lats.T,lats)), np.dot(lats.T,TEC.T))
+ chi2=np.sum(np.square(TEC-lats*slope[:,np.newaxis]),axis=1)
+ chi2select=chi2<np.average(chi2)
+ print "wraps",wraps
+ print "slope",slope[chi2select][0]
+ offsets=-1*(np.average(TEC[chi2select]-lats*slope[chi2select,np.newaxis]),axis=0)*2.*np.pi/steps[0]
print "step",steps[0]
print offsets
- remainingwraps=np.zeros(np.sum(stationIndices),dtype=np.float)
- #remainingwraps[CSstations]=np.round(offsets[CSstations]/(2*np.pi))-np.round(wraps[stationIndices][CSstations])
remainingwraps=np.round(offsets/(2*np.pi))#-np.round(wraps[stationIndices])
print remainingwraps
wraps[stationIndices]+=remainingwraps
@@ -616,7 +608,7 @@ def getAll(ionmodel,refstIdx=0,doClockTEC=True,doRM=False,add_to_h5=True,station
initSol[:,0]=tecarray[timerange[0],stationIndices,pol]+steps[0]*np.round(wraps[stationIndices])
initSol[:,1]=clockarray[timerange[0],stationIndices,pol]+steps[1]*np.round(wraps[stationIndices])
#initSol[:,1]=np.average(clockarray[:,stationIndices,pol]-clockarray[:,[0],pol],axis=0)+steps[1]*np.round(wraps[stationIndices])
- print "wraps",np.round(wraps[stationIndices])
+ print "final wraps",np.round(wraps[stationIndices])
print "prev solutions", clockarray[timerange[0],stationIndices,pol]
print "init Clock with", initSol[:,1]
print "prev solutions TEC", tecarray[timerange[0],stationIndices,pol]