<jupyter_start><jupyter_text>Workflow -- Forward Modeling1. [For a given mantle composition, use perplex to generate material properties for T-P field (convert P to Z)](perplex/Plot_PT.ipynb)2. [Assume a value for heatflow, then use xenoliths to invert for a range of geotherms](geotherm/Invert_for_Geotherm.ipynb)3. Make a velocity model  * Crustal constraints pulled from Crust1.0  * Mantle constraints from perplex calculation4. Run Mineos to compute phase velocities and compare with observations<jupyter_code>#Step 1 - Import mappings for P-T calculation
import warnings
warnings.filterwarnings('ignore')
from matplotlib import pylab as plt
%matplotlib inline
plt.style.use('ggplot')

rocktype = 'dunite'
location = 'ColoradoPlateau'

import pickle
fin = open('perplex/%s_stx11/funs.pickle' % rocktype,'rb')
[frho,fvp,fvs] = pickle.load(fin)
fin.close()

from minos.minos import Minos
#Step 2 - Import range of geotherms
from geotherm import *
fig = plt.figure(1,(6,6))

gts={}

for tag in ['min','max']:

    fin = open('geotherm/geotherm_%s_%s.pickle' % (location, tag) ,'rb')
    gt = pickle.load(fin)
    fin.close()
    
    if tag == 'min':
        endmember = 'Warm'
    else:
        endmember = 'Cool'

    plt.plot(gt.temps,gt.depths,label='%s: %.2f $\mu$W/m$^3$' % (endmember, gt.rhoH*10**6), lw=3)
    
    gts[tag] = gt
    
plt.scatter(gt.xeno_temps,gt.xeno_depths, label = 'xenoliths',zorder=1000, color='black')
plt.title('Heat Flow: %2d mW/m$^2$' % (gt.q0*1e3) )
plt.ylim(max(gt.depths),0)
_,_=plt.ylabel('Depth (km)'), plt.xlabel('Temperature ($^\circ$C)')
_ =plt.legend(loc=3)

plt.savefig('img/%s_geotherm.eps' % location)
#Step 3 - Make velocity model and write to card file
from matplotlib import pylab as plt
%matplotlib inline
plt.style.use('ggplot')

if "Colorado" in location:
    from shen.shen import CrustalModel
    c  = CrustalModel('coloplat')
    zmoho = c.zmoho
    atten_file = 'seismic_JF10_gs10mm_Temp_AD0p40_TP1350_ZG155'
    #lat, lon = 36.,-110.
    lat, lon = 0., 0.
    
elif "Slave" in location:
    from crust1.crust1 import Crust1
    c = Crust1().load('crust1/crustal_model_SlaveCraton.pickle').print()
    zmoho = abs(c.model.Bottom[7])
    atten_file = 'seismic_JF10_gs10mm_Temp_AD0p40_TP1350_ZG205'
    lat, lon = 65.,-110.
    
elif "Africa" in location:
    from crust1.crust1 import Crust1
    c = Crust1().load('crust1/crustal_model_SouthAfrica.pickle').print()
    zmoho = abs(c.model.Bottom[7])
    atten_file = 'seismic_JF10_gs10mm_Temp_AD0p40_TP1350_ZG205'
    lat, lon = -27., 27.
    
elif "Siberia" in location:
    from crust1.crust1 import Crust1
    c = Crust1().load('crust1/crustal_model_Siberia.pickle').print()
    zmoho = abs(c.model.Bottom[7])
    atten_file = 'seismic_JF10_gs10mm_Temp_AD0p40_TP1350_ZG205'
    lat, lon = 62., 96.

def extract_profile(fun, gt):
    from scipy.interpolate import interp1d
    
    vs, zs = [],[]
    for ii, eachz in enumerate(gt.depths):
        
        val = fun(gt.temps[ii], eachz)
            
        vs.append(val[0])
        zs.append(eachz)
    
    return interp1d(zs,vs)
    
    
from minos.minos import CardFile
from attenuation.attenuation import AttenuationProfile

fqmu = AttenuationProfile(atten_file).setup_interpolation_functions()

ls = ['red','blue']

for force_isotropy in [True, False]:

    cf = CardFile().read('minos/STW105.txt')
    cf1 = CardFile().read('minos/STW105.txt')
    cf2 = cf1.regional_profile_from_kustow(lat=lat, lon=lon)

    for kk,tag in enumerate(gts):
        fvp_gt  = extract_profile(fvp, gts[tag])
        fvs_gt  = extract_profile(fvs, gts[tag])
        frho_gt = extract_profile(frho, gts[tag])

        cf0=cf2.set_mantle_region(fvp_gt, fvs_gt, frho_gt, zmoho=zmoho, force_isotropy=force_isotropy)
        cf0=cf0.set_mantle_attenuation(fqmu, zmoho=zmoho)

        #cf0=cf0.set_crustal_region(c.model)
        if "Colorado" in location:
            cf0=cf0.set_crustal_region_shen(c.model, zmoho)
        else:
            cf0=cf0.set_crustal_region(c.model)

        cf0.write(filename='minos/%s.txt' % tag)

        plt.figure(2,figsize=(15,5))
        _ = plt.plot(cf0.df.Vpv, 6371-cf0.df.Radius/1000, lw = 2, label = 'Vpv (m/s)', c=ls[kk])
        _ = plt.plot(cf0.df.Vph, 6371-cf0.df.Radius/1000, lw = 2, label = 'Vph', c=ls[kk], ls ='--')
        #_ = plt.plot(cf.df.Vpv, 6371-cf.df.Radius/1000, '--k', lw = 0.5)
        _ = plt.ylim(500,0)

        _ = plt.plot(cf0.df.Vsv, 6371-cf0.df.Radius/1000,label='Vsv', lw=2, c=ls[kk])
        _ = plt.plot(cf0.df.Vsh, 6371-cf0.df.Radius/1000,label='Vsh', lw=2, c=ls[kk], ls ='--')
        #_ = plt.plot(cf.df.Vsv, 6371-cf.df.Radius/1000, '--k', lw = 0.5)

        _ = plt.plot(cf0.df.Density, 6371-cf0.df.Radius/1000,label='Density (kg/m$^3$)', lw=2, c=ls[kk])
        #_ = plt.plot(cf.df.Density, 6371-cf.df.Radius/1000, '--k', lw = 0.5)

    _ = plt.ylim(500,-10)
    _ = plt.xlim(2000, 10000)
    _ = plt.title("Card File for MINEOS at %s" % location)
    _ = plt.ylabel('Depth (km)')
    _ = plt.xlabel('Density, $V_s$, $V_p$')
    _ = plt.legend(ncol=2,fontsize=8,title='Warm (red) -- Cool (blue)', loc='best', handlelength=3)

    plt.savefig('img/%s_%s_cardfile.eps' % (location, rocktype) )
    plt.show()


    #Step 4 - Compute phase velocities and save to disk

    for tag in gts:
        for ref_period in [-1.0, 1.0]:
            cf = CardFile().read('minos/%s.txt' % tag).write(filename='minos/%s.txt' % tag, ref_period=ref_period)
            minos = Minos(tag, cf).run().read_output()

            if ref_period > 0:
                tag2 = 'anelas'
            else:
                tag2 = ''

            if force_isotropy:
                tag3 = ''
            else:
                tag3 = 'aniso'

            filename='output_minos/minos.%s.%s.%s.%s.%s.pickle' % (tag, tag2, tag3, rocktype, location)
            minos.save(filename)<jupyter_output>Pulling regional profile...
done.<jupyter_text>**Note**: At depths > 350 km, the Card File pulls direct values from STW105, resulting in artificial discontinuities around that depth.<jupyter_code>#Import observations and plot against predictions
from observations.observations import PhaseVelocities, Region

fig = plt.figure(1,(16,16))

tag3 = ''

for iplt, location in enumerate(["ColoradoPlateau","SlaveCraton","SouthAfrica","Siberia"]):
    if "Colorado" in location:
        sources = ["Ekstrom","Ma","Cloud"]
        rgs=[Region(35,39,-111,-109,source) for source in sources]
        
    elif "Africa" in location:
        sources = ["Ekstrom","Ma"]
        rgs=[Region(-28,-26,26,28,source) for source in sources]
        
    elif "Siberia" in location:
        sources = ["Ekstrom","Ma"]
        rgs=[Region(61,63,95,97,source) for source in sources]

    elif "Slave" in location:
        sources = ["Ekstrom","Ma"]
        rgs=[Region(63,66,-112,-108,source) for source in sources]
        
    plt.subplot(2,2,iplt+1)

    colors = ['magenta','darkgreen']
    kk=0

    for rt in ['peridotite','dunite']:
        for ref_period in [-1.0, 1.0]: 
            if ref_period > 0:
                tag2 = ''
            else:
                tag2 = 'no anelas'
                continue

            filename1='output_minos/minos.%s.%s.%s.%s.%s.pickle' % ('max',tag2,tag3,rt,location)
            filename2='output_minos/minos.%s.%s.%s.%s.%s.pickle' % ('min',tag2,tag3,rt,location)

            minos1 = Minos('max').load(filename1)
            minos2 = Minos('min').load(filename2)


            _ = plt.fill_between(minos1.period, minos1.c, minos2.c, label='%s %s' % (rt,tag2), lw=1,
                                 facecolor=colors[kk], edgecolor='darkgray', zorder=1000, alpha=0.25)

            kk += 1


    colors = ['red','black','cyan']
    for isource, rg in enumerate(rgs):
        xs,ys,yerrs = [],[],[]
        for key in rg.mean:
            xs.append(key)
            ys.append(rg.mean[key]/1000.)
            yerrs.append(rg.std[key]/1000. *2.0)

        _ = plt.errorbar(xs, ys, yerr=yerrs,
                         elinewidth=3, label='Obs. from %s' % sources[isource], color=colors[isource], capsize=5, linestyle=' ',capthick=3,
                        zorder = 1)

    plt.xlim(0,200)
    plt.ylim(3.45,4.55)
    plt.ylabel('c (km/s)')
    plt.xlabel('Period (s)')

    _,_ = plt.legend(loc=4), plt.title(location)
    
plt.savefig('img/observations.png')
    
fig = plt.figure(1,(16,8))


for iplt, location in enumerate(["SlaveCraton"]):
    if "Slave" in location:
        sources = ["Ekstrom","Ma"]
        rgs=[Region(63,66,-112,-108,source) for source in sources]
        
    for iplt, tag3 in enumerate(['','aniso']):
        plt.subplot(1,2,iplt+1)

        colors = ['magenta','darkgreen']
        kk=0

        for rt in ['peridotite','dunite']:
            for ref_period in [-1.0, 1.0]: 
                if ref_period > 0:
                    tag2 = ''
                else:
                    tag2 = 'no anelas'
                    continue

                filename1='output_minos/minos.%s.%s.%s.%s.%s.pickle' % ('max',tag2,tag3,rt,location)
                filename2='output_minos/minos.%s.%s.%s.%s.%s.pickle' % ('min',tag2,tag3,rt,location)

                minos1 = Minos('max').load(filename1)
                minos2 = Minos('min').load(filename2)


                _ = plt.fill_between(minos1.period, minos1.c, minos2.c, label='%s %s' % (rt,tag2), lw=1,
                                     facecolor=colors[kk], edgecolor='darkgray', zorder=1000, alpha=0.25)

                kk += 1


        colors = ['red','black','cyan']
        for isource, rg in enumerate(rgs):
            xs,ys,yerrs = [],[],[]
            for key in rg.mean:
                xs.append(key)
                ys.append(rg.mean[key]/1000.)
                yerrs.append(rg.std[key]/1000. *2.0)

            _ = plt.errorbar(xs, ys, yerr=yerrs,
                             elinewidth=3, label='Obs. from %s' % sources[isource], color=colors[isource], capsize=5, linestyle=' ',capthick=3,
                            zorder = 1)

        plt.xlim(0,200)
        plt.ylim(3.45,4.55)
        plt.ylabel('c (km/s)')
        plt.xlabel('Period (s)')

        if 'aniso' in tag3:
            title = '%s with radial aniso. from S362ANI' % location
        else:
            title = '%s without radial aniso.' % location
        
        _,_ = plt.legend(loc=4), plt.title(title)
    
plt.savefig('img/radial_anisotropy_slave.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15514108082257294, 'fraction_numerical': 0.0436154949784792, 'mean_word_length': 2.5432056929854285, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '632810', 'license': 'MIT', 'n_tokens_mistral': 4279, 'n_tokens_neox': 3909, 'n_words': 728, 'path': 'Forward_Modeling_Workflow.ipynb', 'repo_name': 'nmancinelli/ss-precursors', 'size': '10455'}