pro zoom_proj, file_name1, proj_dir, sds_num1, sds_num2, MIN=min_in, MAX=max_in, gif=gif, ps=ps
;
; purpose: produce multiple images of projection data 
;          with progressive zoom-in's
;
; usage: zoom_proj, 'projection_file_name', 'projection_direction', field_num x 2
;    where field_num corresponds to the quantity projected:
;     1 - density           2 - X-ray luminosity (free-free)
;     3 - dark matter density
;     4 - temperature (luminosity-weighted)  (MKRJ  switched 3 & 4)
;     5 - level structure
;     6 - SZ (thermal)      7 - SZ (kinematic)
;     8 - metallicity       9 - stellar density (if star formation on)
;
;   optional parameter: MIN=min_val, MAX=max_val, gif, ps ; added by MKRJ
;
;  (e.g. zoom_proj, 'amr_sim', 'z', 1, 4, MIN=10, MAX=15, /gif)
;  (note: if field_num is negative, then it displays a colour bar)
;

; MKRJ 11/9/06 
;      loop over file_names's and sds_num's
; **********************************
  FOR j=1, 6 DO BEGIN

    if (j eq 1) then begin
        sds_num = sds_num1
        file_name=file_name1+'_'+proj_dir+'.project'
    endif
    if (j eq 2) then file_name=file_name1+'_center_'+proj_dir+'.project'
    if (j eq 3) then file_name=file_name1+'_galaxy_'+proj_dir+'.project'
    if (j eq 4) then begin
        sds_num = sds_num2
        file_name=file_name1+'_'+proj_dir+'.project'
    endif
    if (j eq 5) then file_name=file_name1+'_center_'+proj_dir+'.project'
    if (j eq 6) then file_name=file_name1+'_galaxy_'+proj_dir+'.project'

    icolor_bar = 0
    IF (sds_num lt 0) THEN BEGIN
        sds_num = -sds_num
        icolor_bar = 1
    ENDIF
 
    DFSD_SETINFO, /RESTART
    DFSD_GETINFO, file_name, NSDS=NumSDS

    IF (sds_num gt numSDS) THEN BEGIN
        print, ' read_hdf: Cannot read sds#', sds_num
        print, ' there are only ', numsds,' scientific data sets in ', file_name
        print, ' returning zero valued fltarr'
        data = fltarr(3,3,3)
        RETURN
    ENDIF

    FOR i=1,sds_num DO BEGIN
        r=0.
        l=0.
        HDF_DFSD_GETINFO, file_name, label = l, dims=d, type=t,range=r
;        print,i, l
    ENDFOR
    print, ' read ',l, ' data from ', file_name
    hdf_DFSD_GETDATA, file_name, data

    data_dim = size(data)
    print, 'image info: ',data_dim

    if (sds_num eq 1) then data2 = ALOG10(data)
    if (sds_num eq 2) then data2 = ALOG10(data)
    if (sds_num eq 3) then data2 = ALOG10(data+1.0e12)
    if (sds_num eq 4) then data2 = ALOG10(data)
    if (sds_num eq 5) then data2 = data
    if (sds_num eq 6) then data2 = data
    if (sds_num eq 7) then data2 = data
    if (sds_num eq 8) then data2 = ALOG10(data)
    if (sds_num eq 9) then data2 = ALOG10(data+1.0e9)

    min1 = MIN(data2)
    max1 = MAX(data2)

    n1 = data_dim(1)
    n2 = data_dim(2)

    print, 'min, max =', MIN(data2), MAX(data2)

    IF (icolor_bar eq 1) THEN BEGIN
      nadd = 90    ; MKRJ  11/9/06 initially nadd=40
      nsep = 150   ; MKRJ  added part about nsep
      n3 = n2+nadd+nsep
      data3 = FLTARR(n1, n3)
      data3(*, *) = min1
      print, size(data3)
      data3(0:n1-1, (nadd+nsep):n3-1) = data2

      FOR i=1, n1 DO BEGIN
        data3(i-1, 0:nadd-5) = min1 + (i-0.5)*(max1-min1)/n1
      ENDFOR

      data2 = data3
      n2 = n3
  
    ENDIF

    if (j eq 1) then begin
      data2_1 = data2
      min_a = min1
      max_a = max1
    endif
    if (j eq 2) then data2_2 = data2
    if (j eq 3) then data2_3 = data2
    if (j eq 4) then begin
      data2_4 = data2
      min_b = min1
      max_b = max1
    endif
    if (j eq 5) then data2_5 = data2
    if (j eq 6) then data2_6 = data2

  ENDFOR
; **********************************

  n1 = 2048
  n2 = 2048
  ngap = 200

  data4 = fltarr(3*n1+2*ngap, 2*n2+ngap)
  data4(*, *) = 0.0
  print, size(data4)

  data4(0:n1-1, 0:n2-1) = (data2_4(0:n1-1, 0:n2-1)-min_b)/(max_b-min_b)
  data4(n1+ngap:2*n1+ngap-1, 0:n2-1) = (data2_5(0:n1-1, 0:n2-1)-min_b)/(max_b-min_b)
  data4(2*n1+2*ngap:3*n1+2*ngap-1, 0:n2-1) = (data2_6(0:n1-1, 0:n2-1)-min_b)/(max_b-min_b)

  data4(0:n1-1, n2+ngap:2*n2+ngap-1) = (data2_1(0:n1-1, 0:n2-1)-min_a)/(max_a-min_a)
  data4(n1+ngap:2*n1+ngap-1, n2+ngap:2*n2+ngap-1) = (data2_2(0:n1-1, 0:n2-1)-min_a)/(max_a-min_a)
  data4(2*n1+2*ngap:3*n1+2*ngap-1, n2+ngap:2*n2+ngap-1) = (data2_3(0:n1-1, 0:n2-1)-min_a)/(max_a-min_a)

  IF KEYWORD_SET(min_in) THEN BEGIN
    min1 = min_in
  ENDIF
  IF KEYWORD_SET(max_in) THEN BEGIN
    max1 = max_in
  ENDIF

  byte_image = BYTSCL(data4, MIN=0.0, MAX=1.0, TOP=!D.TABLE_SIZE)

  WINDOW, XSIZE=3*n1+2*ngap, YSIZE=2*n2+ngap
  LOADCT, 4, NCOLORS=256, BOTTOM=0
  TV, byte_image

;  show3, data2
;  surface, data2
;  contour, data2

; MKRJ  produce gif image file, if requested
  IF KEYWORD_SET(gif) then begin
    gif_file = 'zoom_proj_'+file_name1+'_'+proj_dir+'.gif'
    write_gif, gif_file, byte_image
  ENDIF

; MKRJ  produce ps image file, if requested
  IF KEYWORD_SET(ps) then begin
    set_plot, 'ps'
    ps_file = 'zoom_proj_'+file_name1+'_'+proj_dir+'.ps'
    device, bits_per_pixel=8, /color, filename=ps_file
;, $
;            xsize=3*n1+2*ngap, ysize=2*n1+ngap
    TV, byte_image
    device, /close
    set_plot, 'x'
  ENDIF

;  if (sds_num eq 3) then $
;    TVSCL, ALOG10(REBIN(data>1.e12, data_dim(1)*2,data_dim(2)*2, /SAMPLE))+12. $
;  else $
;    TVSCL, ALOG10(REBIN(data, data_dim(1)*2,data_dim(2)*2, /SAMPLE))

END