;+
; $Id: read_prf.pro,v 1.8 2011/04/07 01:41:09 andersl Exp $

 FUNCTION read_prf, file, region=region
;
; return,{nwaves:nwaves, waves:waves, nprf:nprf, nstokes:nstokes, nx:nx, ny:ny, prf:prf}
;-

  ;if (size(file))[1] ne 7 then file='snapshot_'+string(name,format='(i5.5)')+'.prf'
  openr,u,/get,file,/f77
  
  key=0L
  readu,u,key
  print,'key:',key

  if (key gt 0) then begin
    ntot=key
    nregion=0L
    i0=0L
    i1=0L
  end else if (key eq -1) then begin
    nregion=0L
    readu,u,nregion
    print,'nregion:',nregion
    nwaves=lonarr(nregion)
    i0=lonarr(nregion)
    i1=lonarr(nregion)
    ntot=0
    for ir=0,nregion-1 do begin
      nw=0L
      readu,u,nw
      print,'nw:',nw
      nwaves[ir]=nw
      i0[ir]=ntot
      i1[ir]=ntot+nw-1
      ntot=ntot+nw
    end
    ntot=long(total(nwaves,/integer))
  end

  nwaves=ntot
  waves=fltarr(ntot)
  readu,u,waves

  nstokes=0L & nprf=0L & nx=0L & ny=0L
  readu,u,nstokes,nprf,nx,ny
  prf=fltarr(nprf,nstokes,nx,ny)
  prof=fltarr(nstokes,nprf,nx)
  
  print,'nprf, nx, ny =', nprf, nx, ny
  for iy=0,ny-1 do begin
    readu,u,prof
    prf[*,*,*,iy] = transpose(prof,[1,0,2])
  end

  if n_elements(region) gt 0 then begin
    nregion=1
    nwaves=nwaves[region-1]
    nprf=nwaves
    i0=i0[region-1]
    i1=i1[region-1]
    waves=waves[i0:i1]
    prf=prf[i0:i1,*,*,*]
  end

  free_lun, u
  return,{nregion:nregion, i0:i0, i1:i1, nwaves:nwaves, waves:waves, nprf:nprf, nstokes:nstokes, nx:nx, ny:ny, prf:prf}
END