plot_temperature.pro

pic 模拟程序！面向对象

;------------------------------------------------------------------------
;
;  File:    plot_temperature.pro
;
;  Purpose:
;    IDL procedure to read in particle data from HDF5 files generated by
;      a parallel OOPIC run.  One then can plot the temperature.
;
;  Overview:
;
;    Argument list:
;
;    base_name  -> a string that specifies the first input
;                  data file (base_name.txt) and is also used to name 
;                  output files.
;    num_procs  -> number of processors
;    XLabel     -> a string that specifies the coordinate component, e.g. 'x'.
;    XLabel2    -> a string that specifies 2nd coordinate component, e.g. 'y'.
;    ScaleFlag  -> a flag to specify the scaling for length, i.e. of q:
;                        0 -> no scaling
;                        1 -> length scale = c/omega_0 (the laser frequency)
;                        2 -> length scale = c/omega_p (the plasma frequency)
;    ScaleInput -> this is irrelevant for ScaleFlag = 0,
;                          or omega_0 for ScaleFlag = 1,
;          plasma density in [1/cm^3] for ScaleFlag = 2
;
;  Copyright (c) 2003 by Tech-X Corporation.  All rights reserved.
;
;  Change log:
;  1.00 (Bruhwiler 10-24-2003) Initial code.
;
;------------------------------------------------------------------------


pro plot_temperature, base_name, num_procs, XLabel, YLabel, ScaleFlag, ScaleInput

; *****************************************************************
; Names for files and labels
; *****************************************************************

; Specify the text, restore and postscript file names
restore_file  = base_name  + '.dat'

xcontour_file = base_name  + '_vx_rms.ps'
ycontour_file = base_name  + '_vy_rms.ps'
vcontour_file  = base_name  + '_v_rms.ps'

xlineout_file = base_name  + '_lineout_vx_avg.ps'
ylineout_file = base_name  + '_lineout_vy_avg.ps'

x2lineout_file = base_name  + '_lineout_vx_rms.ps'
y2lineout_file = base_name  + '_lineout_vy_rms.ps'

txlineout_file = base_name  + '_lineout_tx.ps'
tylineout_file = base_name  + '_lineout_ty.ps'

tlineout_file  = base_name  + '_lineout_temp.ps'
vlineout_file  = base_name  + '_lineout_v_rms.ps'


; *****************************************************************
; Specify which plots you want (1) or don't want (0):
; *****************************************************************

doTxContour = 1
doTyContour = 1
doTContour  = 1

; Specify non-default window sizes for contour, surface, show3d plots
; Dimension is pixels for the window command
x_win_ccon = 1000
y_win_ccon =  350
; Dimension is cm for the device command
x_device = 25.0
y_device = 12.0

; *****************************************************************
; Set a factor (between 0 and 1) specifying which row of data
;    (r=constant or y=constant) that you want for the lineout plot.
r_factor = 0.5

; *****************************************************************
; Parse the ascii data file or restore from the binary IDL file
; *****************************************************************

; If the restore file exists, then use it.
; Otherwise, parse the text file:

spawn, "ls " + restore_file, check_file, /sh
;print, 'The "check"   file is: ', check_file
;print, 'The "restore" file is: ', restore_file

if (check_file(0) eq restore_file) then begin

  print,'Reading from the restore file: ' + restore_file + ' ...'
  restore, filename = restore_file

endif else begin

; First, we have to determine the total number of particles to be parsed.

  print, ' '
  print, 'Calculating total number of particles to be parsed from the hdf5 files named'
  print, '   ' + base_name + '_p*.h5  ...'
  print, ' '
  total_num_ptcls = 0
  for i_loop=0, num_procs-1 do begin           ; beginning of major for loop

    hdf5_file = base_name + '_p' + strtrim(i_loop,2) + '.h5'
    file_id = h5f_open(hdf5_file)

    group_id     = h5g_open(file_id, '/electrons')
    attribute_id = h5a_open_name(group_id, 'nPtclGrp')
    num_groups = h5a_read(attribute_id)
    h5a_close, attribute_id
    h5g_close, group_id

    for j=0, num_groups-1 do begin           ; loop over the particle groups

      dataset = '/electrons/pGrp' + strtrim(j,2)
      dataset_id   = h5d_open(file_id, dataset+'/ptcls')
      dataspace_id = h5d_get_space(dataset_id)
      npoints = h5s_get_simple_extent_npoints(dataspace_id)
      h5d_close, dataset_id

      total_num_ptcls = total_num_ptcls + npoints/5

    endfor      ; end of sub-loop over particle groups

; Close the file
    h5f_close, file_id

  endfor        ; end of major loop over num_procs
  print, 'The total number of particles to be parsed is: ', total_num_ptcls
  print, ' '




; *****************************************************************
; Define the 1D particle data arrays

  x  = dindgen(total_num_ptcls)
  y  = dindgen(total_num_ptcls)
  vx = dindgen(total_num_ptcls)
  vy = dindgen(total_num_ptcls)
  vz = dindgen(total_num_ptcls)
  wt = dindgen(total_num_ptcls)

; Here we actually read the particle data out of the HDF5 files
  num_parsed = 0
  for i_loop=0, num_procs-1 do begin           ; beginning of major for loop

    hdf5_file = base_name + '_p' + strtrim(i_loop,2) + '.h5'

    print,'Parsing the hdf5 file: ' + hdf5_file + ' ...'
    file_id = h5f_open(hdf5_file)

    group_id     = h5g_open(file_id, '/electrons')
    attribute_id = h5a_open_name(group_id, 'nPtclGrp')
    num_groups = h5a_read(attribute_id)
    h5a_close, attribute_id
    h5g_close, group_id

    for j=0, num_groups-1 do begin           ; loop over the particle groups

      dataset = '/electrons/pGrp' + strtrim(j,2)

      group_id     = h5g_open(file_id, dataset)
      attribute_id = h5a_open_name(group_id, 'np2c')
      ptcl_weight  = h5a_read(attribute_id)
      h5a_close, attribute_id
      h5g_close, group_id

      dataset_id   = h5d_open(file_id, dataset+'/ptcls')
      ptcl = h5d_read(dataset_id)
      h5d_close, dataset_id

      num_ptcl   = n_elements(ptcl) / 5
      num_prev   = num_parsed
      num_parsed = num_parsed + num_ptcl

; Load the raw data into individual arrays.
       x(num_prev:num_parsed-1) = ptcl(0,0:num_ptcl-1)
       y(num_prev:num_parsed-1) = ptcl(1,0:num_ptcl-1)
      vx(num_prev:num_parsed-1) = ptcl(2,0:num_ptcl-1)
      vy(num_prev:num_parsed-1) = ptcl(3,0:num_ptcl-1)
      vz(num_prev:num_parsed-1) = ptcl(4,0:num_ptcl-1)
      wt(num_prev:num_parsed-1) = ptcl_weight

    endfor      ; end of sub-loop over particle groups

; Close the file
    h5f_close, file_id

  endfor        ; end of major loop over num_procs

  print, ' '
  print, 'Number of parsed particles is: ', num_parsed
  print, 'Number of elements in x() is:  ', n_elements(x)

; Save so IDL doesn't have to repeatedly parse the HDF5 file
  print, ' '
  print, 'Now generating an IDL restore file with the name: ' + restore_file + ' ...'
  save, x, y, vx, vy, vz, wt, filename = restore_file

endelse


; *****************************************************************
; Scale the data so it corresponds to the desired units
; *****************************************************************
;
q_label  = '!3' + XLabel
q_label2 = '!3' + YLabel
if ( ScaleFlag eq 1 ) then begin
  ;
  ; scale the length using as a c/omega_0 as a length scale
  ; 
  omega_0 = ScaleInput
  lengthScale = 3.0e8/omega_0
  x = x / lengthScale
  y = y / lengthScale
endif else if (ScaleFlag eq 2) then begin
  ;
  ; scale lengths to the plasma wavenumber k_p
  ; 
  density = ScaleInput
  omega_p = 2. * 3.14 * 9000. * sqrt(density)
  x = x * omega_p / 2.998e+08
  y = y * omega_p / 2.998e+08
endif else begin
  ;
  ; this is the default case of no scaling 
  ;
endelse

; Remove the relativistic factor gamma from the velocities:
 c     = 2.998e+08
 c_sq  = c*c
 v_sq  = vx*vx + vy*vy
 gamma = sqrt(1. + v_sq / c_sq)
 vx = vx / gamma
 vy = vy / gamma

;
; *****************************************************************
; set the min and max values of the data
; *****************************************************************
;
autoscale = 1
if (autoscale eq 1 ) then begin 
  x_min = min(x)
  x_max = max(x)
  y_min = min(y)
  y_max = max(y)
endif else begin
  x_min  =  0.0
  x_max  =  5.3
  y_min =  0.0
  y_max =  2.0
endelse

x_min  =  2.0
y_min =  0.5 * y_max - 5.
y_max =  y_min + 10.


;
; *****************************************************************
; create the spatial grid
; *****************************************************************

numGridsX  = 512
numGridsY =  96

delta_x = (x_max - x_min) / numGridsX
delta_y = (y_max - y_min) / numGridsY

x_grid = dblarr(numGridsX)
y_grid = dblarr(numGridsY)
for i = 0, numGridsX-1  do begin
  x_grid[i] = (i+0.5) * delta_x
end
for i = 0, numGridsY-1  do begin
  y_grid[i] = (i+0.5) * delta_y
end


; *****************************************************************
; Loop for rendering 2-D color contour plot on screen and to a file
; *****************************************************************

if (doTxContour eq 1) then begin
 
; calculate gridded temperature for the 1st velocity component
  print, ' '
  print, 'Now gridding up the temperature to load rms velocity v_x_rms...'
  v_x_rms = temperature_2d(x, y, vx, wt, numGridsX, numGridsY, $
                           XAXIS=xax, YAXIS=yax,               $
                           VMEAN=v_x_avg, DV_MAX=dv_x_max,     $
                           N_MACRO=n_macro)

; Now rotate by 180 degrees to compare with Brad Shadwick's results
  v_x_rms = reverse(v_x_rms, 1)
  v_x_avg = reverse(v_x_avg, 1)
  v_x_ratio = reverse(v_x_ratio, 1)

  dv_x_ratio = dv_x_max
  ind = where(dv_x_max GE 2.e+06)
  dv_x_ratio(ind) = v_x_rms(ind) / dv_x_max(ind)
  ind = where(dv_x_max LT 2.e+06) 
  dv_x_ratio(ind) = 0.

; let's smooth the results a little
;  v_x_rms = smooth(v_x_rms, 3, /EDGE_TRUNCATE)
;  v_x_avg = smooth(v_x_avg, 3, /EDGE_TRUNCATE)

; Calculate the corresponding temperature
  T_x  = 0.511e+06 * (v_x_rms/2.998e+08)^2
;  T_x  = smooth(T_x,     3, /EDGE_TRUNCATE)

; sometimes it's helpful to have a log scale to see what is going on
  T_x_log = alog10( T_x + 1.0e-10 * max(T_x) )

; Get a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Color contour plot will appear in window ', window_number
  window, window_number, xsize=x_win_ccon, ysize=ywin_ccon, RETAIN=2
  xcontour_i = 0
  xcontour_jump:

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    if ( ScaleFlag eq 1 ) then begin
      x_label = q_label  + ' (c/!9w!3!d0!N)'
      y_label = q_label2 + ' (c/!9w!3!d0!N)'
    endif else if ( ScaleFlag eq 2 ) then begin
      x_label = q_label  + ' (c/!9w!3!dp!N)'
      y_label = q_label2 + ' (c/!9w!3!dp!N)'
    endif 
    print, ' '
    print, 'Writing the 2D color contour plot to file ' + xcontour_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0
    if ( ScaleFlag eq 1 ) then begin
      x_label = q_label  + ' (c/!4x!3!d0!N)'
      y_label = q_label2 + ' (c/!4x!3!d0!N)'
    endif else if ( ScaleFlag eq 2 ) then begin
      x_label = q_label  + ' (c/!4x!3!dp!N)'
      y_label = q_label2 + ' (c/!4x!3!dp!N)'
    endif 
    device,decomposed=0
    print, ' '
    print, 'Rendering the 2D color contour plot to the screen...'
  endelse

; Load in the RAINBOW color table
  loadct, 13