docs: def set_plot_tcrhoc(self,linestyle=['-'],burn_limit=0.997,marker=['o'],markevery=500,end_model=[-1],deg_line=True): 
  		'''
  			Plots HRDs
  			end_model - array, control how far in models a run is plottet, if -1 till end
  			symbs_1  - set symbols of runs
  		'''
  		if len(linestyle)==0:
  			linestyle=200*['-']
  
                  for i in range(len(self.runs_H5_surf)):
                          m1p65_last=se(self.runs_H5_out[i])
  			t1_model=-1
  			if end_model[i] != -1:
  				t1_model=end_model[i]
  
  			model=m1p65_last.se.cycles
  			model_list=[]
  			for k in range(0,len(model),5):
  				model_list.append(model[k])
                          rho1=m1p65_last.get(model_list,'rho')   #[:(t1_model-t0_model)]
                          T1=m1p65_last.get(model_list,'temperature')#[:(t1_model-t0_model)]
  			T_unit=m1p65_last.get('temperature_unit')
  			mass=m1p65_last.get(model_list,'mass')
  			mini=m1p65_last.get('mini')
  			zini=m1p65_last.get('zini')
  			#info=m1p65_last.burnstage_upgrade()
  			#burn_info=[]
  			#burn_cycle=[]
  			#for k in range(len(info[0])):
  			#    if 'start' in info[3][k]:
  			#	burn_info.append( info[3][k])
  			#	burn_cycle.append(info[0][k])
  			#print burn_info
  			#print burn_cycle
  			'''
  			#H
  			#Maybe use get_elemental_abunds(), too slow if individually taken
  			h=m1p65_last.get(model_list,'H-1')
  			h_ini=mini*h[0][-1]
  			h_limit=h_ini*burn_limit
  			print h_limit
  			#He
  			he=m1p65_last.get(model_list,'He-4')
  			he+=m1p65_last.get(model_list,'He-3')
  			he_ini=mini*he[0][-1]
  			he_limit=he_ini*burn_limit
  			print he_limit
  			if mini>5:
  				c=m1p65_last.get(model_list,'C-12')
  				c_ini=mini*c[0][-1]
  				c_limit=c_ini*burn_limit
  			if mini>10:
  				ne=m1p65_last.get(model_list,'Ne-20')
  				ne_ini=	mini*ne[0][-1]
  				ne_limit=ne_ini*burn_limit
  				o=m1p65_last.get(model_list,'O-16')
  				o_ini=mini*o[0][-1]
  				o_limit=o_ini*burn_limit
  				si=m1p65_last.get(model_list,'Si-28')
  				si_ini=mini*si[0][-1]
  				si_limit=si_ini*burn_limit
  			print '#############################3'
  			print h_ini,he_ini
  			rho=[]
  			T=[]
  			he_depl_idx=-1
  			h_depl_idx=-1
  			c_depl_idx=-1
  			ne_depl_idx=-1
  			o_depl_idx=-1
  			si_depl_idx=-1
  			'''
  			rho=[]
  			T=[]
  			for k in range(len(model_list)):
  				rho_center=rho1[k][0]
  				T_center=T1[k][0]
  				rho.append(rho_center)
  				T.append(T_center)
  				'''
  				mass11=np.array([0]+list(mass[k]))
  				delta_mass=mass11[1:]-mass11[:-1]
  				#for the low-mass + massive AGB
  				if (sum(np.array(he[k])*np.array(delta_mass))<he_limit) and  (he_depl_idx==-1):
  					he_depl_idx=k			
                                  if (sum(np.array(h[k])*np.array(delta_mass))<h_limit) and  (h_depl_idx==-1):
                                          h_depl_idx=k
  				#for the SAGB + massive
  				if mini>5:
  					if (sum(np.array(c[k])*np.array(delta_mass))<c_limit) and  (c_depl_idx==-1):
                                          	c_depl_idx=k
  				#for the massive stars
  				if mini>10:
  					if (sum(np.array(ne[k])*np.array(delta_mass))<ne_limit) and  (ne_depl_idx==-1):
                                          	ne_depl_idx=k
  					if (sum(np.array(o[k])*np.array(delta_mass))<o_limit) and  (o_depl_idx==-1):
                                          	o_depl_idx=k
  					if (sum(np.array(si[k])*np.array(delta_mass))<si_limit) and  (si_depl_idx==-1):
                                          	si_depl_idx=k
  				'''
  
  			T=np.log10(np.array(T)*T_unit)
  			#T_degeneracy=np.log10(np.array(rho)**(2./3.))
  			rho=np.log10(np.array(rho))
                          figure(1)
                          #plot(logTeff,logL,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
  			#pl.plot(rho,T,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
  			label=str(mini)+'$M_{\odot}$, Z='+str(zini)
  			plt.plot(rho,T,label=label,linestyle=linestyle[i],marker=marker[i],markevery=markevery)
  			'''
  			#plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
  			print burn_info
  			print burn_cycle
  			burn_cycle=np.array(burn_cycle)/100.
  			print 'changed cycles:', burn_cycle
  			print len(rho)
  			#For different burning stages:
  			if 'H_start' in burn_info:
  				plt.plot(rho[burn_cycle[0]],T[burn_cycle[0]],marker='o',color='b')
  			if 'He_start' in burn_info:
  				plt.plot(rho[burn_cycle[1]],T[burn_cycle[1]],marker='o',color='r')
  			if 'C_start' in burn_info:
  				plt.plot(rho[burn_cycle[2]],T[burn_cycle[2]],marker='o',color='g')
  			if 'Ne_start' in burn_info:
  				plt.plot(rho[burn_cycle[3]],T[burn_cycle[3]],marker='D',color='b')
  			if 'O_start' in burn_info:
  				plt.plot(rho[burn_cycle[4]],T[burn_cycle[4]],marker='D',color='r')
  			if 'Si_start' in burn_info:
  				plt.plot(rho[burn_cycle[5]],T[burn_cycle[5]],marker='D',color='g')
  			'''
  		
  			#ax = plt.gca()
  			#ax.invert_xaxis()
  			plt.rcParams.update({'font.size': 16})
  			plt.rc('xtick', labelsize=16)
  			plt.rc('ytick', labelsize=16)
  			legend(loc=4)
  		        plt.xlabel('log $\\rho_{\\rm c}$',fontsize=18)
          		plt.ylabel('log $T_{\\rm c}$',fontsize=18)
  			'''
  			#plt.gca().invert_xaxis()	
  			figure(0)
  			#pl.plot(rho,T,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
                          plt.plot(rho,T,label=self.extra_label[i],linestyle=linestyle[i])
  			plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
  			#For different burning stages:
  			if not h_depl_idx ==-1:
  				plt.plot(rho[h_depl_idx],T[h_depl_idx],marker='o',color='b')
  			if not he_depl_idx ==-1:
  				plt.plot(rho[he_depl_idx],T[he_depl_idx],marker='o',color='r')
  			if not c_depl_idx ==-1:
  				plt.plot(rho[c_depl_idx],T[c_depl_idx],marker='o',color='g')
  			if not ne_depl_idx ==-1:
  				plt.plot(rho[ne_depl_idx],T[ne_depl_idx],marker='D',color='b')
  			if not o_depl_idx ==-1:
  				plt.plot(rho[o_depl_idx],T[o_depl_idx],marker='D',color='r')
  			if not si_depl_idx ==-1:
  				plt.plot(rho[si_depl_idx],T[si_depl_idx],marker='D',color='g')
  			
  
  			#ax = plt.gca()
  			#ax.invert_xaxis()
  			plt.rcParams.update({'font.size': 16})
                  	plt.rc('xtick', labelsize=16)
                  	plt.rc('ytick', labelsize=16)
                  	legend(loc=4)
  		        plt.xlabel('log $\\rho_{\\rm c}$',fontsize=18)
          		plt.ylabel('log $T_{\\rm c}$',fontsize=18)
  			#plt.gca().invert_xaxis()
  			'''
  			i+=1	
  		if deg_line==True:
  			rho=np.arange(0,9,0.01)
  			T_degeneracy=2./3. *rho +np.log10(1.207e5 * 1.8/(2.**(5./3.)))
  			#T_degeneracy=np.log10( (10**np.array(rho))**(2./3.))
  			plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
  		plt.legend(loc=2)