from yt.mods import * from yt.analysis_modules.halo_finding.api import * from matplotlib import use; use('Agg') import matplotlib.pyplot as plt import numpy as np import sys import glob import re #try: # from mpi4py import MPI # comm = MPI.COMM_WORLD # nprocs = comm.size # rank = comm.rank #except: # nprocs=1 # rank=0 rank=0 fn, = glob.glob('RedshiftOutput????') output = int(re.sub(r'\D',"",fn)) pf = load(fn) #haloes = LoadHaloes(pf, "RD%04i"%(output)) f=open("../FOF/RedshiftOutput%04i-halos.dat"%output,'r') line=f.readline();line=f.readline() f.close() c=np.array([float(line.split()[0]),float(line.split()[1]),float(line.split()[2])]) rvir=float(line.split()[8]) #c=haloes[0].center_of_mass() #vl=pf["RefineRegionLeftEdge"] #vr=pf["RefineRegionRightEdge"] #sv = pf.h.region((vl+vr)/2, vl, vr) #width=(vr-vl).max() #c=(vl+vr)/2 z=pf["CosmologyCurrentRedshift"] wkpc = 15.0 #in proper kpc tkpc = 0.6 # in proper kpc width =wkpc*(1.0+z)/40000.0 #in kpc thick =tkpc*(1.0+z)/40000.0 #in kpc #wkpc=width*40000/(1.0+z) #width=wkpc*(z+1.0)/40.0/1000 #c=np.array([0.4341635423, 0.4348677931, 0.5726667823]) #c=[0.5,0.5,0.5] #RL=c-width/2 #RL=[0.425, 0.425, 0.425] #RR=c+width/2 #RR=[0.575, 0.575, 0.575] #sv = pf.h.region(c, RL, RR) RL=np.zeros(3) RR=np.zeros(3) #haloes = LoadHaloes(pf, "RD%04i"%output) if rank == 0: RL[0] = c[0]-width/2 RR[0] = c[0]+width/2 RL[1] = c[1]-width/2 RR[1] = c[1]+width/2 RL[2] = c[2]-thick/2 RR[2] = c[2]+thick/2 sv = pf.h.region(c, RL, RR) p = ProjectionPlot(pf, 2, ["Density","flux","Temperature"], c,weight_field="Density",width=(wkpc, 'kpc'),data_source=sv) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,stars_only=True) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,col='w',ptype=7) #p.annotate_hop_circles(haloes,min_size=200,width=wkpc/pf['kpc']) #p.set_zlim('Density',5e-4,1e-2) p.set_zlim('flux',1e-23,1e-15) p.set_zlim('Density',1e-28,5e-21) p.set_zlim('Temperature',1e3,1e8) p.set_cmap(field="flux", cmap='idl05') p.annotate_sphere(c,radius=rvir/pf['kpc']) for i in ['Density','flux','Temperature']: p.plots[i].hide_colorbar();p.plots[i].hide_axes() p.save('%3.1fkpc_%4.2f'%(wkpc,z)) if rank == 0: RL[0] = c[0]-width/2 RR[0] = c[0]+width/2 RL[2] = c[2]-width/2 RR[2] = c[2]+width/2 RL[1] = c[1]-thick/2 RR[1] = c[1]+thick/2 sv = pf.h.region(c, RL, RR) p = ProjectionPlot(pf, 1, ["Density","flux","Temperature"], c,weight_field="Density", width=(wkpc, 'kpc'),data_source=sv) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,stars_only=True) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,col='w',ptype=7) #p.annotate_hop_circles(haloes,min_size=200,width=wkpc/pf['kpc']) #p.set_zlim('all',5e-4,1e-2) p.set_zlim('flux',1e-23,1e-15) p.set_zlim('Density',1e-28,5e-21) p.set_zlim('Temperature',1e3,1e8) p.set_cmap(field="flux", cmap='idl05') p.annotate_sphere(c,radius=rvir/pf['kpc']) for i in ['Density','flux','Temperature']: p.plots[i].hide_colorbar();p.plots[i].hide_axes() p.save('%3.1fkpc_%4.2f'%(wkpc,z)) if rank == 0: RL[2] = c[2]-width/2 RR[2] = c[2]+width/2 RL[1] = c[1]-width/2 RR[1] = c[1]+width/2 RL[0] = c[0]-thick/2 RR[0] = c[0]+thick/2 sv = pf.h.region(c, RL, RR) p = ProjectionPlot(pf, 0, ["Density","flux","Temperature"], c,weight_field="Density", width=(wkpc, 'kpc'),data_source=sv) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,stars_only=True) #p.annotate_particles(wkpc/pf['kpc'],p_size=0.5,col='w',ptype=7) #p.annotate_hop_circles(haloes,min_size=200,width=wkpc/pf['kpc']) #p.set_zlim('all',5e-4,1e-2) p.set_zlim('Density',1e-28,5e-21) p.set_zlim('Temperature',1e3,1e8) p.set_zlim('flux',1e-23,1e-15) p.set_cmap(field="flux", cmap='idl05') p.annotate_sphere(c,radius=rvir/pf['kpc']) for i in ['Density','flux','Temperature']: p.plots[i].hide_colorbar();p.plots[i].hide_axes() p.save('%3.1fkpc_%4.2f'%(wkpc,z)) #p = ProjectionPlot(pf, 1, ["Density","TotalMetallicity"], c, width=(2*wkpc, 'kpc'),data_source=sv) #p.annotate_particles(2*wkpc/pf['kpc'],ptype=5)#,stars_only=True) #p.annotate_hop_circles(haloes,min_size=200,width=wkpc/pf['kpc']) #p.save('%3.1fkpc_%3.1f'%(2*wkpc,z)) #p = ProjectionPlot(pf, 0, ["Density","TotalMetallicity"], c, width=(2*wkpc, 'kpc'),data_source=sv) #p.annotate_particles(2*wkpc/pf['kpc'],ptype=5)#,stars_only=True) #p.annotate_hop_circles(haloes,min_size=200,width=wkpc/pf['kpc']) #p.save('%3.1fkpc_%3.1f'%(2*wkpc,z)) #pc.add_projection('TotalMetallicity',2,data_source=sv) #pc.add_projection('Temperature',2,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',2,weight_field='Density',data_source=sv) #pc.set_width(2*wkpc,'kpc') #pc.set_zlim(1e-1,1e-3) #pc.save('%3.1fkpc_%3.1f'%(2*wkpc,z)) #del pc #pc = PlotCollection(pf, center=c) #pc.add_projection('Density',2,data_source=sv) #pc.add_projection('TotalMetallicity',2,data_source=sv) #pc.add_projection('Temperature',2,weight_field='Density',data_source=sv) #pc.add_projection('Temperature',1,weight_field='Density',data_source=sv) #pc.add_projection('Temperature',0,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',2,weight_field='Density',data_source=sv) #pc.set_width(2*width,'1') #pc.set_zlim(1e-1,1e-3) #pc.save("Halo") #del pc #pc = PlotCollection(pf, center=c) #pc.add_projection('Density',2,data_source=sv) #pc.add_projection('TotalMetallicity',2,data_source=sv) #pc.add_projection('Temperature',2,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',2,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',1,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',0,weight_field='Density',data_source=sv) #pc.set_width(2*width,'1') #pc.set_zlim(1e-1,1e-6) #pc.save("Halo") #del pc #c=na.array([0.4882009076, 0.4725569205, 0.4654072348]) #pc = PlotCollection(pf, center=c) #sv = pf.h.region(c,c-5.0/pf['kpc'],c+5.0/pf['kpc']) #pc.add_projection('Density',2,data_source=sv) #pc.add_projection('Temperature',2,weight_field='Density',data_source=sv) #pc.add_projection('TotalMetallicity',2,weight_field='Density',data_source=sv) #pc.set_width(10,'kpc') #pc.set_zlim(2e-1,1e-4) #pc.save("most_massive_halo") #del pc