PRO vikram_images, VIK_NAME=vik_name, PSELECT=pselect, $
	NPERRADIUS=nperradius
;+
; PRO vikram_images, VIK_NAME=vik_name, PSELECT=pselect, $
;	NPERRADIUS=nperradius
;
; This routine reads in a Vikram 5-column data file and
; outputs 'thin.png, thick.png and a fits-image file
; of the 3D and 2D projection of the scalar value...
;
; PSELECT value chooses from available scalars:
;  1: Mass Density:
;  3: Velocity:
;  5: Pressure ~ density * temp
;  7: Electron temperture
;  9: Xray flux approx ~ (Density)^2 where Temp > 1.e6
;
;
; - - -
; Can loop over all the files by:
;  IDL> for iv=1000, 1099 do vikram_images, psel=9, vik_name='wrbub.'+strcompress(iv,/remove)
;
; 2/16+/06 - dd
;-

; Get v3d_common variables available...
@v3d_common

; - - - Default values...
; the filename:
if n_elements(vik_name) EQ 0 then vik_name = 'wrbub.1067'
; which parameter selection
if n_elements(PSELECT) EQ 0 then pselect = 1  ; mass density
; default NPERRADIUS
if n_elements(nperradius) EQ 0 then nperradius=47   ; 37
; - - -


; - - - Get the data file values...
;
; Read in Vikram's data file:
; The columns are:
;   Radius (cms), Density (g/cc), Pressure, Velocity (cms/sec), Temperature
;   0.86983E+18  0.74111E-23  0.61268E-11  0.38385E+08  0.10000E+05
;   0.87374E+18  0.74111E-23  0.61268E-11  0.38558E+08  0.10000E+05

; on linux machine
vik_dir = "/nfs/spectra/d1/dd/Science/Vikram/wrbub_files/"
vik_out = "/nfs/spectra/d1/dd/Science/Vikram/wrbub_output/"
;
; on my mac PB
;;vik_dir = './wrbub_files/'
;;vik_out = './'


readcol, vik_dir+vik_name, vik_rad, vik_dens, $
  vik_pres, vik_vel, vik_temp

; Modify the units a bit...
amu = 1.660e-24
parsec = 3.1e18  ; cm

; Radius, pc
vik_rad = vik_rad/parsec

; Density, amu/cm^3
vik_dens = vik_dens/amu

; Velocity, km/s
vik_vel = vik_vel*1.e-5

; - - -


; - - - Select the scalar "intensity" to fill 3D space with...
;

; The Radius
; Scale the radius to 1.0:
rad_scl = vik_rad/MAX(vik_rad)
; add a first point at 0:
rad_scl = [0.0, rad_scl]
;

; Select a value for the scalar in the 3D array:


; Single values from the file:
if n_elements(pselect) GT 0 then begin
; - - - - - - -
case pselect of
1: begin
  ; Mass Density:
  val_scl = vik_dens
  val_name='mass'
end
3: begin
  ;  Velocity:
  val_scl = vik_vel
  val_name = 'velocity'
  end
5: begin
  ; Pressure
  val_scl = vik_pres
  val_name = 'pressure'
  end
7: begin
  ; Electron temperture:
  val_scl = vik_temp
  val_name = 't_e'
  end
9: begin
  ; X-ray flux approx
  val_scl = vik_dens * vik_dens * (vik_temp GT 1.e6)
  val_name='xray'
end
ELSE: print, ' oooopsies'
end ;case
; - - - - - - -
end


; - - - - - - - 
;
; Scale the value array to its maximum:
val_scl = val_scl/MAX(val_scl)
; add a point at 0.0, with the same value as cell(1):
val_scl = [val_scl(0), val_scl]
; 
; and plot it...
window, 1, XSIZE=300, YSIZE=220, YPOS=450, XPOS=0
plot_io, rad_scl, (val_scl > 0.001), YRANGE=[0.01,1.4], /YSTY, $
	XRANGE=[0.0,1.2], /xsty, $
	TITLE= vik_name+': '+val_name
wait, 0.2  ; to let plot appear...
; - - - 


; - - - Create the 3D array filled with this value vs radius...
;
; Since the 3D array will generally have coaser resolution than
; the input scalar-of-r, need to do some (approximate) integrating
; or averaging over the cell size...  This is done in the
; v3d_sphere routine.
;
; Define the v3d range and resolution
v3d_setup, 1.2, nperradius


; Generate a 3D array filled with Vikram-model values!
val3d = v3d_sphere(0.0, 2.0, RLUP_RAD=rad_scl, RLUP_VAL=val_scl, /SILENT)
val3dhalf = 0.5*v3d_sphere(0.0, 2.0, -2.0, 0.0+2.0*v3d_cube_radius/v3d_nz, $
		RLUP_RAD=rad_scl, RLUP_VAL=val_scl, /SILENT)

; and view it...
;
; halfed, thick, off-axis 3D view:
v3d_project, val3dhalf, RXYZ=[15.0, -50.0, 0.0], /BLACK, $
	/PNG_OUT, name_prefix = vik_out + vik_name+'_'+val_name
;
; full, thin, on-axis view:
v3d_project, val3d, /thin, RXYZ=[0.,0.,0.], $
	/PNG_OUT, name_prefix = vik_out + vik_name+'_'+val_name   ;;;;;;, $
	;;;;  /FITS_OUT, /NORMALIZE_FITS
;
; Make a plot of the "radial" distribution along the x axis:
; Use window 4 to keep it on top of screen to avoid weird Mac problem...
; Set the XYpos values as desired...
iwind = 4
window, iwind, XSIZE=300, YSIZE=220, YPOS=450, XPOS=0
plot_io, v3d_vx, val3d(*, v3d_ny, v3d_nz), YRANGE=[0.01,1.4], /YSTY, /xsty, $
	TITLE= vik_name+': '+val_name
oplot, v3d_vx, (val3d(*, v3d_ny, v3d_nz)) > 0.001, linestyle=1, PSYM=-4
oplot, [0.0,0.0], [0.01, 1.0], linestyle=1
; and compare this to the input radial values:
oplot, rad_scl, (val_scl > 0.001)
oplot, -1.0*rad_scl, (val_scl > 0.001)
; write a png file of this plot
wshow, iwind, iconic=0
image = tvrd()
tvlct, red, green, blue, /GET
png_file = vik_out + vik_name+'_'+val_name+'_profile.png'
wrimg, png_file, image, red, green, blue

RETURN
END