1
2
3
4 import g_pipe
5 import math
6
7
8 IB = 0.7
9
10 -def addcurve( xyee, x_range, y_range, bexp, brms, imopfd, xsz, ysz ):
11 if len(xyee) == 0:
12 return
13 xmin, xmax = x_range
14 sx = xsz/float(xmax - xmin)
15 ymin, ymax = y_range
16 sy = ysz/float(ymax - ymin)
17 bs2 = 0.0
18 for (x, y, ex, ey, wt) in xyee:
19 iex = math.hypot(ex*sx, IB)
20 iey = math.hypot(ey*sy, IB)
21 area = iex*iey
22 b = wt * area**(-bexp)
23 bs2 += b*b
24 if not ( bs2 > 0.0):
25 return
26 bfac = brms / math.sqrt(bs2/len(xyee))
27 for (x, y, ex, ey, wt) in xyee:
28 iex = math.hypot(ex*sx, IB)
29 iey = math.hypot(ey*sy, IB)
30 area = iex*iey
31 b = bfac * wt * area**(-bexp)
32 xscaled = (x-xmin)*sx
33 yscaled = (y-ymin)*sy
34 r = math.sqrt(iex * iey)
35 ecc = math.sqrt(1 - min(iex, iey)/max(iex, iey))
36 pa = 90
37 if abs(iey) < abs(iex):
38 pa = 0
39 imopfd.writelines('%.2f %.2f %.3f %g %g %d addGaussGal\n'
40 % (xscaled, yscaled, b, r, ecc, pa))
41
42
43
44
45
46 if __name__ == '__main__':
47 H = 600
48 W = 900
49 si, so = g_pipe.popen2("imop", ['imop', '-write16', '-stdin' ])
50 x_range = (0, 1)
51 y_range = (0, 1)
52 xyee = []
53 for i in range(10):
54 x = 0.1*i
55 y = x - 2*x*x*x + 2*x*x*x*x*x
56 ex = 0.05 - 0.02*x
57 ey = 0.01 + 0.1*y
58 xyee.append( (x, y, ex, ey, 1.0) )
59 si.writelines('%d %d alloc\n' % (H, W))
60 addcurve(xyee, x_range, y_range, 0.5, 200, si, W, H)
61 si.writelines('= foo.fits')
62 si.close()
63 so.close()
64