eso-midas-testdata 17.02pl1.2-2build1 / usr / lib / eso-midas / 17FEB / test / prim / verify16.prg

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  MIDAS procedure verify16.prg  to verify MIDAS commands
!  for big data objects
! 
!  K. Banse     090429, 091228, 100906, 111205
!
!  use as @@ verify16 ffffffffffffffffffff         with f = 1 or 0 (on/off)
!
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
define/par p1 111111111111111 n "Enter control flags for entries: "
!
define/local loop/i/1/1 0
define/local rval/r/1/1 0. ? +lower
define/local ival/i/1/5 0 all +lower
define/local seconds/i/1/2 0,0 ? +lower
define/local rcorr/r/1/20 0. all +lower
define/local icorr/i/1/20 0 all +lower
define/local errsum/i/1/1 0 ? +lower
define/local ccc/c/1/15 000000000000000
define/local secs/i/1/2 0,0 ? +lower
define/local myvals/i/1/2 0,0 ? +lower
define/local longdir1/c/1/60 verylongdirectoryforfilestotestmidas ? +lower
define/local longdir2/c/1/60 anotherverylongdirectoryforfilestotestmidas ? +lower
!
delete/temp                             !get rid of old temporary files
! 
write/key sizez/i/1/5 600,600,0,0,0
write/key dispyes/i/1/2 0,0
!
write/out +------------------------------------------+
write/out Start of procedure verify16.prg
write/out +------------------------------------------+
!
seconds(1) = m$secs()
write/key ccc {p1}
set/format i1
do loop = 1 14
   if ccc({loop}:{loop}) .eq. "1" @@ verify16,000{loop}
   if errsum .gt. 0 then
      write/out "We got problems with entry 000{loop} in verify16.prg!"
      return 1
   endif
enddo
seconds(2) = m$secs()
mid$info(8) = seconds(2)-seconds(1)
! 
!
write/out +------------------------------------------+
write/out procedure verify16.prg:
write/out Total time elapsed = {mid$info(8)} seconds.
write/out All tests o.k. - you deserve a coffee now...
write/out +------------------------------------------+
return 0
!
!  here the different sub-procedures
!
entry 0001
! 
write/out test of CREATE/GRAPHICS + CREATE/DISPLAY 
write/out "----------------------------------------"
! 
reset/display
create/gra 3 600,400,0,380
create/display 7 512,512,616,300
dispyes(1) = 1                          !mark that we have display + graphic
dispyes(2) = 1 
load/lut heat
display/lut
clear/chan overlay
! 
entry 0002
!
write/out test of averaging images
write/out "------------------------"
! 
write/out > create/image a1 2,256,256 ? ? 200.
create/image a1 2,256,256 ? ? 200.
write/out > create/image a2 2,256,256 ? ? 400.
create/image a2 2,256,256 ? ? 400.
write/out > create/image a3 2,256,256 ? ? 100.
create/image a3 2,256,256 ? ? 100.
write/out > create/image a4 2,256,256 ? ? 500.
create/image a4 2,256,256 ? ? 500.
!
write/out > write/descr a1 weight/r/1/1 2.
write/descr a1 weight/r/1/1 2.
write/out > write/descr a2 weight/r/1/1 1.
write/descr a2 weight/r/1/1 1.
write/out > write/descr a3 weight/r/1/1 3.
write/descr a3 weight/r/1/1 3.
! for a4 weight is defaulted to 1.0
!
write/out > average/weights out1 = a1,a2,a3,a4
average/weights out1 = a1,a2,a3,a4
read/image out1
write/out > average/weights out2 = a1,a2,a3,a4 m
average/weights out2 = a1,a2,a3,a4 m
read/image out2
if {out1[@100,@100]} .ne. 228.571 then
   write/out error1
   errsum = errsum + 1
   return
endif
if {out1[@100,@100]} .ne. {out2[@100,@100]} then
   write/out error2
   errsum = errsum + 1
   return
endif
!
write/out > write/descr a2 weight/r/1/1 0.
write/descr a2 weight/r/1/1 0.
write/out > average/weights out1 = a1,a2,a3,a4
average/weights out1 = a1,a2,a3,a4
read/image out1
if {out1[@100,@100]} .ne. 200.0 then
   write/out error3
   errsum = errsum + 1
   return
endif
write/out > write/descr a2 weight/r/1/1 1.5
write/descr a2 weight/r/1/1 1.5
write/out > average/weights out1 = a1,a2,a3,a4
average/weights out1 = a1,a2,a3,a4
read/image out1
if {out1[@100,@100]} .ne. 240.0 then
   write/out error4
   errsum = errsum + 1
   return
endif
write/out > average/images out2 = a1,a2,a3,a4
average/images out2 = a1,a2,a3,a4
read/image out2
if {out2[@100,@100]} .ne. 300.0 then
   write/out error5
   errsum = errsum + 1
   return
endif
write/out > write/descr a1 weight/r/1/1 1.
write/descr a1 weight/r/1/1 1.
write/out > write/descr a2 weight/r/1/1 1.
write/descr a2 weight/r/1/1 1.
write/out > write/descr a3 weight/r/1/1 1.
write/descr a3 weight/r/1/1 1.
write/out > write/descr a3 weight/r/1/1 1.
write/descr a3 weight/r/1/1 1.
write/out > average/weights out1 = a1,a2,a3,a4
average/weights out1 = a1,a2,a3,a4
read/image out1
if {out1[@100,@100]} .ne. {out1[@100,@100]} then
   write/out error6
   errsum = errsum + 1
   return
endif
!
entry 0003
! 
write/out 
write/out test of extending descriptor directory
write/out "--------------------------------------"
! 
define/local nn/i/1/2 0,0
write/out > create/image muchos_dsc 2,1000,1000 ? poly 0.0,0.001,1.0
create/image muchos_dsc 2,1000,1000 ? poly 0.0,0.001,1.0
! 
if mid$sys(1:3) .eq. "MAC" then		!different ls syntax ...
   inputi(20) = 7856
   $ export BLOCKSIZE=512; $ ls -ls muchos_dsc.bdf  | write/keyw out_a
   inputi = {out_a(1:5)}
else
   inputi(20) = 7854
   $ ls -l --block-size=512 muchos_dsc.bdf  | write/keyw out_a
   nn = m$index(out_a," ")
   nn = nn + 1
   out_b = out_a({nn}:)
   nn = m$index(out_b," ")
   nn = nn + 1
   out_a = out_b({nn}:)
   nn = m$index(out_a," ")
   nn = nn + 1
   out_b = out_a({nn}:)
   nn = m$index(out_b," ")
   nn = nn + 1
   nn(2) = nn + 4
   inputi = {out_b({nn}:{nn(2)})}
endif
if inputi .ne. inputi(20) then
   write/out "initial no. blocks should be {inputi(20)} => error"
   errsum = errsum + 1
   return
endif
! 
! write/out > info/frame muchos_dsc debug
! info/frame muchos_dsc debug
! 
write/descr  muchos_dsc histogram/i/1/103 1 all
define/local auxkey/c/1/500 " " all
do nn = 1 20
   write/keyw auxkey "{nn},{nn},"  all
   write/descr muchos_dsc stressdsc{nn}/c/1/{nn} {auxkey}
enddo
delete/descr  muchos_dsc stressdsc0008
do nn = 21 109
   nn(2) = nn - 10
   write/dhelp muchos_dsc stressdsc{nn(2)} "help for descr. stressdsc{nn(2)}"
   write/keyw auxkey "{nn},{nn},"  all
   write/descr muchos_dsc stressdsc{nn}/c/1/{nn} {auxkey}
enddo
write/descr  muchos_dsc stressdsc0110/i/1/120 0 all
write/descr  muchos_dsc stressdsc0109/c/110/150 x all
! 
do nn = 111 500
   nn(2) = nn - 10
   write/dhelp muchos_dsc stressdsc{nn(2)} "help for descr. stressdsc{nn(2)}"
   write/keyw auxkey "{nn},{nn},"  all
   write/descr muchos_dsc stressdsc{nn}/c/1/{nn} {auxkey}
enddo
read/descr muchos_dsc stressdsc000*
! 
write/out > statistics/image muchos_dsc
statistics/image muchos_dsc
read/descr  muchos_dsc histogram f
! 
! write/out > info/frame muchos_dsc debug
! info/frame muchos_dsc debug
write/descr muchos_dsc stressdsc0333/c/330/4 AbcD
read/descr muchos_dsc stressdsc0333 f
read/descr  muchos_dsc stressdsc0109  f
! 
if mid$sys(1:3) .eq. "MAC" then
   inputi(20) = 8192
   $ export BLOCKSIZE=512; $ ls -ls muchos_dsc.bdf  | write/keyw out_a
   inputi = {out_a(1:5)}
else
   inputi(20) = 8190
   $ ls -l --block-size=512 muchos_dsc.bdf  | write/keyw out_a
   nn = m$index(out_a," ")
   nn = nn + 1
   out_b = out_a({nn}:)
   nn = m$index(out_b," ")
   nn = nn + 1
   out_a = out_b({nn}:)
   nn = m$index(out_a," ")
   nn = nn + 1
      out_b = out_a({nn}:)
   nn = m$index(out_b," ")
   nn = nn + 1
   nn(2) = nn + 4
   inputi = {out_b({nn}:{nn(2)})}
endif
if inputi .ne. inputi(20) then
   write/out "image muchos_dsc.bdf has size = {inputi} blocks"
   write/out "final no. blocks should be {inputi(20)} => error"
   errsum = errsum + 1
   return
endif
! 
entry 0004
!
write/out 
write/out test of creating image from ASCII file
write/out "--------------------------------------"
! 
write/out > indisk/ascii matlab.dat matlab 2,500,500
indisk/ascii matlab.dat matlab 2,500,500
write/out > read/descr matlab
read/descr matlab
write/out > statistics/image matlab
statistics/image matlab
write/key rcorr/r/1/12  0.00129361,394.247,1.08023,6.88132,29.1115,1107.43,-
270058.0,0.547258,0.774324,256,1.54606,0.774324
@@ kcompare rcorr outputr 1,1 0.000001 2,2 0.01 3,6 0.005
@@ kcompare rcorr outputr 7,7 0.5 8,9 0.0001 11,12 0.0001
!
entry 0005
!
return			!this takes too long on systems with less memory...
!
write/out 
write/out test of creating big images 
write/out "---------------------------"
! 
$rm -f xbig.big
! 
if mid$sys(31:32) .eq. "64" then
   define/local inni/i/1/1 500
else
   define/local inni/i/1/1 250
endif
! 
write/keyw bigno/s/1/2 1234,1234
bigno(1) = 1000000 * inni
bigno(2) = bigno(1) * 4
set/format i1
write/out > create/image xbig.big 2,1000000,{inni} ? poly 27,1.,1.
create/image xbig.big 2,1000000,{inni} ? poly 27,1.,1. 
$ls -hl xbig.big
write/out > read/descr xbig.big
read/descr xbig.big
write/out > read/image xbig.big @100000,@100,30
read/image xbig.big @100000,@100,30
write/out > statistics/image xbig.big
statistics/image xbig.big
! $rm -f xbig.big
!
entry 0006
!
write/out 
write/out test of TBERDC (center/gauss)
write/out "---------------------------"
! 
if dispyes(1) .eq. 1 reset/display
! 
! the following test procedure was kindly contributed by Cees Bassa
! from the Jodrell Bank Center for Astrophysics, Manchester (July 2010)
! 
! Read image and reset wcs
indisk/fits ISAAC.2006-04-13T06:32:38.944.fits isaac
write/desc isaac start/d/1/2 1.,1.
write/desc isaac step/d/1/2 1.,1.
write/desc isaac cunit/c/1/60 " "

! Create display and load image
if dispyes(1) .eq. 1 then
   create/disp 0 520,520
   load/lut heat
   load isaac scale=-2 cuts=F,1sigma
endif

! Load astrometric standards and position estimates
$echo "define/field 01 12 R F12.6 :xcen" >in.fmt
$echo "define/field 13 25 R F12.6 :ycen" >>in.fmt
$echo "define/field 26 34 R F12.6 :rx" >>in.fmt
$echo "define/field 35 43 R F12.6 :ry" >>in.fmt
$echo "end" >>in.fmt

create/tab astrometry * * astrometry.dat in.fmt

! Load standards
if dispyes(1) .eq. 1 then
   load/tab astrometry :xcen :ycen ? 1 4 4
endif

! Create IDENT column
crea/col astrometry :no I*4 I4
comp/tab astrometry :no = seq
crea/col astrometry :ident C*8 A8
comp/tab astrometry :ident = CONCAT("ID",TOCHAR(:no))

! Create center/gauss search boxes
comp/tab astrometry :xstart = :xcen-10.
comp/tab astrometry :xend = :xcen+10.
comp/tab astrometry :ystart = :ycen-10.
comp/tab astrometry :yend = :ycen+10.

! Run center/gauss
center/gauss isaac,astrometry astrometry

! Offset positions to image center
comp/tab astrometry :x_off = :xcen-512.0
comp/tab astrometry :y_off = :ycen-512.0

! Copy astrometry table because align/ima wants them different
copy/tab astrometry temp

! Load fitted positions
if dispyes(1) .eq. 1 then
   load/tab astrometry :xcen :ycen ? 1 2 5
endif

! Align image
align/ima astrometry,:x_off,:y_off temp,:rx,:ry f ? y

! Plot residuals
if dispyes(1) .eq. 1 then
   create/gra
   plot/tab astrometry :xresidual :yresidual ? 1
endif

! Residual statistics
stat/tab astrometry :xresidual
write/key rcorr/r/1/4 -1.5899,0.467388,5.19838e-07,0.250642  
@@ kcompare rcorr outputr 1,2 0.0001  4,4 0.0001
! 
stat/tab astrometry :yresidual
write/key rcorr/r/1/4  -0.986995,0.577721,-4.47603e-07,0.223052 
@@ kcompare rcorr outputr 1,2 0.0001  4,4 0.0001
!
entry 0007
!
write/out 
write/out test of COPY/IMT3, COPY/T3IMA 
write/out "-----------------------------"
! 
extract/image line_ima = rosat.mt[<,@100:>,@100]
find/minmax line_ima
if dispyes(1) .eq. 1 then
   dele/displ all
   create/display
   load/lut rainbow
   clear/chan ov
   create/graph
   set/graph color=blue
   plot line_ima
endif
! 
write/out > copy/imt3 line_ima tab3d :cola,1 ? 104
copy/imt3 line_ima tab3d :cola,1 ? 104
show/table tab3d
write/out > copy/t3ima tab3d :cola 1 ? -11475.5,90.0 line_tab
copy/t3ima tab3d :cola 1 ? -11475.5,90.0 line_tab
compute/image &d = line_ima - line_tab
find/minmax &d >Null
write/key rcorr/r/1/2  0.,0.
@@ kcompare rcorr outputr 1,2 0.00001  
! 
extract/image sub_ima = rosat.mt[<,@81:>,@180]
if dispyes(1) .eq. 1 then
   load/imag sub_ima
   set/graph color=read
   overplot line_tab
endif
write/out > copy/imt3 sub_ima tab3d :colb,@2 ? update
copy/imt3 sub_ima tab3d :colb,@2 ? update
show/table tab3d
write/key icorr/i/1/8 2,101,0,0,513,104,0,101
@@ kcompare icorr outputi 1,8
! 
write/out > copy/t3ima tab3d :colb 2,100 ? -11475.5,90.0,4274.5,-90.0 sub_tab
copy/t3ima tab3d :colb 2,100 ? -11475.5,90.0,4274.5,-90.0 sub_tab
compute/image &d = sub_ima - sub_tab
find/minmax &d >Null
write/key rcorr/r/1/2  0.,0.
@@ kcompare rcorr outputr 1,2 0.00001  
!
crea/ima rk1 2,100,20 ? ? 0.0,1.0
copy/ii rk1 ik2 i4
read/image ik2
write/out > copy/imt3 ik2 lola
copy/imt3 ik2 lola
write/out > show/tab lola
show/tab lola
write/out > copy/t3ima lola #1 1,20 ? ? ik3
copy/t3ima lola #1 1,20 ? ? ik3
write/out > compute/pix &a = ik3 - ik2 
compute/pix &a = ik3 - ik2 
find/min &a
write/key rcorr/r/1/2  0.,0.
@@ kcompare rcorr outputr 1,2 0.00001
! 
write/out > comp/ima rk4 = 99 - ik2
comp/ima rk4 = 99 - ik2
write/out > copy/imt3 ik2 lola :reverse,@1 ? upda
copy/imt3 ik2 lola :reverse,@1 ? upda
show/table lola
write/out > copy/imt3 rk4 lola :reverse,@1 ? upda
copy/imt3 rk4 lola :reverse,@1 ? upda
show/table lola
! 
entry 0008
! 
write/out 
write/out more tests of 3dim tables
write/out "-------------------------"
!
!convert all extensions to Midas tables
write/out > indisk/fits xamber.fits ROOT=amber
indisk/fits xamber.fits ROOT=amber
! 
!the VIS2DATA is extension #4, so it goes to amber0004.tbl
!copy this table and work on table work.tbl,then
! 
write/out > "$ cp amber0004.tbl work.tbl"
$ cp amber0004.tbl work.tbl
! 
!3 systems are stored in that table, using 3 rows per given time (:time)
!let's work on the systems individually, 
!so we need rows 1,4,... or 2,5,... or 3,6,...
! 
write/out > create/colu work :isi I*4
create/colu work :isi I*4
write/out > compute/table work :isi = MOD(seq,3)
compute/table work :isi = MOD(seq,3)		!holds 0,1,2,0,1,2,...
! 
define/local counter/i/1/2 0,0
set/format i1
! 
!loop through different systems
do counter = 1 3
   inputi = counter
   if counter .eq. 3 inputi = 0
   write/out > select/table work :isi .eq. {inputi}	
   select/table work :isi .eq. {inputi}	
 
   !extract the spectra from the selected rows to get a 510x15 2-dim file 
   write/out > copy/t3ima work :vis2data 1,45 ? ? outi
   copy/t3ima work :vis2data 1,45 ? ? outi >Null

   !get statistics of all lines, and put into table vis1, vis2, vis3
   write/out > statistics/image outi row ? ? R vis{counter} 
   statistics/image outi row ? ? R vis{counter} >Null
   select/table work all
enddo
!
show/table vis2
write/out > statistics/table vis2 :mean
statistics/table vis2 :mean
write/key rcorr   0.0694125,0.886639,0.776187,0.208711 
@@ kcompare rcorr outputr 1,4 0.00005
!  
if dispyes(1) .ne. 1 return
!
!now, plot the stand. dev. of the visibility data of all systems
delete/graph
set/gra
create/graph 0 800,400
set/grap ltype=1 stype=0
plot/table vis2 ? :std
set/gra color=blue
overplot/table vis1 ? :std
set/gra color=red
overplot/table vis3 ? :std
! 
entry 0009
!
entry 00010
!
entry 00011
!
entry 00012
!
! 
entry 00013
! 
entry 00014
!