;-----------------------------------------------------------------------
; $Id: stagger_jacobi.pro,v 1.5 2005/12/08 22:34:57 aake Exp $
;  6th order staggered derivatives, 5th order half zone staggering, on
;  a mesh with non-equidistant spacing.
;-----------------------------------------------------------------------
PRO init_stagger
  mesh_init
END
;-----------------------------------------------------------------------
PRO init_derivs
END
;-----------------------------------------------------------------------
PRO scale_dxdn, f, k
@common_cdata
@common_cmesh
  for i=1,nx-1 do begin
    f(i,*,*) = f(i,*,*)*(xm(i)-xm(i-1))^k
  end
  f(1,*,*) = f(1,*,*)*(xm(2)-xm(1))^k
END
;-----------------------------------------------------------------------
PRO scale_dydn, f, k
@common_cdata
@common_cmesh
  for i=1,ny-1 do begin
    f(*,i,*) = f(*,i,*)*(ym(i)-ym(i-1))^k
  end
  f(*,1,*) = f(*,1,*)*(ym(2)-ym(1))^k
END
;-----------------------------------------------------------------------
PRO scale_dzdn, f, k
@common_cdata
@common_cmesh
  for i=1,nz-1 do begin
    f(*,*,i) = f(*,*,i)*(zm(i)-zm(i-1))^k
  end
  f(*,*,1) = f(*,*,1)*(zm(2)-zm(1))^k
END
;-----------------------------------------------------------------------
FUNCTION xup, f 
;
;  f is centered on (i-.5,j,k), fifth order stagger
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=1,shift=1)) $
  + b*(cshift(f,dim=1,shift=-1) + cshift(f,dim=1,shift=2)) $
  + c*(cshift(f,dim=1,shift=-2) + cshift(f,dim=1,shift=3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION yup, f 
;
;  f is centered on (i,j-.5,k)
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=2,shift=1)) $
  + b*(cshift(f,dim=2,shift=-1) + cshift(f,dim=2,shift=2)) $
  + c*(cshift(f,dim=2,shift=-2) + cshift(f,dim=2,shift=3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION zup, f 
;
;  f is centered on (i,j,k-.5)
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=3,shift=1)) $
  + b*(cshift(f,dim=3,shift=-1) + cshift(f,dim=3,shift=2)) $
  + c*(cshift(f,dim=3,shift=-2) + cshift(f,dim=3,shift=3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION xdn, f 
;
;  f is centered on (i-.5,j,k)
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=1,shift=-1)) $
  + b*(cshift(f,dim=1,shift=1) + cshift(f,dim=1,shift=-2)) $
  + c*(cshift(f,dim=1,shift=2) + cshift(f,dim=1,shift=-3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ydn, f 
;
;  f is centered on (i,j,k-.5)
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=2,shift=-1)) $
  + b*(cshift(f,dim=2,shift=1) + cshift(f,dim=2,shift=-2)) $
  + c*(cshift(f,dim=2,shift=2) + cshift(f,dim=2,shift=-3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION zdn, f 
;
;  f is centered on (i,j,k-.5)
;
@common_cdata
;-----------------------------------------------------------------------
  c = 3./256.
  b = -25./256.
  a = .5-b-c
  val = a*(f + cshift(f,dim=3,shift=-1)) $
  + b*(cshift(f,dim=3,shift=1) + cshift(f,dim=3,shift=-2)) $
  + c*(cshift(f,dim=3,shift=2) + cshift(f,dim=3,shift=-3))
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION xdn1, f 
;
;  f is centered on (i-.5,j,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=1,shift=-1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ydn1, f 
;
;  f is centered on (i,j-.5,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=2,shift=-1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION zdn1, f 
;
;  f is centered on (i,j-.5,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=3,shift=-1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION xup1, f 
;
;  f is centered on (i-.5,j,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=1,shift=+1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION yup1, f 
;
;  f is centered on (i,j-.5,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=2,shift=+1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION zup1, f 
;
;  f is centered on (i,j-.5,k)
;
@common_cdata
;-----------------------------------------------------------------------
  val = 0.5*(f + cshift(f,dim=3,shift=+1))
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddxup,f
;
;  X partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
        a *(cshift(f,dim=1,shift=1) - f) $
      + b *(cshift(f,dim=1,shift=2) - cshift(f,dim=1,shift=-1)) $
      + c *(cshift(f,dim=1,shift=3) - cshift(f,dim=1,shift=-2)))
  nx = (size(f))[1]
  for i=0,nx-1 do val[i,*,*]=dxidxup[i]*val[i,*,*]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddyup, f 
;
;  Y partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
    a *(cshift(f,dim=2,shift=1) - f) $
  + b *(cshift(f,dim=2,shift=2) - cshift(f,dim=2,shift=-1)) $
  + c *(cshift(f,dim=2,shift=3) - cshift(f,dim=2,shift=-2)))
  ny = (size(f))[2]
  for i=0,ny-1 do val[*,i,*]=dyidyup[i]*val[*,i,*]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddzup, f 
;
;  Z partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
     a *(cshift(f,dim=3,shift=1) - f) $
   + b *(cshift(f,dim=3,shift=2) - cshift(f,dim=3,shift=-1)) $
   + c *(cshift(f,dim=3,shift=3) - cshift(f,dim=3,shift=-2)))
  nz = (size(f))[3]
  for i=0,nz-1 do val[*,*,i]=dzidzup[i]*val[*,*,i]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddxdn, f 
;
;  X partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
      a*(f - cshift(f,dim=1,shift=-1)) $
    + b*(cshift(f,dim=1,shift=1) - cshift(f,dim=1,shift=-2)) $
    + c*(cshift(f,dim=1,shift=2) - cshift(f,dim=1,shift=-3)))
  nx = (size(f))[1]
  for i=0,nx-1 do val[i,*,*]=dxidxdn[i]*val[i,*,*]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddydn, f 
;
;  Y partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
      a*(f - cshift(f,dim=2,shift=-1)) $
    + b*(cshift(f,dim=2,shift=1) - cshift(f,dim=2,shift=-2)) $
    + c*(cshift(f,dim=2,shift=2) - cshift(f,dim=2,shift=-3)))
  ny = (size(f))[2]
  for i=0,ny-1 do val[*,i,*]=dyidydn[i]*val[*,i,*]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddzdn, f 
;
;  Z partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  c = (-1.+(3.^5-3.)/(3.^3-3.))/(5.^5-5.-5.*(3.^5-3))
  b = (-1.-120.*c)/24.
  a = (1.-3.*b-5.*c)
  val = ( $
      a*(f - cshift(f,dim=3,shift=-1)) $
    + b*(cshift(f,dim=3,shift=1) - cshift(f,dim=3,shift=-2)) $
    + c*(cshift(f,dim=3,shift=2) - cshift(f,dim=3,shift=-3)))
  nz = (size(f))[3]
  for i=0,nz-1 do val[*,*,i]=dzidzdn[i]*val[*,*,i]
  flop_count, 8
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddxup1, f 
;
;  X partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (cshift(f,dim=1,shift=1) - f)
  nx = (size(f))[1]
  for i=0,nx-1 do val[i,*,*]=dxidxup[i]*val[i,*,*]
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddyup1, f 
;
;  Y partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (cshift(f,dim=2,shift=1) - f)
  ny = (size(f))[2]
  for i=0,ny-1 do val[*,i,*]=dyidyup[i]*val[*,i,*]
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddzup1, f 
;
;  Z partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (cshift(f,dim=3,shift=1) - f)
  nz = (size(f))[3]
  for i=0,nz-1 do val[*,*,i]=dzidzup[i]*val[*,*,i]
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddxdn1, f 
;
;  X partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (f - cshift(f,dim=1,shift=-1))
  nx = (size(f))[1]
  for i=0,nx-1 do val[i,*,*]=dxidxdn[i]*val[i,*,*]
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddydn1, f 
;
;  Y partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (f - cshift(f,dim=2,shift=-1))
  ny = (size(f))[2]
  for i=0,ny-1 do val[*,i,*]=dyidydn[i]*val[*,i,*]
  flop_count, 2
  return, val
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION ddzdn1, f 
;
;  Z partial derivative
;
@common_cdata
@common_cmesh
;-----------------------------------------------------------------------
  val = (f - cshift(f,dim=3,shift=-1))
  nz = (size(f))[3]
  for i=0,nz-1 do val[*,*,i]=dzidzdn[i]*val[*,*,i]
  flop_count, 2
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difx3, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=1,shift=-1) - cshift(f,dim=1,shift=2)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION dify3, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=2,shift=-1) - cshift(f,dim=2,shift=2)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difz3, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=3,shift=-1) - cshift(f,dim=3,shift=2)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difx2, f
;
;  Finite difference one zone left to one zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=1,shift=-1) - cshift(f,dim=1,shift=1)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION dify2, f
;
;  Finite difference one zone left to one zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=2,shift=-1) - cshift(f,dim=2,shift=1)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difz2, f
;
;  Finite difference one zone left to one zones right.
;
;-----------------------------------------------------------------------
  val = cshift(f,dim=3,shift=-1) - cshift(f,dim=3,shift=1)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difx1, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = f - cshift(f,dim=1,shift=1)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION dify1, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = f - cshift(f,dim=2,shift=1)
  return, val
END
;-----------------------------------------------------------------------
FUNCTION difz1, f
;
;  Finite difference one zone left to two zones right.
;
;-----------------------------------------------------------------------
  val = f - cshift(f,dim=3,shift=1)
  return, val
END
;-----------------------------------------------------------------------
PRO flop_count, n
  common cflop, nflop
  default, nflop, 0
  nflop = nflop + n
END
;-----------------------------------------------------------------------
PRO flop_print 
  common cflop, nflop
  print,'number of flops:',nflop
  nflop=0
END
;-----------------------------------------------------------------------
PRO mesh_test,n,pack,dpack
@common_cmesh
@common_cdata

  default,n,64
  default,pack,0.7
  default,dpack,0.2
  mesh_init,pack,dpack,n

  x=findgen(nx)*dx
  xd=x-.5*dx
  k=2.*2.*!PI/(nx*dx)

  ws,2
  !p.multi=[0,1,2]
  f =sin(k*x )
  fd=sin(k*xd)
  fdn=xdn(f)
  plot,xd,fd,psy=-1,title='xdn(f) [equidistant]'
  oplot,xd,fdn,psy=4
  plot,xd,fd-fdn,psy=-1

  ws,0
  !p.multi=[0,1,2]
  f =sin(k*xm )
  fd=sin(k*xmdn)
  fdn=xdn(f)
  plot,xmdn,fd,psy=-1,title='xdn(f) [stretched]'
  oplot,xmdn,fdn,psy=4
  plot,xmdn,fd-fdn,psy=-1

  ws,3
  !p.multi=[0,1,2]
  f =cos(k*xm )+sin(k*xm )
  fd=k*(cos(k*xmdn)-sin(k*xmdn))
  fdn=ddxdn(f)
  plot,xmdn,fd,psy=-1,title='ddxdn(f) [stretched]'
  oplot,xmdn,fdn,psy=4
  plot,xmdn,fd-fdn,psy=-1

  ws,1
  !p.multi=[0,1,2]
  f =   sin(k*xm)+cos(k*xm)
  fd=k*(cos(k*xm)-sin(k*xm))
  fdn=ddxup(xdn(f))
  plot,xm,fd,psy=-1,title='ddxup(xdn(f)) [stretched]'
  oplot,xm,fdn,psy=4
  plot,xm,fd-fdn,psy=-1
  
  ws,4
  !p.multi=[0,1,2]
  plot,xm,psym=-1,title='stretched mesh'
  plot,deriv(xm),yst=0,psym=-1,title='mesh derivative'

END