def update_model(self):
        if not self._model or self._model[0] != self.chosen_model:
            if self.chosen_model.startswith('Cake builtin: '):
                load_model = cake.load_model(
                    self.chosen_model.split(': ', 1)[1])

            elif self.chosen_model.startswith('GF Store: '):
                store_id = self.chosen_model.split(': ', 1)[1]

                load_model = self._engine.get_store(store_id)\
                    .config.earthmodel_1d
            else:
                load_model = cake.load_model(self.chosen_model)

            self._model = (self.chosen_model, load_model)

def init(store_dir, variant):
    if variant is None:
        variant = '2006'

    if variant not in ('2006', '2006a'):
        raise gf.store.StoreError('unsupported variant: %s' % variant)

    modelling_code_id = 'qseis.%s' % variant

    qseis = QSeisConfig(qseis_version=variant)
    qseis.time_region = (
        gf.meta.Timing('begin-50'),
        gf.meta.Timing('end+100'))

    qseis.cut = (
        gf.meta.Timing('begin-50'),
        gf.meta.Timing('end+100'))

    qseis.wavelet_duration_samples = 0.001
    qseis.sw_flat_earth_transform = 1

    store_id = os.path.basename(os.path.realpath(store_dir))

    config = gf.meta.ConfigTypeA(

        id=store_id,
        ncomponents=10,
        sample_rate=0.2,
        receiver_depth=0*km,
        source_depth_min=10*km,
        source_depth_max=20*km,
        source_depth_delta=10*km,
        distance_min=100*km,
        distance_max=1000*km,
        distance_delta=10*km,
        earthmodel_1d=cake.load_model().extract(depth_max='cmb'),
        modelling_code_id=modelling_code_id,
        tabulated_phases=[
            gf.meta.TPDef(
                id='begin',
                definition='p,P,p\\,P\\,Pv_(cmb)p'),
            gf.meta.TPDef(
                id='end',
                definition='2.5'),
            gf.meta.TPDef(
                id='P',
                definition='!P'),
            gf.meta.TPDef(
                id='S',
                definition='!S'),
            gf.meta.TPDef(
                id='p',
                definition='!p'),
            gf.meta.TPDef(
                id='s',
                definition='!s')])

    config.validate()
    return gf.store.Store.create_editables(
        store_dir, config=config, extra={'qseis': qseis})

    def _create_regional_ttt_store(self):

        conf = gf.ConfigTypeA(
            id='empty_regional',
            source_depth_min=0.,
            source_depth_max=20*km,
            source_depth_delta=10*km,
            distance_min=1000*km,
            distance_max=2000*km,
            distance_delta=10*km,
            sample_rate=2.0,
            ncomponents=10,
            earthmodel_1d=cake.load_model(),
            tabulated_phases=[
                gf.TPDef(id=id, definition=defi) for (id, defi) in [
                    ('depthp', 'p'),
                    ('pS', 'pS'),
                    ('P', 'P'),
                    ('S', 'S')
                ]
            ])

        store_dir = mkdtemp(prefix='gfstore')
        self.tempdirs.append(store_dir)

        gf.Store.create(store_dir, config=conf)
        store = gf.Store(store_dir)
        store.make_ttt()

        store.close()
        return store_dir

    def test_angles(self):
        mod = cake.load_model()
        data = [
            [1.0*km, 1.0*km, 1.0*km, 90., 90., 'P'],
            [1.0*km, 2.0*km, 1.0*km, 45., 135., 'P\\'],
            [2.0*km, 1.0*km, 1.0*km, 135., 45., 'p'],
            [1.0*km, 2.0*km, math.sqrt(3.)*km, 60., 120., 'P\\'],
            [2.0*km, 1.0*km, math.sqrt(3.)*km, 120., 60., 'p']]

        for (zstart, zstop, dist, takeoff_want, incidence_want, pdef_want) \
                in data:

            rays = mod.arrivals(
                zstart=zstart,
                zstop=zstop,
                phases=[cake.PhaseDef(sphase)
                        for sphase in 'P,p,P\\,p\\'.split(',')],
                distances=[dist*cake.m2d])

            for ray in rays:
                takeoff = round(ray.takeoff_angle())
                incidence = round(ray.incidence_angle())
                pdef = ray.used_phase().definition()

                assert takeoff == takeoff_want
                assert incidence == incidence_want
                assert pdef == pdef_want

 def update_model(self):
     if not self._model or self._model[0] != self.chosen_model:
         if self.chosen_model in cake.builtin_models():
             load_model = cake.load_model(self.chosen_model)
         else:
             load_model = self._engine.get_store(self.chosen_model).config.earthmodel_1d 
         self._model = (self.chosen_model, load_model)

    def runParallel(inmodel):
        
        print 'does not work properly. EXIT'
        sys.exit(0)
        ProBar = progressbar.ProgressBar(maxval=iterations).start()
        misfits={}
        misfits['pMF']=[]; misfits['sMF']=[]
        misfits['ScsMF']=[]; misfits['ScssMF']=[]

        loadmod=cake.load_model(inmodel)
     
        for latindx, lat in enumerate(_lats):
            for lonindx, lon in enumerate(_lons):
                for zindex, z in enumerate(_depths):
                    #iteration+=1
                    # Start prozess with one event (depth), and one model:
                    eve=model.Event(lat,lon,str_to_time("2010-04-11 22:08:15.500"), 
                       "Spain_Durcal" , z,6.3)
                    [ttpdiff, ttsdiff, ttScsdiff, ttScssdiff] = depthfinder.startup(loadmod, eve, maxdist) 
                   
                    pMF, sMF, ScsMF, ScssMF= map(
                               lambda x: calculateMisfit(x,maxdist), [ttpdiff,ttsdiff,ttScsdiff,ttScssdiff])
                                               
                    resultArray[latindx][lonindx][zindex] = [pMF, sMF, ScsMF, ScssMF]
                    # update progressbar
                    ProBar.update(iteration)
                    identifierstring = inmodel+'.%s.%s.%s'%(lat, lon, z)
                   
                    results[identifierstring]=misfits
        
        try:
            output = open('results.p','w')
            pickle.dump(results, output)
        finally:
            output.close()

 def setTshiftFile(self, reference_event, latindx, lonindx, results, modelstr, phase, depths, maxdist=None):
     '''
     Creates a time shift file relative to a given reference_event.
     
     :param reference_event: time shift will be analyzed relative to this event
     :param latindx:         index of latitude of the reference_events' origin within the results array
     :param lonindx:         index of longitude of the reference_events' origin within the results array
     :param results:         Results as numpy nd array
     :param modelstr:        string modelname
     :param phase:           string phasename
     :param depths:          list of depths which were used in the processing
     :param maxdist:         maximum epicentral distance up to which the time shift is to be analyzed
     '''
     if (type(modelstr)==str):
         modelload=cake.load_model(modelstr)
         modelname=(modelstr.split('/')[1]).split('.')[0]
     if (type(phase)==str):
         phaseload=cake.PhaseDef(phase)
     # get best fit index:
     minimum, minimumIndex = decomposer.getMinimumMisfitForLatindxLonindx(results, 
             modelstr, phase, latindx, lonindx)        
     # set depth of reference event to depth of best fit:
     reference_event.depth=depths[minimumIndex]
     # calculate traveltimes for best fit:
     for card in self._myStationBox.getStationCards:
         theophases = self.calculateTheoreticalArrivals(modelload, reference_event, card)
         if theophases[phase] and card._Phases[phase]: # and card.getDistance2Event(reference_event)<=maxdist:           
             tshift = card._Phases[phase]-theophases[phase].t
             self._myStationBox.setUsedStations(card, tshift)
         else:
             pass
     # create input file for shift plot
     self._myStationBox.writeUsedStationsFile(val='{0}{1}'.format(phase,modelname))   

 def example():
     conf = QSeisRConfigFull()
     conf.source = QSeis2dSource(lat=-80.5, lon=90.1)
     conf.receiver_location = QSeisRReceiver(lat=13.4, lon=240.5, depth=0.0)
     conf.time_reduction = 10.0
     conf.earthmodel_receiver_1d = cake.load_model().extract(
         depth_max='moho')
     return conf

    def runSerial(models=None):
        '''
        Execute serial processing (1 CPU). 
        :param models: list of models to investigate
                if no models are given, all models will be investigated.
        '''

        if models==None:
            models2use = _models
        else:
            models2use = [models]
        
        iteration=0
        iterations = len(models2use)*len(_depths)*len(_lats)*len(_lons)
        sys.stdout.write('calculating misfits... ')
        ProBar = progressbar.ProgressBar(maxval=iterations).start()

        # instantiate result array as numpy nd array:
        resultArray = np.ndarray(shape=(len(_lats), len(_lons), len(_depths), 4)) 

        
        for mod in models2use:
            misfits={}
            misfits['pMF']=[]; misfits['sMF']=[]
            misfits['ScsMF']=[]; misfits['ScssMF']=[]
 
            loadmod=cake.load_model(mod)
            
            for latindx, lat in enumerate(_lats):
                for lonindx, lon in enumerate(_lons):
                    for zindex, z in enumerate(_depths):
                        iteration+=1
                        eve=model.Event(lat,lon,str_to_time("2010-04-11 22:08:15.500"), 
                           "Spain_Durcal" , z,6.3)
                        [ttpdiff, ttsdiff, ttScsdiff, ttScssdiff] = depthfinder.startup(loadmod, eve, maxdist) 
                        
                        [pMF, sMF, ScsMF, ScssMF]= map(
                                    lambda x: calculateMisfit(x,maxdist), [ttpdiff,ttsdiff,ttScsdiff,ttScssdiff])
                        # update progressbar
                        ProBar.update(iteration)

                        # write data to numpy array:
                        resultArray[latindx][lonindx][zindex] = [pMF, sMF, ScsMF, ScssMF]
                        
            results[mod]=resultArray
            depthfinder.storeStationBox()
        
        # finish progressbar:
        ProBar.finish()
        # write dict to pickled data:
        try:
            output = open('numpy_results.p','w')
            pickle.dump(results, output)
        finally:
            output.close()
       
        # write used stations file:
        depthfinder._myStationBox.writeUsedStationsFile()

 def example():
     conf = QSeisConfigFull()
     conf.receiver_distances = [ 2000. ]
     conf.receiver_azimuths = [ 0. ]
     conf.time_start = -10.0
     conf.time_reduction_velocity = 15.0
     conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb')
     conf.sw_flat_earth_transform = 1
     return conf

    def test_earthmodel(self):
        from pyrocko import cake
        mod = cake.load_model()

        nx, ny = 500, 500
        delta = (cake.earthradius * 2.0) / (nx-1)

        x = num.arange(nx) * delta - cake.earthradius
        y = num.arange(ny) * delta - cake.earthradius

        x2 = x[num.newaxis, :]
        y2 = y[:, num.newaxis]

        z = cake.earthradius - num.sqrt(x2**2 + y2**2)

        vp_pro = mod.profile('vp')
        z_pro = mod.profile('z')

        vp = num.interp(z, z_pro, vp_pro)

        inside = z > 0.0

        speeds = num.ones((ny, nx))
        speeds[:, :] = 300.
        speeds[inside] = vp[inside]

        times = num.zeros((ny, nx)) - 1.0

        iy = ny - int(round(600*km / delta))
        ix = nx//2

        times[iy, ix] = 0.0

        @benchmark.labeled('test_earthmodel')
        def run():
            eikonal_ext.eikonal_solver_fmm_cartesian(speeds, times, delta)

        run()

        times[num.logical_not(inside)] = num.nan
        self.compare_with_reference(times, 'test_earthmodel.npy')

        if show_plot:
            from matplotlib import pyplot as plt

            plt.gcf().add_subplot(1, 1, 1, aspect=1.0)
            plt.pcolormesh(x, y, speeds, cmap='gray', edgecolor='none')
            plt.contour(x, y, times, levels=num.linspace(0., 1200, 20))
            plt.gca().axis('off')
            plt.show()

def init(store_dir):
    qssp = QSSPConfig()
    qssp.time_region = (
            gf.meta.Timing('begin-50'),
            gf.meta.Timing('end+100'))

    qssp.cut = (
            gf.meta.Timing('begin-50'),
            gf.meta.Timing('end+100'))

    store_id = os.path.basename(os.path.realpath(store_dir))
    
    config = gf.meta.ConfigTypeA(
            id = store_id,
            ncomponents = 10,
            sample_rate = 0.2,
            receiver_depth = 0*km,
            source_depth_min = 10*km,
            source_depth_max = 20*km,
            source_depth_delta = 10*km,
            distance_min = 100*km,
            distance_max = 1000*km,
            distance_delta = 10*km,
            earthmodel_1d = cake.load_model(),
            modelling_code_id = 'qssp',
            tabulated_phases = [
                    gf.meta.TPDef(
                        id = 'begin',
                        definition = 'p,P,p\\,P\\,Pv_(cmb)p'),
                    gf.meta.TPDef(
                        id = 'end',
                        definition = '2.5'),
                    gf.meta.TPDef(
                        id = 'P',
                        definition = '!P'),
                    gf.meta.TPDef(
                        id = 'S',
                        definition = '!S'),
                    gf.meta.TPDef(
                        id = 'p',
                        definition = '!p'),
                    gf.meta.TPDef(
                        id = 's',
                        definition = '!s')
            ])

    config.validate()
    return gf.store.Store.create_editables(store_dir, config=config, extra={'qssp': qssp})

    def call(self):
        '''Main work routine of the snuffling.'''
        
        self.cleanup()

        self.wanted = []
        for iphase, phase in enumerate(self._phases):
            if getattr(self, 'wantphase_%s' % iphase):
                self.wanted.append(cake.PhaseDef(phase))
        
        if not self.wanted:
            return

        viewer = self.get_viewer()
        pile = self.get_pile
        self.event = viewer.get_active_event()
        if self.event is None:
            self.fail('No active event is marked.')
        
        self.stations = [ s for s in viewer.stations.values() if s.station in pile.stations ]
        
        if not self.stations:
            self.fail('No station information available.')

        self.model = cake.load_model(self.earth_model)

        self.phases_string = []
        for s in self.stations:
            rays = self.model.arrivals([s.dist_deg], phases=self.wanted, zstart=self.event.depth)
            for ray in rays:
                incidence_angle = ray.incidence_angle()
                takeoff_angle = ray.takeoff_angle()
                m = PhaseMarker(nslc_ids = [(s.network,s.station,'*', '*')],
                                tmin = self.event.time+ray.t,
                                tmax=self.event.time+ray.t,
                                kind=1,
                                event=self.event,
                                incidence_angle=incidence_angle,
                                takeoff_angle=takeoff_angle,
                                phasename=ray.given_phase().definition())
                self.add_marker(m)           
                self.phases_string.append(ray.__str__()+'\n')

        self.called = True

    def setup(self):
        '''Customization of the snuffling.'''
        
        self.set_name('Cake')
        self._phases = ('p P pP pP\ pPv3pP\ pPv3pPv3pP\ P(moho)p S s'.split())
        for iphase, phase in enumerate(self._phases):
            self.add_parameter(Switch(phase, 'wantphase_%s' % iphase, iphase==0)) 

        self.choice = Choice('Model','earth_model',cake.builtin_models()[0],(cake.builtin_models()))
        self.add_parameter(self.choice)
        self.add_trigger('Plot Model', self.plot_model)
        self.add_trigger('Plot Rays', self.plot_rays)
        self.add_trigger('Add Model', self.add_model_to_choice)
        self.add_trigger('Add Phase', self.add_phase_definition)
        self.add_trigger('Print Arrivals', self.print_arrivals)
        self.set_live_update(False) 

        self.model = cake.load_model(self.earth_model)
        self.called = False
        self.stations = None

def GetPSArrivalRayTracing(sta_coords = np.array([0,0,0.0]), eq_coords =np.array([0,0,3900]),model_name = 'VpVs.nd'):
    # play witht the Pyrocko modules
    from pyrocko import cake
    import matplotlib
    matplotlib.style.use('ggplot')
    from LocationsOnGrid import LocationsOnGridSmall
    eq_depth = eq_coords[2]
    so_offset = np.linalg.norm(sta_coords[:2] - eq_coords[:2])
    model =cake.load_model('VpVs.nd')
    _,_,_,stCoords = LocationsOnGridSmall(receiver_name='receiver.dat',NX=1,NY = 1,NZ =1) # Get the receiver locations
    
    Distance = so_offset*cake.m2d
    p_transmission_paths = model.arrivals(distances = [Distance],phases = [cake.PhaseDef('p')],zstart = eq_depth)
    s_transmission_paths = model.arrivals(distances = [Distance],phases = [cake.PhaseDef('s')],zstart = eq_depth)
    for rayP,rayS in zip(p_transmission_paths,s_transmission_paths):
        p_arrival  = rayP.t
        print p_arrival
        s_arrival  = rayS.t
        print s_arrival
    return p_arrival,s_arrival,so_offset,model

    def dummy_store(self):
        if self._dummy_store is None:

            conf = gf.ConfigTypeA(
                id='empty_regional',
                source_depth_min=0.,
                source_depth_max=20*km,
                source_depth_delta=10*km,
                distance_min=1000*km,
                distance_max=2000*km,
                distance_delta=10*km,
                sample_rate=2.0,
                ncomponents=10,
                earthmodel_1d=cake.load_model(crust2_profile=(50., 10.)))

            store_dir = mkdtemp(prefix='gfstore')
            self.tempdirs.append(store_dir)

            gf.Store.create(store_dir, config=conf)
            self._dummy_store = gf.Store(store_dir)

        return self._dummy_store

Log.loc[Log.shape[0]-1,'VpMean']=Log.loc[Log.shape[0]-2,'VpMean']


layer_strBot = "%.3f %.3f %.3f %.3f\n" % (
                                           z_layers[i+1]*1.05,Log.loc[mask,'Vp'].mean(),
                                           Log.loc[mask,'Vs'].mean(),
                                           Log.loc[mask,'Rho'].mean()
                                           )

NDFile = NDFile + layer_strBot
#################################################################################################################
# Load the model
with open(('RepsolHighRes.nd'),'w') as f:
    f.write(NDFile)

model =cake.load_model(('VpVs.nd'))



#model =cake.load_model('VpVs.nd')
MC = 100000
nLayers = 6
pertVp = np.zeros((MC,nLayers))
pertVs = np.zeros((MC,nLayers))
tops= np.zeros(nLayers)
bots = np.zeros(nLayers)
k=0.00001
i=0
MCDf=pd.DataFrame()
for l in model.layers():
    

 def test_dump(self):
     f = BytesIO()
     mod = cake.load_model()
     mod.profile('vp').dump(f)
     f.close()

from pyrocko import cake
import collections
import numpy as num
import matplotlib.pyplot as plt

"""
soll klarstellen, warum ich 4 sekunden als static length benutze.
"""

km = 1000

ttt = collections.defaultdict(dict)

m = cake.load_model("earthmodel_castor.txt")
phases = ["Pv12.5p", "P", "p", "pP"]
phases = [cake.PhaseDef(p) for p in phases]

distances = num.linspace(10, 200, 40) * km * cake.m2d
print distances
toplot = []
for d in distances:
    allts = []
    print d
    for p in phases:
        for arrival in m.arrivals([d], phases=p, zstart=2000):
            print arrival
            allts.append(arrival.t)

    toplot.append([d, max(allts) - min(allts)])

import os.path as op
import logging

import numpy as num

from pyrocko import spit, cake

logging.basicConfig(level=logging.INFO)

km = 1000.

model = cake.load_model('ak135-f-average.f')

phases = [cake.PhaseDef('P'), cake.PhaseDef('p'), cake.PhaseDef('Pv_(cmb)p')]


def f(x):
    source_depth, distance = x
    rays = model.arrivals(phases=phases, distances=[distance],
                          zstart=source_depth)
    if rays:
        return rays[0].takeoff_angle()

    return None


fn = 'takeoffangles.sptree'
if not op.exists(fn):
    tree = spit.SPTree(f, 0.01, [[0., 750*km], [0., 120.]], [1.*km, 0.1])
    tree.dump(fn)