def film1D(base, ext, debut, fin, vmin, vmax, dt = 0):
  """ displays the evolution in time of a 1-D solution
  Usage : film1D(base, ext, debut, fin, vmin, vmax)
          film1D(base, ext, debut, fin, vmin, vmax, dt)

  base : root of output files
  ext : extension of output files
  debut : start index
  fin : final index
  vmin, vmax : y-range
  
  Example : film1D("Un", ".dat", 0, 20, -0.1, 1.1)
         will display frames contained in Un0.dat, Un1.dat, ...
         until Un19.dat with y in [-0.1, 1.1]
  
  Put a dt different from 0, if the animation is too fast, so that the system
  will wait dt seconds between each display """
  V = pylab.loadtxt(base + str(debut) + ext)
  line, = pylab.plot(V[:,0], V[:,1]);
  pylab.axis(pylab.array([V[0, 0], V[V.shape[0]-1,0], vmin, vmax]));
  os.system('sleep ' +str(dt))
  for i in range(debut+1, fin):
    V = pylab.loadtxt(base + str(i) + ext)
    line.set_ydata(V[:,1])
    pylab.draw()
    os.system('sleep ' +str(dt))

def main():  

   FILE1 = "./TO_PLOT.dat" 
   FILE2 = "./test.dat" 
   FILE3 = "./train.dat" 
   PATH_TO_OUTPUT = "./plot.png"


   f1 = pylab.loadtxt(FILE1,dtype = str, unpack=True)
   f2 = pylab.loadtxt(FILE2,dtype = str, unpack=True)
   f3 = pylab.loadtxt(FILE3,dtype = str, unpack=True)

   pylab.clf()
   pylab.cla()

   pylab.scatter(f1[0], f1[1], c=color[1], label='10-fold KNN')
   pylab.plot(f2[0], f2[1],  c=color[2], label='Test data')
   pylab.scatter(f3[0], f3[1],  c=color[3], label='Train data')

   pylab.ylabel('Acccuracy', fontsize=14)  
   pylab.xlabel('k', fontsize=14)       
   pylab.legend(loc=1, prop={'size':10})
   pylab.xlim(0.1,100)

		    
   pylab.savefig(PATH_TO_OUTPUT)

def load_act(sym, source):
    ''' Load atomic collision transition file from mechanical quantum computations
        3 ascii input files:
        ael.dat: atomic energy levels in the mael.py format
        temp.dat: list of temperatures
        ups.dat: effective collision strengths
    '''
    path = '/home/tmerle/development/formato2/ad/ct/e/qm/'+SYM+'/'+source+'/'

    dumb = pl.loadtxt(path+'ael.dat', dtype='str')
    ael_list = [State(e=e, g=g, cfg=cfg, term=term, p=p, ref=ref) for e, g, cfg, term, p, ref in dumb]
    #[ael.print() for ael in ael_list]

    t_list = pl.loadtxt(path+'temp.dat', dtype='float')
    nt = len(t_list)

    dumb = pl.loadtxt(path+'ups.dat', dtype='float')
    
    l_list, u_list, ups_list = [], [], []

    for rec in dumb:
        l_list.append(int(rec[0]))
        u_list.append(int(rec[1]))
        ups_list.append(map(float, rec[2:]))            

    col_list = []

    for i, (l, u) in enumerate(zip(l_list, u_list)):
        #print(l, u, ups_list[i])
        col = Collision(type='QM', lower=ael_list[l-1], upper=ael_list[u-1], t_list=t_list, ups_list=ups_list[i], ref=source, kw='UPS_E')
        col_list.append(col)

    return col_list

def p0_1s2():
    pl.clf()
    ourdat  = pl.loadtxt('1p0_1s2_J1.dat')
    ax1 = pl.axes()
    pl.figure(1, figsize=(8, 6), facecolor='white')
    # pl.plot(20.0, 4.7, 'gs', ms=7)
    # pl.plot(20.0, 1.4, 'bo', ms=7)
    label = '1p0' r'$\rightarrow$'  '1s2 ' '($\Delta J=1)$'
    # pl.plot(ourdat[:,0], ourdat[:,1], 'r--', ms=7,  label=label)
    pl.plot(ourdat[:,0], ourdat[:,2], 'b', ms=7, label=label)
    # pl.plot(ourdat[:,0], ourdat[:,3], 'g-.', ms=7)
    zatsdat = pl.loadtxt('../zats_figs/zats_1s2.dat')
    pl.plot(zatsdat[:,0], zatsdat[:,1], 'ro',)
    pl.plot(zatsdat[:,2], zatsdat[:,3], 'y*', ms=9, lw=5) # Hoshino
    pl.plot(zatsdat[:,4], zatsdat[:,5], 'mp', ms=7, lw=5) # Khakoo
    pl.plot(zatsdat[:,6], zatsdat[:,7], 'gs', ms=7, lw=5) # Filopivic
    pl.plot(zatsdat[:,8], zatsdat[:,9], 'cD', ms=9, lw=5) # Chutjian
    pl.legend(loc=1, frameon=False, scatterpoints=0)
    pl.xlim([12,35])
    # pl.ylim([0,5])
    pl.xlabel('Energy (eV)', fontsize=16)
    pl.ylabel('Cross section (MB)', fontsize=16)
    pl.savefig('1p0-1s2_J1_plot.eps')
    pl.show()
    return

	def load_csv(self,f):
		"""
		Loading data from a csv file. Uses pylab's load function. Seems much faster
		than scipy.io.read_array.
		"""
		varnm = f.readline().split(',')

		# what is the date variable's key if any, based on index passed as argument
		if self.date_key != '':
			try:
				rawdata = pylab.loadtxt(f, delimiter=',',converters={self.date_key:pylab.datestr2num})		# don't need to 'skiprows' here
			except ValueError:																				# if loading via pylab doesn't work use csv
				rawdata = self.load_csv_nf(f)	

				# converting the dates column to a date-number
				rawdata[self.date_key] = pylab.datestr2num(rawdata[self.date_key])

			self.date_key = varnm[self.date_key]
		else:
			try:
				rawdata = pylab.loadtxt(f, delimiter=',')														# don't need to 'skiprow' here
			except ValueError:																				# if loading via pylab doesn't work use csv
				rawdata = self.load_csv_nf(f)	

		# making sure that the variable names contain no leading or trailing spaces
		varnm = [i.strip() for i in varnm]

		# transforming the data into a dictionary
		if type(rawdata) == list:
			# if the csv module was used
			self.data = dict(zip(varnm,rawdata))
		else:
			# if the pylab.load module was used
			self.data = dict(zip(varnm,rawdata.T))

def p0_2p1():
    pl.clf()
    ourdat  = pl.loadtxt('1p0_2p1_J0.dat')
    ax1 = pl.axes()
    pl.figure(1, figsize=(8, 6), facecolor='white')
    # pl.plot(20.0, 5.0, 'gs', ms=7)
    # pl.plot(20.0, 5.2, 'bv', ms=7)
    # pl.plot(20.0, 2.7, 'bo', ms=7)
    label = '1p0' r'$\rightarrow$'  '2p1 ' '($\Delta J=0)$'
    # pl.plot(ourdat[:,0], ourdat[:,1], 'r--', ms=7,  label=label)
    pl.plot(ourdat[:,0], ourdat[:,2], 'b', ms=7, label=label)
    # pl.plot(ourdat[:,0], ourdat[:,3], 'g-.', ms=7)
    zatsdat = pl.loadtxt('../zats_figs/zats_2p1.dat')
    pl.plot(zatsdat[:,0], zatsdat[:,1]/10, 'ro',)
    pl.plot(zatsdat[:,2], zatsdat[:,3]/10, 'k^', ms=7, lw=5)
    pl.plot(zatsdat[:,4], zatsdat[:,5]/10, 'gs', ms=7, lw=5)
    pl.plot(zatsdat[:,6], zatsdat[:,7]/10, 'cD', ms=7, lw=5)
    pl.legend(loc=1, frameon=False, scatterpoints=0)
    pl.xlim([12,42])
    pl.ylim([0,15])
    pl.xlabel('Energy (eV)', fontsize=16)
    pl.ylabel('Cross section (MB)', fontsize=16)
    pl.savefig('1p0-2p1_J0_plot.eps')
    pl.show()
    return

def plot1(cursor):
    dage,dew,dmf,dM=gettable(cursor,cols='age,Ha_w,Massfrac,Mr',where='agn=0 AND Mr>-20 and Massfrac NOTNULL',table='sball')
    gage,gew,gmf,gM=gettable(cursor,cols='age,Ha_w,Massfrac,Mr',where='agn=0 AND Mr<-20 and Massfrac NOTNULL',table='sball')
    limit1x,limit1y=P.loadtxt('/home/tom/projekte/sdss/ages/mixred0',unpack=True)[:2]
    limit2x,limit2y=P.loadtxt('/home/tom/projekte/sdss/ages/mixblue0',unpack=True)[:2]

    ax1=P.subplot(1,2,1)
    P.title(r'$\mathrm{Dwarfs,}\,\, M_r > -20$')
    P.loglog(limit1x,limit1y,'--k',linewidth=2)
    P.loglog(limit2x,limit2y,'-k',linewidth=2)
    P.scatter(dage,dew*dcf,c=dmf,vmin=0,vmax=0.05,label='Age, M>-20',alpha=0.2,cmap=P.cm.OrRd,lw=1)
    P.xlabel(r'$\mathrm{Burst\, age}$')
    P.ylabel(r'$\mathrm{EW}(H\alpha)$')
    P.plot([1E6,1E10],[150,150],'k-.',lw=2)
    P.axis([4E6,2E9,55,9E2])

    P.subplot(1,2,2)
    P.title(r'$M_r < -20$')
    P.loglog(limit1x,limit1y,'--k',linewidth=2)
    P.loglog(limit2x,limit2y,'-k',linewidth=2)
    P.scatter(gage,gew*gcf,c=gmf,vmin=0,vmax=0.05,label='Age, M<-20',alpha=0.2,cmap=P.cm.OrRd,lw=1)
    P.xlabel(r'$\mathrm{Burst\, age}$')
    P.plot([1E6,1E10],[150,150],'k-.',lw=2)
    P.axis([4E6,2E9,55,9E2])
    P.setp(P.gca(),'yticklabels',[])
    P.subplots_adjust(wspace=0)

def runAll(ScaleDensity, b, f, Particle, Scale):
    import pylab as p
    import shutil, subprocess
    
    path = "./outData/"
    savepath = "./plots/compare/"
    inifilename = "mybody.ini"

    filenameMine = "mass.dat"
    filenameRockstar = "halos_0.ascii"
    filenameAmiga = "amiga.dat"

    
    inifile = open(inifilename, 'r')
    tmpinifile = open(inifilename + ".tmp", "w")

    for line in inifile:
    
        if line[0:15] == "ScaleDensity = ":
            line = "ScaleDensity = "+ str(ScaleDensity)+"\n "
        if line[0:4] == "b = ":
            line = "b = " + str(b)+"\n " 
        if line[0:4] == "f = ":
            line = "f = " + str(f)+"\n "
        if line[0:21] == "LinkingLenghtScale = ":
            line = "LinkingLenghtScale = "+ str(Scale)+"\n "
        if line[0:20] == "NrParticlesDouble = ":
            line = "NrParticlesDouble = " + str(Particle)+"\n " 
            
        tmpinifile.write(line)
    
    inifile.close()
    tmpinifile.close()
    shutil.move(inifilename + ".tmp", inifilename)

    #Run my halofinder
    #print "Running for: " +"ScaleDensity="+ str(ScaleDensity) + " b=" + str(b) + " f="+str(f)
    subprocess.call(["mpirun","-n","2", "./main"])
    
    #Read data from files
    mineData = p.loadtxt(path + filenameMine)
    rockstarData = p.loadtxt(path + filenameRockstar)
    amigaData = p.loadtxt(path + filenameAmiga)

    sortedDataR = p.sort(rockstarData[:,2])[::-1]
    sortedData = p.sort(mineData)[::-1]
    sortedDataA = p.sort(amigaData[:,3])[::-1]


    p.figure()
    p.loglog(sortedData,range(1,len(sortedData)+1))
    p.loglog(sortedDataR,range(1,len(sortedDataR)+1))
    p.loglog(sortedDataA,range(1,len(sortedDataA)+1))
    p.xlabel("log(Mass) Msun")
    p.ylabel("log(Nr of halos with mass >= Mass)")
    p.legend(("mine","Rockstar","AMIGA"))
    name = "MassFunction_ScaleDensity="+ str(ScaleDensity) + "_b=" + str(b) + "_f="+str(f) + "_Scale=" +str(Scale) + "_Double="+str(double)
    p.savefig(savepath+name+".png")

    def importfile(self,fname,params):
        # if even more sophisticated things are needed, just inherit THzTdData class
        #and override the importfile method
        #try to load the file with name fname
        #it should be possible to write this shorter
        try:
            #if no Y_col is specified            
            if params.has_key('Y_col'):
                #import it right away
                if params['dec_sep']=='.':
                    data=py.loadtxt(fname,
                                usecols=(params['time_col'],
                                         params['X_col'],
                                         params['Y_col']),
                                skiprows=params['skiprows'])
                                
                elif params['dec_sep']==',':
                    #if the decimal separator is , do a replacement
                    str2float=lambda val: float(val.replace(',','.'))
                    data=py.loadtxt(fname,
                                converters={params['time_col']:str2float,
                                            params['X_col']:str2float,
                                            params['Y_col']:str2float},
                                usecols=(params['time_col'],
                                         params['X_col'],
                                         params['Y_col']),
                                skiprows=params['skiprows'])                
            else:
                #import it right away
                if params['dec_sep']=='.':
                    data=py.loadtxt(fname,
                                usecols=(params['time_col'],
                                         params['X_col']),
                                skiprows=params['skiprows'])
                                
                elif params['dec_sep']==',':
                    #if the decimal separator is , do a replacement
                    str2float=lambda val: float(val.replace(',','.'))
                    data=py.loadtxt(fname,
                                converters={params['time_col']:str2float,
                                            params['X_col']:str2float},
                                usecols=(params['time_col'],
                                         params['X_col']),
                                skiprows=params['skiprows'])
                dummy_Y=py.zeros((data.shape[0],1))
                data=py.column_stack((data,dummy_Y))
        except IOError:
            print "File " + fname + " could not be loaded"
            sys.exit()

        #scale the timaaxis
        data[:,0]*=params['time_factor']
        
        #if the measurement was taken in negative time direction, flip the data
        if data[1,0]-data[0,0]<0:
            data=py.flipud(data)
     
        return data

    def on_bl_button(self,event):
        self.selectq_control.SetString(string='Streamwise (u) Velocity')
        self.selectq_control.num = 2
        self.blcall = True 
        self.bldraw = True
        self.on_update_button(event)
        self.blcall = False
        self.jcur = self.Bl_jplane.GetValue()
        gnum = self.qfile_control.manual_text2.GetValue()
        if (self.multgrd):      
            new_str1 = ''.join([self.qfile, '\n','qvalbl.txt', '\n', str(gnum), '\n', str(self.jcur) ,' ', 
                                str(self.jcur) ,' ', str(1), '\n' ,  str(1),' ',
                                str(self.nkind),' ',str(1), '\n' ,str(1) ,
                                ' ',str(1) ,' ', str(1), '\n',
                                str(2), '\n','n'])  
        else:
            new_str1 = ''.join([self.qfile, '\n','qvalbl.txt', '\n', str(self.jcur) ,' ', 
                                str(self.jcur) ,' ', str(1), '\n' ,  str(1),' ',
                                str(self.nkind),' ',str(1), '\n' ,str(1) ,
                                ' ',str(1) ,' ', str(1), '\n',
                                str(2), '\n','n'])  
            
        new_str2 = ''.join([self.grid, '\n',  str(gnum), '\n', 'cordsbl.txt', '\n', str(self.jcur) ,' ', 
                            str(self.jcur) ,' ', str(1), '\n' ,  str(1),' ',
                            str(self.nkind),' ', str(1), '\n' ,str(1) ,
                            ' ',str(1) ,' ', str(1)]) 
        p2 = subprocess.Popen(['getgridcords'], stdin=PIPE, stdout=PIPE)
        o2,e2 = p2.communicate(input=new_str2) 
        p1 = subprocess.Popen(['listplotvar'], stdin=PIPE, stdout=PIPE)
        o1,e1 = p1.communicate(input=new_str1)
        p3 = subprocess.Popen(['tail', '-n', '+7'], stdin=PIPE, stdout=PIPE)
        o3,e3 = p3.communicate(input=o1) 
        p4 = subprocess.Popen(['head', '-n', '-1'], stdin=PIPE, stdout=PIPE)
        o4,e4 = p4.communicate(input=o3) 

        self.file1 = pylab.loadtxt('qvalbl.txt')
        self.file2 = pylab.loadtxt('cordsbl.txt')
        
        os.remove('qvalbl.txt')
        os.remove('cordsbl.txt')
        
        tot = len(self.file1)
        u = [.5 for i in range(tot)]
        x = [.5 for i in range(tot)]
        y = [.5 for i in range(tot)]
        for k in range(0,tot):
            u[k] = float(self.file1[k][3])
            x[k] = float(self.file2[k][3])
            y[k] = float(self.file2[k][4])
        xoverc = round(x[0],5)
        self.x_over_c.SetLabel( ''.join(['x/c =   ', str(xoverc)]) )
        if (self.mom_norm.IsChecked()):
            self.update_bl_mom(x,y,u,(self.nkind-1),(self.jcur))
        else:
            self.update_bl_mom(x,y,u,(self.nkind-1),(self.jcur))
            self.update_bl(x,y,u,(self.nkind-1),(self.jcur))

def getData():
    
    # load the fitting data for X and y and return as elements of a tuple
    # X is a 100 by 10 matrix and y is a vector of length 100
    # Each corresponding row for X and y represents a single data sample

    X = pl.loadtxt('fittingdatap1_x.txt')
    y = pl.loadtxt('fittingdatap1_y.txt')

    return (X,y) 

def plot(experimento = None):
    path = '/home/rafaelbeirigo/ql/experiments/'
    
    ql = pl.loadtxt(path + experimento + "/QL/w.out")
    prql = pl.loadtxt(path + experimento + "/PRQL/w.out")

    pl.xlabel("Episodes")
    pl.ylabel("W")

    pl.plot(ql, label = 'Q-Learning')
    pl.plot(prql, label = 'PRQ-Learning prob')

    pl.legend(loc = 0)

def get_consts(kind):
    if kind == "hs":
        g1, g2 = 1.270042427, 1.922284066
        k_0, d_1, = -1.2540, 2.4001
        eta, Y_0, Y_1, Omega_1, Theta_1, H_A, H_B, Y_a4 = py.loadtxt('../tables/kn-layer-hs.txt').T
        Y_1 = 2*Y_1
        return g1, g2, k_0, d_1, eta, Y_0, H_A, Theta_1, 0.43, 0.078, 0.7
    elif kind == "bkw":
        g1, g2 = 1., 1.
        k_0, d_1 = -1.01619, 1.30272
        eta, Y_0 = py.loadtxt('../tables/kn-layer-bkw.txt').T
        H_A = Y_0 / 2
        return g1, g2, k_0, d_1, eta, Y_0, H_A, Theta_1, 0.345, 0.172
    else:
        raise "Invalid potential!"

 def read_data(cls, path):
     """
     AuxiliaryFiles中のデータ収集
     dictのリスト配列として収集する
     """
     comp_path = os.path.join(path, cls.base_dirc, cls.comp_file)
     comp = pylab.loadtxt(comp_path, comments='#')
     total = np.array([[sum(x) for x in comp]])
     frac = comp / total.T
     energy_path = os.path.join(path, cls.base_dirc, cls.energy_file)
     energy = pylab.loadtxt(energy_path, comments="#")
     energy = energy / total[:][0]
     data = [{'frac_atoms': list(x), 'energy': y}
             for x, y in zip(frac, energy)]
     return(data)

    def getData(self):
        # check if a separator is a white space
        if is_empty(self.params.separator):
            _data = pl.loadtxt(self._file,
                              dtype=(str),
                              skiprows=self.headers_count,
                              unpack=True)
        else:
            _data = pl.loadtxt(self._file,
                              dtype=(str),
                              skiprows=self.headers_count,
                              unpack=True,
                              delimiter=self.params.separator)

        #some issues with preparing data, details in prepare_data_arrays func.
        return prepare_data_arrays(self._file, self.headers_count, _data)

def get_n(material,E=8):
    """
    "by LW 07/04/2011 function get the index of refraction from stored data file,
    index of refraction is a .dat file from http://henke.lbl.gov/optical_constants/getdb2.html 
    (energy range: 2-30keV,delete the header lines, name the file n_material.dat) 
    calling sequence: n=get_n(material,E) where n is the complex refractive index detlta-i*beta, E: X-ray energy in keV"
    """
    #get list_of supported materials from data file directory:
    xdatafiles = [ f for f in listdir(datapath) if isfile(join(datapath,f)) ]    
    name=[]
    for i in range(0, np.size(xdatafiles)):
        m=re.search('(?<=n_)\w+', xdatafiles[i])
        if m is not None:
            name.append(m.group(0))             
    
    E=np.array(E)
    if material in name:
        loadn=datapath+'n_'+material+'.dat'
        n=pl.loadtxt(loadn,comments='%')
        if np.min(E)>=np.min(n[:,0]/1000) and np.max(E)<=np.max(n[:,0]/1000):
            d=np.interp(E*1000,n[:,0],n[:,1])
            b=np.interp(E*1000,n[:,0],n[:,2])
            return d-1j*b
        else: print 'error: energy '+"%3.4f" %E +'[keV] out of range ('+"%3.4f" % np.min(n[:,0]/1000)+'=<E<='+"%3.4f" % np.max(n[:,0]/1000)+'keV)'
    elif material=='material?':
        print 'list of supported materials (based on data files in directory '+datapath+':'
        print name
    else: print 'error: non recognized material, please create index of refraction file first. Type "get_n?" for instructions; type get_n("material?") for list of supported materials'

def get_ac(material,E=8):
    """
    by LW 10/03/2010
    function calculates the critical angle for total external reflection as a function of
    the material and the X-ray energy according to ac=sqrt(2*delta)
    index of refraction is a .dat file from http://henke.lbl.gov/optical_constants/getdb2.html
    (energy range: 2-30keV,delete the header % lines, name the file n_material.dat) % 
    calling sequence: ac=get_ac(material,E) where ac: critial angle in degrees, E [keV] (default: 8keV)
    type get_ac(\'materilal?\') to show list of supported materials"
    """
    
    #get list_of supported materials from data file directory:
    xdatafiles = [ f for f in listdir(datapath) if isfile(join(datapath,f)) ]    
    name=[]
    for i in range(0, np.size(xdatafiles)):
        m=re.search('(?<=n_)\w+', xdatafiles[i])
        if m is not None:
            name.append(m.group(0))             
    
    E=np.array(E)
    if material in name:
        loadn=datapath+'n_'+material+'.dat'
        n=pl.loadtxt(loadn,comments='%')
        if np.min(E)>=np.min(n[:,0]/1000) and np.max(E)<=np.max(n[:,0]/1000):
            d=np.interp(E*1000,n[:,0],n[:,1])
            return np.degrees(np.sqrt(2*d))
        else: print 'error: energy '+"%3.4f" %E +'[keV] out of range ('+"%3.4f" % np.min(n[:,0]/1000)+'=<E<='+"%3.4f" % np.max(n[:,0]/1000)+'keV)'
    elif material=='material?':
        print 'list of supported materials (based on data files in directory '+datapath+':'
        print name       
    else: print 'error: non recognized material, please create index of refraction file first. Type "get_ac?" for instructions; type get_ac("material?") for list of supported materials' 

def meanOfMeans():
    initialExperiment = 25
    finalExperiment = 35
    suffix = "piConc"
    experimentFolder = "170"

    rootDir = os.path.join("/", "home", "rafaelbeirigo", "ql", "experiments", experimentFolder)

    means = []
    meanFileName = "W_avg_list_mean.out"
    for experiment in range(initialExperiment, finalExperiment + 1):
        meanFilePath = os.path.join(rootDir, str(experiment), "PRQL", meanFileName)
        # open experiment the file that contains the "individual" mean
        # into a list
        mean = pl.loadtxt(meanFilePath)

        # this must be done because the function meanError.meanMultiDim()
        # needs each item in the mean list to be a list itself.
        mean = [[item] for item in mean]

        # append it to the list of means
        means.append(mean)

    # obtain the global mean considering the list of individual means
    globalMean = meanError.meanMultiDim(means)
    ## print globalMean

    # means2spreadsheet([globalMean])

    pl.savetxt(os.path.join(rootDir, "W_avg_list_mean." + suffix + ".out"), globalMean, fmt="%1.6f")

    return globalMean

def plot(filename, column, label):
    data = py.loadtxt(filename).T
    X, Y = data[0], data[column]
    mask = (X >= xmin) * (X <= xmax)
    X, Y = X[mask], corrector(Y[mask])
    aY, bY = np.min(Y), np.max(Y)
    pad = 7 if aY < 0 and -bY/aY < 10 else 0
    py.ylabel(r'$' + label + r'$', y=y_coord, labelpad=8-pad, rotation=0)
    py.plot(X, Y, '-', lw=1)
    py.xlabel(r'$\zeta$', labelpad=-5)
    py.xlim(xmin, xmax)
    ax = py.axes()
    ax.axhline(lw=.5, c='k', ls=':')
    specify_tics(np.min(Y), np.max(Y))
    if inside:
        print "Plot inside plot"
        ax = py.axes([.25, .45, .4, .4])
        mask = X <= float(sys.argv[7])
        py.plot(X[mask], Y[mask], '-', lw=1)
        ax.tick_params(axis='both', which='major', labelsize=8)
        ax.axhline(lw=.5, c='k', ls=':')
        ax.xaxis.set_major_locator(MultipleLocator(0.5))
        ax.yaxis.set_major_locator(LinearLocator(3))
        ymin, _, ymax = ax.get_yticks()
        if (ymax + ymin) < (ymax-ymin)/5:
            y = max(ymax, -ymin)
            py.ylim(-y, y)

    def __init__(self, id, offline=False):
        self.id = id
        self.magSampleList = []
        self.count = 0
        self.filename = "./calibration/" + str(self.id) + ".cal"
        if offline:
            # load data from file
            #file = open("filename", mode='r')
            data = loadtxt(self.filename, skiprows=3, delimiter=',')

            for i in range(len(data)):
                x = data[i,0]
                y = data[i,1]
                z = data[i,2]

                mag = (x,y,z)

                self.magSampleList.append(mag)

        else:
            self.collectData = True
            self.c = CaptureThread(self.handleData)
            self.c.start()
            print "rotate device in all axes for 30 seconds..."
            #raw_input("Rotate device. Press ENTER when done...")
            sleep(10)
            print "20" + " (" + str(len(self.magSampleList)) + ")"
            sleep(10)
            print "10" + " (" + str(len(self.magSampleList)) + ")"
            sleep(10)
            print "Done" + " (" + str(len(self.magSampleList)) + ")"
            self.collectData = False
            self.c.stop = True

        self.processData()