def _kip_vline(self, modstart, modstop, sparse, outfile,
                  xlims=[0.,0.], ylims=[0.,0.], ixaxis='log_time_left',
                  mix_zones=5, burn_zones=50):
        """
        *** DEPRECIATED and hence UNSUPPORTED ***
        This function creates a Kippenhahn plot with energy flux using
        vertical lines, better thermal pulse resolution.

        For a more comprehensive plot, your history.data or star.log
        file should contain columns called "mix_type_n", "mix_qtop_n",
        "burn_type_n" and "burn_qtop_n".  The number of columns
        (i.e. the bbiggest 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.

        Parameters
        ----------
        modstart : integer
            Model from which you want to plot (be careful if your
            history.data or star.log output is sparse...).
        modstop : integer
            Model to which you wish to plot.
        sparse : integer
            x-axis sparsity.
        outfile : string
            'filename + extension' where you want to save the figure.
        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.  The default is [0., 0.,].
        ixaxis : string, optional
            Either 'log_time_left', 'age', or 'model_number'.  The
            default "log_time_left".
        mix_zones, burn_zones : integer
            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 default for
            mix_zones is 5, the defalut for burn_zones is 50.

        """


        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])

        Msol=1.98892E+33

        engenstyle = 'full'

        dx = sparse
        x = np.arange(0, nmodels, dx)

        btypemax = 20
        btypemin = -20
        btypealpha=0.

        ########################################################################
        #----------------------------------plot--------------------------------#
        fig = pl.figure()
#       fig.set_size_inches(16,9)
        fsize=15
        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('$\mathrm{log}_{10}(t^*) \, \mathrm{(yr)}$',fontsize=fsize)
        elif ixaxis =='model_number':
            xxx= self.get('model_number')[modstart:modstop]
            print('plot versus model number')
            ax.set_xlabel('Model number',fontsize=fsize)
        elif ixaxis =='age':
            xxx= old_div(self.get('star_age')[modstart:modstop],1.e6)
            print('plot versus age')
            ax.set_xlabel('Age [Myr]',fontsize=fsize)
        else:
            print('ixaxis must be one of: log_time_left, age or model_number')
            sys.exit()

        if xlims == [0.,0.]:
            xlims[0] = xxx[0]
            xlims[1] = xxx[-1]
        if ylims == [0.,0.]:
            ylims[0] = 0.
            ylims[1] = mup


        print('plotting patches')
        ax.plot(xxx[::dx],self.get('star_mass')[modstart:modstop][::dx],'k-')

        print('plotting abund boundaries')
        ax.plot(xxx,self.get('h1_boundary_mass')[modstart:modstop],label='H boundary')
        ax.plot(xxx,self.get('he4_boundary_mass')[modstart:modstop],label='He boundary')
#       ax.plot(xxx,self.get('c12_boundary_mass')[modstart:modstop],label='C boundary')

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

        ax.set_ylabel('Mass [M$_\odot$]')

        ########################################################################

        try:
            self.get('burn_qtop_1')
        except:
            engenstyle = 'twozone'
        old_percent = 0
        if engenstyle == 'full':
            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

                for j in range(1,burn_zones+1):
                    ulimit=self.get('burn_qtop_'+str(j))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                    if j==1:
                        llimit=0.0
                    else:
                        llimit=self.get('burn_qtop_'+str(j-1))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                    btype=float(self.get('burn_type_'+str(j))[modstart:modstop][i*dx])
                    if llimit!=ulimit:
                        if btype>0.:
                            #btypealpha = btype/btypemax
                            #ax.axvline(xxx[i*dx],ymin=(llimit-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimit-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=btypealpha)
                            pass
                        if btype<0.:
                            #btypealpha = (btype/btypemin)/5
                            #ax.axvline(xxx[i*dx],ymin=(llimit-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimit-ylims[0])/(ylims[1]-ylims[0]),color='r',alpha=btypealpha)
                            pass

            print(' \n')

        old_percent = 0

        if engenstyle == 'twozone':
            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:
                    #ax.axvline(xxx[i*dx],ymin=(llimitl1-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimitl1-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=1.)
                #if llimitl2!=ulimitl2:
                    #ax.axvline(xxx[i*dx],ymin=(llimitl2-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimitl2-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=1.)
                #if llimith1!=ulimith1:
                    #ax.axvline(xxx[i*dx],ymin=(llimith1-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimith1-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=4.)
                #if llimith2!=ulimith2:
                    #ax.axvline(xxx[i*dx],ymin=(llimith2-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimith2-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=4.)

            print(' \n')

        mixstyle = 'full'
        try:
            self.get('mix_qtop_1')
        except:
            mixstyle = 'twozone'

        old_percent = 0

        if mixstyle == 'full':
            for i in range(len(x)):
            # writing reading status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map2 " + "...%d%%" % percent)
                    old_percent = percent

                for j in range(1,mix_zones+1):
                    ulimit=self.get('mix_qtop_'+str(j))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                    if j==1:
                        llimit=0.0
                    else:
                        llimit=self.get('mix_qtop_'+str(j-1))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                    mtype=self.get('mix_type_'+str(j))[modstart:modstop][i*dx]
                    if llimit!=ulimit:
                        if mtype == 1:
                            ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=3., linewidth=.5)

            print(' \n')

        old_percent = 0

        if mixstyle == 'twozone':
            for i in range(len(x)):
            # writing reading status
                percent = int(i*100/(len(x) - 1))
                if percent >= old_percent + 5:
                    sys.stdout.flush()
                    sys.stdout.write("\r creating color map2 " + "...%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:
                    ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=5.,linewidth=.5)
                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:
                    ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=3.,linewidth=.5)

            print(' \n')

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

        #fig.savefig(outfile)
        pl.show()
def _kip_vline ( self , modstart , modstop , sparse , outfile , xlims = [ 0. , 0. ] , ylims = [ 0. , 0. ] , ixaxis = 'log_time_left' , mix_zones = 5 , burn_zones = 50 ) : 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 ] ) Msol = 1.98892E+33 engenstyle = 'full' dx = sparse x = np . arange ( 0 , nmodels , dx ) btypemax = 20 btypemin = - 20 btypealpha = 0. fig = pl . figure ( ) fsize = 15 ax = pl . axes ( ) if ixaxis == 'log_time_left' : 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 ( '$\mathrm{log}_{10}(t^*) \, \mathrm{(yr)}$' , fontsize = fsize ) elif ixaxis == 'model_number' : xxx = self . get ( 'model_number' ) [ modstart : modstop ] print ( 'plot versus model number' ) ax . set_xlabel ( 'Model number' , fontsize = fsize ) elif ixaxis == 'age' : xxx = old_div ( self . get ( 'star_age' ) [ modstart : modstop ] , 1.e6 ) print ( 'plot versus age' ) ax . set_xlabel ( 'Age [Myr]' , fontsize = fsize ) else : print ( 'ixaxis must be one of: log_time_left, age or model_number' ) sys . exit ( ) if xlims == [ 0. , 0. ] : xlims [ 0 ] = xxx [ 0 ] xlims [ 1 ] = xxx [ - 1 ] if ylims == [ 0. , 0. ] : ylims [ 0 ] = 0. ylims [ 1 ] = mup print ( 'plotting patches' ) ax . plot ( xxx [ : : dx ] , self . get ( 'star_mass' ) [ modstart : modstop ] [ : : dx ] , 'k-' ) print ( 'plotting abund boundaries' ) ax . plot ( xxx , self . get ( 'h1_boundary_mass' ) [ modstart : modstop ] , label = 'H boundary' ) ax . plot ( xxx , self . get ( 'he4_boundary_mass' ) [ modstart : modstop ] , label = 'He boundary' ) ax . axis ( [ xlims [ 0 ] , xlims [ 1 ] , ylims [ 0 ] , ylims [ 1 ] ] ) ax . set_ylabel ( 'Mass [M$_\odot$]' ) try : self . get ( 'burn_qtop_1' ) except : engenstyle = 'twozone' old_percent = 0 if engenstyle == 'full' : for i in range ( len ( x ) ) : 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 for j in range ( 1 , burn_zones + 1 ) : ulimit = self . get ( 'burn_qtop_' + str ( j ) ) [ modstart : modstop ] [ i * dx ] * self . get ( 'star_mass' ) [ modstart : modstop ] [ i * dx ] if j == 1 : llimit = 0.0 else : llimit = self . get ( 'burn_qtop_' + str ( j - 1 ) ) [ modstart : modstop ] [ i * dx ] * self . get ( 'star_mass' ) [ modstart : modstop ] [ i * dx ] btype = float ( self . get ( 'burn_type_' + str ( j ) ) [ modstart : modstop ] [ i * dx ] ) if llimit != ulimit : if btype > 0. : pass if btype < 0. : pass print ( ' \n' ) old_percent = 0 if engenstyle == 'twozone' : for i in range ( len ( x ) ) : 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 ) print ( ' \n' ) mixstyle = 'full' try : self . get ( 'mix_qtop_1' ) except : mixstyle = 'twozone' old_percent = 0 if mixstyle == 'full' : for i in range ( len ( x ) ) : percent = int ( i * 100 / ( len ( x ) - 1 ) ) if percent >= old_percent + 5 : sys . stdout . flush ( ) sys . stdout . write ( "\r creating color map2 " + "...%d%%" % percent ) old_percent = percent for j in range ( 1 , mix_zones + 1 ) : ulimit = self . get ( 'mix_qtop_' + str ( j ) ) [ modstart : modstop ] [ i * dx ] * self . get ( 'star_mass' ) [ modstart : modstop ] [ i * dx ] if j == 1 : llimit = 0.0 else : llimit = self . get ( 'mix_qtop_' + str ( j - 1 ) ) [ modstart : modstop ] [ i * dx ] * self . get ( 'star_mass' ) [ modstart : modstop ] [ i * dx ] mtype = self . get ( 'mix_type_' + str ( j ) ) [ modstart : modstop ] [ i * dx ] if llimit != ulimit : if mtype == 1 : ax . axvline ( xxx [ i * dx ] , ymin = old_div ( ( llimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , ymax = old_div ( ( ulimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , color = 'k' , alpha = 3. , linewidth = .5 ) print ( ' \n' ) old_percent = 0 if mixstyle == 'twozone' : for i in range ( len ( x ) ) : percent = int ( i * 100 / ( len ( x ) - 1 ) ) if percent >= old_percent + 5 : sys . stdout . flush ( ) sys . stdout . write ( "\r creating color map2 " + "...%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 : ax . axvline ( xxx [ i * dx ] , ymin = old_div ( ( llimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , ymax = old_div ( ( ulimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , color = 'k' , alpha = 5. , linewidth = .5 ) 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 : ax . axvline ( xxx [ i * dx ] , ymin = old_div ( ( llimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , ymax = old_div ( ( ulimit - ylims [ 0 ] ) , ( ylims [ 1 ] - ylims [ 0 ] ) ) , color = 'k' , alpha = 3. , linewidth = .5 ) print ( ' \n' ) print ( 'engenstyle was ' , engenstyle ) print ( 'mixstyle was ' , mixstyle ) print ( '\n finished preparing color map' ) pl . show ( )
*** DEPRECIATED and hence UNSUPPORTED ***
        This function creates a Kippenhahn plot with energy flux using
        vertical lines, better thermal pulse resolution.

        For a more comprehensive plot, your history.data or star.log
        file should contain columns called "mix_type_n", "mix_qtop_n",
        "burn_type_n" and "burn_qtop_n".  The number of columns
        (i.e. the bbiggest 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.

        Parameters
        ----------
        modstart : integer
            Model from which you want to plot (be careful if your
            history.data or star.log output is sparse...).
        modstop : integer
            Model to which you wish to plot.
        sparse : integer
            x-axis sparsity.
        outfile : string
            'filename + extension' where you want to save the figure.
        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.  The default is [0., 0.,].
        ixaxis : string, optional
            Either 'log_time_left', 'age', or 'model_number'.  The
            default "log_time_left".
        mix_zones, burn_zones : integer
            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 default for
            mix_zones is 5, the defalut for burn_zones is 50.