#
# differential operators demo
#
from rsf.proj import *
import sys,wplot
sys.path.append('../../PYUTIL')
from difLOP import *
from DPTEST import *

# ------------------------------------------------------------
do3D = 'y'
doINV = 'n'
# ------------------------------------------------------------

par = dict(
    nx=25,  ox=0, dx=1.0,  lx='x', ux='',
    ny=11,  oy=0, dy=1.0,  ly='y', uy='',
    nz=15,  oz=0, dz=1.0,  lz='z', uz='',
    otype = "c",
    N = 4
    )
wplot.param(par)

# ------------------------------------------------------------
Flow('scl2D',None,
    '''
    spike nsp=1 mag=1 k1=8 k2=13
    n1=%(nz)d o1=%(oz)g d1=%(dz)g label1="" unit1=""
    n2=%(nx)d o2=%(ox)g d2=%(dx)g label2="" unit2=""
    '''%par)
Result('scl2D',wplot.igrey2d('xll=1.9 wantscalebar=y',par))

Flow('vec2D',['scl2D','scl2D'],
    'cat axis=3 space=n ${SOURCES[1:2]}')

Flow('ten2D',['scl2D','scl2D','scl2D'],
    'cat axis=3 space=n ${SOURCES[1:3]}')

if(do3D=='y'):
    Flow('scl3D',None,
    '''
    spike nsp=1 mag=1 k1=8 k2=13 k3=6
    n1=%(nz)d o1=%(oz)g d1=%(dz)g label1="" unit1=""
    n2=%(nx)d o2=%(ox)g d2=%(dx)g label2="" unit2=""
    n3=%(ny)d o3=%(oy)g d3=%(dy)g label3="" unit3=""
    '''%par)
    Result('scl3D',
        'byte gainpanel=a pclip=100|'+wplot.igrey3d('',par))

    Flow('vec3D',['scl3D','scl3D','scl3D'],
        'cat axis=4 space=n ${SOURCES[1:3]}')

    Flow('ten3D',['scl3D','scl3D','scl3D','scl3D','scl3D','scl3D'],
        'cat axis=4 space=n ${SOURCES[1:6]}')
        
# ------------------------------------------------------------
# ------------------------------------------------------------
# ------------------------------------------------------------
# ------------------------------------------------------------

for j in (['x','z','l','g','d','c','v','h','i','k']):
#for j in (['i']):

    if(j=='l' or j=='x' or j=='z' or j=='g' or j=='v'):
        Flow(j+'2D','scl2D','window')
    elif (j=='i'):
        Flow(j+'2D','ten2D','window')
    else:
        Flow(j+'2D','vec2D','window')
    
    O2D = dif2d(j,par['otype'],par['N'])
    O2D.FORW(j+'2D',   j+'2Dfor')
    O2D.ADJT(j+'2Dadj',j+'2Dfor')

    for i in (['for','adj']):
        Result(j+'2D'+i,
            wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))
        
    D2D = DPTEST(O2D,[j+'2Dadj'],[j+'2Dfor'])           
    D2D.RUN()

if(do3D=='y'):
    for j in (['x','y','z','l','g','d','c','v','h','i','k']):
#    for j in (['i']):
            
        if(j=='l' or j=='x' or j=='y' or j=='z' or j=='g' or j=='v'):
            Flow(j+'3D','scl3D','window')
        elif(j=='i'):
            Flow(j+'3D','ten3D','window')
        else:
            Flow(j+'3D','vec3D','window')
        
        O3D = dif3d(j,par['otype'],par['N'])
        O3D.FORW(j+'3D',   j+'3Dfor')
        O3D.ADJT(j+'3Dadj',j+'3Dfor')

        for i in (['for','adj']):
            Result(j+'3D'+i,
                'byte gainpanel=a pclip=100|'
                + wplot.igrey4d('color=E',par))
                    
        D3D = DPTEST(O3D,[j+'3Dadj'],[j+'3Dfor'])           
        D3D.RUN()

# ------------------------------------------------------------
# Laplacian
# ------------------------------------------------------------
L2D = lapl2d(par['otype'],par['N'])
L2D.FORW('scl2D', 'L2Dfor')
L2D.ADJT('L2Dadj','L2Dfor')
Result('L2Dfor',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))
Result('L2Dadj',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))

DL2D = DPTEST(L2D,['L2Dadj'],['L2Dfor'])           
DL2D.RUN()

if(do3D=='y'):
    L3D = lapl3d(par['otype'],par['N'])
    L3D.FORW('scl3D', 'L3Dfor')
    L3D.ADJT('L3Dadj','L3Dfor')
    Result('L3Dfor','byte gainpanel=a pclip=100|'+wplot.igrey3d('color=E',par))
    Result('L3Dadj','byte gainpanel=a pclip=100|'+wplot.igrey3d('color=E',par))
    
    DL3D = DPTEST(L3D,['L3Dadj'],['L3Dfor'])           
    DL3D.RUN()

# ------------------------------------------------------------
# gradient
# ------------------------------------------------------------
G2D = grad2d(par['otype'],par['N'])
G2D.FORW('scl2D', 'G2Dfor')
G2D.ADJT('G2Dadj','G2Dfor')
Result('G2Dfor',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))
Result('G2Dadj',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))

DG2D = DPTEST(G2D,['G2Dadj'],['G2Dfor'])           
DG2D.RUN()

if(do3D=='y'):
    G3D = grad3d(par['otype'],par['N'])
    G3D.FORW('scl3D', 'G3Dfor')
    G3D.ADJT('G3Dadj','G3Dfor')
    Result('G3Dfor','byte gainpanel=a pclip=100|'+wplot.igrey4d('color=E',par))
    Result('G3Dadj','byte gainpanel=a pclip=100|'+wplot.igrey3d('color=E',par))
    
    DG3D = DPTEST(G3D,['G3Dadj'],['G3Dfor'])           
    DG3D.RUN()
 
# ------------------------------------------------------------
# divergence
# ------------------------------------------------------------
D2D = div2d(par['otype'],par['N'])
D2D.FORW('vec2D', 'D2Dfor')
D2D.ADJT('D2Dadj','D2Dfor')
Result('D2Dfor',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))
Result('D2Dadj',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))

DD2D = DPTEST(D2D,['D2Dadj'],['D2Dfor'])           
DD2D.RUN()

if(do3D=='y'):
    D3D = div3d(par['otype'],par['N'])
    D3D.FORW('vec3D', 'D3Dfor')
    D3D.ADJT('D3Dadj','D3Dfor')
    Result('D3Dfor','byte gainpanel=a pclip=100|'+wplot.igrey3d('color=E',par))
    Result('D3Dadj','byte gainpanel=a pclip=100|'+wplot.igrey4d('color=E',par))
    
    DD3D = DPTEST(D3D,['D3Dadj'],['D3Dfor'])           
    DD3D.RUN()

# ------------------------------------------------------------
# curl
# ------------------------------------------------------------
C2D = curl2d(par['otype'],par['N'])
C2D.FORW('vec2D', 'C2Dfor')
C2D.ADJT('C2Dadj','C2Dfor')
Result('C2Dfor',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))
Result('C2Dadj',wplot.igrey2d('xll=1.9 wantscalebar=y color=E',par))

DC2D = DPTEST(C2D,['C2Dadj'],['C2Dfor'])           
DC2D.RUN()

if(do3D=='y'):
    C3D = curl3d(par['otype'],par['N'])
    C3D.FORW('vec3D', 'C3Dfor')
    C3D.ADJT('C3Dadj','C3Dfor')
    Result('C3Dfor','byte gainpanel=a pclip=100|'+wplot.igrey4d('color=E',par))
    Result('C3Dadj','byte gainpanel=a pclip=100|'+wplot.igrey4d('color=E',par))
    
    DC3D = DPTEST(C3D,['C3Dadj'],['C3Dfor'])           
    DC3D.RUN()
    
End()
