def kip_cont(self, modstart, modstop, dcoeff_thresh=1.e12,
                 xres=1000, ylims=[0., 0.], xlims=[0., 0.], yres=2000,
                 ixaxis='log_time_left', outfile='',
                 landscape_plot=False, codev='KEP', kepswitch=12555,
                 outlines=False, write_preproc=False, hatching=False):
        """
        This function creates a Kippenhahn diagram as a contour plot of
        the se output using a convection flag or diffusion coefficient
        threshold.

        Parameters
        ----------
        modstart, modstop : integer
            First and last cycle numbers to plot.
        dcoeff_thresh : float, optional
            Diffusion coefficient threshold, above which assume
            convection is present (only for MESA, where
            convection_indicator is not there)).  The default is 1.e12.
        xres, yres : integer, optional
            x and y resolution.  The default for xres is 1000, the
            default for yres is 2000.
        xlims, ylims : list, optional
            x and y plot limits.  The default is [0., 0.].
        ixaxis : string, optional
            Choose one of 'age', 'model_number' or 'log_time_left'.  The
            default is 'log_time_left'.
        outfile : string, optional
            Name of output file including extension, which will be saved
            in 300 dpi.  The default is ''.
        landscape_plot : boolean, optional
            For a landscape plot.  The default is False.
        codev : string, optional
            Coose one of 'KEP', 'MES' or 'GNV'.  The default is 'KEP'.
        kepswitch : integer, optional
            The default is 12555.
        outlines : boolean, optional
            The default is False.
        write_preproc : boolean, optional
            The default is False.
        hatching : boolean, optional
            The default if False.

        Examples
        --------
        >>> pt.kip_cont(0, -1, xres=1000, yres=1000, ixaxis='log_time_left')

        Notes
        -----
        Unfortunately, although the implicit do loops are faster, I
        cannot implement a % bar, but this example call should not take
        more than 1 minute.

        """

        tmp_cycles=self.se.cycles
        original_cyclelist = [int(tmp_cycles[i]) for i in range(len(tmp_cycles))]
        nmodels=len(self.se.cycles[modstart:modstop])
        print('nmodels           = ', nmodels)
        sparse = int(max(1,old_div(nmodels,xres)))
        cyclelist = original_cyclelist[modstart:modstop:sparse]

        age_unit=self.get('age_unit')
        if codev == 'KEP':
            oneyear = 365.*24.*3600.
        elif codev == 'MES':
            oneyear=self.get('one_year')
        elif codev == 'GNV':
            oneyear=1.e0


        m_min=ylims[0]
        m_max=ylims[1]
        if m_min==0.:
            m_min=0.
            ylims[0] = m_min
        if m_max==0.:
            m_max=float(self.se.get('mini'))
            ylims[1] = m_max


        dy = old_div((m_max-m_min),float(yres))
        y = np.arange(m_min, m_max, dy)

        Z = np.zeros([len(y),len(cyclelist)],float)

        datatype='convection'
        if 'convection_indicator' not in self.se.dcols:
            datatype='dcoeff'

        def list_to_string(array):
            string=''
            for el in array:
                string+='    '+str(el)+'    '
            string+='\n'
            return string

        if write_preproc==True:
            if os.path.exists(self.sedir+'/conv_data_preproc.txt'):
                print('Preprocessor already exists (conv_data_preproc.txt).')
                print('Please remove this file if you want to create a new preprocessor file')
                sys.exit()
            print('write_preproc=True')
            print('getting complete dataset')
            print('getting mass coordinates')
            allcycs=[int(self.se.cycles[i]) for i in range(len(self.se.cycles))]
            mass=self.se.get(allcycs,'mass')
            print('getting conv')
            if datatype=='dcoeff':
                conv=se.f.se.get(allcycs,'dcoeff')
            else:
                conv=self.se.get(allcycs,'convection_indicator')
            print('getting cycles')
            models=self.se.cycles
            print('getting ages')
            if codev=='KEP':
                deltat=self.se.get(allcycs,'deltat')
                ages=np.cumsum(deltat)
            else:
                ages=self.se.ages
            print('writing preprocessor file')
            f=open(self.sedir+'/conv_data_preproc.txt','w')
            #f=open('./conv_data_preproc.txt','w')
            f.write('FORMAT OF THIS FILE:\n')
            f.write('model_number\n')
            f.write('star_age/yr\n')
            f.write('mass/Mo\n')
            f.write('convection_indicator or dcoeff if unavailable \n')
            f.write('*************************************************************\n')

            for i in range(len(models)):
                percent = int(i*100/len(self.se.cycles))
                sys.stdout.flush()
                sys.stdout.write("\rprogress " + "...%d%%" % percent)
                f.write(str(models[i])+'\n')
                f.write(str(ages[i])+'\n')
                f.write(list_to_string(mass[i]))
                f.write(list_to_string(conv[i]))

            f.close()
            print('preprocessor file written')
            sys.exit()

        def realarray(array):
            """
            Parameters
            ----------
            array : list
                Takes array of strings

            Returns
            --------
            list
                array of floats

            """
            floatarray=np.array([float(el) for el in array])
            return floatarray

        print('getting data...')

        if os.path.exists(self.sedir+'/conv_data_preproc.txt'):
            print('... from preprocessor file ...')
            lines=open(self.sedir+'/conv_data_preproc.txt','r').readlines()
            #lines=open('./conv_data_preproc.txt','r').readlines()
            lines=lines[6:]
            mod=np.array([float(line.split()[0]) for line in lines[::4]])
            age=old_div(np.array([float(line.split()[0]) for line in lines[1::4]]),oneyear)
            # masses and mix types will be read just as an array of strings (each string is a profile)
            # to be extracted later
            masses=lines[2::4]
            conv=lines[3::4]

            tmp=[realarray(massco.split()) for massco in masses]
            masses=tmp

            tmp=[realarray(mixvec.split()) for mixvec in conv]
            conv=tmp

            print(len(age),len(mod),len(masses),len(conv))

        else:
            print('... from HDF5 data ...')
            masses=self.se.get(cyclelist,'mass')
            mod=cyclelist
            if self.deltat == None:
                print('getting deltat since not initialised')
                self.deltat = self.se.get(self.se.cycles,'deltat')
            age=old_div(np.cumsum(self.deltat), (3600.*24.*365.))

            if datatype=='convection':
                conv=self.se.get(cyclelist,'convection_indicator')
            else:
                conv=self.se.get(cyclelist,'dcoeff')

        print('calculating convection matrix... ')

        print(len(masses))
        conv_i_vec_matrix=np.zeros([len(y),len(masses)],float)

        for i in range(len(masses)):
            percent = int(i*100/len(masses))
            sys.stdout.flush()
            sys.stdout.write("\rcreating color map mix " + "...%d%%" % percent)

            mass = masses[i]
            mix = conv[i]
#            print i, len(mass),len(mix)
            # rearrange centre to surface:
            if mass[-1]<mass[0]:
                mass=np.append(mass,[0.])[::-1]
                mix=np.append(mix,[0.])[::-1]
            else:
                mass=np.append(mass[::-1],[0.])[::-1]
                mix=np.append(mix[::-1],[0.])[::-1]
            f = interpolate.interp1d(mass, mix, bounds_error=False, fill_value=0.)
            mixnew = f(y)

            for j in range(len(mixnew)):
                conv_i_vec_matrix[j,i] = mixnew[j]


        print('getting stellar mass...')
        mtot = np.array([np.max(m) for m in masses])

        #########################################################################
        # PLOT:
        fsize=20

        params = {'axes.labelsize':  fsize,
          'text.fontsize':   fsize,
          'legend.fontsize': fsize,
          'xtick.labelsize': fsize*0.8,
          'ytick.labelsize': fsize*0.8,
          'xtick.major.pad': 8,          # distance to major tick label in points
          'ytick.major.pad': 6,          # distance to major tick label in points
          'text.usetex': False}
        pl.rcParams.update(params)
        fig = pl.figure(1)
        ax = pl.axes()
        if landscape_plot == True:
            fig.set_size_inches(9,4)
            fsize=20
            pl.gcf().subplots_adjust(bottom=0.15)
            pl.gcf().subplots_adjust(right=0.85)

        ax.set_ylabel('$\mathrm{Mass }(M_\odot)$',fontsize=fsize)

        xxx=[]
        age_array = age

        if ixaxis == 'log_time_left':
        # log of time left until core collapse
            gage= np.array(age_array)
            lage=np.zeros(len(gage))
            agemin = max(old_div(abs(gage[-1]-gage[-2]),5.),1.e-10)
            for i in np.arange(len(gage)):
                if gage[-1]-gage[i]>agemin:
                    lage[i]=np.log10(gage[-1]-gage[i]+agemin)
                else :
                    lage[i]=np.log10(agemin)
            xxx = np.array([lage[int(i)] for i in cyclelist]) # np.array([lage[i] for i in cyclelist])
            print(len(xxx))
            print('plot versus time left')
            ax.set_xlabel('$\mathrm{log}_{10}(t^*/\mathrm{yr})$',fontsize=fsize)
            if xlims == [0.,0.]:
                xlims = [xxx[0],xxx[-1]]
        elif ixaxis =='model_number':
            cyclelist=mod[modstart:modstop:sparse] # np.array([int(cycle) for cycle in cyclelist])
            xxx= cyclelist
            print('plot versus model number')
            ax.set_xlabel('Model number',fontsize=fsize)
            if xlims == [0.,0.]:
                xlims = [cyclelist[0],cyclelist[-1]]
        elif ixaxis =='age':
            xxx = old_div(np.array([age_array[int(i)] for i in cyclelist]),1.e6)
            print('plot versus age')
            ax.set_xlabel('Age [Myr]',fontsize=fsize)
            if xlims == [0.,0.]:
                xlims = [xxx[0],xxx[-1]]

        #cmapMIX=mpl.colors.ListedColormap(['w','#8B8386']) # rose grey
        cmapMIX = mpl.colors.ListedColormap(['#C4C4C4']) # grey77
        #cmap = mpl.colors.ListedColormap(['w','b'])

        # some stuff for rasterizing only the contour part of the plot, for nice, but light, eps:
        class ListCollection(Collection):
            def __init__(self, collections, **kwargs):
                Collection.__init__(self, **kwargs)
                self.set_collections(collections)
            def set_collections(self, collections):
                self._collections = collections
            def get_collections(self):
                return self._collections
            @allow_rasterization
            def draw(self, renderer):
                for _c in self._collections:
                    _c.draw(renderer)

        def insert_rasterized_contour_plot(c):
            collections = c.collections
            for _c in collections:
                _c.remove()
            cc = ListCollection(collections, rasterized=True)
            ax = pl.gca()
            ax.add_artist(cc)
            return cc


        print('plotting contours')
        ax.autoscale(False)
        low_level = dcoeff_thresh - old_div(dcoeff_thresh,10)
        high_level= dcoeff_thresh + old_div(dcoeff_thresh,10)
        if datatype=='dcoeff':
            low_level = dcoeff_thresh - old_div(dcoeff_thresh,10)
            high_level= dcoeff_thresh + old_div(dcoeff_thresh,10)
            CMIX=ax.contourf(xxx,y,conv_i_vec_matrix, cmap=cmapMIX, alpha=1.0,levels=[low_level,1.e99])
        else:
            CMIX=ax.contourf(xxx,y,conv_i_vec_matrix, cmap=cmapMIX, alpha=1.0,levels=[0.5,1.5])

        if outlines == True:
            if datatype=='dcoeff':
                CMIXoutlines=ax.contour(xxx,y,conv_i_vec_matrix, cmap=cmapMIX,levels=[low_level,high_level])
            else:
                CMIX_outlines=ax.contour(xxx,y,conv_i_vec_matrix, cmap=cmapMIX,levels=[0.5,1.5])
        ax.plot(xxx,mtot,color='k')
        ax.axis([xlims[0],xlims[-1],ylims[0],ylims[1]])

        def contour_to_hatched_patches(cntrset, hatch_colors, hatch_patterns,
                               remove_contour=True):
            from itertools import cycle

            patches_list = []
            for pathcollection in cntrset.collections:
                patches_list.append([PathPatch(p1) for p1 in  pathcollection.get_paths()])
                if remove_contour:
                    pathcollection.remove()

            for patches, hc, hp  in zip(patches_list,
                                        cycle(hatch_colors), cycle(hatch_patterns)):
                for p in patches:
                    p.set_fc("none")
                    p.set_ec("k")
                    p.set_hatch(hp)
                    ax.add_patch(p)

        hatch_colors = "k"
        hatch_patterns = "/\|-+xoO.*"
        if hatching==True:
            contour_to_hatched_patches(CMIX, hatch_colors, hatch_patterns,
                    remove_contour=True)

        # uncomment or rasterizing the contour - for making nice eps figures that aren't huge files:
        #insert_rasterized_contour_plot(CMIX)
        #if outlines==True:
        #    insert_rasterized_contour_plot(CMIXoutlines)

        if outfile!='':
            if outfile[-4:] == '.png':
                fig.savefig(outfile,dpi=300)
            elif outfile[-4:] == '.eps':
                fig.savefig(outfile,format='eps')
            elif outfile[-3:] == '.ps':
                fig.savefig(outfile,format='ps')
            elif outfile[-4:] == '.pdf':
                fig.savefig(outfile,format='pdf')
            else:
                print('file format not specified as .eps or .png, saving as eps:')
                print(outfile+'.eps')
                fig.savefig(outfile+'.eps',format='eps')

        pl.show()
This function creates a Kippenhahn diagram as a contour plot of
        the se output using a convection flag or diffusion coefficient
        threshold.

        Parameters
        ----------
        modstart, modstop : integer
            First and last cycle numbers to plot.
        dcoeff_thresh : float, optional
            Diffusion coefficient threshold, above which assume
            convection is present (only for MESA, where
            convection_indicator is not there)).  The default is 1.e12.
        xres, yres : integer, optional
            x and y resolution.  The default for xres is 1000, the
            default for yres is 2000.
        xlims, ylims : list, optional
            x and y plot limits.  The default is [0., 0.].
        ixaxis : string, optional
            Choose one of 'age', 'model_number' or 'log_time_left'.  The
            default is 'log_time_left'.
        outfile : string, optional
            Name of output file including extension, which will be saved
            in 300 dpi.  The default is ''.
        landscape_plot : boolean, optional
            For a landscape plot.  The default is False.
        codev : string, optional
            Coose one of 'KEP', 'MES' or 'GNV'.  The default is 'KEP'.
        kepswitch : integer, optional
            The default is 12555.
        outlines : boolean, optional
            The default is False.
        write_preproc : boolean, optional
            The default is False.
        hatching : boolean, optional
            The default if False.

        Examples
        --------
        >>> pt.kip_cont(0, -1, xres=1000, yres=1000, ixaxis='log_time_left')

        Notes
        -----
        Unfortunately, although the implicit do loops are faster, I
        cannot implement a % bar, but this example call should not take
        more than 1 minute.