<jupyter_start><jupyter_text>2D Isostatic gravity inversion - True Model 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>True model of magma poor margin.<jupyter_output><empty_output><jupyter_text>%matplotlib inlineimport numpy as npfrom scipy.misc import derivativeimport scipy as spyfrom scipy import interpolateimport matplotlibmatplotlib.use('TkAgg', force=True)import matplotlib.pyplot as pltimport mathimport cPickle as pickleimport datetimeimport string as stfrom scipy.misc import imreadfrom __future__ import divisionfrom fatiando import gravmag, mesher, utils, gridderfrom fatiando.mesher import Prism, Polygonfrom fatiando.gravmag import prismfrom fatiando.utils import ang2vec, si2nt, contaminatefrom fatiando.gridder import regular, profilefrom fatiando.vis import mplfrom numpy.testing import assert_almost_equalfrom numpy.testing import assert_array_almost_equalfrom pytest import raisesplt.rc('font', size=16) import functions as fc<jupyter_code>## Observation coordinates.<jupyter_output><empty_output><jupyter_text>Model`s limitsymin = 0.0ymax = 383000.0zmin = -1000.0zmax = 45000.0xmin = -100000.0xmax = 100000.0area = [ymin, ymax, zmax, zmin] ny = 150  number of observation datas and number of prisms along the profile coordinates defining the horizontal boundaries of the adjacent columns along the profiley = np.linspace(ymin, ymax, ny)  coordinates of the center of the columns forming the interpretation modeln = ny - 1dy = (ymax - ymin)/nycmin = ymin + 0.5*dyycmax = ymax - 0.5*dyyc = 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 = 2.0edge = sigma*dy*n<jupyter_code>## Model parameters<jupyter_output><empty_output><jupyter_text>Model densities Indices and polygons relationship: cc = continental crust layer oc = ocean crust layer w = water layer s = sediment layer m = mantle layerdw = np.array([1030.0])ds0 = np.array([2350.0])ds1 = np.array([2600.0])dcc = np.array([2870.0])doc = np.array([2885.0])dm = np.array([3300.0])dc = dcc coordinate defining the horizontal boundaries of the continent-ocean boundaryCOT = 350000.0 list defining crust density variancedc = np.zeros_like(yc)aux = yc <= COTfor i in range(len(yc[aux])):    dc[i] = dccfor i in range(len(yc[aux]),n):    dc[i] = doc     defining sediments layers density matrixds = np.vstack((np.reshape(np.repeat(ds0,n),(1,n)),np.reshape(np.repeat(ds1,n),(1,n)))) S0 => isostatic compensation surface (Airy's model) SR = S0+dS0 => reference Moho (Forward modeling)S0 = np.array([44000.0]) originaldS0 = np.array([2200.0]) original<jupyter_code>## For display<jupyter_output><empty_output><jupyter_text>img = imread('../data/fig/zalan[2015]_fig2_cropped.png') shape=(499,1320)extent = [ymin, ymax, zmax, 0.0]def calculate_aspect(shape,extent):    dx = (abs(extent[1] - extent[0]))/float(shape[1])    dy = (abs(extent[3] - extent[2]))/float(shape[0])    return dx/dy<jupyter_code>## Water bottom<jupyter_output><empty_output><jupyter_text>mpl.close('all')axes = mpl.subplot(1,1,1)mpl.ylim(zmax, zmin)mpl.xlim(ymin, ymax)mpl.xticks(fontsize=12)mpl.yticks(fontsize=12)mpl.xlabel('y (km)')mpl.ylabel('z (km)')mpl.imshow(img, extent=[ymin, ymax, zmax, 0.0], alpha=0.3)mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=0.3)mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='-k', linewidth=0)mpl.m2km()bath_picks = mpl.draw_polygon(area, axes, marker='.', color='r') bath_picks bath_picks = np.array([[  1.15826613e+03,   3.08550569e+02],       [  1.57524194e+05,   9.08540872e+02],       [  1.62929435e+05,   1.02853893e+03],       [  1.70265121e+05,   1.38853311e+03],       [  1.77214718e+05,   1.86852536e+03],       [  1.88411290e+05,   2.10852148e+03],       [  2.06943548e+05,   2.58851372e+03],       [  2.54432460e+05,   2.82850984e+03],       [  2.75667339e+05,   2.94850790e+03],       [  2.94199597e+05,   3.18850402e+03],       [  3.10801411e+05,   3.30850208e+03],       [  3.32422379e+05,   3.78849433e+03],       [  3.81455645e+05,   3.90849239e+03]])  change the coordinates of the extremum points in order to avoid problems for constructing the interpolatorbath_picks[0,0] = yminbath_picks[-1,0] = ymax bathymetry = fc.surface_interpolate_function(bath_picks,yc)print bathymetry.shape np.savetxt('../data/volcanic-margin-bathymetry.txt', bathymetry, fmt='%.18f')  layer water thicknesstw = np.reshape(bathymetry,(n,1))<jupyter_code>bathymetry = np.reshape(np.loadtxt('../data/etopo1-pelotas.txt'),(n,1))
tw = 0.0 - bathymetry<jupyter_output><empty_output><jupyter_text>Sediments surfaces<jupyter_code>mpl.close('all')

axes = mpl.subplot(1,1,1)
mpl.ylim(zmax, zmin)
mpl.xlim(ymin, ymax)
mpl.xticks(fontsize=12)
mpl.yticks(fontsize=12)
mpl.xlabel('y (km)')
mpl.ylabel('z (km)')
#mpl.imshow(img, extent=[ymin, ymax, zmax, 0.0], alpha=0.3)
mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=0.3)
mpl.plot(yc, tw, '-b', linewidth=1)
mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='og', linewidth=0)
mpl.m2km()

toi_picks = mpl.draw_polygon(area, axes, marker='.', color='r')
toi_picks<jupyter_output><empty_output><jupyter_text>toi_picks = np.array([[    772.17741935,    1148.53699334],       [  16601.81451613,    1148.53699334],       [  25095.76612903,    1148.53699334],       [  29728.83064516,    1388.53311462],       [  42855.84677419,    2228.51953908],       [  63318.5483871 ,    3668.49626673],       [  71426.41129032,    4388.48463055],       [  77603.83064516,    4748.47881246],       [  84167.33870968,    4628.48075183],       [  88028.22580645,    4628.48075183],       [  91116.93548387,    4748.47881246],       [  93433.46774194,    5108.47299438],       [  98066.53225806,    5948.45941884],       [  98452.62096774,    4988.47493374],       [ 100769.15322581,    4988.47493374],       [ 103857.86290323,    4868.4768731 ],       [ 110035.28225806,    4988.47493374],       [ 120459.67741935,    4988.47493374],       [ 126251.00806452,    5828.4613582 ],       [ 137447.58064516,    7388.43614649],       [ 142852.82258065,    7868.42838904],       [ 145169.35483871,    7868.42838904],       [ 148258.06451613,    7268.43808585],       [ 151346.77419355,    7748.4303284 ],       [ 165632.05645161,    7508.43420712],       [ 173353.83064516,    7388.43614649],       [ 176056.4516129 ,    7028.44196457],       [ 179145.16129032,    6788.4458433 ],       [ 182619.95967742,    6908.44390394],       [ 186480.84677419,    7148.44002521],       [ 189569.55645161,    7148.44002521],       [ 193430.44354839,    6788.4458433 ],       [ 213507.05645161,    6788.4458433 ],       [ 239761.08870968,    6908.44390394],       [ 245938.50806452,    7028.44196457],       [ 253660.28225806,    7028.44196457],       [ 263312.5       ,    6908.44390394],       [ 297674.39516129,    7028.44196457],       [ 306940.52419355,    7028.44196457],       [ 318909.27419355,    7148.44002521],       [ 324700.60483871,    7148.44002521],       [ 327789.31451613,    6908.44390394],       [ 332036.29032258,    7268.43808585],       [ 335897.17741935,    7028.44196457],       [ 337055.44354839,    6788.4458433 ],       [ 341302.41935484,    6908.44390394],       [ 345163.30645161,    6548.44972202],       [ 347479.83870968,    6548.44972202],       [ 354815.52419355,    7028.44196457],       [ 359062.5       ,    6908.44390394],       [ 363309.47580645,    6668.44778266],       [ 365239.91935484,    6548.44972202],       [ 367942.54032258,    6548.44972202],       [ 372189.51612903,    6668.44778266],       [ 376050.40322581,    6668.44778266],       [ 381455.64516129,    6548.44972202]])  change the coordinates of the extremum points in order to avoid problems for constructing the interpolatortoi_picks[0,0] = ymintoi_picks[-1,0] = ymax toi = fc.surface_interpolate_function(toi_picks,yc)for i in range(len(toi)):    if toi[i] < tw[i]:        toi[i] = tw[i]<jupyter_code>np.savetxt('../data/volcanic-margin-true-toi-surface.txt', toi, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>layer sediments thicknessts0 = toi - tw<jupyter_code>mpl.close('all')

axes = mpl.subplot(1,1,1)
mpl.ylim(zmax, zmin)
mpl.xlim(ymin, ymax)
mpl.xticks(fontsize=12)
mpl.yticks(fontsize=12)
mpl.xlabel('y (km)')
mpl.ylabel('z (km)')
#mpl.imshow(img, extent=[ymin, ymax, zmax, 0.0], alpha=0.3)
mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=0.3)
mpl.plot(yc, toi, '-b', linewidth=1)
mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='og', linewidth=0)
mpl.m2km()

basement_picks = mpl.draw_polygon(area, axes, marker='.', color='r')
basement_picks<jupyter_output><empty_output><jupyter_text>basement_picks = np.array([[   1158.26612903,    3000.00],       [   8880.04032258,    2588.51372099],       [  14671.37096774,    1748.52729653],       [  18918.34677419,    1988.5234178 ],       [  26640.12096774,    3548.49820609],       [  43241.93548387,   13628.33529963],       [  55982.86290323,   20348.22669533],       [  63318.5483871 ,   23228.18015062],       [  70654.23387097,   25388.1452421 ],       [  76059.47580645,   24428.160757  ],       [  81850.80645161,   25988.13554528],       [  86483.87096774,   26348.1297272 ],       [  89958.66935484,   25628.14136337],       [  99996.97580645,   18908.24996768],       [ 103085.68548387,   19028.24802831],       [ 106560.48387097,   20108.23057405],       [ 110421.37096774,   20108.23057405],       [ 113896.16935484,   19028.24802831],       [ 122004.03225806,   16868.28293684],       [ 129725.80645161,   12548.35275389],       [ 132814.51612903,   12788.34887517],       [ 134744.95967742,   13028.34499644],       [ 137447.58064516,   12308.35663262],       [ 141694.55645161,   12548.35275389],       [ 146713.70967742,   13388.33917836],       [ 149030.24193548,   13508.33723899],       [ 152891.12903226,   12788.34887517],       [ 155979.83870968,   12668.35081453],       [ 157910.28225806,   12668.35081453],       [ 162929.43548387,   13508.33723899],       [ 167562.5       ,   13748.33336027],       [ 169879.03225806,   13148.34305708],       [ 176828.62903226,   13148.34305708],       [ 181075.60483871,   13148.34305708],       [ 188025.2016129 ,   14108.32754218],       [ 191113.91129032,   14108.32754218],       [ 195360.88709677,   14348.32366346],       [ 200380.04032258,   14828.31590601],       [ 203854.83870968,   15548.30426983],       [ 205399.19354839,   15068.31202728],       [ 207715.72580645,   16268.29263365],       [ 210418.34677419,   17228.27711875],       [ 213120.96774194,   17348.27517939],       [ 217367.94354839,   17228.27711875],       [ 221614.91935484,   15548.30426983],       [ 223545.36290323,   15548.30426983],       [ 228950.60483871,   16148.29457302],       [ 232039.31451613,   16388.29069429],       [ 234355.84677419,   16268.29263365],       [ 237444.55645161,   17108.27905812],       [ 239375.        ,   17108.27905812],       [ 245166.33064516,   16028.29651238],       [ 248641.12903226,   15068.31202728],       [ 250957.66129032,   13868.33142091],       [ 252888.10483871,   13988.32948154],       [ 255976.81451613,   14948.31396664],       [ 262154.23387097,   13508.33723899],       [ 265629.03225806,   15068.31202728],       [ 268331.65322581,   15668.30233047],       [ 274122.98387097,   16148.29457302],       [ 276825.60483871,   15428.30620919],       [ 280686.49193548,   13748.33336027],       [ 283389.11290323,   12188.35857198],       [ 286863.91129032,   14468.32172409],       [ 293041.33064516,   15908.29845174],       [ 300763.10483871,   16508.28875493],       [ 309643.14516129,   16268.29263365],       [ 321997.98387097,   15188.31008792],       [ 328947.58064516,   13148.34305708],       [ 332036.29032258,   10988.37796561],       [ 333966.73387097,    7628.43226776],       [ 334352.82258065,    4988.47493374],       [ 382613.91129032,    5108.47299438]])  change the coordinates of the extremum points in order to avoid problems for constructing the interpolatorbasement_picks[0,0] = yminbasement_picks[-1,0] = ymax basement = fc.surface_interpolate_function(basement_picks,yc)for i in range(len(basement)):    if basement[i] < toi[i]:        basement[i] = toi[i]<jupyter_code>np.savetxt('../data/volcanic-margin-true-basement-surface.txt', basement, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>layer igneous thicknessts1 = basement - toi  thickness sediments vectorts = np.vstack((np.reshape(ts0,(1,n)),np.reshape(ts1,(1,n))))<jupyter_code>## Moho surface<jupyter_output><empty_output><jupyter_text>mpl.close('all')axes = mpl.subplot(1,1,1)mpl.ylim(zmax+10000., zmin)mpl.xlim(ymin, ymax)mpl.xticks(fontsize=12)mpl.yticks(fontsize=12)mpl.xlabel('y (km)')mpl.ylabel('z (km)')mpl.imshow(img, extent=[ymin, ymax, zmax, 0.0], alpha=0.3)mpl.imshow(img, extent=extent, aspect=calculate_aspect(shape,extent), alpha=0.3)mpl.plot(yc, toi, '-b', linewidth=1)mpl.plot(yc, basement, '-g', linewidth=1)mpl.plot(yc, moho, '-g', linewidth=1)mpl.paths([[ymin, 0.0]], [[ymax, 0.0]], style='og', linewidth=0)mpl.m2km()area2 = [ymin, ymax, zmax+10000., zmin]moho_picks = mpl.draw_polygon(area2, axes, marker='.', color='r') moho_picks<jupyter_code>moho_picks = np.array([[   1518.33556011,   37872.3690103 ],
       [   9985.48495019,   37967.53246753],
       [  18981.83117715,   37586.8786386 ],
       [  34857.73628355,   36349.75369458],
       [  56554.80659562,   34160.99417824],
       [  75076.69588642,   32543.21540528],
       [  88306.61680842,   31686.74429019],
       [ 104182.52191482,   30830.2731751 ],
       [ 118470.83651057,   29402.82131661],
       [ 127996.37957441,   27309.2252575 ],
       [ 138580.31631201,   24073.6677116 ],
       [ 148105.85937585,   21313.92745186],
#       [ 151810.23723401,   20647.78325123],
#       [ 153927.02458153,   20362.29287953],
#       [ 159748.18978721,   20076.80250784],
#       [ 163981.76448225,   20457.45633677],
#       [ 151810.23723401,   20647.78325123],
#       [ 153927.02458153,   20362.29287953],
#       [ 159748.18978721,   20076.80250784],
#       [ 163981.76448225,   20457.45633677],
       [ 173507.30754608,   21884.90819525],
       [ 181445.26009928,   22931.70622481],
       [ 192558.39367376,   23122.03313927],
       [ 201554.73990072,   23217.19659651],
       [ 207375.9051064 ,   23026.86968204],
       [ 213726.26714896,   22360.72548142],
       [ 221135.02286527,   21028.43708016],
       [ 226426.99123407,   20267.1294223 ],
       [ 231189.76276599,   20267.1294223 ],
       [ 236481.73113479,   20647.78325123],
       [ 242832.09317735,   21409.09090909],
       [ 249182.45521991,   21409.09090909],
       [ 253945.22675183,   21789.74473802],
       [ 259237.19512063,   22265.56202418],
       [ 271408.72236887,   22265.56202418],
       [ 275642.2970639 ,   21789.74473802],
       [ 278817.47808518,   21028.43708016],
       [ 283580.2496171 ,   20933.27362293],
       [ 289401.41482278,   19791.31213614],
       [ 295751.77686534,   18459.02373489],
       [ 299456.1547235 ,   18268.69682042],
       [ 304748.1230923 ,   18839.67756382],
       [ 312686.0756455 ,   19600.98522167],
       [ 321682.42187246,   19886.47559337],
       [ 332795.55544693,   19791.31213614],
       [ 337029.13014197,   19600.98522167],
       [ 342321.09851077,   18363.86027765],
       [ 345496.27953205,   17412.22570533],
       [ 380423.27076613,   17031.5718764 ]])
# change the coordinates of the extremum points in order to
# avoid problems for constructing the interpolator
moho_picks[0,0] = ymin
moho_picks[-1,0] = ymax
moho = fc.surface_interpolate_function(moho_picks,yc)
for i in range(len(moho)):
    if moho[i] < basement[i]:
        moho[i] = basement[i]
# layer mantle thickness
tm = S0 - moho

# layer crust thickness
toc = moho - tw - ts0 - ts1<jupyter_output><empty_output><jupyter_text>np.savetxt('../data/volcanic-margin-true-moho-surface.txt', moho, fmt='%.18f')<jupyter_code>## Known depths (basement)<jupyter_output><empty_output><jupyter_text>print '[', np.reshape(yc,(n,))[3], ',' , np.reshape(basement,(n,))[3], '],'print '[', np.reshape(yc,(n,))[9], ',' , np.reshape(basement,(n,))[9], '],'print '[', np.reshape(yc,(n,))[42], ',' , np.reshape(basement,(n,))[42], '],'print '[', np.reshape(yc,(n,))[59], ',' , np.reshape(basement,(n,))[59], '],'print '[', np.reshape(yc,(n,))[67], ',' , np.reshape(basement,(n,))[67], '],'print '[', np.reshape(yc,(n,))[91], ',' , np.reshape(basement,(n,))[91], '],'print '[', np.reshape(yc,(n,))[112], ',' , np.reshape(basement,(n,))[112], '],'print '[', np.reshape(yc,(n,))[132], ',' , np.reshape(basement,(n,))[132], ']'<jupyter_code># Basement surface known values
basement_known = np.array([[ 8996.6442953 , 2571.60124265 ],
[ 24419.4630872 , 3099.87458409 ],
[ 109244.966443 , 20108.2305741 ],
[ 235197.986577 , 16497.3157779 ],
[ 289177.852349 , 15007.7089799 ],
[ 340587.248322 , 6888.23678696 ]])

#np.savetxt('../data/volcanic-margin-basement-known-depths.txt', basement_known, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>print '[', np.reshape(yc,(n,))[3], ',' , np.reshape(basement,(n,))[3], '],'print '[', np.reshape(yc,(n,))[9], ',' , np.reshape(basement,(n,))[9], '],'print '[', np.reshape(yc,(n,))[15], ',' , np.reshape(basement,(n,))[15], '],'print '[', np.reshape(yc,(n,))[22], ',' , np.reshape(basement,(n,))[22], '],'print '[', np.reshape(yc,(n,))[42], ',' , np.reshape(basement,(n,))[42], '],'print '[', np.reshape(yc,(n,))[51], ',' , np.reshape(basement,(n,))[51], '],'print '[', np.reshape(yc,(n,))[59], ',' , np.reshape(basement,(n,))[59], '],'print '[', np.reshape(yc,(n,))[67], ',' , np.reshape(basement,(n,))[67], '],'print '[', np.reshape(yc,(n,))[91], ',' , np.reshape(basement,(n,))[91], '],'print '[', np.reshape(yc,(n,))[101], ',' , np.reshape(basement,(n,))[101], '],'print '[', np.reshape(yc,(n,))[112], ',' , np.reshape(basement,(n,))[112], '],'print '[', np.reshape(yc,(n,))[121], ',' , np.reshape(basement,(n,))[121], '],'print '[', np.reshape(yc,(n,))[132], ',' , np.reshape(basement,(n,))[132], ']'  Basement surface known valuesbasement_known_more = np.array([[ 8996.6442953 , 2571.60124265 ],[ 24419.4630872 , 3099.87458409 ],[ 39842.2818792 , 11564.2263147 ],[ 57835.5704698 , 21075.5904256 ],[ 109244.966443 , 20108.2305741 ],[ 132379.194631 , 12754.5239185 ],[ 152942.95302 , 12786.3354849 ],[ 173506.711409 , 13148.3430571 ],[ 235197.986577 , 16497.3157779 ],[ 260902.684564 , 13800.0774904 ],[ 289177.852349 , 15007.7089799 ],[ 312312.080537 , 16034.9910351 ],[ 340587.248322 , 6888.23678696 ]])np.savetxt('../data/volcanic-margin-basement-more-known-depths.txt', basement_known_more, fmt='%.18f') Basement surface known valuesbasement_known_new = np.array([[ 39842.2818792 , 11564.2263147 ],[ 57835.5704698 , 21075.5904256 ],[ 132379.194631 , 12754.5239185 ],[ 152942.95302 , 12786.3354849 ],[ 173506.711409 , 13148.3430571 ],[ 260902.684564 , 13800.0774904 ],[ 312312.080537 , 16034.9910351 ]])np.savetxt('../data/volcanic-margin-basement-new-known-depths.txt', base_known_new, fmt='%.18f')  Basement surface known valuesbasement_known_more = np.array([[ 8996.6442953 , 2571.60124265 ],[ 24419.4630872 , 3099.87458409 ],[ 57835.5704698 , 21075.5904256 ],[ 109244.966443 , 20108.2305741 ],[ 152942.95302 , 12786.3354849 ],[ 235197.986577 , 16497.3157779 ],[ 260902.684564 , 13800.0774904 ],[ 289177.852349 , 15007.7089799 ],[ 340587.248322 , 6888.23678696 ]])np.savetxt('../data/volcanic-margin-basement-few-more-known-depths.txt', basement_known_more, fmt='%.18f')basement_known_new = np.array([[ 57835.5704698 , 21075.5904256 ],[ 152942.95302 , 12786.3354849 ],[ 260902.684564 , 13800.0774904 ]])np.savetxt('../data/volcanic-margin-basement-few-new-known-depths.txt', basement_known_new, fmt='%.18f')<jupyter_code>## Known depths (moho)<jupyter_output><empty_output><jupyter_text>print '[', np.reshape(yc,(n,))[59], ',' , np.reshape(moho,(n,))[59], '],'print '[', np.reshape(yc,(n,))[85], ',' , np.reshape(moho,(n,))[85], '],'print '[', np.reshape(yc,(n,))[133], ',' , np.reshape(moho,(n,))[133], ']'<jupyter_code># Moho surface known values
moho_known = np.array([[ 152942.95302 , 21422.6569741 ],
[ 219775.167785 , 21272.9746543 ],
[ 343157.718121 , 18113.1167328 ]])

#np.savetxt('../data/volcanic-margin-moho-known-depths.txt', moho_known, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>Synthetic data<jupyter_code># parameters vector
psyn = []
psyn = np.vstack((ts1, 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,ts0,two_layers=True)
prism_c_syn = fc.prism_c_function(xmax,xmin,dy,edge,S0,dcc,dc,tw,psyn,yc,ts0,two_layers=True)
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)
gzs0_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_s_syn[1])
gzs1_syn = prism.gz(np.reshape(x,(n,)),np.reshape(yc,(n,)),np.reshape(z,(n,)),prism_s_syn[2])
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 + gzs0_syn + gzs1_syn + gzc_syn + gzm_syn
#gsyn = fc.g_function(x,yc,z,gzw_syn,prism_s_syn,prism_c_syn,prism_m_syn)
np.savetxt('../data/magma-poor-margin-synthetic-gravity-data.txt', gsyn, fmt='%.18f')<jupyter_output><empty_output><jupyter_text>Model plot<jupyter_code>polygons_water = []
for (yi, twi) in zip(yc, tw):
    y1 = yi - 0.5*dy
    y2 = yi + 0.5*dy

    polygons_water.append(Polygon(np.array([[y1, y2, y2, y1], 
                                            [0.0, 0.0, twi, twi]]).T,
                                               props={'density': dw - dcc}))
polygons_sediments0 = []
for (yi, twi, s0i) in zip(yc, np.reshape(tw,(n,)), np.reshape(toi,(n,))):
    y1 = yi - 0.5*dy
    y2 = yi + 0.5*dy

    polygons_sediments0.append(Polygon(np.array([[y1, y2, y2, y1], 
                                                [twi, twi, s0i, s0i]]).T,
                                                    props={'density': ds0 - dcc}))
    
polygons_sediments1 = []
for (yi, s0i, s1i) in zip(yc, np.reshape(toi,(n,)), np.reshape(basement,(n,))):
    y1 = yi - 0.5*dy
    y2 = yi + 0.5*dy

    polygons_sediments1.append(Polygon(np.array([[y1, y2, y2, y1], 
                                                [s0i, s0i, s1i, s1i]]).T,
                                                   props={'density': ds1 - dcc}))
    
polygons_crust = []
for (yi, si, Si, dci) in zip(yc, np.reshape(basement,(n,)), np.reshape(moho,(n,)), dc):
    y1 = yi - 0.5*dy
    y2 = yi + 0.5*dy

    polygons_crust.append(Polygon(np.array([[y1, y2, y2, y1], 
                                            [si, si, Si, Si]]).T,
                                                props={'density': dci - dcc})) 

polygons_mantle = []
for (yi, Si) in zip(yc, np.reshape(moho,(n,))):
    y1 = yi - 0.5*dy
    y2 = yi + 0.5*dy

    polygons_mantle.append(Polygon(np.array([[y1, y2, y2, y1], 
                                                [Si, Si, S0+dS0, S0+dS0]]).T,
                                                    props={'density': dm - dcc}))<jupyter_output><empty_output><jupyter_text>aux2 = yc <= COTaux3 = yc >= COT-(dy)ybottom = ((moho[n-len(yc[aux3])]/zmax) - 1.)*(-1) - 0.02ytop = ((basement[n-len(yc[aux3])]/zmax) - 1.)*(-1) - 0.03print ybottomprint ytop<jupyter_code>%matplotlib inline
plt.close('all')
fig = plt.figure(figsize=(12,7))

import matplotlib.gridspec as gridspec
heights = [8, 1]
gs = gridspec.GridSpec(2, 1, height_ratios=heights)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

ax1.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=1)
aux = yc <= COT
for (pwi) in (polygons_water):
    tmpx = [x for x in pwi.x]
    tmpx.append(pwi.x[0])
    tmpy = [y for y in pwi.y]
    tmpy.append(pwi.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='lightskyblue')
for (ps0i) in (polygons_sediments0):
    tmpx = [x for x in ps0i.x]
    tmpx.append(ps0i.x[0])
    tmpy = [y for y in ps0i.y]
    tmpy.append(ps0i.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='tan')
for (ps1i) in (polygons_sediments1):
    tmpx = [x for x in ps1i.x]
    tmpx.append(ps1i.x[0])
    tmpy = [y for y in ps1i.y]
    tmpy.append(ps1i.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='rosybrown')
for (pci) in (polygons_crust[:len(yc[aux])]):
    tmpx = [x for x in pci.x]
    tmpx.append(pci.x[0])
    tmpy = [y for y in pci.y]
    tmpy.append(pci.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='orange')
for (pcoi) in (polygons_crust[len(yc[aux]):n]):
    tmpx = [x for x in pcoi.x]
    tmpx.append(pcoi.x[0])
    tmpy = [y for y in pcoi.y]
    tmpy.append(pcoi.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='olive')
for (pmi) in (polygons_mantle):
    tmpx = [x for x in pmi.x]
    tmpx.append(pmi.x[0])
    tmpy = [y for y in pmi.y]
    tmpy.append(pmi.y[0])
    ax1.plot(tmpx, tmpy, linestyle='None')
    ax1.fill(tmpx, tmpy, color='pink')

#ax1.axhline(y=S0, xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=3)
#ax1.axvline(x=COT, ymin=ybottom, ymax=ytop, linestyle=':', color='w', linewidth=3, label='COT')
ax1.plot(yc, tw, '-k', linewidth=3)
ax1.plot(yc, toi, '-k', linewidth=3)
ax1.plot(yc, basement, '-k', linewidth=3, label='true surfaces')
ax1.plot(yc, moho, '-k', linewidth=3)
ax1.axhline(y=S0+dS0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=3)
#ax1.set_ylim((S0+dS0), zmin)
ax1.set_ylim((56000.0), zmin)
ax1.set_xlim(ymin, ymax)
ax1.set_xlabel('y (km)', fontsize=16)
ax1.set_ylabel('z (km)', fontsize=16)
ax1.set_xticklabels(['%g'% (0.001*l) for l in ax1.get_xticks()], fontsize=14)
ax1.set_yticklabels(['%g'% (0.001*l) for l in ax1.get_yticks()], fontsize=14)
#ax1.legend(loc='lower right', fontsize=14, facecolor='silver')    

X, Y = fig.get_dpi()*fig.get_size_inches()
plt.title('Density (kg/m$^{3}$)', fontsize=18)
ax2.axis('off')

layers_list1 = ['water', 'sediment 1', 'sediment 2', 'continental', 'oceanic', 'mantle']
layers_list2 = ['', '', '', 'crust', 'crust', '']
colors_list = ['lightskyblue', 'tan', 'rosybrown', 'orange', 'olive', 'pink']
density_list = ['1030', '2350', '2600', '2870', '2885', '3300'] #original
ncols = len(colors_list)
nrows = 1
h = Y / nrows
w = X / (ncols + 1)

i=ncols-1
for color, density, layers1, layers2 in zip(colors_list, density_list, layers_list1, layers_list2):
    col = i // nrows
    row = i % nrows
    x = X - (col*w) - w
    
    yi_line = Y
    yf_line = Y - Y*0.15
    yi_text1 = Y - Y*0.2
    yi_text2 = Y - Y*0.27
    yi_text3 = Y - Y*0.08
    
    i-=1
    poly = Polygon(np.array([[x, x+w*0.75, x+w*0.75, x], [yi_line, yi_line, yf_line, yf_line]]).T)
    tmpx = [x for x in poly.x]
    tmpx.append(poly.x[0])
    tmpy = [y for y in poly.y]
    tmpy.append(poly.y[0])
    ax2.plot(tmpx, tmpy, linestyle='-', color='k', linewidth=1)
    ax2.fill(tmpx, tmpy, color=color)
    ax2.text(x+w*0.375, yi_text1, layers1, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax2.text(x+w*0.375, yi_text2, layers2, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax2.text(x+w*0.375, yi_text3, density, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='center')

plt.tight_layout()    
    
#mpl.savefig('../manuscript/figures/magma-poor-margin-true-model.png', dpi='figure', bbox_inches='tight')
plt.show()
%matplotlib inline
plt.close('all')
fig = plt.figure(figsize=(12,13))

import matplotlib.gridspec as gridspec
heights = [8, 8, 1]
gs = gridspec.GridSpec(3, 1, height_ratios=heights)
ax3 = plt.subplot(gs[0])
ax4 = plt.subplot(gs[1])
ax5 = plt.subplot(gs[2])

ax3.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='--', linewidth=1)
#ax3.plot(0.001*yc, gobs, 'or', mfc='none', markersize=8, label='observed data')
ax3.plot(0.001*yc, gsyn, '-g', linewidth=2, label='simulated data')
ax3.set_xlim(0.001*ymin, 0.001*ymax)
ax3.set_ylabel('gravity disturbance (mGal)', fontsize=16)
ax3.set_xticklabels(['%g'% (l) for l in ax3.get_xticks()], fontsize=14)
ax3.set_yticklabels(['%g'% (l) for l in ax3.get_yticks()], fontsize=14)
ax3.legend(loc='best', fontsize=14, facecolor='silver')

ax4.axhline(y=0.0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=1)
aux = yc <= COT
for (pwi) in (polygons_water):
    tmpx = [x for x in pwi.x]
    tmpx.append(pwi.x[0])
    tmpy = [y for y in pwi.y]
    tmpy.append(pwi.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='lightskyblue')
for (ps0i) in (polygons_sediments0):
    tmpx = [x for x in ps0i.x]
    tmpx.append(ps0i.x[0])
    tmpy = [y for y in ps0i.y]
    tmpy.append(ps0i.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='tan')
for (ps1i) in (polygons_sediments1):
    tmpx = [x for x in ps1i.x]
    tmpx.append(ps1i.x[0])
    tmpy = [y for y in ps1i.y]
    tmpy.append(ps1i.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='rosybrown')
for (pci) in (polygons_crust[:len(yc[aux])]):
    tmpx = [x for x in pci.x]
    tmpx.append(pci.x[0])
    tmpy = [y for y in pci.y]
    tmpy.append(pci.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='orange')
for (pcoi) in (polygons_crust[len(yc[aux]):n]):
    tmpx = [x for x in pcoi.x]
    tmpx.append(pcoi.x[0])
    tmpy = [y for y in pcoi.y]
    tmpy.append(pcoi.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='olive')
for (pmi) in (polygons_mantle):
    tmpx = [x for x in pmi.x]
    tmpx.append(pmi.x[0])
    tmpy = [y for y in pmi.y]
    tmpy.append(pmi.y[0])
    ax4.plot(tmpx, tmpy, linestyle='None')
    ax4.fill(tmpx, tmpy, color='pink')

#ax4.axhline(y=S0, xmin=ymin, xmax=ymax, color='w', linestyle='--', linewidth=3)
#ax4.axvline(x=COT, ymin=ybottom, ymax=ytop, linestyle=':', color='w', linewidth=3, label='COT')
ax4.plot(yc, tw, '-k', linewidth=3)
ax4.plot(yc, toi, '-k', linewidth=3)
ax4.plot(yc, basement, '-k', linewidth=3, label='true surfaces')
ax4.plot(yc, moho, '-k', linewidth=3)
ax4.axhline(y=S0+dS0, xmin=ymin, xmax=ymax, color='k', linestyle='-', linewidth=3)
#ax4.set_ylim((S0+dS0), zmin)
ax4.set_ylim((56000.0), zmin)
ax4.set_xlim(ymin, ymax)
ax4.set_xlabel('y (km)', fontsize=16)
ax4.set_ylabel('z (km)', fontsize=16)
ax4.set_xticklabels(['%g'% (0.001*l) for l in ax4.get_xticks()], fontsize=14)
ax4.set_yticklabels(['%g'% (0.001*l) for l in ax4.get_yticks()], fontsize=14)
#ax4.legend(loc='lower right', fontsize=14, facecolor='silver')    

X, Y = fig.get_dpi()*fig.get_size_inches()
plt.title('Density (kg/m$^{3}$)', fontsize=18)
ax5.axis('off')

layers_list1 = ['water', 'sediment 1', 'sediment 2', 'continental', 'oceanic', 'mantle']
layers_list2 = ['', '', '', 'crust', 'crust', '']
colors_list = ['lightskyblue', 'tan', 'rosybrown', 'orange', 'olive', 'pink']
density_list = ['1030', '2350', '2600', '2870', '2885', '3300'] #original
ncols = len(colors_list)
nrows = 1
h = Y / nrows
w = X / (ncols + 1)

i=ncols-1
for color, density, layers1, layers2 in zip(colors_list, density_list, layers_list1, layers_list2):
    col = i // nrows
    row = i % nrows
    x = X - (col*w) - w
    
    yi_line = Y
    yf_line = Y - Y*0.15
    yi_text1 = Y - Y*0.2
    yi_text2 = Y - Y*0.27
    yi_text3 = Y - Y*0.08
    
    i-=1
    poly = Polygon(np.array([[x, x+w*0.75, x+w*0.75, x], [yi_line, yi_line, yf_line, yf_line]]).T)
    tmpx = [x for x in poly.x]
    tmpx.append(poly.x[0])
    tmpy = [y for y in poly.y]
    tmpy.append(poly.y[0])
    ax5.plot(tmpx, tmpy, linestyle='-', color='k', linewidth=1)
    ax5.fill(tmpx, tmpy, color=color)
    ax5.text(x+w*0.375, yi_text1, layers1, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax5.text(x+w*0.375, yi_text2, layers2, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='top')
    ax5.text(x+w*0.375, yi_text3, density, fontsize=(w*0.14), horizontalalignment='center', verticalalignment='center')

plt.tight_layout()    
    
#mpl.savefig('../manuscript/figures/magma-poor-margin-true-model-grafics.png', dpi='figure', bbox_inches='tight')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.184731327737271, 'fraction_numerical': 0.20881073463635808, 'mean_word_length': 3.759148563930602, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '317191', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 18044, 'n_tokens_neox': 13927, 'n_words': 2133, 'path': 'code/magma-poor-margin-model.ipynb', 'repo_name': 'pinga-lab/paper-isostatic-grav2D', 'size': '32662'}