[yt-users] Overplotting equipotential contours from total potential

[trobolo dinni]

Dear yt users,

I am running a binary star simulation with enzo, and in my box there is a
star made up of a gas envelope and a particle as core, and a secondary star
represented by a single particle (no gas).

I would like to plot the equipotential surfaces of the system on a slice
plot, hence, since the default Grav_Potential field does not include the
particles potential, I created a custom made field as follows:

*  def _Total_Grav_Potential(field,data): *
* #G in cgs*
* grav_const_cgs = 6.674e-8*
* *
* #data for gas*
* gas_grav_potential_cgs = data['Grav_Potential']  *
(length_unit1/time_unit1)**2.0 #in units of velocity^2*
* gas_position_x_cgs = data['x'] * length_unit1*
* gas_position_y_cgs = data['y'] * length_unit1*
* gas_position_z_cgs = data['z'] * length_unit1*
* *
* #initialize the array for the total grav potential to the gas grav
potential*
* total_grav_potential_cgs = gas_grav_potential_cgs*
* *
* #loops on the particles*
* for i in range(len(data['ParticleMassMsun'])):*
* *
* #data of the current particle*
* current_particle_mass_cgs= data['ParticleMassMsun'][i] * mass_unit*
* current_particle_position_x_cgs = data['particle_position_x'][i] *
length_unit1*
* current_particle_position_y_cgs = data['particle_position_y'][i] *
length_unit1*
* current_particle_position_z_cgs = data['particle_position_z'][i] *
length_unit1*
* *
* #computes the array of grav potential for the current particle*
* current_particle_grav_potential_cgs = -(grav_const_cgs *
current_particle_mass_cgs)/(((gas_position_x_cgs -
current_particle_position_x_cgs)**2.0 + (gas_position_y_cgs -
current_particle_position_y_cgs)**2.0 + (gas_position_z_cgs -
current_particle_position_z_cgs)**2.0)**0.5)*
* *
* #adds the current particle potential to the total potential*
* total_grav_potential_cgs = total_grav_potential_cgs +
current_particle_grav_potential_cgs*
* *
* return total_grav_potential_cgs*
* #adds the field to the available yt fields*
* add_field("Total_Grav_Potential", function=_Total_Grav_Potential,
take_log=True ,units=r"\rm{erg}/\rm{g}")*
*
*
where the length unit and time unit variables have the function to convert
my quantities in cgs.
If I then call my "Total_Grav_Potential" field, I get reasonable values,
but when I try to inspect the shape of the equipotential surfaces on the
xy-plane by plotting the following:

* sp = SlicePlot(pf, 'z', "Density", width = 1.0)*
* sp.annotate_velocity(factor=16, normalize=True)*
*
sp.annotate_contour('Total_Grav_Potential',clim=(max_total_pot,min_total_pot))
*
* sp.annotate_particles(1.0, p_size=50.0, marker='o', col='black')*
* sp.annotate_text((0.7,1.05), "time = "+str(pf.current_time)+"yr") *
* sp.save(plot_dir+"density_slice_z_"+str(index)+".png")*

nothing happens, and the plot does not show any additional contour; while
if I do the same with the default "Grav_Potential" field, the gas
equipotential contours get plotted without problems.
I also tried to play with the "annotate_contour" options "ncont", "factor"
and "clim" to refine my plot but the result is always the same.
Additionally, if I do not set the "clim" option the following error shows
up:

*ValueError: zero-size array to minimum.reduce without identity*
*
*

I would like to ask if my custom definition of the Total_Grav_Potential
field is correct from a yt coding point of view and if I should modify it
(or if there is any improvement I can do).
Also, I would like to ask if you have any idea of what could be the problem
with the plotting.


Thanks for the help,
                                 Roberto
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20130924/1264c635/attachment.htm>