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