<jupyter_start><jupyter_text>Synthetic Model  Notebook to generate a set of models.  Import libraries<jupyter_code>from IPython.display import Markdown as md
from IPython.display import display as dp
import string as st
import sys
import numpy as np
import matplotlib.pyplot as plt
import cPickle as pickle
import datetime

from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
from fatiando.utils import ang2vec, vec2ang
from fatiando.mesher import Sphere, Prism, PolygonalPrism
from fatiando.gravmag import sphere, prism, polyprism
notebook_name = 'synthetic_model.ipynb'<jupyter_output><empty_output><jupyter_text>Importing auxiliary functions<jupyter_code>dir_modules = '../../../mypackage'
sys.path.append(dir_modules)
import auxiliary_functions as func<jupyter_output><empty_output><jupyter_text>Loading 2D grid properties<jupyter_code>with open('data/regular_grid.pickle') as f:
        regular = pickle.load(f)
with open('data/airborne_survey.pickle') as f:
        airborne = pickle.load(f)<jupyter_output><empty_output><jupyter_text>List of saved files<jupyter_code>saved_files = []<jupyter_output><empty_output><jupyter_text>Observation area<jupyter_code>print 'Area limits: \n x_max = %.1f m \n x_min = %.1f m \n y_max = %.1f m \n y_min = %.1f m' % (regular['area'][1], regular['area'][0],regular['area'][3],regular['area'][2])<jupyter_output>Area limits: 
 x_max = 5500.0 m 
 x_min = -6500.0 m 
 y_max = 6500.0 m 
 y_min = -5500.0 m<jupyter_text>Regional Field<jupyter_code>inc_gf , dec_gf = (-40.,-22.) ## https://geomag.bgs.ac.uk/data_service/models_compass/wmm_calc.html no RJ<jupyter_output><empty_output><jupyter_text>Create a model w/ the presence of Remanence Magnetization<jupyter_code>model_multi = dict()<jupyter_output><empty_output><jupyter_text>main field<jupyter_code>model_multi['main_field'] = (inc_gf,dec_gf)<jupyter_output><empty_output><jupyter_text>Setting bounds of a polygonal cross-section<jupyter_code>% matplotlib notebook
fig = plt.figure(figsize=(7,7), tight_layout=True)

ax = fig.add_subplot(111)
ax.set_title('Click on figure to set any point')
ax.axis([regular['y'].min(), regular['y'].max(), regular['x'].min(), regular['x'].max()])
ax.set_ylabel('x(m)')
ax.set_xlabel('y(m)')


line, = ax.plot([], [])
line.figure.canvas.draw()
x = []
y = []
plotx = []
ploty = []
def pick(event):
    if event.inaxes != ax.axes:
        return 'It must be any point over de area '
    x.append(event.ydata)
    y.append(event.xdata)
    plotx.append(event.xdata)
    ploty.append(event.ydata)
    line.set_color('r')
    line.set_marker('o')
    line.set_linestyle('')
    line.set_data(plotx,ploty)
    line.figure.canvas.draw()

line.figure.canvas.mpl_connect('button_press_event', pick)

plt.show()

print x
print y<jupyter_output>[]<jupyter_text>Magnetization Polygonal prism<jupyter_code>model_multi['m_Rpp'] = 4.
model_multi['inc_R'] = -25.
model_multi['dec_R'] = 30.

mag_tot_Rpp = ang2vec(model_multi['m_Rpp'],
                  model_multi['inc_R'],
                  model_multi['dec_R'])

model_multi['magnetization_pp'] = mag_tot_Rpp<jupyter_output><empty_output><jupyter_text>Vertices of a Polygonal prism-1<jupyter_code>model_multi['x1_verts'] = [-120.86173387456893, 217.88504245497916, -139.68099922621059, -1005.3672054017161, -2002.788269038715, -3019.0285980273529, -3884.7148042028603, -4204.6423151807649, -3809.4377427962945, -2360.3543107199021, -1494.6681045443956, -666.62042907216983, -120.86173387456893]
model_multi['y1_verts'] = [-1475.5336403581159, -822.01437309183802, -227.90594830431655, -386.33486158098731, -623.97823149599935, -604.17461733641267, -227.90594830431655, -1020.0505146876794, -1594.3553253156206, -1574.5517111560348, -1297.3011129218594, -1396.3191837197774, -1475.5336403581159]
model_multi['z1_top'] = 450.
model_multi['z1_bottom'] = 3150.
model_multi['verts1'] = zip(model_multi['x1_verts'],
                                   model_multi['y1_verts'] )<jupyter_output><empty_output><jupyter_text>Creating a polyprism model<jupyter_code>model_multi['polygons'] = [PolygonalPrism(model_multi['verts1'],
                                              model_multi['z1_top'],
                                              model_multi['z1_bottom'],
                                              {'magnetization':mag_tot_Rpp})]<jupyter_output><empty_output><jupyter_text>Magnetization spheres<jupyter_code>model_multi['m_Rs'] = 3.
model_multi['inc_R'] = -25.
model_multi['dec_R'] = 30.


mag_tot_Rs = ang2vec(model_multi['m_Rs'],
                  model_multi['inc_R'],
                  model_multi['dec_R'])

model_multi['magnetization_s'] = mag_tot_Rs<jupyter_output><empty_output><jupyter_text>Sphere-1<jupyter_code>model_multi['xc1'] = 1800.
model_multi['yc1'] = -1800.
model_multi['zc1'] = 1000.
model_multi['radius1'] = 500.<jupyter_output><empty_output><jupyter_text>Sphere-2<jupyter_code>model_multi['xc2'] = 800.
model_multi['yc2'] = 800.
model_multi['zc2'] = 1000.
model_multi['radius2'] = 500.<jupyter_output><empty_output><jupyter_text>Creating a spheres model<jupyter_code>model_multi['spheres'] = [Sphere(model_multi['xc1'],
                                         model_multi['yc1'],
                                         model_multi['zc1'],
                                         model_multi['radius1'],
                                         {'magnetization':mag_tot_Rs}),
                            Sphere(model_multi['xc2'],
                                         model_multi['yc2'],
                                         model_multi['zc2'],
                                         model_multi['radius2'],
                                         {'magnetization':mag_tot_Rs})]
theta = np.linspace(0, 2 * np.pi, 100)<jupyter_output><empty_output><jupyter_text>Magnetization of rectangular prisms<jupyter_code>model_multi['m_Rp'] = 2.5
model_multi['inc_R'] = -25.
model_multi['dec_R'] = 30.


mag_tot_Rp = ang2vec(model_multi['m_Rp'],
                  model_multi['inc_R'],
                  model_multi['dec_R'])

model_multi['magnetization_p'] = mag_tot_Rp<jupyter_output><empty_output><jupyter_text>Rectangular prisms  Prism-1<jupyter_code>model_multi['x1_max'] = 3200.
model_multi['x1_min'] = 2200.
model_multi['y1_max'] = 3000.
model_multi['y1_min'] = 2300
model_multi['z1_min'] = 450.
model_multi['z1_max'] = 950. 
model_multi['x_verts1'] = [model_multi['x1_max'],model_multi['x1_min'],model_multi['x1_min'],model_multi['x1_max'],model_multi['x1_max']]
model_multi['y_verts1'] = [model_multi['y1_max'],model_multi['y1_max'],model_multi['y1_min'],model_multi['y1_min'],model_multi['y1_max']]<jupyter_output><empty_output><jupyter_text>Prism-2<jupyter_code>model_multi['x2_max'] = -2000.
model_multi['x2_min'] = -3000.
model_multi['y2_max'] = 4000.
model_multi['y2_min'] = 2000.
model_multi['z2_min'] = 500.
model_multi['z2_max'] = 2050. 
model_multi['x_verts2'] = [model_multi['x2_max'],model_multi['x2_min'],model_multi['x2_min'],model_multi['x2_max'],model_multi['x2_max']]
model_multi['y_verts2'] = [model_multi['y2_max'],model_multi['y2_max'],model_multi['y2_min'],model_multi['y2_min'],model_multi['y2_max']]<jupyter_output><empty_output><jupyter_text>Creating model<jupyter_code>model_multi['prisms'] = [Prism(model_multi['x1_min'],model_multi['x1_max'],
                               model_multi['y1_min'],model_multi['y1_max'],
                               model_multi['z1_min'],model_multi['z1_max'],
                               {'magnetization':mag_tot_Rp}),
                         Prism(model_multi['x2_min'],model_multi['x2_max'],
                               model_multi['y2_min'],model_multi['y2_max'],
                               model_multi['z2_min'],model_multi['z2_max'],
                               {'magnetization':mag_tot_Rp})]<jupyter_output><empty_output><jupyter_text>Generating .pickle file<jupyter_code>now = datetime.datetime.utcnow().strftime('%d %B %Y %H:%M:%S UTC')
model_multi['metadata'] = 'Generated by {name} on {date}'.format(date=now, name=notebook_name)
file_name = 'data/model_multi.pickle'
with open(file_name, 'w') as f:
    pickle.dump(model_multi, f)
    
saved_files.append(file_name)<jupyter_output><empty_output><jupyter_text>Visualization of model projection<jupyter_code>title_font = 20
bottom_font = 18
saturation_factor = 1.
plt.close('all')
plt.figure(figsize=(9,9), tight_layout=True)

plt.title('Model projection',fontsize=title_font)

plt.plot(model_multi['y1_verts'],model_multi['x1_verts'], 
         color='k',linestyle='-',linewidth=2)

plt.plot(model_multi['y_verts1'],model_multi['x_verts1'], 
         color='k',linestyle='-',linewidth=2)
plt.plot(model_multi['y_verts2'],model_multi['x_verts2'], 
         color='k',linestyle='-',linewidth=2)


plt.plot(model_multi['radius1']*np.sin(theta)+ model_multi['yc1'] ,
         model_multi['radius1']*np.cos(theta)+ model_multi['xc1'], 
         color='k',linestyle='-',linewidth=2)
plt.plot(model_multi['radius2']*np.sin(theta)+ model_multi['yc2'] ,
         model_multi['radius2']*np.cos(theta)+ model_multi['xc2'], 
         color='k',linestyle='-',linewidth=2)

plt.xlabel('y (km)', fontsize = title_font)
plt.ylabel('x (km)', fontsize = title_font)
plt.ylim(np.min(airborne['x']),np.max(airborne['x']))
plt.xlim(np.min(airborne['y']),np.max(airborne['y']))
plt.tick_params(labelsize=15)

file_name = 'figs/model_projection'
plt.savefig(file_name+'.png',dpi=600)
saved_files.append(file_name+'.png')

plt.savefig(file_name+'.eps',dpi=600)
saved_files.append(file_name+'.eps')


plt.show()<jupyter_output><empty_output><jupyter_text>Saved files<jupyter_code>with open('reports/report_%s.md' % notebook_name[:st.index(notebook_name, '.')], 'w') as q:
    q.write('# Saved files \n')
    now = datetime.datetime.utcnow().strftime('%d %B %Y %H:%M:%S UTC')
    header = 'Generated by {name} on {date}'.format(date=now, name=notebook_name)
    q.write('\n\n'+header+'\n\n')
    for i, sf in enumerate(saved_files):
        print '%d  %s' % (i+1,sf)
        q.write('*  `%s` \n' % (sf))<jupyter_output>1  data/model_multi.pickle
2  figs/model_projection.png
3  figs/model_projection.eps
{'fraction_non_alphanumeric': 0.1570572635301051, 'fraction_numerical': 0.06954130242608782, 'mean_word_length': 4.353312302839116, 'pattern_counts': {'":': 0, '<': 90, '<?xml version=': 0, '>': 90, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '247338', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4222, 'n_tokens_neox': 3787, 'n_words': 480, 'path': 'code/notebooks/synthetic_tests/model_multibody/.ipynb_checkpoints/synthetic_model-checkpoint.ipynb', 'repo_name': 'pinga-lab/eqlayer-magnetization-direction', 'size': '10181'}