Changeset 30 for trunk/TOOLS/render_particle_data.pro

Timestamp:

07/30/10 13:16:33 (22 months ago)

Author:

tabel

Message:

Added radial and tangential velocity as well as cylindrical radius for the particles as a derived field.

File:

• trunk/TOOLS/render_particle_data.pro (modified) (8 diffs)

trunk/TOOLS/render_particle_data.pro

@@ -77 +77 @@
 @common_blocks.inc
 if N_elements(tcenter) le 1 then tcenter = center
-  Nderived = 7
+  Nderived = 20
+  id = 0L
   dfields = PTRARR(Nderived) 
-  dfields[0] = PTR_NEW(["abs(part. velocity)", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
-  dfields[1] = PTR_NEW(["particle radius", "particle_position_x", "particle_position_y", "particle_position_z"])
-  dfields[2] = PTR_NEW(["particle age", "creation_time"])
-  dfields[3] = PTR_NEW(["enclosed mass in particles","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
-  dfields[4] = PTR_NEW(["particle voronoi density","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
-  dfields[5] = PTR_NEW(["particle voronoi radius","particle_position_x", "particle_position_y", "particle_position_z"])
-  dfields[6] = PTR_NEW(["particle voronoi volume","particle_position_x", "particle_position_y", "particle_position_z"])
+  dfields[id++] = PTR_NEW(["abs(part. velocity)", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
+  dfields[id++] = PTR_NEW(["particle radial velocity","particle_position_x", "particle_position_y", "particle_position_z", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
+  dfields[id++] = PTR_NEW(["particle tangential velocity", "particle_position_x", "particle_position_y", "particle_position_z", "particle_velocity_x", "particle_velocity_y","particle_velocity_z"])
+  dfields[id++] = PTR_NEW(["particle radius", "particle_position_x", "particle_position_y", "particle_position_z"])
+  dfields[id++] = PTR_NEW(["particle cyl. radius", "particle_position_x", "particle_position_y", "particle_position_z"])
+  dfields[id++] = PTR_NEW(["particle age", "creation_time"])
+  dfields[id++] = PTR_NEW(["enclosed mass in particles","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
+  dfields[id++] = PTR_NEW(["particle voronoi density","particle_position_x", "particle_position_y", "particle_position_z", "particle massdensity"])
+  dfields[id++] = PTR_NEW(["particle voronoi radius","particle_position_x", "particle_position_y", "particle_position_z"])
+  dfields[id++] = PTR_NEW(["particle voronoi volume","particle_position_x", "particle_position_y", "particle_position_z"])
 ; the voronoi calculations are done in read_particle after all
 ; particles have been read.
+  dfields = dfields[0:id-1]
 
   if N_elements(stride) eq 1 then stride = LONARR(N_elements(usegi))+stride
@@ -94 +99 @@
      Nvar = N_ELEMENTS(VariableNames)
      data_format = determine_data_format_from_file_name(usegi[0].baryon_file)
-     if data_format lt 0 then return,0
+     if data_format lt 0 then return, 0
 ; determine all data fields to read and read data
      field_names = identify_unique_field_names(VariableNames, dfields)
@@ -106 +111 @@
         pack_particles, a
      
+     ax = double(slice_ori eq 0.)
      print, field_names
      for i=0, Nvar-1 do begin   ; if derived field was requested compute it here
@@ -115 +121 @@
                        -1: junk = ''
                        0: a[j].data[i] = ptr_new(sqrt(fn(a[j], "particle_velocity_x")^2+ $
-                                               fn(a[j], "particle_velocity_y")^2+ $
-                                               fn(a[j], "particle_velocity_z")^2))
-                       1: a[j].data[i] = ptr_new(sqrt((tcenter[0]-fn(a[j], "particle_position_x"))^2 + $
-                                               (tcenter[1]-fn(a[j], "particle_position_y"))^2 + $
-                                               (tcenter[2]-fn(a[j], "particle_position_z"))^2))
-                       2:  a[j].data[i] = ptr_new(times[time_index] - fn(a[j], "creation_time"))
-                       3: a[j].data[i] = $  ; allocate array
+                                                      fn(a[j], "particle_velocity_y")^2+ $
+                                                      fn(a[j], "particle_velocity_z")^2))
+                       1: a[j].data[i] = ptr_new(radial_velocity(tcenter, $
+                           fn(a[j], "particle_position_x"),fn(a[j], "particle_position_y"),fn(a[j], "particle_position_z"), $
+                           fn(a[j], "particle_velocity_x"),fn(a[j], "particle_velocity_y"),fn(a[j], "particle_velocity_z"),/cons_vel))
+
+                       2: a[j].data[i] = ptr_new(tangential_velocity(tcenter, $
+                           fn(a[j], "particle_position_x"),fn(a[j], "particle_position_y"),fn(a[j], "particle_position_z"), $
+                           fn(a[j], "particle_velocity_x"),fn(a[j], "particle_velocity_y"),fn(a[j], "particle_velocity_z")))
+
+                       3: a[j].data[i] = ptr_new(sqrt((tcenter[0]-fn(a[j], "particle_position_x"))^2 + $
+                                                      (tcenter[1]-fn(a[j], "particle_position_y"))^2 + $
+                                                      (tcenter[2]-fn(a[j], "particle_position_z"))^2))
+
+                       4: a[j].data[i] = ptr_new(sqrt(ax[0]*(tcenter[0]-fn(a[j], "particle_position_x"))^2 + $
+                                                      ax[1]*(tcenter[1]-fn(a[j], "particle_position_y"))^2 + $
+                                                      ax[2]*(tcenter[2]-fn(a[j], "particle_position_z"))^2))
+
+                       5: a[j].data[i] = ptr_new(times[time_index] - fn(a[j], "creation_time"))
+                       6: a[j].data[i] = $  ; allocate array
                           ptr_new(fn(a[j], "particle massdensity")) ; 
-                       4: a[j].data[i] = $
+                       7: a[j].data[i] = $
                           ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
-                       5: a[j].data[i] = $
+                       8: a[j].data[i] = $
                           ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
-                       6: a[j].data[i] = $
+                       9: a[j].data[i] = $
                           ptr_new(fn(a[j], "particle_position_x")) ; voronoi calculations done in read_particle_data
                    endcase
-                   print, max(*a[j].data[i])
-                   a[j].data_fields[i] = (*dfields[match])[0]
+;                   print, max(*a[j].data[i])
+;                   a[j].data_fields[i] = (*dfields[match])[0]
                endif
            endfor               ; loop over grids
@@ -280 +299 @@
 if npoints lt 10 then print, devicePart1, devicePart2
 
-pointimage = hist2d(devicePart1, devicePart2, $
-                      min1=1.,min2=1.,max1=xy_sl_size,max2=xy_sl_size, $
-                    binsize1=parti.size,binsize2=parti.size, density=rho)
+;pointimage = hist2d(devicePart1, devicePart2, $
+;                      min1=1.,min2=1.,max1=xy_sl_size,max2=xy_sl_size, $
+;                    binsize1=parti.size,binsize2=parti.size, density=rho)
 ;pointimage = rho
 
@@ -288 +307 @@
                 min1=1.,min2=1.,max1=xy_sl_size,max2=xy_sl_size, density=rho) ;/pointimage
 
-ind  = where(rho ge 0, count)
+print, size(himage)
+ind  = where(rho gt 0, count)
 if not shouldbe_histogram(names[parti.index]) then begin
    if count gt 0 then himage[ind] /= float(rho[ind])
@@ -294 +314 @@
 endif
 
-if parti.index le 2 then himage = pointimage
-;breaki
-;if parti.size ge 2 then pimage = congrid(pimage, xy_sl_size,
-;xy_sl_size, /center, cubic=(interpolate_i < 1), /minus_one)
-
-;print, 'minimum values in hist2d used:', om1, om2
-;breaki
-
 secmin = min(himage)
 ind  = where(himage gt secmin, count)
-
-print, secmin
-
-if ind[0] ge 0 then secmin = min(himage[ind])
-print, secmin
-
-
-if shouldbe_logged(names[parti.index]) then begin
-   oind = where(himage le 0. )
-   himage = himage > 0.5*secmin
-   himage = alog10(himage )
-endif
-
-
-print, min(himage)
-
-s = (size(himage))[1:2]
-dix = s[0]
-diy = s[1]
+;print, secmin
+
+
+if shouldbe_logged(names[parti.index]) or parti.index le 2 then begin
+;   himage = himage > 0.5*secmin
+   himage = alog10(himage > secmin)
+endif
+
+s = size(himage)
+dix = s[1]
+diy = s[2]
 locs = lindgen(dix*diy)
 nind = ARRAY_INDICES([dix,diy], locs, /DIMENSIONS)
@@ -332 +336 @@
   z_stuetz = himage[*]
   n_stuetz = N_elements(z_stuetz)
-
-  if parti.size lt 4 then $
-     pimage = (rebin(himage, dix*parti.size,diy*parti.size,  sample=1-(interpolate_i < 1)))[0:xy_sl_size-1, 0:xy_sl_size-1] $
-  else  pimage = render_image()
+if N_elements(z_stuetz) ne dix*diy then breaki
+;  if parti.size lt 4 then $
+;     pimage = (rebin(himage, dix*parti.size,diy*parti.size,  sample=1-(interpolate_i < 1)))[0:xy_sl_size-1, 0:xy_sl_size-1] $
+;  else  
+  pimage = render_image()
 
 return, pimage