-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtesttillpressure.py
103 lines (85 loc) · 1.8 KB
/
testtillpressure.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
"""
This script compares the pressure obtained from tilleos.c with the one from tillotson.c.
obtained from gendatotipsy.
"""
from matplotlib import *
rcParams['font.family'] = 'serif'
from numpy import *
from matplotlib.pyplot import *
"""
data1 = loadtxt('out.old.dat',skiprows=3)
data2 = loadtxt('out.dat',skiprows=3)
"""
# Tillotson material
data1 = loadtxt('testtillpressure.txt')
data2 = loadtxt('testtillpressure2.txt')
"""
The region where the pressure is zero.
"""
np = loadtxt('till_granite_np.txt')
rho_np = np[:,0]
u_np = np[:,1]
print size(np)
"""
# The cold curve for basalt.
coldcurve = loadtxt('/home/ics/creinh/code/ballic-test/tools/cold_basalt.txt')
# Values for basalt (in code units!)
us = 4.72084
us2 = 18.2032
rho0 = 7.32745
"""
# The cold curve for granite.
coldcurve = loadtxt('/home/ics/creinh/code/ballic-test/tools/cold_granite.txt')
# Values for granite (in code units!)
us = 3.5
us2 = 18.0
rho0 = 7.33
rhocold = coldcurve[:,0]
ucold = coldcurve[:,1]
"""
data3 = loadtxt('testtillpressure3.txt')
rho = data3[:,0]
u = data3[:,1]
"""
print where(fabs(data1-data2) > 1e-30)
print
index = where(fabs(data1-data2) > 1e-30)
print data1[index]
print
print data2[index]
figure(1)
imshow(data1-data2)
xlim(0,20)
ylim(0,7.33)
xlabel('Int. energy')
ylabel('Density')
colorbar()
figure(2)
imshow(data1)
xlabel('Int. energy')
ylabel('Density')
colorbar()
figure(3)
imshow(data2)
xlabel('Int. energy')
ylabel('Density')
colorbar()
"""
#
# Plot the points where the pressure differ in the u(rho) plane.
#
figure(4)
scatter(rho,u,color='blue')
plot(rhocold,ucold,color='red')
xmax = max(rho)*1.1
ymax = max(u)*1.1
plot(rho_np,u_np,color='red')
xlim(0,xmax)
ylim(0,ymax)
plot([0,rho0],[us,us],'r--')
plot([0,rho0],[us2,us2],'r--')
plot([rho0,rho0],[0,ymax],'r--')
xlabel('Density')
ylabel('Int. energy')
"""
show()