Is anybody here an expert on the way enzo_anyl calculates spin parameter? I was hoping to convince myself that:
a) I understand how it does it in enzo_anyl b) That is being replicated in DerivedQuantities.py
I have gone back and forth on this with Brian and Britton, but that was almost a year ago. I think it needs another look. Here is how yt works right now:
am = data["SpecificAngularMomentum"]*data["CellMassMsun"] j_mag = am.sum(axis=1) m_enc = data["CellMassMsun"].sum() + data["ParticleMassMsun"].sum() e_term_pre = na.sum(data["CellMassMsun"]*data["VelocityMagnitude"]**2.0) weight=data["CellMassMsun"].sum()
so we get the sum of Lx, Ly, Lz, the total mass, the total kinetic energy in the baryons, and the total *baryon* mass.
then during the combination step:
W = weight.sum() M = m_enc.sum() J = na.sqrt(((j_mag.sum(axis=0))**2.0).sum())/W E = na.sqrt(e_term_pre.sum()/W) G = 6.67e-8 # cm^3 g^-1 s^-2 spin = J * E / (M*1.989e33*G)
What this does it combine all the weights from the individual grid or processors to get the total baryon mass in the entire region, the entire *enclosed mass* (which includes the particles), and then the magnitude of the angular momentum vector for all the enclosed baryons, which gets divided by the enclosed *baryon mass* to get the average specific angular momentum for the region. E is then the total kinetic energy divided by the total enclosed mass, which gives a characteristic baryon velocity. Finally, we take the average specific angular momentum, multiply that by the characteristic velocity, and divide by the total enclosed (baryon+particle) mass.
Does this make sense? Should the characteristic velocity and angular momentum include the particles? Or does this not matter?
-Matt
Okay, after some tests, it turns out that I think I was doing this all correctly. So I apologize for the noise.
On my current run, enzo_anyl gives:
spin gas: 0.0251921 for the same radius, yt gives: 0.0251090534601
This is not a huge deal, I think. I'm not sure *why* it doesn't give the same result past four decimal places, but I am not sure that it is a meaningful difference.
I will put in the DM spin parameter later today...
-Matt
On Sat, Jan 24, 2009 at 9:20 AM, Matthew Turk matthewturk@gmail.com wrote:
Is anybody here an expert on the way enzo_anyl calculates spin parameter? I was hoping to convince myself that:
a) I understand how it does it in enzo_anyl b) That is being replicated in DerivedQuantities.py
I have gone back and forth on this with Brian and Britton, but that was almost a year ago. I think it needs another look. Here is how yt works right now:
am = data["SpecificAngularMomentum"]*data["CellMassMsun"] j_mag = am.sum(axis=1) m_enc = data["CellMassMsun"].sum() + data["ParticleMassMsun"].sum() e_term_pre = na.sum(data["CellMassMsun"]*data["VelocityMagnitude"]**2.0) weight=data["CellMassMsun"].sum()
so we get the sum of Lx, Ly, Lz, the total mass, the total kinetic energy in the baryons, and the total *baryon* mass.
then during the combination step:
W = weight.sum() M = m_enc.sum() J = na.sqrt(((j_mag.sum(axis=0))**2.0).sum())/W E = na.sqrt(e_term_pre.sum()/W) G = 6.67e-8 # cm^3 g^-1 s^-2 spin = J * E / (M*1.989e33*G)
What this does it combine all the weights from the individual grid or processors to get the total baryon mass in the entire region, the entire *enclosed mass* (which includes the particles), and then the magnitude of the angular momentum vector for all the enclosed baryons, which gets divided by the enclosed *baryon mass* to get the average specific angular momentum for the region. E is then the total kinetic energy divided by the total enclosed mass, which gives a characteristic baryon velocity. Finally, we take the average specific angular momentum, multiply that by the characteristic velocity, and divide by the total enclosed (baryon+particle) mass.
Does this make sense? Should the characteristic velocity and angular momentum include the particles? Or does this not matter?
-Matt
Done with particle spin parameter. Caught some other bugs -- particle velocities not getting turned to CGS, mainly.
As a note, I think that all particle fields using velocities should carefully choose between using bulk_velocity and particle_bulk_velocity or something similarly named. That's another new quantity that should be added. Furthermore, I believe the minor differences between spin parameters in YT and in enzo_anyl is due to using actual enclosed mass versus rho_vir * 4/3 pi r_vir^3.
This raises a new point, though -- perhaps it'd be nice if we had a routine to set up bulk_velocity and whatnot. It'd be good if we did this so that you could specify the sub-radius for a smaller bulk_velocity calculation.
-Matt
On Sat, Jan 24, 2009 at 10:41 AM, Matthew Turk matthewturk@gmail.com wrote:
Okay, after some tests, it turns out that I think I was doing this all correctly. So I apologize for the noise.
On my current run, enzo_anyl gives:
spin gas: 0.0251921 for the same radius, yt gives: 0.0251090534601
This is not a huge deal, I think. I'm not sure *why* it doesn't give the same result past four decimal places, but I am not sure that it is a meaningful difference.
I will put in the DM spin parameter later today...
-Matt
On Sat, Jan 24, 2009 at 9:20 AM, Matthew Turk matthewturk@gmail.com wrote:
Is anybody here an expert on the way enzo_anyl calculates spin parameter? I was hoping to convince myself that:
a) I understand how it does it in enzo_anyl b) That is being replicated in DerivedQuantities.py
I have gone back and forth on this with Brian and Britton, but that was almost a year ago. I think it needs another look. Here is how yt works right now:
am = data["SpecificAngularMomentum"]*data["CellMassMsun"] j_mag = am.sum(axis=1) m_enc = data["CellMassMsun"].sum() + data["ParticleMassMsun"].sum() e_term_pre = na.sum(data["CellMassMsun"]*data["VelocityMagnitude"]**2.0) weight=data["CellMassMsun"].sum()
so we get the sum of Lx, Ly, Lz, the total mass, the total kinetic energy in the baryons, and the total *baryon* mass.
then during the combination step:
W = weight.sum() M = m_enc.sum() J = na.sqrt(((j_mag.sum(axis=0))**2.0).sum())/W E = na.sqrt(e_term_pre.sum()/W) G = 6.67e-8 # cm^3 g^-1 s^-2 spin = J * E / (M*1.989e33*G)
What this does it combine all the weights from the individual grid or processors to get the total baryon mass in the entire region, the entire *enclosed mass* (which includes the particles), and then the magnitude of the angular momentum vector for all the enclosed baryons, which gets divided by the enclosed *baryon mass* to get the average specific angular momentum for the region. E is then the total kinetic energy divided by the total enclosed mass, which gives a characteristic baryon velocity. Finally, we take the average specific angular momentum, multiply that by the characteristic velocity, and divide by the total enclosed (baryon+particle) mass.
Does this make sense? Should the characteristic velocity and angular momentum include the particles? Or does this not matter?
-Matt
-
Matthew Turk