In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from numpy import poly1d, polyfit, power
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd
import PIL as pil
import heapq
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 18, 14
# import seaborn as sns
# sns.set_palette("deep", desat=.6)
try:
from PIL import Image
except:
import Image
incl = 30.
sin_i2 = np.sin(incl*np.pi/180.)**2
cos_i2 = np.cos(incl*np.pi/180.)**2
In [2]:
os.chdir("C:\\science\\2FInstability\\data\\ngc1068")
plt.imshow(np.asarray(Image.open("shapiro_fg3a.png")))
plt.plot()
Out[2]:
In [3]:
plt.imshow(np.asarray(Image.open("shapiro_fg4a.png")))
plt.plot()
Out[3]:
In [4]:
pylab.rcParams['figure.figsize'] = 15, 5
#Velocity
lu = (215, 67)
lu_val = (30, 200)
rd = (1217, 463)
rd_val = (80, 25)
data = [(389,259), (411,257), (467,241), (513,249), (575,263), (589,245), (693,275), (699,261), (809,269), (871,293), (953,293)]
def extend_values(data, lu, lu_val, rd, rd_val):
xscale = 1.0*(rd_val[0]-lu_val[0])/(rd[0]-lu[0])
yscale = 1.0*(rd_val[1]-lu_val[1])/(lu[1]-rd[1])
extended = []
for d in data:
extended.append((xscale*(d[0]-lu[0])+lu_val[0], yscale*(rd[1] - d[1])+rd_val[1]))
return extended
e_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*e_data)[0],zip(*e_data)[1], 'or')
vel_fit = lambda l: 356.*np.power(l, -0.21)
plt.plot(np.linspace(30., 90., 100), map(vel_fit, np.linspace(30., 90., 100)), '--')
plt.ylim(25, 225)
plt.xlim(30, 95)
plt.show()
In [5]:
#Sig_maj
lu = (219, 509)
lu_val = (30, 120)
rd = (1215, 913)
rd_val = (80, 30)
data = [(387,689), (409,699), (463,721), (471,753), (515,731), (575,765), (587,783), (691,763), (699,813), (809,797), (873,819), (953,775)]
sig_maj_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_maj_data)[0],zip(*sig_maj_data)[1], 'or')
sigR = lambda l: 213.*np.exp(-l/72.)
sigZ = lambda l: 124.*np.exp(-l/72.)
phi_to_R = lambda l: 0.5*(1 - 0.21)
sigR_min = lambda l: 193.*np.exp(-l/66.)
sigZ_min = lambda l: 115.*np.exp(-l/66.)
phi_to_R_min = lambda l: 0.5*(1 - 0.24)
sig_maj = lambda l: sqrt(sigR(l)**2 * (phi_to_R(l) * sin_i2 + sigZ(l)**2 * cos_i2/sigR(l)**2))
sig_maj_min = lambda l: sqrt(sigR_min(l)**2 * (phi_to_R_min(l) * sin_i2 + sigZ_min(l)**2 * cos_i2/sigR_min(l)**2))
plt.plot(np.linspace(30., 90., 100), map(sig_maj, np.linspace(30., 90., 100)), '--')
plt.plot(np.linspace(30., 90., 100), map(sig_maj_min, np.linspace(30., 90., 100)), '-')
plt.ylim(30, 130)
plt.xlim(30, 95)
plt.show()
In [6]:
#Sig_min
lu = (219, 957)
lu_val = (30, 120)
rd = (1217, 1359)
rd_val = (80, 30)
data = [(387,1071), (449,1083), (467,1085), (539,1103), (547,1133), (609,1149), (667,1187), (689,1191), (779,1193), (813,1215),
(875,1205), (963,1229), (971,1237), (1123,1225), (1131,1249), (1331,1269), (1365,1227)]
sig_min_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_min_data)[0], zip(*sig_min_data)[1], 'or')
sig_min = lambda l: sqrt(sigR(l)**2 * sin_i2 + sigZ(l)**2 * cos_i2)
sig_min_min = lambda l: sqrt(sigR_min(l)**2 * sin_i2 + sigZ_min(l)**2 * cos_i2)
plt.plot(np.linspace(30., 90., 100), map(sig_min_min, np.linspace(30., 90., 100)), '-')
plt.plot(np.linspace(30., 90., 100), map(sig_min, np.linspace(30., 90., 100)), '--')
plt.ylim(30, 130)
plt.xlim(30, 95)
plt.show()