pro get_SN_shock,field,tmin=tmin,tmax=tmax,min=amin,max=amax, $
njump=njump,datadir=datadir,OLDFILE=OLDFILE,$
proc=proc,exp=exp,mesh=mesh, NOPLOT=NOPLOT, transform=transform, $
QUIET=QUIET, object=shockinfo, POWERLAW=POWERLAW, tolerance=tolerance, $
midpnt=midpnt
;
; $Id$
;
;  Reads in 4 slices as they are generated by the pencil code.
;  The variable "field" can be changed. Default is 'lnrho'.
;
;  if the keyword /mpeg is given, the file movie.mpg is written.
;  tmin is the time after which data are written
;  nrepeat is the number of repeated images (to slow down movie)
;
;  Typical calling sequence
;  rvid_box,'bz',tmin=190,tmax=200,min=-.35,max=.35,/mpeg
;
default,extension,'.xz'
default,amax,.05
default,amin,-amax
default,field,'lnrho'
default,datadir,'data'
default,nrepeat,0
default,njump,0
default,tmin,0.
default,tmax,1e38
default,wait,.03
;
if keyword_set(proc) then begin
file_slice=datadir+'/'+proc+'/slice_'+field+extension
endif else begin
file_slice=datadir+'/slice_'+field+extension
endelse
;
;  Read the dimensions and precision (single or double) from dim.dat
;
mx=0L & my=0L & mz=0L & 
nghostx=0L & nghosty=0L & nghostz=0L
;

pc_read_dim,object=dim,proc=proc,QUIET=QUIET
mx=dim.mx
my=dim.my
mz=dim.mz
nghostx=dim.nghostx
nghosty=dim.nghosty
nghostz=dim.nghostz
;
nx=dim.nx
ny=dim.ny
nz=dim.nz
;
;
  if keyword_set(mesh) then begin
zaxisscale=(findgen(nz)-(nz/2))
  endif else begin
  pc_read_grid,obj=gridobj,proc=proc,QUIET=QUIET,/TRIMXYZ
  zaxisscale=gridobj.z
  ;  zaxisscale=zaxisscale[nghostz:nghostz+nz-1]
  endelse
  ;
  ;
  ;
  ;
  t=0. & islice=0
  axz=fltarr(nx,nz)
  slice_xpos=0.
  slice_ypos=0.
  slice_zpos=0.
  slice_z2pos=0.
  ;
  ;  open MPEG file, if keyword is set
  ;
  dev='x' ;(default)
  ;
  ;  allow for jumping over njump time slices
  ;  initialize counter
  ;
  ijump=njump ;(make sure the first one is written)
  ;
  it=0
  map=fltarr(nx,10000) 
  timescale=fltarr(10000) 
shockinfo = fltarr(10000,3)
  print,'      it        t      left shock   right shock '
  close,99
  openw,99,'shock.txt'

  close,1 & openr,1,file_slice,/f77
  while not eof(1) do begin
  if keyword_set(OLDFILE) then begin ; For files without position
  readu,1,axz,t
  endif else begin
  readu,1,axz,t,slice_z2pos
  endelse

  if n_elements(transform) ne 0 then begin
res= execute(transform)
  endif

  default,midpnt,((size(axz))[2])/2
  if it eq 0 then tt=t else tt=[tt,t]
  map[*,it]=axz[*,midpnt]
  timescale[it]=t
  it=it+1L
;
  if t ge tmin and t le tmax then begin
    if ijump eq njump then begin
      find_SN_shock,axz[*,midpnt], zaxisscale,leftpnt=shock_left, rightpnt=shock_right, $
                             tol=tolerance
      shockinfo[it,0]=t
      shockinfo[it,1]=shock_left
      shockinfo[it,2]=shock_right
      if not keyword_set(noplot) then begin
        if keyword_set(exp) then begin
          plot,zaxisscale,exp(axz[*,midpnt]),yr=[amin,amax] 
          print,islice,t,exp(min([axz])),exp(max([axz])),shock_left,shock_right
        endif else begin
          plot,zaxisscale,axz[*,midpnt],yr=[amin,amax]
          print,islice,t,min([axz]),max([axz]),shock_left,shock_right
        endelse
      endif
      printf,99,islice,shock_left,shock_right 
;      ijump=0
    endif else begin
      ijump=ijump+1
    endelse
  endif
  islice=islice+1
endwhile
close,1
close,99
;
;  reform map appropriately
;
;nxz=n_elements(axz) & nt=it
map=map[*,0:it-1]
timescale=timescale[0:it-1]
shockinfo=shockinfo[0:it-1,*]
;map=reform(map,nx,it)

contour,transpose(map),timescale,zaxisscale,/fil,nlev=60

if keyword_set(powerlaw) then begin
default,tail,5
notail1=min([tail,it-1])
notail2=max([it-tail,notail1+1])
plot_oo,shockinfo[notail1+1:notail2,0],abs(shockinfo[notail1+1:notail2,2])

powerlaw=(alog(shockinfo[notail1+1:notail2,2])-alog(shockinfo[notail1,2])) $ 
     /(alog(shockinfo[notail1+1:notail2,0])-alog(shockinfo[notail1,0]))

print,"Using ", (notail2-notail1), " points :-"
print,"Mean for power law: ", mean(powerlaw)
print,"Variance power law: ", mean(powerlaw^2)-mean(powerlaw)^2
endif
END