pythonFLAC/2.dip_angle_loop.py at main · chingchenn/pythonFLAC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#!/usr/bin/env python
import math
import flac
import os,sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
from matplotlib import cm
import matplotlib.pyplot as plt
import function_for_flac as f2
model_list=['Chi01','chih0606','chih0607','chih0608','chih0609']
name_list=['3','5','7','9','11']
newcolors = ['#2F4F4F','#4682B4','#CD5C5C','#708090','#AE6378','#282130','#7E9680','#24788F','#849DAB','#EA5E51','#35838D','#4198B9','#414F67','#97795D','#6B0D47','#A80359','#52254F']
# depth1=-5
# depth2=-150
# case =1
if len(sys.argv) <= 1:
    print('''
          input [1] = plot style : 1 = plot dip only
                                  : 2 = plot dip with smoothing
          input [2] [3] = dip depth range from [2] to [3]
                          Note that [2] [3] is negative
          ''')
    sys.exit()
case =sys.argv[1]
depth1=sys.argv[2]
depth2=sys.argv[3]
if int(case)==1:
    print(11111111111111111)
    fig, (ax2)= plt.subplots(1,1,figsize=(13,8))
if int(case)==2:
    print(22222222222222222)
    fig, (ax,ax2)= plt.subplots(2,1,figsize=(10,10))
for qq,model in enumerate(model_list):
    path = '/scratch2/jiching/22winter/'+model+'/'
    path = '/home/jiching/geoflac/'+model+'/'
#    path = '/scratch2/jiching/03model/'+model+'/'
    #path = 'D:/model/'+model+'/'
    #path = '/scratch2/jiching/sem02model/'+model+'/'
    #path = '/scratch/jiching/summer2021/week11/'+model+'/'
    #path = '/scratch2/jiching/'+model+'/'
    # path = '/Volumes/My Book/model/'+model+'/'
    #path = '/Volumes/SSD500/model/'+model+'/'
    os.chdir(path)
    fl = flac.Flac();end = fl.nrec
    nex = fl.nx - 1;nez = fl.nz - 1

    phase_oceanic = 3
    phase_ecolgite = 13
    phase_oceanic_1 = 17
    phase_ecolgite_1 = 18
    angle = np.zeros(end)


    for i in range(1,end):
        x, z = fl.read_mesh(i)
        ele_x = (x[:fl.nx-1,:fl.nz-1] + x[1:,:fl.nz-1] + x[1:,1:] + x[:fl.nx-1,1:]) / 4.
        ele_z = (z[:fl.nx-1,:fl.nz-1] + z[1:,:fl.nz-1] + z[1:,1:] + z[:fl.nx-1,1:]) / 4.
        phase = fl.read_phase(i)
        trench_ind = np.argmin(z[:,0])
        crust_x = np.zeros(nex)
        crust_z = np.zeros(nex)
        for j in range(trench_ind,nex):
            ind_oceanic = (phase[j,:] == phase_oceanic) + (phase[j,:] == phase_ecolgite)+(phase[j,:] == phase_oceanic_1) + (phase[j,:] == phase_ecolgite_1)
            if True in ind_oceanic:
                crust_x[j] = np.average(ele_x[j,ind_oceanic])
                crust_z[j] = np.average(ele_z[j,ind_oceanic])

        ind_within_80km = (crust_z >= int(depth2)) * (crust_z < int(depth1))
        if not True in (crust_z < int(depth2)):
            continue

        crust_xmin = np.amin(crust_x[ind_within_80km])
        crust_xmax = np.amax(crust_x[ind_within_80km])
        crust_zmin = np.amin(crust_z[ind_within_80km])
        crust_zmax = np.amax(crust_z[ind_within_80km])
        dx = crust_xmax - crust_xmin
        dz = crust_zmax - crust_zmin
        angle[i] = math.degrees(math.atan(dz/dx))
    if int(case)==2:
        nnewangle = f2.moving_window_smooth(angle[angle>0],3)
        ax.plot(fl.time[angle>0],nnewangle,c=newcolors[qq],lw=3,label=name_list[qq])
        # ax.legend(title = "convergent velocity mm/year")
        ax.set_xlim(0,fl.time[-1])
        ax.set_title('Smoothing Angle Variation')
        ax.set_ylabel('Angel ($^\circ$)')
    ax2.plot(fl.time[angle>0],angle[angle>0],c=newcolors[qq],lw=3,label=name_list[qq])
    ax2.legend(title = "geology zone",fontsize = 16)
    ax2.set_xlim(0,fl.time[-1])
    ax2.grid(axis='y')
    ax2.set_xlabel('Time (Myr)')
    ax2.set_ylabel('Angel ($^\circ$)')
    ax2.set_title('Angle Variation')
plt.savefig('/home/jiching/geoflac/figure/'+'dip_'+model_list[0]+'_'+model_list[1]+'_'+str(depth2)+'km'+'_.png')