[yt-users] Surface density profile

[Jun-Hwan Choi]

Hi yt user,

I try to compute the surface density profile of the gas disk from my 
enzo simulation.
I perform this by using profile model after define disk object.
In order to check the output, I compute the surface density by making 
the annulus object after making two disk object and using boolean object.
The resulting surface density profiles are quite different.
I think that it may result from the interpolation of intersecting grid 
cell with the object.
Is there any explanation that the interpolation?
And, does anyone can suggest which is the more accurate approach?

FYI, I will past my yt script here
======================================
from yt.mods import *
import numpy as na
import matplotlib.pyplot as plt
import pylab as pl
import math as math
from itertools import cycle

lines = ["-","--"]
linecycler = cycle(lines)
colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k')
colorcycler = cycle(colors)


def _RadiuspcXY(field, data):
     center = data.get_field_parameter("center")
     DW = data.pf.domain_right_edge - data.pf.domain_left_edge
     radius = na.zeros(data["x"].shape, dtype='float64')
     for i, ax in enumerate('xy'):
         r = na.abs(data[ax] - center[i])
         radius += na.minimum(r, na.abs(DW[i]-r))**2.0
     na.sqrt(radius, radius)
     return radius
def _ConvertRadius(data):
     return data.convert("pc")
add_field("RadiuspcXY", function=_RadiuspcXY,
           validators=[ValidateParameter("center")],
           convert_function = _ConvertRadius, units=r"\rm{pc}",
           display_name="RadiusXY")

Msun_in_cgs = 2e33
rmin = 2e-4
rmax  = 20
rdisk = rmax
hdisk = 1
Nbin = 75
dr = (na.log(rmax)-na.log(rmin))/Nbin

index = 153
pf = load("DD0153/DD0153")
center = pf.h.find_max("Density")[1]

radius = []
Sigma  = []
for i in range(0,Nbin):
     rlow  = na.exp(na.log(rmin)+dr*i)
     rhigh = na.exp(na.log(rmin)+dr*(i+1))
     rcenter = 0.5*(rlow+rhigh)

     disk1 = pf.h.disk(center, [0,0,1], rhigh/pf.units['pc'], 
hdisk/pf.units['pc'] )
     disk2 = pf.h.disk(center, [0,0,1], rlow/pf.units['pc'], 
hdisk/pf.units['pc'] )
     disk3 = pf.h.boolean([disk1, "NOT", disk2])

     Area = math.pi*(rhigh**2 - rlow**2)
     radius.append(rcenter)
     Sigma.append(disk3["CellMassMsun"].sum()/Area)


plt.loglog(radius, Sigma, label='Profile1',
            linestyle=next(linecycler),color=next(colorcycler))

# Profile
disk = pf.h.disk(center, [0,0,1], rdisk/pf.units['pc'], 
hdisk/pf.units['pc'] )
profile = BinnedProfile1D(disk, Nbin, "RadiuspcXY", rmin, rdisk, 
log_space=True, lazy_reader=True, end_collect=False)
profile.add_fields('Density',weight='CellMassMsun')

n_bins = profile['RadiuspcXY'].size
Sigma2 = na.zeros(n_bins)
Sigma2 = 
profile['Density']*((pf.units['cm']/pf.units['pc'])**3/Msun_in_cgs)*hdisk

plt.loglog(profile['RadiuspcXY'], Sigma2, label='Profile2',
            linestyle=next(linecycler),color=next(colorcycler))

plt.xlabel(r'$R_{xy}$ [pc]')
plt.ylabel(r'$\Sigma$')
plt.legend(loc=1)

plt.savefig("SurfaceDensity_comp")
========================================================

Thank you in advance,
Junhwan

P.S.: I think there is a way to post my script in web, but I forget where.



-- 
--------------------------------------------------------------
Jun-Hwan Choi, Ph.D.
Department of Physics and Astronomy, University of Kentucky
Tel: (859) 897-6737        Fax: (859) 323-2846
Email: jhchoi at pa.uky.edu   URL: http://www.pa.uky.edu/~jhchoi
--------------------------------------------------------------