Skip to content

Commit 2259625

Browse files
author
bigobs
committedMar 28, 2018
Changes to scripts/*focus* and added .gitignore file.
1 parent 065b724 commit 2259625

File tree

5 files changed

+841
-0
lines changed

5 files changed

+841
-0
lines changed
 

‎.gitignore

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
gui/flask_rts2/build/
2+
mont4k_pipeline/build/
3+
**/*.swp
4+
*.swp
5+
junk
6+
**/*.pyc
7+
8+
9+
10+

‎scripts/focusing.py

+347
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,347 @@
1+
#!/usr/bin/python
2+
#
3+
# Autofocosing routines.
4+
#
5+
# You will need: scipy matplotlib sextractor
6+
# This should work on Debian/ubuntu:
7+
# sudo apt-get install python-matplotlib python-scipy python-pyfits sextractor
8+
#
9+
# If you would like to see sextractor results, get DS9 and pyds9:
10+
#
11+
# http://hea-www.harvard.edu/saord/ds9/
12+
#
13+
# Please be aware that current sextractor Ubuntu packages does not work
14+
# properly. The best workaround is to install package, and the overwrite
15+
# sextractor binary with one compiled from sources (so you will have access
16+
# to sextractor configuration files, which program assumes).
17+
#
18+
# (C) 2002-2008 Stanislav Vitek
19+
# (C) 2002-2010 Martin Jelinek
20+
# (C) 2009-2010 Markus Wildi
21+
# (C) 2010-2014 Petr Kubanek, Institute of Physics <kubanek@fzu.cz>
22+
# (C) 2010 Francisco Forster Buron, Universidad de Chile
23+
#
24+
# This program is free software; you can redistribute it and/or
25+
# modify it under the terms of the GNU General Public License
26+
# as published by the Free Software Foundation; either version 2
27+
# of the License, or (at your option) any later version.
28+
#
29+
# This program is distributed in the hope that it will be useful,
30+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
31+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32+
# GNU General Public License for more details.
33+
#
34+
# You should have received a copy of the GNU General Public License
35+
# along with this program; if not, write to the Free Software
36+
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
37+
38+
from rts2 import scriptcomm
39+
from rts2 import sextractor
40+
from scottSock import scottSock
41+
sepPresent = False
42+
try:
43+
import sep
44+
sepPresent = True
45+
except Exception as ex:
46+
pass
47+
48+
from pylab import *
49+
from scipy import *
50+
from scipy import optimize
51+
import numpy
52+
import pickle
53+
54+
LINEAR = 0
55+
"""Linear fit"""
56+
P2 = 1
57+
"""Fit using 2 power polynomial"""
58+
P4 = 2
59+
"""Fit using 4 power polynomial"""
60+
H3 = 3
61+
"""Fit using general Hyperbola (three free parameters)"""
62+
H2 = 4
63+
"""Fit using Hyperbola with fixed slope at infinity (two free parameters)"""
64+
65+
class Focusing (scriptcomm.Rts2Comm):
66+
"""Take and process focussing data."""
67+
68+
def __init__(self,exptime = 10,step=20,attempts=10,filterGalaxies=False):
69+
scriptcomm.Rts2Comm.__init__(self)
70+
self.log('I', 'This is a test')
71+
self.exptime = exptime
72+
self.step = step
73+
self.focuser = "F0"
74+
self.attempts = attempts
75+
76+
# if |offset| is above this value, try linear fit
77+
self.linear_fit = self.step * self.attempts / 2.0
78+
# target FWHM for linear fit
79+
self.linear_fit_fwhm = 3.5
80+
self.filterGalaxies = filterGalaxies
81+
82+
def doFit(self,fit):
83+
b = None
84+
errfunc = None
85+
fitfunc_r = None
86+
p0 = None
87+
88+
# try to fit..
89+
# this function is for flux..
90+
#fitfunc = lambda p, x: p[0] * p[4] / (p[4] + p[3] * (abs(x - p[1])) ** (p[2]))
91+
92+
# prepare fit based on its type..
93+
if fit == LINEAR:
94+
fitfunc = lambda p, x: p[0] + p[1] * x
95+
errfunc = lambda p, x, y: fitfunc(p, x) - y # LINEAR - distance to the target function
96+
p0 = [1, 1]
97+
fitfunc_r = lambda x, p0, p1: p0 + p1 * x
98+
elif fit == P2:
99+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2)
100+
errfunc = lambda p, x, y: fitfunc(p, x) - y # P2 - distance to the target function
101+
p0 = [1, 1, 1]
102+
fitfunc_r = lambda x, p0, p1, p2 : p0 + p1 * x + p2 * (x ** 2)
103+
elif fit == P4:
104+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2) + p[3] * (x ** 3) + p[4] * (x ** 4)
105+
errfunc = lambda p, x, y: fitfunc(p, x) - y # P4 - distance to the target function
106+
p0 = [1, 1, 1, 1, 1]
107+
fitfunc_r = lambda x, p0, p1: p0 + p1 * x + p2 * (x ** 2) + p3 * (x ** 3) + p4 * (x ** 4)
108+
elif fit == H3:
109+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + p[1] ** 2 * (x - p[2])**2)
110+
errfunc = lambda p, x, y: fitfunc(p, x) - y # H3 - distance to the target function
111+
p0 = [400., 3.46407715307, self.fwhm_MinimumX] # initial guess based on real data
112+
fitfunc_r = lambda x, p0, p1, p2 : sqrt(p0 ** 2 + p1 ** 2 * (x - p2) ** 2)
113+
elif fit == H2:
114+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + 3.46407715307 ** 2 * (x - p[1])**2) # 3.46 based on H3 fits
115+
errfunc = lambda p, x, y: fitfunc(p, x) - y # H2 - distance to the target function
116+
p0 = [400., self.fwhm_MinimumX] # initial guess based on real data
117+
fitfunc_r = lambda x, p0, p1 : sqrt(p0 ** 2 + 3.46407715307 ** 2 * (x - p1) ** 2)
118+
else:
119+
raise Exception('Unknow fit type {0}'.format(fit))
120+
121+
self.fwhm_poly, success = optimize.leastsq(errfunc, p0[:], args=(self.focpos, self.fwhm))
122+
123+
b = None
124+
125+
if fit == LINEAR:
126+
b = (self.linear_fit_fwhm - self.fwhm_poly[0]) / self.fwhm_poly[1]
127+
elif fit == H3:
128+
b = self.fwhm_poly[2]
129+
self.log('I', 'found minimum FWHM: {0}'.format(abs(self.fwhm_poly[0])))
130+
self.log('I', 'found slope at infinity: {0}'.format(abs(self.fwhm_poly[1])))
131+
elif fit == H2:
132+
b = self.fwhm_poly[1]
133+
self.log('I', 'found minimum FWHM: {0}'.format(abs(self.fwhm_poly[0])))
134+
else:
135+
b = optimize.fmin(fitfunc_r,self.fwhm_MinimumX,args=(self.fwhm_poly), disp=0)[0]
136+
self.log('I', 'found FWHM minimum at offset {0}'.format(b))
137+
return b
138+
139+
def tryFit(self,defaultFit):
140+
"""Try fit, change to linear fit if outside allowed range."""
141+
b = self.doFit(defaultFit)
142+
if (abs(b - numpy.average(self.focpos)) >= self.linear_fit):
143+
self.log('W','cannot do find best FWHM inside limits, trying H2 fit - best fit is {0}, average focuser position is {1}'.format(b, numpy.average(self.focpos)))
144+
b = self.doFit(H2)
145+
if (abs(b - numpy.average(self.focpos)) >= self.linear_fit):
146+
self.log('W','cannot do find best FWHM inside limits, trying linear fit - best fit is {0}, average focuser position is {1}'.format(b, numpy.average(self.focpos)))
147+
b = self.doFit(LINEAR)
148+
return b,LINEAR
149+
return b,H2
150+
return b,defaultFit
151+
152+
153+
def doFitOnArrays(self,fwhm,focpos,defaultFit):
154+
self.fwhm = array(fwhm)
155+
self.focpos = array(focpos)
156+
self.fwhm_MinimumX = 0
157+
min_fwhm=fwhm[0]
158+
for x in range(0,len(fwhm)):
159+
if fwhm[x] < min_fwhm:
160+
self.fwhm_MinimumX = x
161+
min_fwhm = fwhm[x]
162+
return self.tryFit(defaultFit)
163+
164+
def findBestFWHM(self,tries,defaultFit=P2,min_stars=95,ds9display=False,threshold=2.7,deblendmin=0.03):
165+
# X is FWHM, Y is offset value
166+
self.focpos=[]
167+
self.fwhm=[]
168+
fwhm_min = None
169+
self.fwhm_MinimumX = None
170+
keys = list(tries.keys())
171+
keys.sort()
172+
sextr = sextractor.Sextractor(threshold=threshold,deblendmin=deblendmin)
173+
for k in keys:
174+
try:
175+
sextr.runSExtractor(tries[k])
176+
fwhm,fwhms,nstars = sextr.calculate_FWHM(min_stars,self.filterGalaxies)
177+
except Exception as ex:
178+
self.log('W','offset {0}: {1}'.format(k,ex))
179+
continue
180+
self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,fwhm,nstars))
181+
focpos.append(k)
182+
fwhm.append(fwhm)
183+
if (fwhm_min is None or fwhm < fwhm_min):
184+
fwhm_MinimumX = k
185+
fwhm_min = fwhm
186+
return focpos,fwhm,fwhm_min,fwhm_MinimumX
187+
188+
def __sepFindFWHM(self,tries):
189+
from astropy.io import fits
190+
import math
191+
import traceback
192+
focpos=[]
193+
fwhm=[]
194+
fwhm_min=None
195+
fwhm_MinimumX=None
196+
keys = list(tries.keys())
197+
keys.sort()
198+
ln2=math.log(2)
199+
for k in keys:
200+
try:
201+
fwhms=[]
202+
ff=fits.open(tries[k])
203+
# loop on images..
204+
for i in range(1,len(ff)-1):
205+
data=ff[i].data
206+
bkg=sep.Background(numpy.array(data,numpy.float))
207+
sources=sep.extract(data-bkg, 5.0 * bkg.globalrms)
208+
self.log('I','bkg gobalrms {}'.format(bkg.globalrms))
209+
210+
for s in sources:
211+
fwhms.append(2 * math.sqrt(ln2 * (s[15]**2 + s[16]**2)))
212+
213+
214+
im_fwhm=numpy.median(fwhms)
215+
# find median from fwhms measurements..
216+
217+
self.log('I','median fwhm {}'.format(numpy.median(fwhms)))
218+
self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,im_fwhm,len(fwhms)))
219+
focpos.append(k)
220+
fwhm.append(im_fwhm)
221+
if (fwhm_min is None or im_fwhm < fwhm_min):
222+
fwhm_MinimumX = k
223+
fwhm_min = im_fwhm
224+
except Exception as ex:
225+
self.log('W','offset {0}: {1} {2}'.format(k,ex,traceback.format_exc()))
226+
227+
self.log('I','pickling')
228+
fd = open( "rts2.pkl", 'w' )
229+
pickle.dump(sources, fd)
230+
fd.close()
231+
return focpos,fwhm,fwhm_min,fwhm_MinimumX
232+
233+
234+
def findBestFWHM(self,tries,defaultFit=H3,min_stars=15,ds9display=False,threshold=2.7,deblendmin=0.03):
235+
# X is FWHM, Y is offset value
236+
self.focpos=[]
237+
self.fwhm=[]
238+
self.fwhm_min = None
239+
self.fwhm_MinimumX = None
240+
if sepPresent:
241+
self.focpos,self.fwhm,self.fwhm_min,self.fwhm_MinimumX = self.__sepFindFWHM(tries)
242+
else:
243+
self.focpos,self.fwhm,self.fwhm_min,self.fwhm_MinimumX = self.__sexFindFWHM(tries,threshold,deblendmin)
244+
self.focpos = array(self.focpos)
245+
self.fwhm = array(self.fwhm)
246+
247+
return self.tryFit(defaultFit)
248+
249+
def beforeReadout(self):
250+
self.current_focus = self.getValueFloat('FOC_POS',self.focuser)
251+
if (self.num == self.attempts):
252+
self.setValue('FOC_TOFF',0,self.focuser)
253+
else:
254+
self.off += self.step
255+
self.setValue('FOC_TOFF',self.off,self.focuser)
256+
257+
def takeImages(self):
258+
self.setValue('exposure',self.exptime)
259+
self.setValue('SHUTTER','LIGHT')
260+
self.off = -1 * self.step * (self.attempts / 2)
261+
self.setValue('FOC_TOFF',self.off,self.focuser)
262+
tries = {}
263+
# must be overwritten in beforeReadout
264+
self.current_focus = None
265+
266+
for self.num in range(1,self.attempts+1):
267+
self.log('I','starting {0}s exposure on offset {1}'.format(self.exptime,self.off))
268+
img = self.exposure(self.beforeReadout,'%b/foc_%N_{0}.fits'.format(self.num))
269+
270+
tries[self.current_focus] = img
271+
272+
self.log('I','all focusing exposures finished, processing data')
273+
274+
return self.findBestFWHM(tries)
275+
276+
def run(self):
277+
self.focuser = self.getValue('focuser')
278+
# send to some other coordinates if you wish so, or disable this for target for fixed coordinates
279+
#self.altaz (89,90)
280+
b,fit = self.takeImages()
281+
if fit == LINEAR:
282+
self.setValue('FOC_DEF',b,self.focuser)
283+
b,fit = self.takeImages()
284+
285+
self.setValue('FOC_DEF',b,self.focuser)
286+
287+
def plotFit(self,b,ftype):
288+
"""Plot fit graph."""
289+
fitfunc = None
290+
291+
if ftype == LINEAR:
292+
fitfunc = lambda p, x: p[0] + p[1] * x
293+
elif ftype == P2:
294+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2)
295+
elif ftype == P4:
296+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2) + p[3] * (x ** 3) + p[4] * (x ** 4)
297+
elif ftype == H3:
298+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + p[1] ** 2 * (x - p[2]) ** 2)
299+
elif ftype == H2:
300+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + 3.46407715307 ** 2 * (x - p[1]) ** 2) # 3.46 based on HYPERBOLA fits
301+
else:
302+
raise Exception('Unknow fit type {0}'.format(ftype))
303+
304+
x = linspace(self.focpos.min() - 1, self.focpos.max() + 1)
305+
306+
plot (self.focpos, self.fwhm, "r+", x, fitfunc(self.fwhm_poly, x), "r-")
307+
308+
show()
309+
310+
311+
def to_dataserver( fname, outfile='test.fits', clobber=True ):
312+
313+
fitsfd = fits.open( fname )
314+
315+
316+
width = 0
317+
height = 0
318+
for ext in fitsfd:
319+
if hasattr( ext, 'data' ):
320+
if ext.data is not None:
321+
width+=ext.data.shape[0]
322+
height+=ext.data.shape[1]
323+
324+
fitsfd.close()
325+
fsize = os.stat(fname).st_size
326+
327+
fd = open(fname, 'rb')
328+
329+
330+
if clobber:
331+
clobber_char = '!'
332+
else:
333+
clobber_char = ''
334+
meta = " {} {}{} 1 {} {} 0".format( fsize, clobber_char, '/home/bigobs/data/rts2'+outfile, width, height )
335+
meta = meta + (256-len(meta))*' '
336+
337+
data = meta+fd.read()
338+
lendata = len(data)
339+
soc = scottSock( '10.30.1.1', 6543 )
340+
341+
counter = 0
342+
socsize = 1024
343+
buffsize = 0
344+
while buffsize < len(data):
345+
sent = soc.send( data[buffsize:buffsize+1024] )
346+
buffsize+=sent
347+

‎scripts/focusing_test.py

+356
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,356 @@
1+
#!/usr/bin/python
2+
#
3+
# Autofocosing routines.
4+
#
5+
# You will need: scipy matplotlib sextractor
6+
# This should work on Debian/ubuntu:
7+
# sudo apt-get install python-matplotlib python-scipy python-pyfits sextractor
8+
#
9+
# If you would like to see sextractor results, get DS9 and pyds9:
10+
#
11+
# http://hea-www.harvard.edu/saord/ds9/
12+
#
13+
# Please be aware that current sextractor Ubuntu packages does not work
14+
# properly. The best workaround is to install package, and the overwrite
15+
# sextractor binary with one compiled from sources (so you will have access
16+
# to sextractor configuration files, which program assumes).
17+
#
18+
# (C) 2002-2008 Stanislav Vitek
19+
# (C) 2002-2010 Martin Jelinek
20+
# (C) 2009-2010 Markus Wildi
21+
# (C) 2010-2014 Petr Kubanek, Institute of Physics <kubanek@fzu.cz>
22+
# (C) 2010 Francisco Forster Buron, Universidad de Chile
23+
#
24+
# This program is free software; you can redistribute it and/or
25+
# modify it under the terms of the GNU General Public License
26+
# as published by the Free Software Foundation; either version 2
27+
# of the License, or (at your option) any later version.
28+
#
29+
# This program is distributed in the hope that it will be useful,
30+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
31+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32+
# GNU General Public License for more details.
33+
#
34+
# You should have received a copy of the GNU General Public License
35+
# along with this program; if not, write to the Free Software
36+
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
37+
38+
from rts2 import scriptcomm
39+
from rts2 import sextractor
40+
from scottSock import scottSock
41+
sepPresent = False
42+
try:
43+
import sep
44+
sepPresent = True
45+
except Exception as ex:
46+
pass
47+
48+
import os
49+
from pylab import *
50+
from scipy import *
51+
from scipy import optimize
52+
import numpy
53+
import pickle
54+
55+
LINEAR = 0
56+
"""Linear fit"""
57+
P2 = 1
58+
"""Fit using 2 power polynomial"""
59+
P4 = 2
60+
"""Fit using 4 power polynomial"""
61+
H3 = 3
62+
"""Fit using general Hyperbola (three free parameters)"""
63+
H2 = 4
64+
"""Fit using Hyperbola with fixed slope at infinity (two free parameters)"""
65+
66+
class Focusing (scriptcomm.Rts2Comm):
67+
"""Take and process focussing data."""
68+
def getTriestmp(self, pathway="/home/rts2obs/rts2images"):
69+
focusfiles = [x for x in os.listdir(pathway) if "foc" in x]
70+
tries = {}
71+
for f in focusfiles:
72+
num = f.split('_')[2].split(".")[0]
73+
tries[float(num)] = f
74+
75+
return tries
76+
77+
def __init__(self,exptime = 10,step=20,attempts=10,filterGalaxies=False):
78+
scriptcomm.Rts2Comm.__init__(self)
79+
self.log('I', 'This is a test')
80+
self.exptime = exptime
81+
self.step = step
82+
self.focuser = "F0"
83+
self.attempts = attempts
84+
85+
# if |offset| is above this value, try linear fit
86+
self.linear_fit = self.step * self.attempts / 2.0
87+
# target FWHM for linear fit
88+
self.linear_fit_fwhm = 3.5
89+
self.filterGalaxies = filterGalaxies
90+
91+
def doFit(self,fit):
92+
b = None
93+
errfunc = None
94+
fitfunc_r = None
95+
p0 = None
96+
97+
# try to fit..
98+
# this function is for flux..
99+
#fitfunc = lambda p, x: p[0] * p[4] / (p[4] + p[3] * (abs(x - p[1])) ** (p[2]))
100+
101+
# prepare fit based on its type..
102+
if fit == LINEAR:
103+
fitfunc = lambda p, x: p[0] + p[1] * x
104+
errfunc = lambda p, x, y: fitfunc(p, x) - y # LINEAR - distance to the target function
105+
p0 = [1, 1]
106+
fitfunc_r = lambda x, p0, p1: p0 + p1 * x
107+
elif fit == P2:
108+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2)
109+
errfunc = lambda p, x, y: fitfunc(p, x) - y # P2 - distance to the target function
110+
p0 = [1, 1, 1]
111+
fitfunc_r = lambda x, p0, p1, p2 : p0 + p1 * x + p2 * (x ** 2)
112+
elif fit == P4:
113+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2) + p[3] * (x ** 3) + p[4] * (x ** 4)
114+
errfunc = lambda p, x, y: fitfunc(p, x) - y # P4 - distance to the target function
115+
p0 = [1, 1, 1, 1, 1]
116+
fitfunc_r = lambda x, p0, p1: p0 + p1 * x + p2 * (x ** 2) + p3 * (x ** 3) + p4 * (x ** 4)
117+
elif fit == H3:
118+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + p[1] ** 2 * (x - p[2])**2)
119+
errfunc = lambda p, x, y: fitfunc(p, x) - y # H3 - distance to the target function
120+
p0 = [400., 3.46407715307, self.fwhm_MinimumX] # initial guess based on real data
121+
fitfunc_r = lambda x, p0, p1, p2 : sqrt(p0 ** 2 + p1 ** 2 * (x - p2) ** 2)
122+
elif fit == H2:
123+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + 3.46407715307 ** 2 * (x - p[1])**2) # 3.46 based on H3 fits
124+
errfunc = lambda p, x, y: fitfunc(p, x) - y # H2 - distance to the target function
125+
p0 = [400., self.fwhm_MinimumX] # initial guess based on real data
126+
fitfunc_r = lambda x, p0, p1 : sqrt(p0 ** 2 + 3.46407715307 ** 2 * (x - p1) ** 2)
127+
else:
128+
raise Exception('Unknow fit type {0}'.format(fit))
129+
130+
self.fwhm_poly, success = optimize.leastsq(errfunc, p0[:], args=(self.focpos, self.fwhm))
131+
132+
b = None
133+
134+
if fit == LINEAR:
135+
b = (self.linear_fit_fwhm - self.fwhm_poly[0]) / self.fwhm_poly[1]
136+
elif fit == H3:
137+
b = self.fwhm_poly[2]
138+
self.log('I', 'found minimum FWHM: {0}'.format(abs(self.fwhm_poly[0])))
139+
self.log('I', 'found slope at infinity: {0}'.format(abs(self.fwhm_poly[1])))
140+
elif fit == H2:
141+
b = self.fwhm_poly[1]
142+
self.log('I', 'found minimum FWHM: {0}'.format(abs(self.fwhm_poly[0])))
143+
else:
144+
b = optimize.fmin(fitfunc_r,self.fwhm_MinimumX,args=(self.fwhm_poly), disp=0)[0]
145+
self.log('I', 'found FWHM minimum at offset {0}'.format(b))
146+
return b
147+
148+
def tryFit(self,defaultFit):
149+
"""Try fit, change to linear fit if outside allowed range."""
150+
b = self.doFit(defaultFit)
151+
if (abs(b - numpy.average(self.focpos)) >= self.linear_fit):
152+
self.log('W','cannot do find best FWHM inside limits, trying H2 fit - best fit is {0}, average focuser position is {1}'.format(b, numpy.average(self.focpos)))
153+
b = self.doFit(H2)
154+
if (abs(b - numpy.average(self.focpos)) >= self.linear_fit):
155+
self.log('W','cannot do find best FWHM inside limits, trying linear fit - best fit is {0}, average focuser position is {1}'.format(b, numpy.average(self.focpos)))
156+
b = self.doFit(LINEAR)
157+
return b,LINEAR
158+
return b,H2
159+
return b,defaultFit
160+
161+
162+
def doFitOnArrays(self,fwhm,focpos,defaultFit):
163+
self.fwhm = array(fwhm)
164+
self.focpos = array(focpos)
165+
self.fwhm_MinimumX = 0
166+
min_fwhm=fwhm[0]
167+
for x in range(0,len(fwhm)):
168+
if fwhm[x] < min_fwhm:
169+
self.fwhm_MinimumX = x
170+
min_fwhm = fwhm[x]
171+
return self.tryFit(defaultFit)
172+
173+
def findBestFWHM(self,tries,defaultFit=P2,min_stars=95,ds9display=False,threshold=2.7,deblendmin=0.03):
174+
# X is FWHM, Y is offset value
175+
self.focpos=[]
176+
self.fwhm=[]
177+
fwhm_min = None
178+
self.fwhm_MinimumX = None
179+
keys = list(tries.keys())
180+
keys.sort()
181+
sextr = sextractor.Sextractor(threshold=threshold,deblendmin=deblendmin)
182+
for k in keys:
183+
try:
184+
sextr.runSExtractor(tries[k])
185+
fwhm,fwhms,nstars = sextr.calculate_FWHM(min_stars,self.filterGalaxies)
186+
except Exception as ex:
187+
self.log('W','offset {0}: {1}'.format(k,ex))
188+
continue
189+
self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,fwhm,nstars))
190+
focpos.append(k)
191+
fwhm.append(fwhm)
192+
if (fwhm_min is None or fwhm < fwhm_min):
193+
fwhm_MinimumX = k
194+
fwhm_min = fwhm
195+
return focpos,fwhm,fwhm_min,fwhm_MinimumX
196+
197+
def __sepFindFWHM(self,tries):
198+
from astropy.io import fits
199+
import math
200+
import traceback
201+
focpos=[]
202+
fwhm=[]
203+
fwhm_min=None
204+
fwhm_MinimumX=None
205+
keys = list(tries.keys())
206+
keys.sort()
207+
ln2=math.log(2)
208+
for k in keys:
209+
try:
210+
fwhms=[]
211+
ff=fits.open(tries[k])
212+
# loop on images..
213+
for i in range(1,len(ff)-1):
214+
data=ff[i].data
215+
bkg=sep.Background(numpy.array(data,numpy.float))
216+
sources=sep.extract(data-bkg, 5.0 * bkg.globalrms)
217+
self.log('I','bkg gobalrms {}'.format(bkg.globalrms))
218+
219+
for s in sources:
220+
fwhms.append(2 * math.sqrt(ln2 * (s[15]**2 + s[16]**2)))
221+
222+
223+
im_fwhm=numpy.median(fwhms)
224+
# find median from fwhms measurements..
225+
226+
self.log('I','median fwhm {}'.format(numpy.median(fwhms)))
227+
self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,im_fwhm,len(fwhms)))
228+
focpos.append(k)
229+
fwhm.append(im_fwhm)
230+
if (fwhm_min is None or im_fwhm < fwhm_min):
231+
fwhm_MinimumX = k
232+
fwhm_min = im_fwhm
233+
except Exception as ex:
234+
self.log('W','offset {0}: {1} {2}'.format(k,ex,traceback.format_exc()))
235+
236+
self.log('I','pickling')
237+
fd = open( "rts2.pkl", 'w' )
238+
pickle.dump(sources, fd)
239+
fd.close()
240+
return focpos,fwhm,fwhm_min,fwhm_MinimumX
241+
242+
243+
def findBestFWHM(self,tries,defaultFit=H3,min_stars=15,ds9display=False,threshold=2.7,deblendmin=0.03):
244+
# X is FWHM, Y is offset value
245+
self.focpos=[]
246+
self.fwhm=[]
247+
self.fwhm_min = None
248+
self.fwhm_MinimumX = None
249+
if sepPresent:
250+
self.focpos,self.fwhm,self.fwhm_min,self.fwhm_MinimumX = self.__sepFindFWHM(tries)
251+
else:
252+
self.focpos,self.fwhm,self.fwhm_min,self.fwhm_MinimumX = self.__sexFindFWHM(tries,threshold,deblendmin)
253+
self.focpos = array(self.focpos)
254+
self.fwhm = array(self.fwhm)
255+
256+
return self.tryFit(defaultFit)
257+
258+
def beforeReadout(self):
259+
self.current_focus = self.getValueFloat('FOC_POS',self.focuser)
260+
if (self.num == self.attempts):
261+
self.setValue('FOC_TOFF',0,self.focuser)
262+
else:
263+
self.off += self.step
264+
self.setValue('FOC_TOFF',self.off,self.focuser)
265+
266+
def takeImages(self):
267+
self.setValue('exposure',self.exptime)
268+
self.setValue('SHUTTER','LIGHT')
269+
self.off = -1 * self.step * (self.attempts / 2)
270+
self.setValue('FOC_TOFF',self.off,self.focuser)
271+
tries = {}
272+
# must be overwritten in beforeReadout
273+
self.current_focus = None
274+
275+
for self.num in range(1,self.attempts+1):
276+
self.log('I','starting {0}s exposure on offset {1}'.format(self.exptime,self.off))
277+
img = self.exposure(self.beforeReadout,'%b/foc_%N_{0}.fits'.format(self.num))
278+
279+
tries[self.current_focus] = img
280+
281+
self.log('I','all focusing exposures finished, processing data')
282+
283+
return self.findBestFWHM(tries)
284+
285+
def run(self):
286+
self.focuser = self.getValue('focuser')
287+
# send to some other coordinates if you wish so, or disable this for target for fixed coordinates
288+
#self.altaz (89,90)
289+
b,fit = self.takeImages()
290+
if fit == LINEAR:
291+
self.setValue('FOC_DEF',b,self.focuser)
292+
b,fit = self.takeImages()
293+
294+
self.setValue('FOC_DEF',b,self.focuser)
295+
296+
def plotFit(self,b,ftype):
297+
"""Plot fit graph."""
298+
fitfunc = None
299+
300+
if ftype == LINEAR:
301+
fitfunc = lambda p, x: p[0] + p[1] * x
302+
elif ftype == P2:
303+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2)
304+
elif ftype == P4:
305+
fitfunc = lambda p, x: p[0] + p[1] * x + p[2] * (x ** 2) + p[3] * (x ** 3) + p[4] * (x ** 4)
306+
elif ftype == H3:
307+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + p[1] ** 2 * (x - p[2]) ** 2)
308+
elif ftype == H2:
309+
fitfunc = lambda p, x: sqrt(p[0] ** 2 + 3.46407715307 ** 2 * (x - p[1]) ** 2) # 3.46 based on HYPERBOLA fits
310+
else:
311+
raise Exception('Unknow fit type {0}'.format(ftype))
312+
313+
x = linspace(self.focpos.min() - 1, self.focpos.max() + 1)
314+
315+
plot (self.focpos, self.fwhm, "r+", x, fitfunc(self.fwhm_poly, x), "r-")
316+
317+
show()
318+
319+
320+
def to_dataserver( fname, outfile='test.fits', clobber=True ):
321+
322+
fitsfd = fits.open( fname )
323+
324+
325+
width = 0
326+
height = 0
327+
for ext in fitsfd:
328+
if hasattr( ext, 'data' ):
329+
if ext.data is not None:
330+
width+=ext.data.shape[0]
331+
height+=ext.data.shape[1]
332+
333+
fitsfd.close()
334+
fsize = os.stat(fname).st_size
335+
336+
fd = open(fname, 'rb')
337+
338+
339+
if clobber:
340+
clobber_char = '!'
341+
else:
342+
clobber_char = ''
343+
meta = " {} {}{} 1 {} {} 0".format( fsize, clobber_char, '/home/bigobs/data/rts2'+outfile, width, height )
344+
meta = meta + (256-len(meta))*' '
345+
346+
data = meta+fd.read()
347+
lendata = len(data)
348+
soc = scottSock( '10.30.1.1', 6543 )
349+
350+
counter = 0
351+
socsize = 1024
352+
buffsize = 0
353+
while buffsize < len(data):
354+
sent = soc.send( data[buffsize:buffsize+1024] )
355+
buffsize+=sent
356+

‎scripts/rts2-focusing

+43
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
#!/usr/bin/python
2+
#
3+
# Autofocosing routines.
4+
#
5+
# You will need: scipy matplotlib sextractor
6+
# This should work on Debian/Ubuntu:
7+
# sudo apt-get install python-matplotlib python-scipy python-pyfits sextractor
8+
#
9+
# If you would like to see sextractor results, get DS9 and pyds9:
10+
#
11+
# http://hea-www.harvard.edu/saord/ds9/
12+
#
13+
# Please be aware that current sextractor Ubuntu packages does not work
14+
# properly. The best workaround is to install package, and the overwrite
15+
# sextractor binary with one compiled from sources (so you will have access
16+
# to sextractor configuration files, which program assumes).
17+
#
18+
# Alternatively, use SEP package (pip install sep; http://sep.readthedocs.org).
19+
#
20+
# (C) 2002-2008 Stanislav Vitek
21+
# (C) 2002-2010 Martin Jelinek
22+
# (C) 2009-2010 Markus Wildi
23+
# (C) 2010-2016 Petr Kubanek, Institute of Physics <kubanek@fzu.cz>
24+
# (C) 2010 Francisco Forster Buron, Universidad de Chile
25+
#
26+
# This program is free software; you can redistribute it and/or
27+
# modify it under the terms of the GNU General Public License
28+
# as published by the Free Software Foundation; either version 2
29+
# of the License, or (at your option) any later version.
30+
#
31+
# This program is distributed in the hope that it will be useful,
32+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
33+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
34+
# GNU General Public License for more details.
35+
#
36+
# You should have received a copy of the GNU General Public License
37+
# along with this program; if not, write to the Free Software
38+
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
39+
40+
import focusing
41+
42+
a = focusing.Focusing()
43+
a.run()

‎scripts/test_focus.py

+85
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
from rts2 import scriptcomm
2+
from rts2 import sextractor
3+
from scottSock import scottSock
4+
sepPresent = False
5+
try:
6+
import sep
7+
sepPresent = True
8+
except Exception as ex:
9+
pass
10+
11+
import os
12+
from pylab import *
13+
from scipy import *
14+
from scipy import optimize
15+
import numpy
16+
import pickle
17+
18+
19+
20+
def getTries(pathway="/home/rts2obs/rts2images"):
21+
focusfiles = [x for x in os.listdir(pathway) if "foc_20171204" in x and ".fits" in x]
22+
tries = {}
23+
for f in focusfiles:
24+
#print f
25+
num = f.split('_')[2].split(".")[0]
26+
tries[float(num)] = pathway+"/"+f
27+
28+
return tries
29+
30+
def __sepFindFWHM(tries):
31+
from astropy.io import fits
32+
import math
33+
import traceback
34+
focpos=[]
35+
fwhm=[]
36+
fwhm_min=None
37+
fwhm_MinimumX=None
38+
keys = list(tries.keys())
39+
keys.sort()
40+
ln2=math.log(2)
41+
for k in keys:
42+
try:
43+
fwhms=[]
44+
ff=fits.open(tries[k])
45+
# loop on images..
46+
for i in range(1,len(ff)-1):
47+
data=ff[i].data
48+
bkg=sep.Background(numpy.array(data,numpy.float))
49+
sources=sep.extract(data-bkg, 5.0 * bkg.globalrms)
50+
#self.log
51+
print('I','bkg gobalrms {}'.format(bkg.globalrms))
52+
53+
for s in sources:
54+
fwhms.append(2 * math.sqrt(ln2 * (s[15]**2 + s[16]**2)))
55+
56+
57+
im_fwhm=numpy.median(fwhms)
58+
# find median from fwhms measurements..
59+
60+
#self.log
61+
print('I','median fwhm {}'.format(numpy.median(fwhms)))
62+
#self.log
63+
print('I','offset {0} fwhm {1} with {2} stars'.format(k,im_fwhm,len(fwhms)))
64+
focpos.append(k)
65+
fwhm.append(im_fwhm)
66+
if (fwhm_min is None or im_fwhm < fwhm_min):
67+
fwhm_MinimumX = k
68+
fwhm_min = im_fwhm
69+
except Exception as ex:
70+
#self.log
71+
print('W','offset {0}: {1} {2}'.format(k,ex,traceback.format_exc()))
72+
73+
#self.log('I','pickling')
74+
#fd = open( "rts2.pkl", 'w' )
75+
#pickle.dump(sources, fd)
76+
#fd.close()
77+
return focpos,fwhm,fwhm_min,fwhm_MinimumX
78+
79+
80+
def main():
81+
tries = getTries()
82+
focpos, fwhm, fwhm_min, fwhm_MinimumX = __sepFindFWHM(tries)
83+
print focpos, fwhm, fwhm_min, fwhm_MinimumX
84+
85+
main()

0 commit comments

Comments
 (0)
Please sign in to comment.