def kip_cont(self, ifig=110, modstart=0, modstop=-1,t0_model=0,
                 outfile='out.png', xlims=None, ylims=None,
                 xres=1000, yres=1000, ixaxis='model_number',
                 mix_zones=20, burn_zones=20, plot_radius=False,
                 engenPlus=True, engenMinus=False,
                 landscape_plot=False, rad_lines=False, profiles=[],
                 showfig=True, outlines=True, boundaries=True,
                 c12_boundary=False, rasterise=False, yscale='1.',
                 engenlevels=None,CBM=False,fsize=14):
        """
        This function creates a Kippenhahn plot with energy flux using
        contours.

        This plot uses mixing_regions and burning_regions written to
        your history.data or star.log. Set both variables in the
        log_columns.list file to 20 as a start.

        The output log file should then contain columns called
        "mix_type_n", "mix_qtop_n", "burn_type_n" and "burn_qtop_n".
        The number of columns (i.e. the biggest value of n) is what
        goes in the arguments as mix_zones and burn_zones.  DO NOT
        WORRY! if you do not have these columns, just leave the default
        values alone and the script should recognise that you do not
        have these columns and make the most detailed plot that is
        available to you.

        Defaults are set to get some plot, that may not look great if
        you zoom in interactively. Play with xres and yres as well as
        setting the xlims to ylims to the region you are interested in.

        Parameters
        ----------
        ifig : integer, optional
            Figure frame number.  The default is 110.
        modstart : integer, optional
            Model from which you want to plot (be careful if your
            history.data or star.log output is sparse...).  If it is
            0 then it starts from the beginning, works even if
            log_cnt > 1.
            The default is 0.
        modstop : integer, optional
            Model to which you wish to plot, -1 corresponds to end
            [if log_cnt>1, devide modstart and modstop by log_cnt,
            this needs to be improved! SJ: this should be ficed now].
            The defalut is -1.
        t0_model : integer, optional
            Model number from which to reset the time to 0. Typically,
            if modstart!=0, t0_model=modstart is a good choice, but we
            leave the choice to the user in case the time is wished to
            start from 0 at a different key point of the evolution.
            The default value is 0.
        outfile : sting, optional
            'filename + extension' where you want to save the figure.
            The defalut is "out.png".
        xlims, ylims : list, optional
            Plot limits, however these are somewhat obsolete now that
            we have modstart and modstop.  Leaving them as 0. is
            probably no slower, and you can always zoom in afterwards
            in mpl.  ylims is important for well resolved thermal
            pulse etc plots; it's best to get the upper and lower limits
            of he-intershell using s.kippenhahn_CO(1,'model') first.
            The default is [0., 0.].
        xres, yres : integer, optional
            plot resolution. Needless to say that increasing these
            values will yield a nicer plot with some slow-down in
            plotting time.  You will most commonly change xres.  For a
            prelim plot, try xres~200, then bump it up to anywhere from
            1000-10000 for real nicely resolved, publication quality
            plots.  The default is 1000.
        ixaxis : string, optional
            Either 'log_time_left', 'age', or 'model_number'.  The
            default is "model_number".
        mix_zones, burn_zones : integer, optional
            As described above, if you have more detailed output about
            your convection and energy generation boundaries in columns
            mix_type_n, mix_qtop_n, burn_type_n and burn_qtop_n, you
            need to specify the total number of columns for mixing zones
            and burning zones that you have.  Can't work this out from
            your history.data or star.log file?  Check the
            history_columns.list that you used, it'll be the number
            after "mixing regions" and "burning regions".  Can't see
            these columns?  leave it and 2 conv zones and 2 burn zones
            will be drawn using other data that you certainly should
            have in your history.data or star.log file.  The defalut for
            both is 20.
        plot_radius : boolean, optional
            Whether on a second y-axis you want to plot the radius of
            the surface and the he-free core.  The default is False.
        engenPlus : boolean
            Plot energy generation contours for eps_nuc>0.  The default
            is True.
        endgenMinus : boolean, optional
            Plot energy generation contours for eos_nuc<0.  The default
            is True.
        landscape_plot : boolean, optionla
            The default is False.
        rad_lines : boolean, optional
            The deafault is False.
        profiles : list, optional
            The default is [].
        showfig : boolean, optional
            The default is True.
        outlines : boolean, optional
            Whether or not to plot outlines of conv zones in darker
            colour.
        boundaries : boolean, optional
            Whether or not to plot H-, He- and C-free boundaries.
        c12_boundary : boolean, optional
            The default is False.
        rasterise : boolean, optional
            Whether or not to rasterise the contour regions to make
            smaller vector graphics figures.  The default is False.
        yscale : string, optional
            Re-scale the y-axis by this amount
        engenlevels : list
            Give cusstom levels to the engenPlus contour. If None,
            the levels are chosen automatically.
            The default is None.
        CBM : boolean, optional
            plot contours for where CBM is active?

        fsize : integer
            font size for labels

        Notes
        -----
        The parameter xlims is depricated.

        """
        if ylims is None:
            ylims=[0.,0.]
        if xlims is None:
            xlims=[0.,0.]
        
        # Find correct modstart and modstop:
        mod=np.array([int(i) for i in self.get('model_number')])
        mod1=np.abs(mod-modstart).argmin()
        mod2=np.abs(mod-modstop).argmin()
        if modstart != 0 : modstart=mod1
        if modstop != -1 : modstop=mod2

        xxyy=[self.get('star_age')[modstart:modstop],self.get('star_age')[modstart:modstop]]
        mup = max(float(self.get('star_mass')[0])*1.02,1.0)
        nmodels=len(self.get('model_number')[modstart:modstop])

        if ylims == [0.,0.]:
            mup   = max(float(self.get('star_mass')[0])*1.02,1.0)
            mDOWN = 0.
        else:
            mup = ylims[1]
            mDOWN = ylims[0]

        # y-axis resolution
        ny=yres
        #dy=mup/float(ny)
        dy = old_div((mup-mDOWN),float(ny))

        # x-axis resolution
        maxpoints=xres
        dx=int(max(1,old_div(nmodels,maxpoints)))

        #y = np.arange(0., mup, dy)
        y = np.arange(mDOWN, mup, dy)
        x = np.arange(0, nmodels, dx)
        Msol=1.98892E+33

        engenstyle = 'full'

        B1=np.zeros([len(y),len(x)],float)
        B2=np.zeros([len(y),len(x)],float)
        try:
            self.get('burn_qtop_1')
        except:
            engenstyle = 'twozone'
        if engenstyle == 'full' and (engenPlus == True or engenMinus == True):
            ulimit_array = np.array([self.get('burn_qtop_'+str(j))[modstart:modstop:dx]*\
                           self.get('star_mass')[modstart:modstop:dx] for j in range(1,burn_zones+1)])
            #ulimit_array = np.around(ulimit_array,decimals=len(str(dy))-2)
            llimit_array = np.delete(ulimit_array,-1,0)
            llimit_array = np.insert(ulimit_array,0,0.,0)
            #llimit_array = np.around(llimit_array,decimals=len(str(dy))-2)
            btype_array = np.array([self.get('burn_type_'+str(j))[modstart:modstop:dx] for j in range(1,burn_zones+1)])
            old_percent = 0
            for i in range(len(x)):
                # writing status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map burn " + "...%d%%" % percent)
                    old_percent = percent

                for j in range(burn_zones):
                    if btype_array[j,i] > 0. and abs(btype_array[j,i]) < 99.:
                        B1[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 10.0**(btype_array[j,i])
                    elif btype_array[j,i] < 0. and abs(btype_array[j,i]) < 99.:
                        B2[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 10.0**(abs(btype_array[j,i]))
            print(' \n')

        if engenstyle == 'twozone' and (engenPlus == True or engenMinus == True):
            V=np.zeros([len(y),len(x)],float)
            old_percent = 0
            for i in range(len(x)):
                # writing status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map1 " + "...%d%%" % percent)
                    old_percent = percent
                llimitl1=old_div(self.get('epsnuc_M_1')[modstart:modstop][i*dx],Msol)
                ulimitl1=old_div(self.get('epsnuc_M_4')[modstart:modstop][i*dx],Msol)
                llimitl2=old_div(self.get('epsnuc_M_5')[modstart:modstop][i*dx],Msol)
                ulimitl2=old_div(self.get('epsnuc_M_8')[modstart:modstop][i*dx],Msol)
                llimith1=old_div(self.get('epsnuc_M_2')[modstart:modstop][i*dx],Msol)
                ulimith1=old_div(self.get('epsnuc_M_3')[modstart:modstop][i*dx],Msol)
                llimith2=old_div(self.get('epsnuc_M_6')[modstart:modstop][i*dx],Msol)
                ulimith2=old_div(self.get('epsnuc_M_7')[modstart:modstop][i*dx],Msol)
                # lower thresh first, then upper thresh:
                if llimitl1!=ulimitl1:
                    for k in range(ny):
                        if llimitl1<=y[k] and ulimitl1>y[k]:
                            V[k,i]=10.
                if llimitl2!=ulimitl2:
                    for k in range(ny):
                        if llimitl2<=y[k] and ulimitl2>y[k]:
                            V[k,i]=10.
                if llimith1!=ulimith1:
                    for k in range(ny):
                        if llimith1<=y[k] and ulimith1>y[k]:
                            V[k,i]=30.
                if llimith2!=ulimith2:
                    for k in range(ny):
                        if llimith2<=y[k] and ulimith2>y[k]:
                            V[k,i]=30.
            print(' \n')

        mixstyle = 'full'
        try:
            self.get('mix_qtop_1')
        except:
            mixstyle = 'twozone'
        if mixstyle == 'full':
            old_percent = 0
            Z=np.zeros([len(y),len(x)],float)
            if CBM:
                Zcbm=np.zeros([len(y),len(x)],float)
            ulimit_array = np.array([self.get('mix_qtop_'+str(j))[modstart:modstop:dx]*self.get('star_mass')[modstart:modstop:dx] for j in range(1,mix_zones+1)])
            llimit_array = np.delete(ulimit_array,-1,0)
            llimit_array = np.insert(ulimit_array,0,0.,0)
            mtype_array = np.array([self.get('mix_type_'+str(j))[modstart:modstop:dx] for j in range(1,mix_zones+1)])
            for i in range(len(x)):
                # writing status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map mix " + "...%d%%" % percent)
                    old_percent = percent
                for j in range(mix_zones):
                    if mtype_array[j,i] == 1.:
                        Z[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 1.
                    if CBM:
                        if mtype_array[j,i] == 2.:
                            Zcbm[(np.abs(y-llimit_array[j][i])).argmin():(np.abs(y-ulimit_array[j][i])).argmin()+1,i] = 1.
            print(' \n')

        if mixstyle == 'twozone':
            Z=np.zeros([len(y),len(x)],float)
            old_percent = 0
            for i in range(len(x)):
            # writing reading status
                # writing status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map mix " + "...%d%%" % percent)
                    old_percent = percent

                ulimit=self.get('conv_mx1_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                llimit=self.get('conv_mx1_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                if llimit!=ulimit:
                    for k in range(ny):
                        if llimit<=y[k] and ulimit>y[k]:
                            Z[k,i]=1.
                ulimit=self.get('conv_mx2_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                llimit=self.get('conv_mx2_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                if llimit!=ulimit:
                    for k in range(ny):
                        if llimit<=y[k] and ulimit>y[k]:
                            Z[k,i]=1.
            print(' \n')

        if rad_lines == True:
            masses = np.arange(0.1,1.5,0.1)
            rads=[[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
            modno=[]
            for i in range(len(profiles)):
                p=mesa_profile('./LOGS',profiles[i])
                modno.append(p.header_attr['model_number'])
                for j in range(len(masses)):
                    idx=np.abs(p.get('mass')-masses[j]).argmin()
                    rads[j].append(p.get('radius')[idx])

        print('engenstyle was ', engenstyle)
        print('mixstyle was ', mixstyle)
        print('\n finished preparing color map')

        ########################################################################
        #----------------------------------plot--------------------------------#
        fig = pyl.figure(ifig)
        #fsize=20
        if landscape_plot == True:
            fig.set_size_inches(9,4)
            pl.gcf().subplots_adjust(bottom=0.2)
            pl.gcf().subplots_adjust(right=0.85)

        params = {'axes.labelsize':  fsize,
          'axes.labelsize': fsize,
          'font.size': fsize,
          'legend.fontsize': fsize,
          'xtick.labelsize': fsize,
          'ytick.labelsize': fsize,
                  'text.usetex': False}
        pyl.rcParams.update(params)

        #ax=pl.axes([0.1,0.1,0.9,0.8])

        #fig=pl.figure()
        ax=pl.axes()

        if ixaxis == 'log_time_left':
        # log of time left until core collapse
            gage= self.get('star_age')
            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 = lage[modstart:modstop]
            print('plot versus time left')
            ax.set_xlabel('$ \\log_{10}(t-t_\mathrm{end})\ /\  \mathrm{[yr]}$',fontsize=fsize)
            if xlims[1] == 0.:
                xlims = [xxx[0],xxx[-1]]
        elif ixaxis =='model_number':
            xxx= self.get('model_number')[modstart:modstop]
            print('plot versus model number')
            ax.set_xlabel('Model number',fontsize=fsize)
            if xlims[1] == 0.:
                xlims = [self.get('model_number')[modstart],self.get('model_number')[modstop]]
        elif ixaxis =='age':
            if t0_model != 0:
                t0_mod=np.abs(mod-t0_model).argmin()
                xxx= self.get('star_age')[modstart:modstop] - self.get('star_age')[t0_mod]
                print('plot versus age')
                ax.set_xlabel('t - %.5e / [yr]' %self.get('star_age')[modstart],fontsize=fsize)
            else:
                xxx= old_div(self.get('star_age')[modstart:modstop],1.e6)
                ax.set_xlabel('t [Myr]',fontsize=fsize)
            if xlims[1] == 0.:
                xlims = [xxx[0],xxx[-1]]

        ax.set_ylabel('$\mathrm{enclosed\ mass\ /\ [M_\odot]}$',fontsize=fsize)

        # 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

        cmapMIX = matplotlib.colors.ListedColormap(['w','#8B8386']) # rose grey
        if CBM:
            cmapCBM = matplotlib.colors.ListedColormap(['w','g']) # green
        cmapB1  = pyl.cm.get_cmap('Blues')
        cmapB2  = pl.cm.get_cmap('Reds')

        ylims1=[0.,0.]
        ylims1[0]=ylims[0]
        ylims1[1]=ylims[1]
        if ylims == [0.,0.]:
            ylims[0] = 0.
            ylims[1] = mup
            #print("Setting ylims[1] to mup="+str(mup))
        if ylims[0] != 0.:
            ylab='$(\mathrm{Mass }$ - '+str(ylims[0])
            if yscale!='1.':
                ylab+=') / '+yscale+' $M_\odot$'
            else:
                ylab+=') / $M_\odot$'
            ax.set_ylabel(ylab)
            y = y - ylims[0]
            y = y*float(yscale) # SJONES tweak
            ylims[0] = y[0]
            ylims[1] = y[-1]

        print('plotting contours')
        CMIX    = ax.contourf(xxx[::dx],y,Z, cmap=cmapMIX,alpha=0.6,levels=[0.5,1.5])
        #CMIX    = ax.pcolor(xxx[::dx],y,Z, cmap=cmapMIX,alpha=0.6,vmin=0.5,vmax=1.5)
        if rasterise==True:
            insert_rasterized_contour_plot(CMIX)
        if outlines == True:
            CMIX_outlines    = ax.contour(xxx[::dx],y,Z, cmap=cmapMIX)
            if rasterise==True:
                insert_rasterized_contour_plot(CMIX_outlines)

        if CBM:
            CCBM    = ax.contourf(xxx[::dx],y,Zcbm, cmap=cmapCBM,alpha=0.6,levels=[0.5,1.5])
            if rasterise==True:
                insert_rasterized_contour_plot(CCBM)
            if outlines == True:
                CCBM_outlines    = ax.contour(xxx[::dx],y,Zcbm, cmap=cmapCBM)
                if rasterise==True:
                    insert_rasterized_contour_plot(CCBM_outlines)

        if engenstyle == 'full' and engenPlus == True:
            if engenlevels!= None:
                CBURN1  = ax.contourf(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.5,\
                        locator=matplotlib.ticker.LogLocator(),levels=engenlevels)
                if outlines:
                    CB1_outlines  = ax.contour(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.7, \
                        locator=matplotlib.ticker.LogLocator(),levels=engenlevels)
            else:
                CBURN1  = ax.contourf(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.5, \
                                      locator=matplotlib.ticker.LogLocator())
                if outlines:
                    CB1_outlines  = ax.contour(xxx[::dx],y,B1, cmap=cmapB1, alpha=0.7, \
                                               locator=matplotlib.ticker.LogLocator())
            CBARBURN1 = pyl.colorbar(CBURN1)
            CBARBURN1.set_label('$|\epsilon_\mathrm{nuc}-\epsilon_{\\nu}| \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$',fontsize=fsize)
            if rasterise==True:
                insert_rasterized_contour_plot(CBURN1)
                if outlines:
                    insert_rasterized_contour_plot(CB1_outlines)

        if engenstyle == 'full' and engenMinus == True:
            CBURN2  = ax.contourf(xxx[::dx],y,B2, cmap=cmapB2, alpha=0.5, locator=matplotlib.ticker.LogLocator())
            if outlines:
                CBURN2_outlines  = ax.contour(xxx[::dx],y,B2, cmap=cmapB2, alpha=0.7, locator=matplotlib.ticker.LogLocator())
            CBARBURN2 = pl.colorbar(CBURN2)
            if engenPlus == False:
                CBARBURN2.set_label('$|\epsilon_\mathrm{nuc}-\epsilon_{\\nu}| \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$',fontsize=fsize)
            if rasterise==True:
                insert_rasterized_contour_plot(CBURN2)
                if outlines:
                    insert_rasterized_contour_plot(CB2_outlines)

        if engenstyle == 'twozone' and (engenPlus == True or engenMinus == True):
            ax.contourf(xxx[::dx],y,V, cmap=cmapB1, alpha=0.5)

        print('plotting patches')
        mtot=self.get('star_mass')[modstart:modstop][::dx]
        mtot1=(mtot-ylims1[0])*float(yscale)
        ax.plot(xxx[::dx],mtot1,'k-')

        if boundaries == True:
            print('plotting abund boundaries')
            try:
                bound=self.get('h1_boundary_mass')[modstart:modstop]
                bound1=(bound-ylims1[0])*float(yscale)
                ax.plot(xxx,bound1,label='H boundary',linestyle='-')

                bound=self.get('he4_boundary_mass')[modstart:modstop]
                bound1=(bound-ylims1[0])*float(yscale)
                ax.plot(xxx,bound1,label='He boundary',linestyle='--')

                bound=self.get('c12_boundary_mass')[modstart:modstop]
                bound1=(bound-ylims1[0])*float(yscale)
                ax.plot(xxx,bound1,label='C boundary',linestyle='-.')

            except:
                try:
                    bound=self.get('he_core_mass')[modstart:modstop]
                    bound1=(bound-ylims1[0])*float(yscale)
                    ax.plot(xxx,bound1,label='H boundary',linestyle='-')

                    bound=self.get('c_core_mass')[modstart:modstop]-ylims[0]
                    bound1=(bound-ylims1[0])*float(yscale)
                    ax.plot(xxx,bound1,label='He boundary',linestyle='--')

                    bound=self.get('o_core_mass')[modstart:modstop]-ylims[0]
                    bound1=(bound-ylims1[0])*float(yscale)
                    ax.plot(xxx,bound1,label='C boundary',linestyle='-.')

                    bound=self.get('si_core_mass')[modstart:modstop]-ylims[0]
                    bound1=(bound-ylims1[0])*float(yscale)
                    ax.plot(xxx,bound1,label='C boundary',linestyle='-.')

                    bound=self.get('fe_core_mass')[modstart:modstop]-ylims[0]
                    bound1=(bound-ylims1[0])*float(yscale)
                    ax.plot(xxx,bound1,label='C boundary',linestyle='-.')

                except:
#                    print 'problem to plot boundaries for this plot'
                    pass

        ax.axis([xlims[0],xlims[1],ylims[0],ylims[1]])

        if plot_radius == True:
            ax2=pyl.twinx()
            ax2.plot(xxx,np.log10(self.get('he4_boundary_radius')[modstart:modstop]),label='He boundary radius',color='k',linewidth=1.,linestyle='-.')
            ax2.plot(xxx,self.get('log_R')[modstart:modstop],label='radius',color='k',linewidth=1.,linestyle='-.')
            ax2.set_ylabel('log(radius)')
        if rad_lines == True:
            ax2=pyl.twinx()
            for i in range(len(masses)):
                ax2.plot(modno,np.log10(rads[i]),color='k')

        if outfile[-3:]=='png':
            fig.savefig(outfile,dpi=300)
        elif outfile[-3:]=='eps':
            fig.savefig(outfile,format='eps')
        elif outfile[-3:]=='pdf':
            fig.savefig(outfile,format='pdf')
        if showfig == True:
            pyl.show()