-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmeanqu.py
More file actions
executable file
·117 lines (100 loc) · 4.64 KB
/
meanqu.py
File metadata and controls
executable file
·117 lines (100 loc) · 4.64 KB
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
#!/usr/bin/env python
import sys
if len(sys.argv) < 2:
print 'dude, give me a log file. wtf?'
exit(1)
import asciidata, numpy, pylab
# load files
print 'loading ', sys.argv[1]
f = asciidata.AsciiData(sys.argv[1], comment_char='#')
nu = numpy.array(f.columns[0])
q = numpy.array(f.columns[1])
u = numpy.array(f.columns[2])
err = numpy.array(f.columns[3])
# flux model for 3c286 stokes parameters
# init with 1.0 ghz index
trueq = 0.638 * (0.96/nu)**(0.337)
trueu = 1.433 * (0.96/nu)**(0.337)
ind = numpy.where((nu >= 1.39) & (nu <= 1.47))
trueq[ind] = 0.526 * (1.39/nu[ind])**(0.353)
trueu[ind] = 1.263 * (1.39/nu[ind])**(0.353)
ind = numpy.where((nu >= 1.76) & (nu <= 1.84))
trueq[ind] = 0.515 * (1.76/nu[ind])**(0.353)
trueu[ind] = 1.158 * (1.76/nu[ind])**(0.353)
ind = numpy.where((nu >= 1.97) & (nu <= 2.05))
trueq[ind] = 0.494 * (1.97/nu[ind])**(0.410)
trueu[ind] = 1.110 * (1.97/nu[ind])**(0.410)
# fluxweight = numpy.sum(q/err**2)/numpy.sum(1/err**2)
# print 1/numpy.sqrt(numpy.sum(1/err**2))
ind = numpy.where(((nu >= 0.96) & (nu <= 1.04)) | ((nu >= 1.39) & (nu <= 1.47)) | ((nu >= 1.76) & (nu <= 1.84)) | ((nu >= 1.97) & (nu <= 2.05)))
meanq = (q-trueq).mean()
meanu = (u-trueu).mean()
stdq = (q-trueq).std()/numpy.sqrt(len(ind[0]))
stdu = (u-trueu).std()/numpy.sqrt(len(ind[0]))
print 'Total:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)
ind = numpy.where((nu >= 0.96) & (nu <= 1.04))
meanq = (q[ind]-trueq[ind]).mean()
meanu = (u[ind]-trueu[ind]).mean()
stdq = (q[ind]-trueq[ind]).std()/numpy.sqrt(len(ind[0]))
stdu = (u[ind]-trueu[ind]).std()/numpy.sqrt(len(ind[0]))
print '1000:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)
ind = numpy.where((nu >= 1.39) & (nu <= 1.47))
meanq = (q[ind]-trueq[ind]).mean()
meanu = (u[ind]-trueu[ind]).mean()
stdq = (q[ind]-trueq[ind]).std()/numpy.sqrt(len(ind[0]))
stdu = (u[ind]-trueu[ind]).std()/numpy.sqrt(len(ind[0]))
print '1430:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)
ind = numpy.where((nu >= 1.76) & (nu <= 1.84))
meanq = (q[ind]-trueq[ind]).mean()
meanu = (u[ind]-trueu[ind]).mean()
stdq = (q[ind]-trueq[ind]).std()/numpy.sqrt(len(ind[0]))
stdu = (u[ind]-trueu[ind]).std()/numpy.sqrt(len(ind[0]))
print '1800:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)
ind = numpy.where((nu >= 1.97) & (nu <= 2.05))
meanq = (q[ind]-trueq[ind]).mean()
meanu = (u[ind]-trueu[ind]).mean()
stdq = (q[ind]-trueq[ind]).std()/numpy.sqrt(len(ind[0]))
stdu = (u[ind]-trueu[ind]).std()/numpy.sqrt(len(ind[0]))
print '2010:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)
ind = numpy.where(((nu >= 1.39) & (nu <= 1.47)) | ((nu >= 1.97) & (nu <= 2.05)))
meanq = (q[ind]-trueq[ind]).mean()
meanu = (u[ind]-trueu[ind]).mean()
stdq = (q[ind]-trueq[ind]).std()/numpy.sqrt(len(ind[0]))
stdu = (u[ind]-trueu[ind]).std()/numpy.sqrt(len(ind[0]))
print '1430-2010:'
#print '%2.0f\pm%2.0f & %2.0f\pm%2.0f' % (1000*meanq,1000*stdq,1000*meanu,1000*stdu)
p = numpy.sqrt(trueq[ind].mean()**2 + trueu[ind].mean()**2)
perr = numpy.sqrt(meanq**2 + meanu**2)/(p/0.095)
stdp = numpy.sqrt(stdq**2 + stdu**2)/(p/0.095)
therr = numpy.degrees(0.5*numpy.arctan2(meanu,meanq))
print '$%1.2f\pm%1.2f$ & $%3.0f$ ' % (100*perr,100*stdp,therr)