mksedov.py
No OneTemporary
Actions

Subscribers

None

File Metadata

	#!/usr/bin/env python

	from pNbody import *
	from pNbody import ic
	from pNbody import libutil

	from numpy import *


	import sys
	import time

	import copy

	from PyGear import gadget


	def grid(n=10,m=10,M=1.,name='grid.dat',ftype='binary'):

	dx = 1/float(n)
	dy = 1/float(m)

	xx = arange(0,1+dx,dx)
	yy = arange(0,1+dy,dy)

	x = zeros(n*m)
	y = zeros(n*m)
	z = zeros(n*m)

	k = 0
	for i in xrange(n):
	for j in xrange(m):
	x[k+0] = xx[i+0] - 0.5
	y[k+0] = yy[j+0] - 0.5
	k = k + 1


	n = len(x)
	pos = transpose((x,y,z))
	mass= M*ones(n)/n
	nb = Nbody(status='new',pos=pos,mass=mass,p_name=name,ftype=ftype)
	return nb


	###################################
	# generate model

	n = 500

	'''
	# random
	nb = ic.box(n,0.5,0.5,0.5,ftype='gadget',irand=0)
	nb.pos[:,2] = nb.pos[:,2]*0
	nb.pos = nb.pos + 0.5
	nb.u = ones(nb.nbody)
	'''


	# grid
	n = int(sqrt(500))
	nb = grid(n,n,ftype='gadget')
	nb.pos = nb.pos + [0.5/n,0.5/n, 0]
	nb.pos[:,2] = nb.pos[:,2]*0
	nb.pos = nb.pos + 0.5

	nb.pos = nb.pos + random.random(nb.pos.shape)1./float(n)0.3

	nb.u = ones(nb.nbody)/float(nb.nbody)




	# sedov
	r = nb.rxyz(center=[0.5,0.5,0.5])
	idx = argmin(r)
	nb.u[idx] = 10


	nb.rename('snap.dat')
	nb.write()