# Author: G. Hennenfent
#         Seismic Laboratory for Imaging and Modeling
#         Department of Earch & Ocean Sciences
#         The University of British Columbia
#         
# Date  : July, 07

from rsfproj import *

# define customized plotting function
def mygrey(title,label1="Time",label2="Offset",unit1="s",unit2="m"):
    return '''
    sfgrey title="%s" minval=-.005 maxval=.005 clip=.0008
    label1="%s" unit1="%s" label2="%s" unit2="%s" crowd1=.5 crowd2=.75
    '''%(title,label1,unit1,label2,unit2)
def myfkgrey(title,label1="Time",label2="Offset",unit1="s",unit2="m"):
    return '''
    sfgrey title="%s" minval=-.005 maxval=.005 clip=.0008
    label1="%s" unit1="%s" label2="%s" unit2="%s" crowd1=.5 crowd2=.25
    '''%(title,label1,unit1,label2,unit2)

# define SLIM FTP server information
FTPserver = {
    'server': 'slim.eos.ubc.ca',
    'login':'anonymous',
    'password':'ghennenfent@eos.ubc.ca'}

# define dataset folder on FTP server (/data/hegilles) and dataset
# name
logname = 'synth/delphi'
dataset = 'shotdeep.p.rsf'

# fetch data from FTP server
Fetch(dataset,logname,FTPserver)
Flow('data_12p5m',dataset,
     '''sfdd form=native | sfcostaper nw1=30 nw2=30 |
     sfbandpass fhi=60 | sfwindow j1=2''')
Result('data_12p5m',mygrey("") )

Flow('fkdata_12p5m','data_12p5m','sfspectra2 | sfwindow max1=60')
Result('fkdata_12p5m',myfkgrey("",label1="Frequency",
                       label2="Wavenumber",unit1="Hz",unit2="1/m") )

jump = 3
name = 'data_subREG'
fkname = 'fkdata_subREG'
Flow([name,'mask_subREG2D'],'data_12p5m',
     'sfwindow j2=%d | sflpad mask=${TARGETS[1]} jump=%d | sfwindow n2=481'%(jump,jump) )
mask = 'mask_subREG'
Flow(mask,'mask_subREG2D',
     'sftransp | sfwindow n2=1 n1=481 | sfput n3=1 | sfdd type=int')

Flow(fkname,name,'sfspectra2 | sfwindow max1=60')
Result(name,mygrey("") )
Result(fkname,myfkgrey("",
                     label1="Frequency",label2="Wavenumber",
                     unit1="Hz",unit2="1/m") )
name2 = 'data_subIRREG'
fkname2 = 'fkdata_subIRREG'
mask2 = 'mask_subIRREG'
Flow(mask2,'data_12p5m',
     'sfjitter perc=%f seed=1982 | sfdd type=int'%(1-1./jump) )
Flow(name2,['data_12p5m',mask2],
     'sfheadercut mask=${SOURCES[1]}')
Flow(fkname2,name2,'sfspectra2 | sfwindow max1=60')
Result(name2,mygrey("") )
Result(fkname2,myfkgrey("",
                      label1="Frequency",label2="Wavenumber",
                      unit1="Hz",unit2="1/m") )
name3 = 'data_subRAND'
fkname3 = 'fkdata_subRAND'
mask3 = 'mask_subRAND'
Flow(mask3,'data_12p5m',
     'sfput n3=1 | sfshotholes perc=%f'%(1-1./jump) )
Flow(name3,['data_12p5m',mask3],
     'sfheadercut mask=${SOURCES[1]}')
Flow(fkname3,name3,'sfspectra2 | sfwindow max1=60')
Result(name3,mygrey("") )
Result(fkname3,myfkgrey("",
                      label1="Frequency",label2="Wavenumber",
                      unit1="Hz",unit2="1/m") )

interp = WhereIs('interpolation.py')

# INTERPOLATION BY CRSI
# curvelet transform parameters
nbs = 7         # nb of scales
nba = 64        # nb of angles at 2nd coarsest scale
ac = 1          # all curvelets
# solver parameters
nbout = 50      # nb of optimization subproblems (SPs)
nbin = 5        # nb of itr used to approx. solution of SP
# sparsity constraint parameters
thrmax = .01    # perc of curvelet coeff. kept for 1st SP
thrmin = .99    # perc of curvelet coeff. kept for last SP
# angular constraint parameters
angmaxL = 90    # max dip (degres) for negative wave numbers
angmaxR = 90    # max dip (degres) for positive wave numbers

# look for interpolation.py in your path (see SLIMpy installation
# instructions)
for mode in ['REG','IRREG','RAND']:
    # interpolate
    res = 'crsi_sub%(mode)s'%vars()
    Flow(res,['data_sub%(mode)s'%vars(),'mask_sub%(mode)s'%vars(),interp],
         interp + '''
         -q mask=${SOURCES[1]} data=$SOURCE output=$TARGET
         solverparams=%(nbout)d,%(nbin)d transparams=%(nbs)d,%(nba)d,%(ac)d
         thrparams=%(thrmax)f,%(thrmin)f angconst=%(angmaxL)f,%(angmaxR)f
         '''%vars(),stdin=0,stdout=-1)
    Result(res,mygrey(""))
    fkres = 'fk'+res
    Flow(fkres,res,'sfspectra2 | sfwindow max1=60')
    Result(fkres,myfkgrey("",
                          label1="Frequency",label2="Wavenumber",
                          unit1="Hz",unit2="1/m") )

End()