Astronomska skupina Fakultete za matematiko in fizilo Univerze v Ljubljani
Teleskop Vega

Spremenljive zvezde v razsuti kopici M67 (aperturna fotometrija)

Delovna mapa s podatki je http://astro.ago.fmf.uni-lj.si/podatki/2005/V2005-04-04/

Obdelava meritev

1.korak

vse slike obdelamo za 'zero' 'dark' in 'flat'. Glej navodila.

2.korak

vse slike opremimo z ekvatorialnimi 'wcs' koordinatami. Vzamem eno izmed slik, recimo serija1-20.fts, s programom fitsblink v 'edit parameters' vnesem ročno koordinate kopice, vnesem velikost piksla, ki je 0.914", potem naredim 'detct stars' in 'match'. Dam še 'write to file' in grem iz programa. Nato v irafu wcs koordinate iz slike serija1-20.fts, ki jo vzamem za referenčno, prekopiram v vse ostale slike:

 ecl> wcscopy *.fts serija1-20

in si prihranim goro tipkanja in klikanja, ker so parametri za vse slike skoraj enaki, nato v konzoli uporabim se skriptni programček, ki osveži konstante:

 $ writewcs.pl *.fts

na vsaki sliki znova poračuna konstante polje in natančne wcs koordinate vpiše v zaglavje slike za nadaljno procesiranje. Ce pri kaksni sliki ne uspe najti astrometricnega fita, lahko sicer poizkusim se rocno v fitsblinku, tako da popravljam nastavitev zacetnih parametrov za astrometricni fit. Ce ne najdem matcha je slika prakticno neuporabna za nadaljno obdelavo.

Potem se pocistim mapo slabih slik, pozenem program

 $ checkbad.sh

vse 'slabe' slike, to je tiste kjer ni mogoce najti astrometricnega fita pomete v podmapo bad.

3.korak

sestavimo spisek spremenljivih zvezd, ki jih želimo izmeriti in dodamo se 5-7 primerialnih zvezd. Datoteko imenujem vhodna.wcs

  132.907867   11.849253 AH Cnc 13.31 13.80 EB*WUMa?           
  132.861337   11.936761 ET Cnc 15.96  0    EB*WUMa?        
  132.867375   11.824378 EV Cnc 12.89 13.35 EB*WUMa?       
  132.9201     11.81817  00214  13.85 0.60
  132.9218     11.908008 00215  13.763 0.592
  132.8937     11.730394 02029  14.901 0.733
  132.8776     11.796464 02003  15.77 0.99

Prva in druga kolona je rektascenzija in deklinacija v stopinjah. Ostalo, ki se sicer ne upošteva pri računu je ime zvezde, magnituda in tip spremenljivke, kot za komentar.

Datoteko napišem ročno ali pa si pomagam recimo s Simbadom tako da naredim overlay na našo sliko in potem poklikam izbrane zvede. Uporabi program skycat ali gaia ali jskycat. Če to delam v irafu, je ukaz:

 imexamine wcs=world

tako dobim rektascenzijo in deklinacijo, ce pa imam samo xy koordinate in uporabim še program wcsctran za transformacijo v wcs koordinatni sistem, ki se edini ohranja na vseh slikah. Primerjalne zvezde izberemo iz UBV kataloga kopice M67, moramo pa seveda preveriti če je izbrana zvezda na sliki in ali ni presvetljena (saturirana). Za to je zelo primerna je procedura 'imexamin' v osnovnem paketu iraf.

3a.korak

Preverimo koordinate zvezd tako da narišemo na sliko:

 ecl> epar wcsctran
 PACKAGE = imcoords
   TASK = wcsctran

 input   =           vhodna.wcs  The input coordinate files
 output  =            vhodna.xy  The output coordinate files
 image   =           serija1-15  The input images
 inwcs   =                world  The input coordinate system
 outwcs  =              logical  The output coordinate system
 (columns=        1 2 3 4 5 6 7) List of input file columns
 (units  =                 hour) List of input coordinate units
 (formats=                     ) List of output coordinate formats
 (min_sig=                    7) Minimum precision of output coordinates
 (verbose=                  yes) Write comments to the output file ?
 (mode   =                   ql)

 ecl> display serija1-15 1
 ecl> wcsctran 
 ecl> tvmark 1 vhodna.xy label=yes

in dobimo

kliki na sliko za povečavo

4.korak

Pripravimo spisek slik za obdelavo:

 ecl> files serija1*.fts > seznam_slik

5.korak

V irafu v podpaketu noao.digiphot.apphot nastavi parametre za fotometrijo:

 appgot> photpars
                            Image Reduction and Analysis Facility
 PACKAGE = apphot
   TASK = photpars

 (weighti=             constant) Photometric weighting scheme for wphot
 (apertur=                    2) List of aperture radii in scale units
 (zmag   =                  25.) Zero point of magnitude scale
 (mkapert=                   no) Draw apertures on the display
 (mode   =                   ql)
 ($nargs =                    0)

Nastavim velikost zaslonke, na dva 2 piksla. Sicer priporocajo vzeti 0.8*FWHM, ker naj bi bi tam nekje maksimum S/N. Nato se nastavim parametre za phot:

 apphot> epar phot

 PACKAGE = apphot
   TASK = phot

 image   =         @seznam_slik  The input image(s)
 skyfile =                       The input sky file(s)
 (coords =           vhodna.wcs) The input coordinate files(s) (default: image.coo.?)
 (output =              default) The output photometry file(s) (default: image.mag.?)
 (plotfil=                     ) The output plots metacode file
 (datapar=                     ) Data dependent parameters
 (centerp=                     ) Centering parameters
 (fitskyp=                     ) Sky fitting parameters
 (photpar=                     ) Photometry parameters
 (interac=                   no) Interactive mode ?
 (radplot=                   no) Plot the radial profiles in interactive mode ?
 (icomman=                     ) Image cursor: [x y wcs] key [cmd]
 (gcomman=                     ) Graphics cursor: [x y wcs] key [cmd]
 (wcsin  =                world) The input coordinate system  (logical,tv,physical,world)
 (wcsout =              logical) The output coordinate system (logical,tv,physical)
 (cache  =             )_.cache) Cache the input image pixels in memory ?
 (verify =            )_.verify) Verify critical parameters in non-interactive mode ?
 (update =            )_.update) Update critical parameters in non-interactive mode ?
 (verbose=           )_.verbose) Print messages in non-interactive mode ?
 (graphic=          )_.graphics) Graphics device
 (display=           )_.display) Display device
 (mode   =                   ql)

6.korak

pozenemo fotometrijo na celi seriji slik, medtem gremo lahko na malco:

 apphot> phot

7.korak

iz dobljenih datotek mag.1 izločimo samo nekatere podatke

 apphot> epar pdump

 PACKAGE = apphot
   TASK = pdump

 infiles = @seznam_slik//.mag.1  Input apphot/daophot databases(s)
 fields  = ID,MAG,MERR,FLUX,IMAG  Fields to be extracted
 expr    =     PERROR="NoError?"  Boolean expression
 (headers=                   no) Print field headers?
 (paramet=                  yes) Print parameters?
 (inlist =                     )
 (mode   =                   ql)

in poženem v datoteko:

 apphot> pdump @seznam_slik//.mag.1 > rezultat.dat

ki izgleda takole

 1  13.291  0.005  48268.11  serija1-0
 2  15.347  0.015  7262.556  serija1-0
 3  12.337  0.003  116217.5  serija1-0
 4  13.429  0.005  42515.93  serija1-0
 5  12.525  0.004  97755.13  serija1-0
 6  13.525  0.006  38890.88  serija1-0
 7  14.164  0.008  21588.11  serija1-0
 8  15.168  0.014  8569.697  serija1-0
 1  13.046  0.004  60469.31  serija1-1
 2  15.146  0.013  8740.249  serija1-1
 3  12.194  0.003  132534.1  serija1-1
 4  13.183  0.005  53301.91  serija1-1
 5  12.286  0.003  121841.9  serija1-1
 6  13.932  0.007  26733.48  serija1-1
 7  14.985  0.012  10142.86  serija1-1
 1  13.039  0.004  60853.91  serija1-2
 2  15.199  0.013  8324.11  serija1-2
 se nadaljuje...

8.korak

Podatek o času posameznega posnetka, je najbolje potegniti kar iz glave slike:

 ecl> imhead @seznam_slik lo+ | grep TIME-OBS > glave_slik.txt

in datoteko še malo poeditiram tako da dobim:

 serija1-0   19:19:25
 serija1-1   19:34:09
 serija1-2   19:36:14
 serija1-3   19:38:20
 serija1-4   19:40:27
 serija1-5   19:42:33
 serija1-6   19:44:39
 ... 

Še bolje bi bilo naštimat, da bi iraf phot pobral recimo JD ze kar iz glave, drugič.

9.korak

Datoteko rezultat moramo se preurediti, zato smo napisali programček v pythonu:

#! /usr/bin/env python2
"""
    Program prebere in preuredi rezultate fotometrije
<:vspace>
"""
import os
import glob
import string
import re
from math import *
<:vspace>
def str_to_dec(x):
        """ Funkcija pretvoti zapis ure v decimalno stevilo
        """
	z=1.0
	z1,z2,z3='','',''
	p=re.match(r'(?P<z1>[+-])?(?P<h>\d{1,2})[: ](?P<z2>[+-])?(?P<m>\d{1,2})[: ](?P<z3>[+-])?(?P<s>\d{1,2}(\.\d*)?)',x)
	if p:
		if p.group('z1')=='-' or p.group('z2')=='-' or p.group('z3')=='-' : z=-1.0
		return z*(float(p.group('h'))+float(p.group('m'))/60+float(p.group('s'))/3600)
<:vspace>
<:vspace>
""" 
    Glavni program 
"""
<:vspace>
# prebere datoteko 
d = open('rezultat.dat','r')
line = d.readline()
<:vspace>
No =[]
Mag=[]
Merr=[]
Flux=[]
Image=[]
<:vspace>
while line:
	rec = string.split(string.strip(line))
	if len(rec)==5:
		n, mag,merr,flux,image = rec
         	No.append(int(n))
		Mag.append(float(mag))
		Merr.append(float(merr))
		Flux.append(float(flux))
		Image.append(image)
	line = d.readline()
d.close()
<:vspace>
# prebere informacije o casih posnetka 
Slika=[]
Cas=[]
d = open('seznam_slik_casi.txt','r')
line = d.readline()
while line:
	rec = string.split(string.strip(line))
	slika,cas = rec
	ncas = str_to_dec(cas)
	Slika.append(slika)
	Cas.append(ncas)
	line = d.readline()
d.close()
<:vspace>
<:vspace>
""" 
   analizira podatke
"""
j=0
k=0
for i in No:
	if i==1:
		"""
		print  'primerjalni flux=',Flux[j+3] + Flux[j+4] + Flux[j+5] + Flux[j+6] 
		print  'AHCnc            ',Flux[j+0]
		"""
		mag1 = -2.5 * log10( Flux[j+0]/(Flux[j+3] + Flux[j+4] + Flux[j+5] + Flux[j+6]))
		mag2 = -2.5 * log10( Flux[j+1]/(Flux[j+3] + Flux[j+4] + Flux[j+5] + Flux[j+6]))
		mag3 = -2.5 * log10( Flux[j+2]/(Flux[j+3] + Flux[j+4] + Flux[j+5] + Flux[j+6]))
		print Cas[k],mag1, mag2, mag3
		k = k+1
	j=j+1

Poženem:

 $./analiza.py > obdelano.dat

Datoteka obdelano.dat izgleda takole:

 19.3236111111 3.56322376935 8.29832192225 1.36649291601
 19.5691666667 3.1363633843 7.9718565959 1.1746033663
 19.6038888889 3.1249447242 8.09824460043 1.16537026057
 19.6388888889 3.15129794826 8.20088659039 1.1276509186
 19.6741666667 3.1386160781 8.11418128338 1.28333700942
 19.7091666667 3.14663139041 8.03406194278 1.23571726145
 19.7441666667 3.14174002525 8.20221146788 1.22542396475
 19.7791666667 3.12857319041 7.96959396419 1.25071514042
 se nadaljuje ....

Prvi stolpec je čas UT v urah, drugi in naprej je delta magnitude za AH Cnc, ET Cnc, EV Cnc .

10. korak

Narišem časovni potek sija spremenljive zvezde:

 $ gnuplot

 gnuplot> plot 'obdelano.dat' using 1:2

in dobimo:

bd

Zadnja sprememba 13.06.2013 14:25 CEST
© Copyright 2004-2017
B.Dintinjana
Univerza v Ljubljani, F M F
Pot na Golovec 25
1000 Ljubljana, Slovenija