1# This file is part of the Astrometry.net suite.
2# Licensed under a 3-clause BSD style license - see LICENSE
3from __future__ import print_function
4import matplotlib
5matplotlib.use('Agg')
6import os.path
7from pylab import *
8from numpy import *
9# test_tweak 2>tt.py
10import tt
11
12def savefig(fn):
13    from pylab import savefig as sf
14    print('Saving', fn)
15    sf(fn)
16
17if __name__ == '__main__':
18    #print 'me:', __file__
19    #tt = os.path.join(os.path.dirname(__file__), 'test_tweak')
20
21    for run in [2,3]:
22
23        tanxy = getattr(tt, 'origxy_%i' % run)
24        xy = getattr(tt, 'xy_%i' % run)
25        noisyxy = getattr(tt, 'noisyxy_%i' % run)
26        gridx = getattr(tt, 'gridx_%i' % run)
27        gridy = getattr(tt, 'gridy_%i' % run)
28        truesip_a = getattr(tt, 'truesip_a_%i' % run)
29        truesip_b = getattr(tt, 'truesip_b_%i' % run)
30        sip_a = getattr(tt, 'sip_a_%i' % run)
31        sip_b = getattr(tt, 'sip_b_%i' % run)
32        x0,y0 = tt.x0, tt.y0
33
34        truedxy = xy - tanxy
35        obsdxy  = noisyxy - tanxy
36
37        xlo,xhi = -500, 2500
38        ylo,yhi = -500, 2500
39
40        X1 = linspace(xlo, xhi, 100)
41        Y1 = gridy
42        X1,Y1 = meshgrid(X1,Y1)
43        X1 = X1.T
44        Y1 = Y1.T
45        X2 = gridx
46        Y2 = linspace(ylo, yhi, 100)
47        X2,Y2 = meshgrid(X2,Y2)
48        truesipx_x = zeros_like(X1)
49        truesipy_x = zeros_like(X1)
50        truesipx_y = zeros_like(Y2)
51        truesipy_y = zeros_like(Y2)
52        for xo,yo,c in truesip_a:
53            truesipx_y += c * (X2 - x0)**xo * (Y2 - y0)**yo
54            truesipx_x += c * (X1 - x0)**xo * (Y1 - y0)**yo
55        for xo,yo,c in truesip_b:
56            truesipy_y += c * (X2 - x0)**xo * (Y2 - y0)**yo
57            truesipy_x += c * (X1 - x0)**xo * (Y1 - y0)**yo
58
59        x = xy[:,0]
60        y = xy[:,1]
61        truedx = truedxy[:,0]
62        truedy = truedxy[:,1]
63        obsdx = obsdxy[:,0]
64        obsdy = obsdxy[:,1]
65
66        for order in range(2,6):
67            clf()
68
69            sipx_x = zeros_like(X1)
70            sipy_x = zeros_like(X1)
71            sipx_y = zeros_like(Y2)
72            sipy_y = zeros_like(Y2)
73            for xo,yo,c in sip_a[order]:
74                sipx_y += c * (X2 - x0)**xo * (Y2 - y0)**yo
75                sipx_x += c * (X1 - x0)**xo * (Y1 - y0)**yo
76            for xo,yo,c in sip_b[order]:
77                sipy_y += c * (X2 - x0)**xo * (Y2 - y0)**yo
78                sipy_x += c * (X1 - x0)**xo * (Y1 - y0)**yo
79
80            subplot(2,2,1)
81            plot(x, truedx, 'bs', mec='b', mfc='None')
82            plot(x, obsdx, 'r.')
83            plot(X1, -truesipx_x, 'b-', alpha=0.2)
84            plot(X1, -sipx_x, 'r-', alpha=0.2)
85            xlabel('x')
86            ylabel('dx')
87            xlim(xlo, xhi)
88
89            subplot(2,2,2)
90            plot(x, truedy, 'bs', mec='b', mfc='None')
91            plot(x, obsdy, 'r.')
92            plot(X1, -truesipy_x, 'b-', alpha=0.2)
93            plot(X1, -sipy_x, 'r-', alpha=0.2)
94            xlabel('x')
95            ylabel('dy')
96            xlim(xlo, xhi)
97
98            subplot(2,2,3)
99            plot(y, truedx, 'bs', mec='b', mfc='None')
100            plot(y, obsdx, 'r.')
101            plot(Y2, -truesipx_y, 'b-', alpha=0.2)
102            plot(Y2, -sipx_y, 'r-', alpha=0.2)
103            xlabel('y')
104            ylabel('dx')
105            xlim(xlo, xhi)
106
107            subplot(2,2,4)
108            plot(y, truedy, 'bs', mec='b', mfc='None')
109            plot(y, obsdy, 'r.')
110            plot(Y2, -truesipy_y, 'b-', alpha=0.2)
111            plot(Y2, -sipy_y, 'r-', alpha=0.2)
112            xlabel('y')
113            ylabel('dy')
114            xlim(xlo, xhi)
115
116            savefig('tt%i-%i.png' % (run, order))
117
118            clf()
119            subplot(111)
120
121            plot(tanxy[:,0], tanxy[:,1], 'b.')
122            plot(noisyxy[:,0], noisyxy[:,1], 'r.')
123            plot(xy[:,0], xy[:,1], 'bo', mec='b', mfc='None')
124
125            if False:
126                X3,Y3 = meshgrid(linspace(xlo, xhi, 11),
127                                 linspace(ylo, yhi, 11))
128                truesipx = X3
129                truesipy = Y3
130                for xo,yo,c in truesip_a:
131                    truesipx += c * (X3 - x0)**xo * (Y3 - y0)**yo
132                for xo,yo,c in truesip_b:
133                    truesipy += c * (X3 - x0)**xo * (Y3 - y0)**yo
134                sipx = X3
135                sipy = Y3
136                for xo,yo,c in sip_a[order]:
137                    sipx += c * (X3 - x0)**xo * (Y3 - y0)**yo
138                for xo,yo,c in sip_b[order]:
139                    sipy += c * (X3 - x0)**xo * (Y3 - y0)**yo
140                plot(truesipx, truesipy, 'bs', mec='b', mfc='None')
141                plot(sipx, sipy, 'ms', mec='m', mfc='None')
142
143            plot(X1, Y1, 'g-', alpha=0.25)
144            plot(X2, Y2, 'g-', alpha=0.25)
145
146            plot(truesipx_x + X1, truesipy_x + Y1, 'b-', alpha=0.25)
147            plot(truesipx_y + X2, truesipy_y + Y2, 'b-', alpha=0.25)
148
149            plot(sipx_x + X1, sipy_x + Y1, 'r-', alpha=0.25)
150            plot(sipx_y + X2, sipy_y + Y2, 'r-', alpha=0.25)
151
152            xlim(xlo,xhi)
153            ylim(ylo,yhi)
154            savefig('ttxy%i-%i.png' % (run,order))
155