<jupyter_start><jupyter_text>Sensibility Matrix - G Este [IPython Notebook](http://ipython.org/videos.htmlthe-ipython-notebook) utiliza a biblioteca de código aberto [Fatiando a Terra](http://fatiando.org/)<jupyter_code>#%matplotlib inline
import numpy as np
from scipy.misc import derivative
import scipy as spy
from scipy import interpolate
import matplotlib
matplotlib.use('TkAgg', force=True)
import matplotlib.pyplot as plt
import math
import cPickle as pickle
import datetime
#from IPython.display import Image as img
#from IPython.display import Markdown as md
#from IPython.display import display as dp
import string as st

from scipy.misc import imread

from __future__ import division

from fatiando import gravmag, mesher, utils, gridder
from fatiando.mesher import Prism, Polygon
from fatiando.gravmag import prism
from fatiando.utils import ang2vec, si2nt, contaminate
from fatiando.gridder import regular
from fatiando.vis import mpl

from numpy.testing import assert_almost_equal
from numpy.testing import assert_array_almost_equal
from pytest import raises

plt.rc('font', size=16)
import functions as fc
# Model`s limits
ymin = 0.0
ymax = 100000.0
zmin = -1000.0
zmax = 100000.0
xmin = -100000.0
xmax =  100000.0

area = [ymin, ymax, zmax, zmin]
ny = 10 # number of observation datas and number of prisms along the profile
# coordinates defining the horizontal boundaries of the
# adjacent columns along the profile
y = np.linspace(ymin, ymax, ny)
# coordinates of the center of the columns forming the
# interpretation model
n = ny - 1
dy = (ymax - ymin)/n
ycmin = ymin + 0.5*dy
ycmax = ymax - 0.5*dy
yc = np.reshape(np.linspace(ycmin, ycmax, n),(n,1))
x = np.zeros_like(yc)
z = np.zeros_like(yc)-150.0
## Edge extension (observation coordinates)
sigma = 0.5
edge = sigma*dy*n
# Model densities
# Indices and polygons relationship:
# cc = continental crust layer
# oc = ocean crust layer
# w = water layer
# s = sediment layer
# m = mantle layer
dw = np.array([1030.0])
ds0 = np.array([2550.0])
dcc = np.array([2670.0])
doc = np.array([2840.0])
dm = np.array([3200.0])
#dc = dcc
# coordinate defining the horizontal boundaries of the continent-ocean boundary
COT = 60000.0

# list defining crust density variance
dc = np.zeros_like(yc)
aux = yc <= COT
for i in range(len(yc[aux])):
    dc[i] = dcc
for i in range(len(yc[aux]),n):
    dc[i] = doc
# defining sediments layers density matrix
ds = np.reshape(np.repeat(ds0,n),(n,1))
    
# S0 => isostatic compensation surface (Airy's model)
# SR = S0+dS0 => reference Moho (Forward modeling)
S0 = np.array([40000.0])
dS0 = np.array([5000.0])

pjmin = np.array([-50.0]) 
pjmax = np.array([100000.0])
dp = 1000.
A = 500.
B = 3500.
c = 0.00009
D = 20000.
bathymetry = fc.bathymetry_function(A,B,c,D,yc)
tw = bathymetry - 0.0<jupyter_output><empty_output><jupyter_text>mpl.close('all')axes = mpl.figure().add_subplot(1,1,1)mpl.ylim(zmax, zmin)mpl.xlim(ymin, ymax)mpl.xlabel('y (m)')mpl.ylabel('z (m)')mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=1)mpl.plot(yc, bathymetry, '-b', linewidth=1)mpl.plot(yc, sediments_limit, '--r', linewidth=1)draw = mpl.draw_polygon(area, axes, color='r') draw<jupyter_code>sediments = np.array([[  6250.        ,  19331.16883117],
       [ 21774.19354839,  17418.29004329],
       [ 38911.29032258,  16325.21645022],
       [ 63306.4516129 ,  14685.60606061],
       [ 92741.93548387,  13592.53246753]])
# change the coordinates of the extremum points in order to
# avoid problems for constructing the interpolator
sediments[0,0] = ymin
sediments[-1,0] = ymax

sediments_interpolated = fc.surface_interpolate_function(sediments,yc)

ts = sediments_interpolated - tw
S = fc.moho_function(S0,dw,ds,dcc,dm,dc,tw,ts)

tm = S0 - S

toc = S - tw - ts #thickness of the oceanic crust along the profile
print ts.shape, tm.shape, S0.shape, dS0.shape
psyn = []
psyn = np.vstack((ts, tm, dS0))
# prisms calculation by <NAME>
prism_w_syn = fc.prism_w_function(xmax,xmin,dy,edge,dw,dcc,tw,yc)
prism_s_syn = fc.prism_s_function(xmax,xmin,dy,edge,ds,dcc,tw,psyn,yc)
prism_c_syn = fc.prism_c_function(xmax,xmin,dy,edge,S0,dcc,dc,tw,psyn,yc)
prism_m_syn = fc.prism_m_function(xmax,xmin,dy,edge,S0,dcc,dm,psyn,yc)
# z component of gravity calculation by <NAME>
gzw_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_w_syn)
gzs_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_s_syn[0])
gzc_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_c_syn)
gzm_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_m_syn)

#Observed data calculation:
gsyn = gzw_syn + gzs_syn + gzc_syn + gzm_syn
#gsyn = fc.g_function(x,yc,z,gzw_syn,prism_s_syn,prism_c_syn,prism_m_syn)<jupyter_output><empty_output><jupyter_text>Calcula matriz de sensibilidade com plot da variacao do parametro<jupyter_code>G1 = fc.G_matrix_function(xmax,xmin,dy,edge,1000.,1000.,S0,dw,ds,dm,dcc,dc,tw,psyn,yc)
%matplotlib inline
G2 = fc.G_matrix_function_all(xmax,xmin,dy,edge,1000.,1000.,COT,S0,dw,ds,dm,dcc,dc,tw,psyn,yc,area,11)<jupyter_output><empty_output><jupyter_text>Compara dois metodos de calculo da matriz sensibilidade<jupyter_code>assert_array_almost_equal(G1, G2, decimal=3)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup