home |
|
Sugestão de scripts para redução
|
|
|
|
Instrumentos e Detetores
CAMIV
Informações p/ o usuário |
|
Estes são os scripts usados pelo Francisco Jablonski na redução dos dados
da CamIV. Para facilitar, ele criou um package (chama-se camiv) que deve ser
declarado no começo do login.cl, conforme o modelo abaixo. Edite de tal
forma que a sua estrutura de diretório substitua a dele no exemplo.
Então, modifique seu login.cl para que o caminho para o package camiv seja
conhecido:
------------------- primeiras linhas do login.cl ----------------------
set home = "/export/usuarios/chico/iraf/"
set imdir = "HDR$"
set uparm = "home$uparm/"
set userid = "chico"
task $camiv = "/export/usuarios/chico/iraf/NIR/camiv.cl"
-------------------------continua--------------------------------------
Não esqueça de editar também o imtype no seu login.cl. A linha onde ele
aparece deve ser
set imtype = "fits"
(Se vc ainda usa o formato .imh/.pix para alguma coisa, precisa editar outra
vez o login.cl antes de usar esse tipo de imagem. Eu acho vantajoso usar
direto o formato fits para tudo.)
O arquivo camiv.cl (listagem abaixo) contém os nomes dos tasks criados.
Veja que o diretório /export/usuarios/chico/iraf/NIR contém tanto o
file camiv.cl (que declara o package) como os files correspondentes aos
diversos tasks.
-------------------------file camiv.cl----------------------------------
procedure camiv()
begin
package camiv
task lineariza = "/export/usuarios/chico/iraf/NIR/lineariza.cl"
task badfaz = "/export/usuarios/chico/iraf/NIR/badfaz.cl"
task prepara = "/export/usuarios/chico/iraf/NIR/prepara.cl"
task reduz = "/export/usuarios/chico/iraf/NIR/reduz.cl"
clbye()
end
------------------------- fim -----------------------------------------
A listagem dos quatro scripts esta no final deste texto. Corte e crie
um arquivo texto para cada script, e coloque-os com os nomes corretos
no diretório que escolher para substituir o do exemplo acima.
Com as definições desse jeito, depois de entrar no cl você vera um
package a mais (camiv) entre os outros. Teclando "camiv ?" aparecem
os task definidos aí acima. Note que *não* foram escritos help pages
para esses tasks.
Uma rápida descrição dos tasks e seu funcionamento:
1) lineariza
Corrige as imagens da Camiv dos efeitos de não-linearidade do detetor.
O modelo utilizado na correção é o mais simples possível e pressupõe
que todos os pixels se comportam igualmente:
ADU_corrigido = ADU_bruto + coef * (ADU_bruto)**2
O valor default para coef = 3.6e-6
O F. Jablonski determinou esse valor a partir de muitas noites distintas
de uso da CamIV. Você pode ter uma idéia da qualidade dessa determinação
olhando o gráfico que está em http://old.lna.br/instrum/Camiv/coef_nlin.htm
O uso é como em qualquer outro task IRAF,
lineariza obj00*.fits obj00*.fits
2) badfaz
Esse task produz uma máscara de pixels ruins. Note a opção de deixar
o próprio task normalizar a imagem de flat-field fornecida. O príncipio de
funcionamento é o exame do histograma dos pixels do flat-field normalizado.
Caso queira fornecer os limites de corte no histograma, use o task imhistogram
de preferência com as opçõ logy+ e z1=0.5 z2=1.5
Uma dica prática para escolher os limites de corte no histograma das intensidades
dos pixles é o seguinte: a máscara de pixels ruins não é *muito*
diferente em J e H. Então, "sintonize" os limites de corte de modo a que
as máscaras em J e H fiquem parecidas.
3) prepara
Este task faz a maior parte das operações preparatórias para reduzir
os dados da CamIV, e já foi testado várias vezes em combate.
Dadas N imagens de um campo obtidas com pontilhamento, o task estima
o valor do céu a partir da mediana delas, subtrai esse céu de cada
imagem individual (mas você pode fornecer uma imagem de céu alternativa,
se quizer), divide as imagens individuais pelo flat-field normalizado
(opcional) e interpola sobre os pixels ruins (opcional).
Caso o parâmetro 'registra' seja 'yes', mostra na tela cada imagem individual
e pede para apontar um mesmo objeto de referência. A imagem registrada
pode ser armazenada e o programa entra novamente no modo cursor no qual
você pode apontar estrelas e teclando qualquer coisa, ter o valor dos
magnitudes brutas salvas no file logfile.
Se não optar por registrar as imagens (por exemplo, no caso de estrelas
padrão estamos interessados em saber os valores de mags instrumentais
uma a uma nas imagens individuais, para estimar melhor os erros no
ponto zero), o programa salvará no logfile apenas os resultados da
fotometria para um objeto no campo. Mostra esse valor na tela.
Exemplos:
cl> prepara obj0005 nima=7 flat=flatjn.fits mask=badj.pl registra+
output=teste.fits logf=99ago24.fot
neste caso, 7 imagens a partir da imagem obj0005 serão processadas, e o
resultado registrado será guardado na imagem teste.fits.
caso queira processar imagens que tenham nomes não ordenados, coloque-as
num file, uma por linha e use
cl> prepara @file_texto.txt etc etc
nesse caso, não á necessário dizer quantas imagens são.
4) reduz
Este task reduz duas imagens preparadas pelo task prepara, achando
as fontes pontuais em cada uma, as fontes comuns e as fontes exclusivas.
Digamos que os filtros usados sejam J e H. Os arquivos de saída
terão nomes
campo.J
campo.H
onde "campo" é o keyword 'object' no header das imagens
|
**************************************** os scripts ***************************************************
procedure lineariza(lista_in, lista_out)
# Lineariza uma lista de imagens da CamIV
# a correcao e' ADU(corr) = ADU(raw) + coef * ADU(raw)**2
#
# By Chico Jablonski, 21/05/99
# 13/09/99 - leva em conta nacc, lista_out adicionada
# 04/05/00 - flag flaglin adicionado, keyword lincof adicionado na imagem linearizada
# - trata passo a passo as imagens (mais seguro em termos de espaco em disco)
# 17/08/00 - corrigido bug de ordem de nomes quando ja' havia imagens linearizadas na lista.
string lista_in {prompt = "Lista de imagens de entrada"}
string lista_out {prompt = "Lista de imagens de saida"}
real coef {0.0000035, prompt = "Coeficiente quadratico"}
string *list1
begin
string infile,outfile,nome,nometemp,nomesai
real c
int nacc
cache sections, imarith, imdelete
set imclobber=yes
infile = mktemp("tmp$lin")
sections(lista_in, option="fullname", > infile)
outfile = mktemp("tmp$lin")
sections(lista_out, option="fullname", > outfile)
list = infile
list1 = outfile
while (fscan(list,nome) != EOF)
{
imgets(nome,param="FLAG-LIN", >& "dev$null")
if (imgets.value == '0')
{
j = fscan(list1,nomesai)
nometemp = mktemp("tmp$lin")
imgets(nome,param="NACCUM", >& "dev$null")
if (imgets.value != '0') { nacc = int(imgets.value) } else { nacc = 1 }
c = coef
if (nacc > 1) c = c/real(nacc)
imarith(nome,"*",c,nometemp,pix="r",calc="r",verb-)
imarith(nometemp,"+",1.,nometemp,pix="r",calc="r",verb-)
imarith(nometemp,"*",nome,nomesai,pix="r",calc="r",verb-)
imdelete(nometemp, ver-, >& "dev$null")
hedit(nomesai,fields='FLAG-LIN',value="1",add+,dele-,veri-,show-,up+)
hedit(nomesai,fields='LINCOEF',value=coef,add+,dele-,veri-,show-,up+)
print(nome,' -> ',nomesai)
}
else
{
j = fscan(list1,nomesai)
print ('Imagem ',nome,' ja linearizada...')
}
}
delete(infile, ver-, >& "dev$null")
delete(outfile, ver-, >& "dev$null")
end
*************************************************************************
procedure badfaz(flat,bad)
string flat {prompt="Imagem de flat-field"}
string bad {prompt="Imagem de bad-pixels"}
bool norm {no,prompt="Normaliza o flat-field"}
real z1 {0.85,prompt="Limiar inferior p/ BPM"}
real z2 {1.15,prompt="Limiar superior p/ BPM"}
bool auto {no,prompt="Determina limites automaticamente?"}
real pclip {0.2,prompt="Area do histograma rejeitada (%)"}
bool c513 {no,prompt="Forca marcacao da coluna 513?"}
# 29/08/00 - opcao para forcar coluna 513
begin
string ff,bpm,aux,aux1,tmp,tmp1
real l1,l2,clip,x,y,soma,soma1,x1,x2
bool normaliza,automatico
ff = flat
bpm = bad
normaliza =norm
l1 = z1 ; l2 = z2 ; clip = pclip
automatico = auto
tmp1= mktemp('tmpbad') // '.fits'
if (normaliza) {
unlearn imsurfit
tmp = mktemp('tmpbad') // '.fits'
imsurfit(inpu=ff,outp=tmp,xor=5,yor=5,func="spline3",xmed=8,ymed=8)
imarith(ff,'/',tmp,tmp1,ver-)
imdel(tmp,ver-)
} else { imcopy(ff,tmp1,ver-) }
if (automatico) {
unlearn imhistogram
aux = mktemp('tmpbad')
aux1= mktemp('tmpbad')
imhist(ima=tmp1,z1=0.5,z2=1.5,nbin=1024,list+,logy-, >>aux)
list = aux
soma = 0.0
while (fscan(list,x,y) != EOF) {
if (x>0.7 && x<1.3) { soma="soma" + y } } list # determina limiar inferior no histograma normalizado soma1="0.0" x1="clip/100./2." list="aux" while (fscan(list,x,y) !="EOF)" { if (x>0.7) {
soma1 = soma1 + y/soma
if (soma1 >= x1) { l1 = x; break } }
}
list=""
# determina limiar superior no histograma normalizado
sort(aux,col=1,num+,rev+, >>aux1)
soma1 = 0.0
x2 = x1
list = aux1
while (fscan(list,x,y) != EOF) {
if (x<1.3) { soma1="soma1" + y/soma if (soma1>= x2) { l2 = x; break } }
}
list=""
dele(aux,ver-)
dele(aux1,ver-)
}
print('Limiares p/ bad pixel map: ',l1,' e ',l2)
unlearn imreplace
imreplace(tmp1,value=2,low=l2)
imreplace(tmp1,value=2,upp=l1)
imarith(tmp1,'/',2,tmp1,ver-)
imreplace(tmp1,value=0,upp=0.9)
if (c513)
{
imreplace(tmp1//"[513,*]",value=2,low=INDEF,upp=INDEF)
}
if (stridx(bpm,".pl") == 0) {bpm = bpm // '.pl'}
imcopy(tmp1,bpm,ver-)
print('Codigo p/ FIXPIX interpolar em colunas = 1')
print('Codigo p/ FIXPIX interpolar em linhas = 2')
imdel(tmp1,ver-)
end
***************************************************************************
procedure prepara(ima1)
# F. Jablonski 09/99
# 14/10/99 - Corrigido erro no datapars.epadu
# 04/07/00 - Checagem da consistencia de filtros e objetos
# 13/11/01 - Forcei as modas a serem iguais depois da subtracao do ceu.
string ima1 {prompt="Imagem 1"}
int nimas {5,prompt="Numero de imagens"}
string sky {" ",prompt="Imagem de ceu (nada para calcular das imagens)"}
string metsky {"median",prompt="Metodo p/ combinar ceu"}
string mask {" ",prompt="Mascara de pixels maus"}
string flat {" ",prompt="Imagem de flatfield"}
bool registra {no,prompt="Registra imagens?"}
string metreg {"median",prompt="Metodo p/ combinar em registro"}
string output {" ",prompt="Nome da imagem registrada"}
string regions {" ",prompt="Arquivo de regioes p/ registro"}
string logfile {"teste.log",prompt="Logfile de saida"}
real diaf {5.,prompt="Raio do diafragma p/ fotometria"}
real ceu {7.,prompt="Raio do anel de ceu"}
real larg {20.,prompt="Largura do anel de ceu"}
real RON {15.,prompt="Ruido de leitura (e- rms)"}
real gain {4.5,prompt="Ganho (e-/ADU)"}
real angrot {0.,prompt="Angulo de rotacao da imagem registrada"}
struct *lista
begin
string ima,imav,it,itv,xyp,bad,imaout,versao,overlap
string nome, saida,raiz,aux,aux1,aux2,aux3,aux4,listain,tsky
string dif1,dif2,dif3,dif4,dif5,difs,ff,filtro
string shifts,regimas,regioes,templist,fotlog
real x1,x2,x3,x4,x5,dj,djm,rot
int j,j1,j2,nf,i1,i2,k1,k2,linear
bool reg,flaglista
struct linha
dj=0.0; djm=0.0
x1=0.0;x2=0.0;x3=0.0;x4=0.0;x5=0.0
linear=0
flaglista=no
versao = ""
show version | scan(versao,versao)
versao = substr(versao,1,7)
#cache wphot,txdump,display,imcombine,imarith,fixpix
if (! defpac("astutil") || ! defpac("digiphot") || ! defpac("apphot") ){
print('Astutil ou digiphot ou apphot nao carregados.'); bye}
unlearn imcombine
imcombine.rdnoise = RON
imcombine.gain = gain
unlearn wphot
unlearn photpars
photpars.apertures = str(diaf)
photpars.zmag = 25.
unlearn fitskypars
fitskypars.annulus = ceu
fitskypars.dannulus = larg
fitskypars.smooth = yes
fitskypars.sloclip = 5
fitskypars.shiclip = 5
unlearn centerpars
centerpars.cbox = 5.
unlearn datapars
datapars.readnoise = RON
datapars.epadu = gain
datapars.exposure="itime"
datapars.airmass="airmass"
datapars.filter="filter"
datapars.obstime="UT"
unlearn xregister
unlearn rvcorrect
keywpars.utmiddle="UT"
keywpars.exptime ="ITIME"
unlearn observatory
observatory(command="set",obsid="lna", >& "dev$null")
imaout = output
listain=mktemp('tmp')
difs = mktemp('tmp')
tsky = sky
ff = flat
bad = mask
saida = logfile
nf = nimas
if (nf <2) bye reg="registra" regioes="regions" # just in case... if (access('dif1.fits')){delete('dif?.fits',ver-)} if (access('dif1.fits.omag.1')){delete('dif?.fits.omag.*',ver-)} ima="ima1" rot="angrot" # verifique se e' uma lista do tipo @ if (stridx("@",ima) !="0)" { sections(imag="ima,opt='fullname',">listain)
nf=int(sections.nimages)
if (nf == 0) {
print('Template de entrada invalido.'); bye }
flaglista=yes
list=listain
j=fscan(list,aux)
list=""
j=stridx(".",aux)
i=stridx("0123456789",aux)
raiz = substr(aux,1,i-1)
j1 = int(substr(aux,i,j-1))
j2 = j1 + nf - 1
j1 = 1
j2 = nf
}
else {
# nome na linha de comando, ve se tem extensao .fits
if (stridx(".",ima)==0){ima = ima // '.fits'}
j = stridx(".",ima)
raiz = substr(ima,1,j-5)
# acha imagem final
j1 = int(substr(ima,j-4,j-1))
j2 = j1 + nf - 1
aux=""
for (j=j1; j<=j2; j="j+1)" { aux="str(j)" if (j<10 ) {aux="0" // aux} if (j<100 ) {aux="0" // aux} if (j<1000) {aux="0" // aux} print(raiz // aux // '.fits',>>listain)
}
}
list=listain
j=0
aux=""; aux1=""; aux2=""; aux3=""
while(fscan(list,aux) != EOF)
{
j=j+1
if (j == 1)
{
imgets(aux,'filter')
aux2 = imgets.value
filtro = aux2
imgets(aux,'object')
aux3 = imgets.value
print(' ')
print('Reduzindo ',aux3,' filtro ',aux2)
print(' ')
}
imgets(aux,'filter')
aux1 = imgets.value
if (aux1 != aux2) {beep ; print'Warning: Filtros diferentes...')}
imgets(aux,'object')
aux1 = imgets.value
if (aux1 != aux3) {beep ; print'Warning: Objetos diferentes...')}
imgets(aux,'flag-lin')
linear = linear + int(imgets.value)
}
if (linear != nf) {beep ; print'Warning: imagens nao linearizadas...')}
list=""
# lista das imagens a reduzir
aux=""
for (j=j1; j<=j2; j="j+1)" { aux="dif" // str(j) // '.fits' print(aux,>>difs)
}
# caso de 2 imagens somente e ceu nao fornecido pelo usuario
if ((nf == 2) && (! access(tsky)))
{
print ('Subtraindo ceu...')
list = listain
j = fscan(list,aux1)
j = fscan(list,aux2)
list = ""
aux3 = 'dif' // str(j1) // '.fits'
imarith(aux1,'-',aux2,aux3,pix='r',calc='r')
aux3 = 'dif' // str(j2) // '.fits'
imarith(aux2,'-',aux1,aux3,pix='r',calc='r')
}
else
{
# caso de > 2 imagens
# Se nao for dada imagem de ceu, usa mediana das imagens de entrada
if (! access(tsky))
{
print('Estimando ceu...')
tsky = mktemp('tmp') // '.fits'
if (metsky == 'median' )
{
imcombine(input='@'//listain,outpu=tsky,comb='median',reject='none',zero='mode')
}
else
{
imcombine(input='@'//listain,outpu=tsky,comb='average',reject='avsigclip',zero='mode')
}
}
else
{
imgets(tsky,'filter')
if (filtro != imgets.value) { beep ; print('Imagem de ceu em filtro diferente...')}
}
print ('Subtraindo ceu...')
imarith('@'//listain,'-',tsky,'@'//difs,ver-)
# forca todas a terem moda zero
print ('Forcando modas a zero...')
list=difs
while (fscan(list,aux) != EOF) {
imstat(images=aux,fields="mode", form-) | scan(x1)
imarith(aux,'-',x1,aux)
}
list=""
}
# trata flafield
if (access(ff)) {
imgets(ff,'filter')
if (filtro != imgets.value) { beep ; print('Imagem de flat-field em filtro diferente...')}
print ('Dividindo pelo flatfield...')
imarith('@'//difs,'/',ff,'@'//difs,ver-)}
if (access(bad)) {
imgets(bad,'filter')
if (filtro != imgets.value) { beep ; print('Imagem de pixels maus em filtro diferente...')}
print ('Aplicando mascara de pixels ruins...')
fixpix('@'//difs,masks=bad,lint=2,cint=1)}
# Faz funcionar imagens velhas da CamIV com o IRAF 2.11.3
if (versao == 'V2.11.3') { velho2novo('@'//difs,'@'//difs) }
if (! reg) {
print('Clique mouse sobre o objeto e tecle espaco (em cada imagem)...')
list = difs
j=0; aux=""
while (fscan(list,aux) != EOF) {
display(aux,1, >& "dev$null")
xyp = mktemp('tmp')
print(imcur, >xyp)
fotlog=mktemp('tmp')
wphot(ima=aux,coords=xyp,outp=fotlog,interac-,verify-)
delete(xyp,ver-)
# converte DATE-OBS para o novo formato (Y2K)
#velho2novo(aux,aux)
rvcorrect(imag=aux,header=no,imup=yes, >& "dev$null")
imgets(aux,'OBJECT')
nome = imgets.value
imgets(aux,'DATE-OBS')
aux3 = imgets.value
imgets(aux,param="HJD")
djm = real(imgets.value)-2450000.
aux1 = mktemp('tmp')
txdump(fotlog,'ifilter,xairmass,mag,merr,xcenter,ycenter,otime',yes, >>aux1)
delete(fotlog,ver-)
lista = aux1
i1=fscan(lista,aux2,x1,x2,x3,x4,x5)
lista = ""
delete(aux1,ver-)
printf("%11s %3s %6.3f %7.3f %6.3f %7.2f %7.2f %11.5f %s \n",nome,aux2,x1,x2,x3,x4,x5,djm,aux3)
printf("%11s %3s %6.3f %7.3f %6.3f %7.2f %7.2f %11.5f %s\n",nome,aux2,x1,x2,x3,x4,x5,djm,aux3, >>saida)
}
list=""
}
else {
print('Aponte *o mesmo* objeto de referencia em todas as imagens...')
xyp = mktemp('tmp') # shifts do objeto de referencia
shifts = mktemp('tmp') # log do xreg
regimas = mktemp('tmp') # imagens registradas
for (j=j1; j<=j2; j="j+1)" { aux="tmpr" // str(j) // '.fits' print(aux,>>regimas)
}
j=0
list = difs
while (fscan(list,aux) != EOF) {
display(aux,1, >& "dev$null")
print(imcur, >>xyp)
if(j==0) {type(xyp, >>xyp)
type(xyp) | scan(x1,x2) # p/ uma pequena secao p/ registro
j=1}
#velho2novo(aux,aux)
rvcorrect(imag=aux,header=no,imup=yes, >& "dev$null")
imgets(aux,param="HJD")
djm = djm + real(imgets.value) - 2450000.0
}
list=""
djm = djm / (j2-j1+1)
aux = 'dif' // str(j1) // '.fits' # referencia e' a 1a. da sequencia
print('Registrando...')
aux4 = mktemp('tmp')
if (access(regioes))
{
xregister(input='@'//difs,refer=aux,regions=regioes,shifts=shifts,output='@'//regimas,coords=xyp,background='median', > aux4)
}
else
{
# secao fixa para registro, de 64 x 64 pixels em torno da referencia
i1 = int(x1) - 32 ; i2 = i1 + 64
k1 = int(x2) - 32 ; k2 = k1 + 64
aux1 = '['//str(i1)//':'//str(i2)//','//str(k1)//':'//str(k2)//']'
xregister(input='@'//difs,refer=aux,regions=aux1,shifts=shifts,output='@'//regimas,coords=xyp,background="median", > aux4)
}
list = aux4
type(aux4)
while(fscan(list, linha) != EOF) { }
i = stridx("[",linha)
j = stridx("]",linha)
overlap=substr(linha,i,j)
list=""
delete(aux4,ver-)
if (imaout == "" || imaout == " ") {
imaout = raiz // str(j1) // '-' // str(j2) // '.fits' } # ima combinada
#print ('Forcando modas a zero...')
#list=regimas
#while (fscan(list,aux) != EOF) {
# imstat(images=aux//'[256:768,256:768]',fields="mode", form-) | scan(x1)
#print(x1)
# imarith(aux,'-',x1,aux)
# }
#list=""
print('Combinando...')
if (metreg == 'median')
{
imcombine(input='@'//regimas,outpu=imaout,comb='median',reject='none',zero='mode')
}
else
{
imcombine(input='@'//regimas,outpu=imaout,comb='average',reject='avsigclip',zero='mode')
}
if (rot == 90.)
{
rotate(input=imaout,output=imaout,rotation=rot)
j = stridx(",",overlap)
aux = '[' // substr(overlap,j+1,strlen(overlap)-1) // ',' // substr(overlap,2,j-1) // ']'
overlap = aux
}
hedit(imaout,fields='overlap',value=overlap,add+,dele-,verify-,show-,up+)
display(imaout,1,bpd="none", >& "dev$null")
xyp = mktemp('tmp')
aux1 = mktemp('tmp')
print('Aponte objetos, ctrl-d termina')
rimcur(imaout, >>xyp)
fotlog=mktemp('tmp')
wphot(imaout,coords=xyp,outp=fotlog,interac-,verify-)
imgets(imaout,'OBJECT')
nome = imgets.value
txdump(fotlog,'ifilter,xairmass,mag,merr,xcenter,ycenter,otime',yes, >>aux1)
delete(fotlog,ver-)
imgets(imaout,'DATE-OBS')
aux3 = imgets.value
list = aux1
while (fscan(list,aux2,x1,x2,x3,x4,x5,aux1) != EOF ) {
printf("%11s %3s %6.3f %7.3f %6.3f %7.2f %7.2f %11.5f %s \n",nome,aux2,x1,x2,x3,x4,x5,djm,aux3)
printf("%11s %3s %6.3f %7.3f %6.3f %7.2f %7.2f %11.5f %s\n",nome,aux2,x1,x2,x3,x4,x5,djm,aux3, >>saida) }
list=""
imdel('@'//regimas, ver-)
}
delete ('tmp*', ver-)
delete ('dif*', ver-)
delete ('*omag*',ver-)
end
****************************************************************************
procedure reduz(ima1,ima2)
# Reduz imagens camiv em filtros diferentes (ou epocas diferentes)
# produzindo uma lista de fontes comuns, com os indices de cor ou
# diferencas de magnitude, e as fontes detectadas so' numa ou noutra
# imagem
string ima1 {prompt="Imagem no filtro 1"}
string ima2 {prompt="Imagem no filtro 2"}
real fwhm1 {prompt="FWHM da psf na imagem 1"}
real fwhm2 {prompt="FWHM da psf na imagem 2"}
real limiar1 {5.,prompt="Limiar (em sigmas) p/ imagem 1"}
real limiar2 {5.,prompt="Limiar (em sigmas) p/ imagem 2"}
real width1 {1.5,prompt="Width of conv kernel in 1"}
real width2 {1.5,prompt="Width of conv kernel in 2"}
real Zp1 {25.,prompt = "Ponto zero p/ imagem 1"}
real Zp2 {25.,prompt = "Ponto zero p/ imagem 2"}
real ext1 {.10,prompt = "Extincao filtro 1"}
real ext2 {.10,prompt = "Extincao filtro 2"}
real errofot {0.1,prompt = "Limiar de erro na fotometria"}
real diaf1 {5.,prompt="Raio do diafragma p/ fotometria em 1"}
real diaf2 {5.,prompt="Raio do diafragma p/ fotometria em 2"}
real ceu {7.,prompt="Raio do anel de ceu"}
real larg {20.,prompt="Largura do anel de ceu"}
real RON {15.,prompt="Read-out noise (e- rms)"}
real gain {4.5,prompt="Ganho (e-/ADU)"}
begin
string im1,im2,coo1,coo2,aux1,aux2,aux3,aux4,saida1,saida2,fotmag1,fotmag2
string campo,lixo1,lixo2,lixo3,lixo4,lixo5,lixo6,aux5,aux6
real sig1,sig2,lim1,lim2,w1,w2,x1,x2,y1,y2,fw1,fw2
real xmin, xmax, ymin, ymax
real ruido,ganho,x,y
int k1,k2
im1 = ima1
im2 = ima2
fw1 = fwhm1
fw2 = fwhm2
lim1 = limiar1
lim2 = limiar2
w1 = width1
w2 = width2
coo1 = mktemp('tmpcom'); coo2 = mktemp('tmpcom')
unlearn daofind
unlearn findpars
unlearn wphot
unlearn datapars
unlearn photpars
unlearn centerpars
unlearn fitskypars
unlearn average
datapars.readnoise = RON
datapars.epadu = gain
datapars.exposure="itime"
datapars.airmass ="airmass"
datapars.filter ="filter"
datapars.obstime ="UT"
centerpars.cbox=5.
fitskypars.salgorithm="centroid"
fitskypars.annulus = ceu
fitskypars.dannulus = larg
fitskypars.smooth = yes
fitskypars.sloclip = 5.
fitskypars.shiclip = 5.
# checa se e' o mesmo objeto nas duas imagens
imgets(im1,param="object")
aux1 = imgets.value
imgets(im2,param="object")
aux2 = imgets.value
if (aux1 != aux2) {
print('Verificar se eh o mesmo objeto nas duas imagens...')
bye }
imgets(im1,param="filter")
lixo1 = imgets.value
imgets(im2,param="filter")
lixo2 = imgets.value
if (lixo1 == lixo2)
{
saida1 = aux1 // '.1'
saida2 = aux2 // '.2'
print('Duas imagens de ',aux1,' no filtro ',lixo1)
}
else
{
saida1 = aux1 // '.' // lixo1
saida2 = aux2 // '.' // lixo2
print('Imagens de ',aux1,' nos filtros ',lixo1,' ',lixo2)
}
date | scan(lixo1,lixo2,lixo3,lixo4,lixo5,lixo6) # data da reducao
print('# reducao em ',lixo1,' ',lixo2,' ',lixo3,' ',lixo4,' ',lixo5,' ',lixo6,>saida1)
print('# imagem = ',im1, >>saida1)
imgets(im1,param='object')
print('# objeto = ', imgets.value,>>saida1)
imgets(im1,param='date-obs')
print('# date-obs = ',imgets.value, >>saida1)
imgets(im1,param='airmass')
print('# airmass = ',imgets.value,>>saida1)
imgets(im1,param='HJD')
print('# HJD = ',imgets.value, >>saida1)
print('# fwhm = ',fwhm1, >>saida1)
print('# limiar = ',limiar1, >>saida1)
print('# width = ',width1, >>saida1)
print('# Zp = ',Zp1, >>saida1)
print('# ext = ',ext1, >>saida1)
print('# errofot = ',errofot, >>saida1)
print('# diaf = ',diaf1, >>saida1)
print('# ceu = ',ceu, >>saida1)
print('# larg = ',larg, >>saida1)
print('# RON = ',RON, >>saida1)
print('# gain = ',gain, >>saida1)
print('# reducao em ',lixo1,' ',lixo2,' ',lixo3,' ',lixo4,' ',lixo5,' ',lixo6,>saida2)
print('# imagem = ',im2, >>saida2)
imgets(im2,param='object')
print('# objeto = ', imgets.value,>>saida2)
imgets(im2,param='date-obs')
print('# date-obs = ',imgets.value, >>saida2)
imgets(im2,param='airmass')
print('# airmass = ',imgets.value,>>saida2)
imgets(im2,param='HJD')
print('# HJD = ',imgets.value, >>saida2)
print('# fwhm = ',fwhm2, >>saida2)
print('# limiar = ',limiar2, >>saida2)
print('# width = ',width2, >>saida2)
print('# Zp = ',Zp2, >>saida2)
print('# ext = ',ext2, >>saida2)
print('# errofot = ',errofot, >>saida2)
print('# diaf = ',diaf2, >>saida2)
print('# ceu = ',ceu, >>saida2)
print('# larg = ',larg, >>saida2)
print('# RON = ',RON, >>saida2)
print('# gain = ',gain, >>saida2)
# coordenadas do objeto de referencia (comum) para as saidas
aux1 = mktemp('tmpcom')
aux2 = mktemp('tmpcom')
aux3 = mktemp('tmpcom')
aux4 = mktemp('tmpcom')
display(im1,1, >& "dev$null")
print ('Aponte objeto de referencia na imagem 1...')
print (imcur,>>aux1)
display(im2,2, >& "dev$null")
print ('Aponte objeto de referencia na imagem 2...')
print (imcur,>>aux2)
# centra precisamente
center(im1,coords=aux1,output=aux3,interac-,veri-)
center(im2,coords=aux2,output=aux4,interac-,veri-)
#type(aux1)
#type(aux2)
x1=0.0; x2=0.0; y1=0.0; y2=0.0
txdump(aux3,fields='xcenter, ycenter',expr=yes) | scan(x1,y1)
txdump(aux4,fields='xcenter, ycenter',expr=yes) | scan(x2,y2)
print(x1,' ',y1, >>saida1)
print(x2,' ',y2, >>saida2)
#type(saida1)
#type(saida2)
#estima sigma do ceu (no ceu do objeto de referencia)
aux5 = mktemp('tmpcom')
wphot(im1,coords=saida1,output=aux5,interac-,veri-)
txdump(aux5,fields='stdev',expr=yes) | scan(sig1)
aux6 = mktemp('tmpcom')
wphot(im2,coords=saida2,output=aux6,interac-,veri-)
txdump(aux6,fields='stdev',expr=yes) | scan(sig2)
print('Achando objetos...')
center(im1,coords=aux1,output=aux3,interac-,veri-)
findpars.threshold = limiar1
findpars.nsigma = w1
datapars.sigma = sig1
datapars.fwhmpsf = fw1
daofind(im1,output=coo1,veri-)
#tvmark(1,coords=coo1)
findpars.threshold = limiar2
findpars.nsigma = w2
datapars.sigma = sig2
datapars.fwhmpsf = fw2
daofind(im2,output=coo2,veri-)
#tvmark(2,coords=coo2)
print('Sharpness da PSF na imagem 1:')
txdump(coo1,fields='sharpness',expr='yes') | average(option='new_sample')
print('Roundness da PSF na imagem de 1:')
txdump(coo1,fields='sround',expr='yes') | average(option='new_sample')
print('Sharpness da PSF na imagem 2:')
txdump(coo2,fields='sharpness',expr='yes') | average(option='new_sample')
print('Roundness da PSF na imagem de 2:')
txdump(coo2,fields='sround',expr='yes') | average(option='new_sample')
# realiza a fotometria
print('Realizando a fotometria...')
fotmag1 = mktemp('tmpcom')
photpars.zmag = Zp1
photpars.apertures = str(diaf1)
wphot(im1,coords=coo1,output=fotmag1,interac-,veri-)
fotmag2 = mktemp('tmpcom')
photpars.zmag = Zp2
photpars.apertures = str(diaf2)
wphot(im2,coords=coo2,output=fotmag2,interac-,veri-)
# seleciona os campos de interesse
hselect(im1,fields="naxis[1],naxis[2]",expr=yes) | scan(xmax,ymax)
xmax = xmax - diaf1/2.; ymax = ymax - diaf1/2.
xmin = 1+diaf1/2. ; ymin = 1+diaf1/2.
campo = 'xcenter >= ' // str(xmin) // ' && xcenter <= ' // str(xmax) // ' && ' campo="campo" // 'ycenter>= ' // str(ymin) // ' && ycenter <= ' // str(ymax) campo="campo" // ' && mag !="INDEF" && merr < ' // str(errofot) campo="campo" // ' && PIER="=" 0' txdump(fotmag1,fields="xcenter,ycenter,mag,merr,stdev" ,expr="campo,">>saida1)
hselect(im2,fields="naxis[1],naxis[2]",expr=yes) | scan(xmax,ymax)
xmax = xmax - diaf2/2.; ymax = ymax - diaf2/2.
xmin = 1+diaf2/2. ; ymin = 1+diaf2/2.
campo = 'xcenter >= ' // str(xmin) // ' && xcenter <= ' // str(xmax) // ' && ' campo="campo" // 'ycenter>= ' // str(ymin) // ' && ycenter <= ' // str(ymax) campo="campo" // ' && mag !="INDEF" && merr < ' // str(errofot) campo="campo" // ' && PIER="=" 0' txdump(fotmag2,fields="xcenter,ycenter,mag,merr,stdev" ,expr="campo,">>saida2)
# mostra com tvmark os objetos medidos
#campo = str(diaf) // ',' // str(ceu) // ',' // str(ceu+larg)
#tvmark(1, coords=saida1, radii=campo, color=212)
#tvmark(2, coords=saida2, radii=campo, color=212)
delete(coo1//','//coo2//','//fotmag1//','//fotmag2,ver-)
delete(aux1//','//aux2//','//aux3//','//aux4)
delete(aux5//','//aux6)
end