### Script piapp pour afficher les resultats du programme 
### de reconstruction de beam par la technique de double FFT
###  Reza , Sep 2006
# ------------------------------------------


# 0/ Verification, modif des parametres
#  Modifier si besoin le code dans treccyl.cc et faire ensuite make
csh> make 
# Le programme se trouve dans Objs/treccyl 
# 1/ Faire tourner le programme - reconstruction plan X-freq (en dehors de piapp)
Objs/treccyl X onecyl.ppf


# 2/ Affichage resultat
#  Afficher le plan source reconstruit (les frequences le 
# long de l'axe Y (les lignes de la matrice), et les directions
# des sources (angles) suivant l'axe Ox (colonnes)
# Le zero des angles se trouve au milieu
delobjs *
openppf onecyl.ppf
disp recsrcplane 'lut=sqrt,0.,0.1 invy grey32 win'

# 2.b/ le ntuple ntsrc contient la liste des sources de depart
print ntsrc
# Attention, ce qui est reconstruit n'est pas directement l'angle, mais
# sin(angle)*frequence - pour une matrice 512lignesx256colonnes
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'
addtitle 'PlanSource (Y=freq,X=AngNS)' ' ' 'font=times,bold,18 red'
w2eps recbeam.eps


# 3/ Faire tourner le programme - reconstruction cube X-Y-freq (en dehors de piapp)
Objs/treccyl XY12 n12cyl.ppf
OU 
Objs/treccyl XY5 n5cyl.ppf

# 4/ Affichage resultats
delobjs *
openppf n12cyl.ppf
OU 
openppf n5cyl.ppf

listobjs

## Pour afficher les plans reconstruits (angX-freq) sur les cylindres 0 , 3 avec les sources
## ATTENTION, si on change les parametres, la conversion angle->num pixel doit etre adapte
disp recsrcplane0 'lut=sqrt,0.,0.1 invy grey32 win'
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'
disp recsrcplane3 'lut=sqrt,0.,0.1 invy grey32 win'
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 freqred*512*2 amp 1 'nsta same red marker=circle,7'

## Pour afficher trois plans (angX-angY) avec les sources, autour des trois frequences des sources 
zone 1 3
disp recXYf0 'zoomx2 invy grey32 lut=sqrt,0.,1.5'
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 600.*angyrad/M_PI+10  amp fabs(freqred-0.1)<0.01 'nsta same red marker=circle,7'
disp recXYf1 'zoomx2 invy grey32 lut=sqrt,0.,1.5'
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 600.*angyrad/M_PI+10  amp fabs(freqred-0.27)<0.01 'nsta same red marker=circle,7'
disp recXYf2 'zoomx2 invy grey32 lut=sqrt,0.,1.5'
plot2dw ntsrc sin(angrad)*freq*0.25*256+128 600.*angyrad/M_PI+10  amp fabs(freqred-0.38)<0.01 'nsta same red marker=circle,7'

## On peut aussi extraire une coupe angX-angY de la boite 3D reconstruit
objaoper recsrcbox slicexy 101
set cobj /home/slicexy_101
disp $cobj 'zoomx2 grey32 lut=sqrt,0.,1'