-
Notifications
You must be signed in to change notification settings - Fork 0
/
GH_plots.py
executable file
·128 lines (89 loc) · 3.92 KB
/
GH_plots.py
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
## ==================================================================
## Produce Gauss-Hermite plots
## ==================================================================
## warrenj 20150825 Routine to plot h_3 and h_4 vs v/sigma for all
## bins
## warrenj 20150917 Loop added to plot both without editing plot
import numpy as np # for reading files
import glob # for searching for files
from astropy.io import fits as pyfits # reads fits files (is from astropy)
import matplotlib.pyplot as plt # used for plotting
from scipy.stats import gaussian_kde # for calc plot density
from checkcomp import checkcomp
cc = checkcomp()
#-----------------------------------------------------------------------------
def GH_plots(galaxy, wav_range="", plots=False):
if wav_range:
wav_range_dir = wav_range + "/"
else:
wav_range_dir = ""
wav_range = ""
limit=1 # how many of the lowest/highest errors to clip
scale=1 #scale the errors
tessellation_File = "%s/Data/vimosindi/analysis/%s/" %(cc.base_dir, galaxy) +\
"voronoi_2d_binning_output_kin.txt"
tessellation_File2 = "%s/Data/vimosindi/analysis/%s/" %(cc.base_dir, galaxy) +\
"voronoi_2d_binning_output2_kin.txt"
output_v = "%s/Data/vimosindi/analysis/%s/results/" % (cc.base_dir, galaxy) +\
"%sgal_stellar_vel.dat" % (wav_range_dir)
output_v_uncert = output_v
output_sigma = "%s/Data/vimosindi/analysis/%s/results/" % (cc.base_dir, galaxy) +\
"%sgal_stellar_sigma.dat" % (wav_range_dir)
output_sigma_uncert = output_sigma
output_h3 = "%s/Data/vimosindi/analysis/%s/results/" % (cc.base_dir, galaxy) +\
"%sgal_stellar_h3.dat" % (wav_range_dir)
output_h3_uncert = output_h3
output_h4 = "%s/Data/vimosindi/analysis/%s/results/" % (cc.base_dir, galaxy) +\
"%sgal_stellar_h4.dat" % (wav_range_dir)
output_h4_uncert = output_h4
vel, vel_uncert = np.loadtxt(output_v, unpack=True)
sigma, sigma_uncert = np.loadtxt(output_sigma, unpack=True)
h3, h3_uncert = np.loadtxt(output_h3, unpack=True)
h4, h4_uncert = np.loadtxt(output_h4, unpack=True)
vel += -np.mean(vel)
x = vel/sigma
x_uncert = np.sqrt((vel_uncert/vel)**2 + (sigma_uncert/sigma)**2)*x/scale
x_uncert_sorted = sorted(x_uncert)
x_uncert_min = x_uncert_sorted[limit]
x_uncert_max = x_uncert_sorted[-limit-1]
plot = "h3"
for i in [0,1]:
if plot == "h3":
y = h3
y_uncert = h3_uncert/scale
ytitle = r"$h_3$"
else:
y = h4
y_uncert = h4_uncert/scale
ytitle = r"$h_4$"
# Calculate the point density
xy = np.vstack([x,y])
z = gaussian_kde(xy)(xy)
# Sort the points by density, so that the densest points are plotted last
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
#plt.plot(x, y, 'bx')
plt.title("Local " + ytitle + " - (v/sigma) relation")
plt.xlabel(r"$v/\sigma$")
plt.ylabel(ytitle)
plt.errorbar(x, y, xerr=np.clip(x_uncert,x_uncert_min,x_uncert_max),
yerr=y_uncert, fmt='+',zorder=1)
plt.scatter(x, y, c=z, s=50, edgecolor='',zorder=2)
ax=plt.gca()
plt.text(0.02,0.98, "Galaxy: " + galaxy.upper(),
verticalalignment='top',
transform=ax.transAxes)
#plt.savefig("/home/warrenj/Desktop/" + plot + "-(v-sigma)_" + wav_range + \
#".png", bbox_inches="tight")
plt.savefig("%s/Data/vimosindi/analysis/%s/results/" % (cc.base_dir,
galaxy) + "%s/plots/%s-(v-sigma)_%s.png" % (wav_range_dir,
plot, wav_range), bbox_inches="tight")
if plots:
plt.show()
plot = "h4"
##############################################################################
# Use of plot_results.py
if __name__ == '__main__':
galaxy = "ngc3557"
wav_range="4200-"
GH_plots(galaxy, wav_range, plots=True)