Linux Astronomia

Gnuplot

Creare grafici utilizzando dati del MPC (e non solo)
# grafico inclinazione semiasse maggiore
set xlabel "a (UA)"
set ylabel "incl. "
set grid x2tics
set x2tics ("2:9" 1.9,"1:4" 2.06,"2:7" 2.25, "1:3" 2.5, "2:5" 2.82, "1:2" 3.27, "3:5 "3.7)
set x2label " Lacune di Kirkwood "
set title " Asteroidi San Marcello (104) "
plot 'SM_Numerati.dat' u 2:3 w p ps 1 lt 1  pt 5 ti "Numerati", /
'SM_unnumerati.dat' u 2:3 w p ps 1 lt 3  pt 5 ti "non Numerati"

# come scrivere testo nell'area del  grafico
set termoption enhanced # importante altrimento non riconosce i caratteri
set label 1 at -0.65, 1.0
set label 1 "These examples require no extra fonts:" tc lt 3 # esempio label 1
set label 2 at 0.4, 13  "10^{-2}   A_{j,k} " tc lt 6         # label colore 
set grid  lt 0 lw 1 lc rgb "magenta"                         # colore griglia
set grid ytics lt 0 lw 1 lc rgb "green"                      # colore 


# Per salvare un grafico in una immagine in formato PNG

set termina png
set output "sm_E_Unm-f001.png"

# Per calcolare  la frequanza di dati e creare il relativo grafico

classe(x, a) = a*int(x/a)
a=0.05
Nt=1
plot 'SM_unnumerati.dat' u (classe($1,a)):(1./Nt) smooth frequency t '' w boxes, /
 'SM_unnumerati.dat' u (classe($1,a)):(1./Nt) smooth frequency t '' w impulses 4,/
'SM_unnumerati.dat' u (classe($1,a)):(1./Nt) smooth frequency t '' w points 3

# altro modo 
plot 'V_B.tsv' u (classe($10,a)):(1./Nt) s f t ' ' smooth frequency  w boxes lc 1,  'V_B.tsv' u (classe($10,a)):(1./Nt) smooth frequency  w impulses 4 t '', 'V_B.tsv' u (classe($10,a)):(1./Nt) smooth frequency  w points 3 t '

script per creare grafico delle frequenze

CALCOLARE LA LINE DI DEGRESSIONE E MOSTRARLA IN GRAFICO

f(x) = m*x + b
fit f(x) '888.txt u 1:2 via m,b
plot '888.txt' u 1:2 pt 7 title 'mag B', f(x) title 'fit'

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

m               = 1.01952          +/- 0.004464     (0.4379%)
b               = 0.696638         +/- 0.07518      (10.79%)

#deviazione standar
stddev_y = sqrt(FIT_WSSR / (FIT_NDF + 1 ))
plot '888.txt' u 1:2 pt 7 title 'mag B', f(x) title 'fit', '888.txt' u 1:2:(stddev_y) w ye


#  formula del grafico f(x)= 1.01952(x)+0.696638 

Calcolare Valore medio

stat 'press.dat' u 3 name "A"
stat 'press.dat' u 2 name "B"
plot 'press.dat' u 3 w l, A_mean, '' u 2 w l, B_mean

Deviazione Standar calcolo
stat 'press.dat' u 3 every ::0::1 
plot 'press.dat' u 3:(9.5) w yerrorbar, A_mean, '' u 2:(6.5) w yerrorbars, B_mean

Min Max

 stat 'press.dat' u 2
plot 'press.dat' u 3:(9.5) w yerrorbar, A_mean, '' u 2:(6.5) w yerrorbars, B_mean, STATS_min_y w l



Grafico 3d

set pm3d
set cblabel "colour gradient" tc rgb 'red' # label bel grafico  

set palette model RGB defined (0 "dark-blue", 1 "dark-blue", 1 "blue", 2 "blue", 2 "#00BFFF", 3 "#00BFFF", 3 "#00FFFF", 4 "#00FFFF", 4 "yellow", 5 "yellow", 5 "dark-yellow", 6 "dark-yellow", 6 "orange", 7 "orange", 7 "#FF4500", 8 "#FF4500", 8 "magenta", 9 "magenta", 9 "dark-red", 10 "dark-red" )

# set palette model RGB defined (0 "blue", 1 "dark-blue", 2 "yellow", 3 "orange", 4 "red", 5 "dark-red" )
# oppure
set palette defined (0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1)
splot 'V_B.tsv' u 1:2:10 with points palette pt 5 ps 0.5
 # adesso ho un grafico in 3d  con colori
 # ora voglio impostarlo in piano ma con 3 dati
plot 'V_B.tsv' u 1:2:10 with points palette pt 7 ps 0.8 t '

set cbrange [-1:3]         # imposta in maniera fissa la scala dei colori
gnuplot> set zrange [-1:3]

set datafile separator ',' #cvs file esempio

 # Gnuplot color map 
 # vedi il risultato

colori
#cambiare colore bordo
set border 15 ls 2 lw 1 lc 8
#cambiare colore sfondo
set object 1 rectangle from screen 0,0 to screen 1,1 fillcolor rgb"black" behind
set border 15 ls 2 lw 1 lc 9
set x2label 'Temperature 'tc rgb "green"
set grid x2tics lt 0 lw 1 lc 4
# legenda 
set key tc rgb 'white'
set title 'pippo' tc rgb "white"
set key bmargin   # la mette fuori dal grafico

# garfico aspetto quadrato
set size square


Filtri sui dati
plot  'NGC1245/V_Br.tsv'   u 3:4:($5>1.8 & $5<3 ? $5 :1/0)  w point palette pt 7 ps .5  t ''
plot [>1: ][ :>34] 'nomefile' u 1:2 


Creare una legenda per gli assi  partendo da un file Tabulato
Tipo    Teff    V-B     V-R
O8-9V   31111   0.31    -0.15
O5-7V   33333   0.32    -0.15
O9.5V   29999   0.30    -0.14
B0V     26500   0.30    -0.13
B0.5V   24500   0.28    -0.12

cat cat HR-VBRI.csv | cut -f 1,2 | sed -n '2~2P' | awk '{ print "\""$1"\" " $2 "" }' | tr '\n' ',' 

risultato

"O8-9V" 31111,"O9.5V" 29999,"B0.5V" 24500,

data
set xdata time
set timefmt "%Y-%m-%d"

set yrange[] reverse

Gaussiana 

distribuzione asteroidi in funzione del valore H si crea il grafico con la frequenza dei dati
si esportano i dati con il comando
set table "hist.dat"
abbiamo creato un file con i dati dela frequenza
poi
unset table "hist.dat"
 un po grossolano ma funziona:
PI=3.14; s=1; m=1; a=1;
gauss(x) = a/(2*PI*s**2)**0.5*exp(-(x-m)**2/(2*s**2))
fit gauss(x) 'hist.dat' using 1:2 via s, m, a
plot 'hist.dat'using 1:2 w boxes, gauss(x)
vedi risultato 


Animazione

set xdata time
set timefmt "%H.%M.%S"
set yrange [16.5 : 18 ] reverse
set xrange ["02.50.00" : "03.20.00"]
set title " 104 San Marcello - Haumea 21 Gen 2017 " 
set xlabel " Time UT "
set ylabel " Mag. "
set xtics rotate 90
set term png  
do for [i=2:131] {
set output sprintf('Haumea_frame%03.0f.png',i)
plot  'Haumea_dat2.dat' u 2:3:($1"<"i ? $1 :1/0)  pt 7 t'Haumea' ...
}

MPC data to Xephem

Comandi utili

# awk per determinare quanti decimali stampare 
cat -30 T5_01.fit.dat | sed '/#/d' | awk  -v OFMT=%2.4f '{print $6*1, $7*1}'
# awk "taglaire" i decimali 
cat  T5_01.fit.dat | sed '/#/d' | awk   '{print substr($1,1,2),substr($6,1,5), substr($7,1,5)}'

# confrontare due file
# 
 cat  M67B_\@.tsv| sed '1d' | awk '{ print $1,$2,$3 }' | while read a b c; do cat M67R_@R.tsv| sed '1d' | awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)}  { print $1,$2,$3, '$a','$b', '$c', abs($1/'$a'), abs($2/'$b'), abs(($1/'$a'))-abs(($2/'$b')), $3-'$c' }' ;done | sort | awk '{ if ( $9 <= 0.00005)  print $0}' | tr ' ' '\t'

# rinominare file da xxxx.jpg a 000?.jpg
ls *.jpg | nl -nrz -w4 | while read a b; do mv "$b" $a.jpg; done  





DARK FLAT con Eclipse
# usa eclispe per rimuvere il dark
# 2007 11 02
# backman
ls *D*.[fF]*[tT]*
echo -n "digitare l' immagine prima  sommare  "
read DARK
ls *.[fF]*[tT]*
echo -n "digitare il suffisso delle immagini (es 17P*N.*) : "
read IMAGE
rm out
for var in `ls $IMAGE*$1`
do
ccube "@$var @$DARK  - " out.fit
cp out.fit dark_$var
rm out.fit
done




# usa eclispe per rimuvere il dark
# 2007 11 02
# backman
ls *.[fF]*[tT]*
echo -n "digitare il suffisso delle immagini (es 17P*N.*) : "
read IMAGE
for var in `ls $IMAGE*$1`
do
ccube "@$var @$IMAGE  - " IMAGE.fit
cp out.fit dark_$var
done


SEX e DS9
Utilizzzo di sextracor e DS9

# backman asa
# script che utilizza sex 
# per creare multifile 
# il file copia nella ditectory
# in cui viene lanciato lo script
# i file di configurazione 
# di sex

echo " lo script copia alcuni file necessari al funizonamento"
echo
echo .
# file da copiare
cp /home/backman/SCRIPT_ASA/default.* .
# adesso inzia lo script
echo -n " inserire suffisso immagine (* tutte): "
read IMGSEX
rm asa
for var  in `ls $IMGSEX*$1.[fF]*[tT]*`
do
sex $var 
wc -l test.cat > asa
cp test.cat $var.dat
echo $var.dat >> asa
sort < asa
rm asa

done

# rimuve i file copiati
rm default.*

Ecco cosa crea:

#   1 MAG_ISO         Isophotal magnitude                             [mag]
#   2 MAG_AUTO        Kron-like elliptical aperture magnitude         [mag]
#   3 MAG_BEST        Best of MAG_AUTO and MAG_ISOCOR                 [mag]
#   4 X_IMAGE         Object position along x                         [pixel]
#   5 Y_IMAGE         Object position along y                         [pixel]
#   6 ALPHA_J2000     Right ascension of barycenter (J2000)           [deg]
#   7 DELTA_J2000     Declination of barycenter (J2000)               [deg]
#   8 FWHM_IMAGE      FWHM assuming a gaussian core                   [pixel]
 19.0390  18.9915  18.9915    277.761     29.192 133.0084427 +11.9990342     3.62
 16.2839  16.2785  16.2785    451.813     25.370 132.9322309 +12.0153072     3.57


Adesso che abbiamo estrapolato le stelle dalle immagini

creaiamo per ogni immagine un file di registro per ds9

rm asa
for var in `ls $1*.dat`
do
cat $var | sed -e /#/d |awk '{print ("circle""("$13"," $14","5"\") # color=red")}' > $var.reg
wc -l $var.reg > asa
echo $var.reg >> asa
sort < asa
rm asa

done


ecco il risultato>
circle(506.502,16.647,5") # color=red
circle(817.063,14.446,5") # color=red
circle(934.623,13.380,5") # color=red
circle(873.211,13.812,5") # color=red
circle(187.395,13.995,5") # color=red






Header FITS 

#2010 03 02
# prova scritp per estrapolare dati dai file fits
# questo per i file aquisiti con ccdsoft e maximdl
#(.[fF]*[tT]*)
#
#find /home/backman/FotoAstro/2008/ -name "*.[fF]*[tT]*" > aa.txt
#
#cerca i file fits nelle directory e crea il file aa.txt
# gethead -tpu -n 2  date-obs  exptime filter ccd-temp objctra objctdec IMAGETYP naxis naxis1  @bb.txt | cut -c 24-100
# il parametro p mostra il path completo . sotto l'ho tolto per adesso non mi serve

echo -n " inserire  percorso completo dei file da cercare /home/backman/ppp/123/  :"
read DIR

find /home/backman$DIR -name "*.[fF]*[tT]*" > filehead.txt

gethead -tu -n 2 OBJECT  date-obs  exptime filter ccd-temp objctra objctdec AIRMASS IMAGETYP naxis naxis1  @filehead.txt 
echo "creo il file header $DIR .txt"

gethead -tu -n 2 OBJECT  date-obs  exptime filter ccd-temp objctra objctdec AIRMASS IMAGETYP naxis naxis1  @filehead.txt > header.txt

rm filehead.txt

echo " FATTO !!!!!"

esempio:
 
04.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T20:57:58.687 30.00   Lunar   -2.74   05 53 47.875    +19 41 54.75    1.27    Light Frame     2       382
05.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T20:58:49.75  30.00   Lunar   -4.82   05 53 47.875    +19 41 54.75    1.29    Light Frame     2       382
06.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T20:59:26.968 30.00   Lunar   -5.65   05 53 47.875    +19 41 54.75    1.29    Light Frame     2       382
07.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T21:00:44.843 30.00   Lunar   -4.82   05 53 53.875    +19 41 54.84    1.29    Light Frame     2       382
08.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T21:01:21.921 30.00   Lunar   -4.82   05 53 53.875    +19 41 54.84    1.29    Light Frame     2       382
09.217P_LINEAR.FIT        217P_LINEAR     2010-03-02T21:01:58.718 30.00   Lunar   -4.82   05 53 53.875    +19 41 54.84    1.29    Light Frame     2       382



Inserire dati negli header fits 


#!/bin/bash
# programma inserimento dati immagini fits
# sett 2006
# use wctolls software
# p. bacci 
#
echo -n " inserisci le coordinate RA [00:00:00] = "

read RA

echo -n " inserisci le coordinate DEC [-00:00:00] = "
read DEC

echo " "
echo " "
ls  *.fits
echo " "
echo -n " inserisci il nome dei file da MODIFICARE senza estenzione .fits!! : "
read FILE
echo " " 
echo "============================================= "
echo " parametri x y w e rot settati su 0.0052 187° "
echo " ============================================="
echo " "
echo " "
echo " --------- dati inseriti ----------"
echo " RA: $RA ///// DEC: $DEC  file : $FILE.fits"
echo "-----------------------------------"
echo " "
echo " "

echo -n GO "premi - 1  - per confermare......  -  0  - per inserire di nuovo  "
read GO

if [ $GO = 1 ]
then
echo " inserimentodati  dati"

echo -n " premi -- a  -- per inserire RA DEC  CDELTA e CROT2 == altro  == inserire  CDELTA E CROT2 "
read S


if [ $S = a ]
then
echo "tutti i dati"
echo " $FILE.fits"

sethead -nkhv 'CDELT1=0.00052' 'CDELT2=-0.00052' $FILE.fits

echo " $FILE "
sethead -khv 'CROTA2=187' $FILEe*e.fits
sethead -khv "RA= $RA"  "DEC= $DEC"  $FILE*e.fits

gethead -ah CDELT1 CDELT2 CROTA2 RA DEC OBJECT $FILE*e.fits
  
 else 
echo "solo delta e crot2"

sethead -nkv 'CDELT1=0.00052' 'CDELT2=-0.00052' $FILE.fits
sethead -kv 'CROTA2=187' $FILE*e.fits
gethead -ah CDELT1 CDELT2 CROTA2 RA DEC OBJECT $FILE*e.fits

fi


# scelta forna indietro
else 


echo "sbagliato tutto !!!!!!"
./afitins
fi 
exit 0;




NEOCP 
#!/bin/bash
# 10/05/2009
# questo script crea da due immagini iniserite nella cartella
# un immagine nella uqale si vedono gli oggetti NEOCP
# 

cd /home/backman/ASA_SCRIPT/skyneo

# crea immagini archivio immagini old.
num=`ls old_*.jpg | wc -l`
cp skyneocp.jpg old_skyneocp_$num.jpg

rm skyneocp.jpg
rm 555.txt
rm ToConfirm.html

wget http://www.cfa.harvard.edu/iau/NEO/ToConfirm.html


cat ToConfirm.html  | grep CHECKED | grep UT |cut -c 66-156 |sed s/""//g | sed s%""%%g | sed s%]%%g |awk '{print ( $1, $8, $9, $12, $16)}' | awk '{print ($2, $3/60, $4, $1, $5 )}' | awk '{print (($1+$2)*42.8,$3*2.8, $4, $5 )}'| awk '{print (1061-$1, 261-$2, $3, $4) }' > 555.txt


#cat ToConfirm.html  grep CHECKED | grep UT |cut -c 66-156 |sed -e 's%%%g' |sed -e s%]%%g | awk '{print ( $1, $9, $10, $13, $17)}' |awk '{print ($2, $3/60, $4, $1, $5 )}' | awk '{print (($1+$2)*42.8,$3*2.8, $4, $5 )}'| awk '{print (1061-$1, 261-$2, $3, $4) }' > 555.txt


mm=`cat -n 555.txt | cut -c 5-6`


cp cielo2.jpg skyneocp.jpg


for i in $mm
do


ar=`head -n $i 555.txt  | tail -1 | awk '{print ($1)}'`
dec=`head -n $i 555.txt  | tail -1 | awk '{print ($2)}'`
ldec=`head -n $i 555.txt  | tail -1 | awk '{print ($2 - 5)}'`
neo=`head -n $i 555.txt  | tail -1 | awk '{print ($3)}'`
mag=`head -n $i 555.txt  | tail -1 | awk '{print ($4)}'`
l="${neo:4}_$mag"
k="+$ar+$dec"


#ldec calcolo posizone label
# l = label nome sull'immagine
# k posizione asteroidi
# kl posizone label
# ${neo:3}


kl="+$ar+$ldec"
echo "neocp: $neo - mag: $mag ---- label: $kl -- punto: $k -- declinazione label: $ldec"

convert -background '#00000080' -fill gray -pointsize 12  label:$l miff:- |  composite  -geometry $k  - skyneocp.jpg skyneocp.jpg
composite  -geometry $kl neo_01.gif  skyneocp.jpg skyneocp.jpg

#wget http://scully.cfa.harvard.edu/~cgi/ShowObsOrbs.COM?Obj=$neo&obs=y >> obs_neo.txt
#wget http://scully.cfa.harvard.edu/~cgi/ShowObsOrbs.COM?Obj=$neoD&orb=y >> orb_neo.txt



done



convert -background '#00000080' -fill white label:` date +%H:%M:%S-%d/%m/%Y` miff:- |
 composite -gravity south -geometry +0+3 - skyneocp.jpg skyneocp.jpg

#display skyneocp.jpg

#-------------------------
# crea pagina neocp
rm neocpA.txt
rm neoTOcp.php
rm no_neo.txt

cat ToConfirm.html | grep "None at the moment"  | sed -e 's/None at the moment/Non ci sono Neocp da confermare al momento/g' > no_neo.txt
cat ToConfirm.html | grep CHECKED | grep "name=\"Obj"  > neocpA.txt
cat neocp_u.php no_neo.txt neocpA.txt neocp_d.php > neoTOcp.php


# invio delle immagine al server

userid='user'
passid='password'

ftp -niv << EOF
open ftp.domain.org
user $userid $passid

#
#pwd
#ls

# cartella remota
cd meteo
delete skyneocp.jpg

#cartella locale
lcd /home/backman/ASA_SCRIPT/skyneo

mput skyneocp.jpg
mput neoTOcp.php

bye
EOF