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