diff --git a/CEP/CPA/PSS3/CAL/src/demo.g b/CEP/CPA/PSS3/CAL/src/demo.g index b1fdfaf627230f98814b3fa47999ae344acd0ede..11013dca5bd70c9b57cb8b75573b9c2867e0b4a7 100644 --- a/CEP/CPA/PSS3/CAL/src/demo.g +++ b/CEP/CPA/PSS3/CAL/src/demo.g @@ -1,20 +1,20 @@ -include 'table.g' -include 'meqcalibrater.g' -include 'parmtable.g' -include 'gsm.g' -include 'imgannotator.g' -include 'mkimg.g' +include 'table.g'; +include 'meqcalibrater.g'; +include 'parmtable.g'; +include 'gsm.g'; +include 'imgannotator.g'; +include 'mkimg.g'; # # Demo function showing the predict functionality and creating an image of it. # predict := function(fname='demo', ant=4*[0:20], - modeltype='LOFAR', calcuvw=F, trace=T) + modeltype='LOFAR', calcuvw=F, trace=T) { local mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant, modeltype=modeltype, calcuvw=calcuvw); + ant=ant, modeltype=modeltype, calcuvw=calcuvw); if (is_fail(mc)) { print "meqcalibratertest(): could not instantiate meqcalibrater"; fail; @@ -23,7 +23,7 @@ predict := function(fname='demo', ant=4*[0:20], mc.settimeinterval(3600); # calibrate per 1 hour mc.clearsolvableparms(); - mc.resetiterator() + mc.resetiterator(); while (mc.nextinterval()) { d := mc.getsolvedomain(); @@ -32,23 +32,23 @@ predict := function(fname='demo', ant=4*[0:20], mc.predict('MODEL_DATA'); # mc.saveresidualdata(); } - + mc.done(); - + mssel := ''; if (len(ant) > 0) { - ant +:=1; # msselect adds 1 to ANTENNA1,2 - mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); + ant +:=1; # msselect adds 1 to ANTENNA1,2 + mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } - print mssel + print mssel; mkimg(spaste(fname, '.MS'), spaste(fname, '.img'), msselect=mssel); - + return T; } solve := function(fname='demo', ant=4*[0:20], - modeltype='LOFAR', calcuvw=F, - niter=1, sleep=F, sleeptime=2, wait=F) + modeltype='LOFAR', calcuvw=F, + niter=1, sleep=F, sleeptime=2, wait=F) { annotator := imgannotator(spaste(fname, '.img'), 'raster'); @@ -56,20 +56,20 @@ solve := function(fname='demo', ant=4*[0:20], # First get the values used in the simulation. t := table(spaste(fname,'_gsm.MEP/DEFAULTVALUES')); t1 := t.query('NAME==pattern("{RA,DEC,StokesI}.*")', - sortlist='SRCNR,NAME'); + sortlist='SRCNR,NAME'); t.close(); if (t1.nrows() == 0) fail "No sources in MEP table"; if (t1.nrows()%3 != 0) fail "_gsm.MEP inconsistent for RA,DEC,StokesI"; # MEP is sorted, thus order is DEC,RA,StokesI annotator.hold(); for (i in [1:(t1.nrows()/3)]) { - dec := t1.getcell("SIM_VALUES", 3*i-2)[1,1]; - ra := t1.getcell("SIM_VALUES", 3*i-1)[1,1]; - stk := t1.getcell("SIM_VALUES", 3*i)[1,1]; - print stk; - src_mrk[i] := annotator.add_marker(1, ra, dec, F, F, 1, - as_integer(100*stk), 3); - annotator.add_marker(2, ra, dec, F, F, 1, as_integer(100*stk), 12); + dec := t1.getcell("SIM_VALUES", 3*i-2)[1,1]; + ra := t1.getcell("SIM_VALUES", 3*i-1)[1,1]; + stk := t1.getcell("SIM_VALUES", 3*i)[1,1]; + print stk; + src_mrk[i] := annotator.add_marker(1, ra, dec, F, F, 1, + as_integer(100*stk), 3); + annotator.add_marker(2, ra, dec, F, F, 1, as_integer(100*stk), 12); } annotator.release(); t1.close(); @@ -78,13 +78,13 @@ solve := function(fname='demo', ant=4*[0:20], # Create calibrater object # mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant, - modeltype=modeltype, calcuvw=calcuvw); + ant=ant, + modeltype=modeltype, calcuvw=calcuvw); if (wait) { - print "Press RETURN to continue."; - shell("read"); + print "Press RETURN to continue."; + shell("read"); } if (is_fail(mc)) { @@ -111,91 +111,91 @@ solve := function(fname='demo', ant=4*[0:20], if (sleep) shell(sleep_cmd); mc.resetiterator() - while (mc.nextinterval()) - { - d := mc.getsolvedomain(); - print 'solvedomain = ', d; - - parms := mc.getparms("RA.* DEC.* StokesI.*"); - print parms - print len(parms) - nrpos := len(parms) / 3; - if (nrpos > 0) { - annotator.hold(); - for (i in [1:nrpos]) { - ra := parms[spaste('RA.CP',i)].value[1]; - dec := parms[spaste('DEC.CP',i)].value[1]; - stokesI := parms[spaste('StokesI.CP',i)].value[1]; - print 'src = ', i, ' ra = ', ra, ' dec = ', dec, - ' I = ', stokesI; - - annotator.change_marker_size(src_mrk[i], stokesI*100); - if (is_fail(annotator)) fail; - annotator.add_marker(i, real(ra), real(dec), i==nrpos); - } - annotator.release(); - } - - for (i in [1:niter]) { - print "iteration", i; + while (mc.nextinterval()) + { + d := mc.getsolvedomain(); + print 'solvedomain = ', d; - srec := mc.solve(); - solverec[spaste("iter",i)] := srec; - parms := mc.getparms("RA.* DEC.* StokesI.*"); + print parms; + print len(parms); nrpos := len(parms) / 3; if (nrpos > 0) { - annotator.hold(); - for (j in [1:nrpos]) { - ra := parms[spaste('RA.CP',j)].value[1]; - dec := parms[spaste('DEC.CP',j)].value[1]; - stokesI := parms[spaste('StokesI.CP',j)].value[1]; - print 'src = ', j, ' ra = ', ra, ' dec = ', dec, - ' I = ', stokesI; + annotator.hold(); + for (i in [1:nrpos]) { + ra := parms[spaste('RA.CP',i)].value[1]; + dec := parms[spaste('DEC.CP',i)].value[1]; + stokesI := parms[spaste('StokesI.CP',i)].value[1]; + print 'src = ', i, ' ra = ', ra, ' dec = ', dec, + ' I = ', stokesI; - annotator.change_marker_size(src_mrk[j], stokesI*100); + annotator.change_marker_size(src_mrk[i], stokesI*100); if (is_fail(annotator)) fail; - annotator.add_marker(j, real(ra), real(dec), j==nrpos); + annotator.add_marker(i, real(ra), real(dec), i==nrpos); + } + annotator.release(); + } + + for (i in [1:niter]) { + print "iteration", i; + + srec := mc.solve(); + solverec[spaste("iter",i)] := srec; + + parms := mc.getparms("RA.* DEC.* StokesI.*"); + nrpos := len(parms) / 3; + if (nrpos > 0) { + annotator.hold(); + for (j in [1:nrpos]) { + ra := parms[spaste('RA.CP',j)].value[1]; + dec := parms[spaste('DEC.CP',j)].value[1]; + stokesI := parms[spaste('StokesI.CP',j)].value[1]; + print 'src = ', j, ' ra = ', ra, ' dec = ', dec, + ' I = ', stokesI; + + annotator.change_marker_size(src_mrk[j], stokesI*100); + if (is_fail(annotator)) fail; + annotator.add_marker(j, real(ra), real(dec), j==nrpos); + } + annotator.release(); } - annotator.release(); + sleep_cmd := spaste('sleep ', sleeptime); + if (sleep) shell(sleep_cmd); } - sleep_cmd := spaste('sleep ', sleeptime); - if (sleep) shell(sleep_cmd); + print mc.getstatistics(); + mc.saveresidualdata(); + mc.saveparms(); } - print mc.getstatistics() - mc.saveresidualdata(); - mc.saveparms(); - } mc.done(); mssel := ''; if (len(ant) > 0) { - ant +:=1; # msselect adds 1 to ANTENNA1,2 - mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); + ant +:=1; # msselect adds 1 to ANTENNA1,2 + mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } mkimg(spaste(fname, '.MS'), spaste(fname, '.imgs', src[1]), - msselect=mssel, type='corrected'); + msselect=mssel, type='corrected'); return solverec; # return ref annotator; } solveej := function(fname='demo', ant=4*[0:20], - modeltype='LOFAR', calcuvw=F, - niter=1, sleep=F, sleeptime=2, wait=F) + modeltype='LOFAR', calcuvw=F, + niter=1, sleep=F, sleeptime=2, wait=F) { # # Create calibrater object # mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant, - modeltype=modeltype, calcuvw=calcuvw); + ant=ant, + modeltype=modeltype, calcuvw=calcuvw); if (wait) { - print "Press RETURN to continue."; - shell("read"); + print "Press RETURN to continue."; + shell("read"); } if (is_fail(mc)) { @@ -215,59 +215,59 @@ solveej := function(fname='demo', ant=4*[0:20], if (sleep) shell(sleep_cmd); mc.resetiterator() - while (mc.nextinterval()) - { - d := mc.getsolvedomain(); - print 'solvedomain = ', d; - - parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); - nrp := len(parms); - if (nrp > 0) { - for (nm in field_names(parms)) { - print nm, '=', parms[nm].value; - } - } - - for (i in [1:niter]) { - print "iteration", i; - - srec := mc.solve(); - solverec[spaste("iter",i)] := srec; + while (mc.nextinterval()) + { + d := mc.getsolvedomain(); + print 'solvedomain = ', d; + + parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); + nrp := len(parms); + if (nrp > 0) { + for (nm in field_names(parms)) { + print nm, '=', parms[nm].value; + } + } - parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); - nrp := len(parms); - if (nrp > 0) { - for (nm in field_names(parms)) { - print nm, '=', parms[nm].value; - } - } - - sleep_cmd := spaste('sleep ', sleeptime); - if (sleep) shell(sleep_cmd); + for (i in [1:niter]) { + print "iteration", i; + + srec := mc.solve(); + solverec[spaste("iter",i)] := srec; + + parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); + nrp := len(parms); + if (nrp > 0) { + for (nm in field_names(parms)) { + print nm, '=', parms[nm].value; + } + } + + sleep_cmd := spaste('sleep ', sleeptime); + if (sleep) shell(sleep_cmd); + } + print mc.getstatistics(); + mc.saveresidualdata(); + mc.saveparms(); } - print mc.getstatistics() - mc.saveresidualdata(); - mc.saveparms(); - } mc.done(); mssel := ''; if (len(ant) > 0) { - ant +:=1; # msselect adds 1 to ANTENNA1,2 - mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); + ant +:=1; # msselect adds 1 to ANTENNA1,2 + mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } mkimg(spaste(fname, '.MS'), spaste(fname, '.imgs', src[1]), - msselect=mssel, type='corrected'); + msselect=mssel, type='corrected'); return solverec; # return ref annotator; } peel := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], - datacol='MODEL_DATA', mapnr='', - modeltype='LOFAR', calcuvw=F, - niter=1, sleep=F, sleeptime=2, wait=F) + datacol='MODEL_DATA', mapnr='', + modeltype='LOFAR', calcuvw=F, + niter=1, sleep=F, sleeptime=2, wait=F) { annotator := imgannotator(spaste(fname, '.img', mapnr), 'raster'); @@ -275,20 +275,20 @@ peel := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], # First get the values used in the simulation. t := table(spaste(fname,'_gsm.MEP/DEFAULTVALUES')); t1 := t.query('NAME==pattern("{RA,DEC,StokesI}.*")', - sortlist='SRCNR,NAME'); + sortlist='SRCNR,NAME'); t.close(); if (t1.nrows() == 0) fail "No sources in MEP table"; if (t1.nrows()%3 != 0) fail "_gsm.MEP inconsistent for RA,DEC,StokesI"; # MEP is sorted, thus order is DEC,RA,StokesI annotator.hold(); for (i in [1:(t1.nrows()/3)]) { - dec := t1.getcell("SIM_VALUES", 3*i-2)[1,1]; - ra := t1.getcell("SIM_VALUES", 3*i-1)[1,1]; - stk := t1.getcell("SIM_VALUES", 3*i)[1,1]; - print stk; - src_mrk[i] := annotator.add_marker(1, ra, dec, F, F, 1, - as_integer(100*stk), 3); - annotator.add_marker(2, ra, dec, F, F, 1, as_integer(100*stk), 12); + dec := t1.getcell("SIM_VALUES", 3*i-2)[1,1]; + ra := t1.getcell("SIM_VALUES", 3*i-1)[1,1]; + stk := t1.getcell("SIM_VALUES", 3*i)[1,1]; + print stk; + src_mrk[i] := annotator.add_marker(1, ra, dec, F, F, 1, + as_integer(100*stk), 3); + annotator.add_marker(2, ra, dec, F, F, 1, as_integer(100*stk), 12); } annotator.release(); t1.close(); @@ -297,13 +297,13 @@ peel := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], # Create calibrater object # mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant, datacolname=datacol, - modeltype=modeltype, calcuvw=calcuvw); + ant=ant, datacolname=datacol, + modeltype=modeltype, calcuvw=calcuvw); if (wait) { - print "Press RETURN to continue."; - shell("read"); + print "Press RETURN to continue."; + shell("read"); } if (is_fail(mc)) { @@ -316,9 +316,9 @@ peel := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], mc.settimeinterval(3600); # calibrate per 3600 seconds mc.clearsolvableparms(); for (i in src) { - mc.setsolvableparms(spaste('StokesI.CP', i)); - mc.setsolvableparms(spaste('RA.CP', i)); - mc.setsolvableparms(spaste('DEC.CP', i)); + mc.setsolvableparms(spaste('StokesI.CP', i)); + mc.setsolvableparms(spaste('RA.CP', i)); + mc.setsolvableparms(spaste('DEC.CP', i)); } solverec := [=]; @@ -327,92 +327,92 @@ peel := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], if (sleep) shell(sleep_cmd); mc.resetiterator() - while (mc.nextinterval()) - { - mc.peel (src-1, predsrc-1); - - d := mc.getsolvedomain(); - print 'solvedomain = ', d; + while (mc.nextinterval()) + { + mc.peel (src-1, predsrc-1); - parms := mc.getparms("RA.* DEC.* StokesI.*"); - print parms - print len(parms) - nrpos := len(parms) / 3; - if (nrpos > 0) { - annotator.hold(); - for (i in [1:nrpos]) { - ra := parms[spaste('RA.CP',i)].value[1]; - dec := parms[spaste('DEC.CP',i)].value[1]; - stokesI := parms[spaste('StokesI.CP',i)].value[1]; - print 'src = ', i, ' ra = ', ra, ' dec = ', dec, - ' I = ', stokesI; - - annotator.change_marker_size(src_mrk[i], stokesI*100); - if (is_fail(annotator)) fail; - annotator.add_marker(i, real(ra), real(dec), i==nrpos); - } - annotator.release(); - } - - for (i in [1:niter]) { - print "iteration", i; + d := mc.getsolvedomain(); + print 'solvedomain = ', d; - srec := mc.solve(); - solverec[spaste("iter",i)] := srec; - parms := mc.getparms("RA.* DEC.* StokesI.*"); - nrpos := len(parms) / 3; + print parms + print len(parms) + nrpos := len(parms) / 3; if (nrpos > 0) { - annotator.hold(); - for (j in [1:nrpos]) { - ra := parms[spaste('RA.CP',j)].value[1]; - dec := parms[spaste('DEC.CP',j)].value[1]; - stokesI := parms[spaste('StokesI.CP',j)].value[1]; - print 'src = ', j, ' ra = ', ra, ' dec = ', dec, - ' I = ', stokesI; - - annotator.change_marker_size(src_mrk[j], stokesI*100); + annotator.hold(); + for (i in [1:nrpos]) { + ra := parms[spaste('RA.CP',i)].value[1]; + dec := parms[spaste('DEC.CP',i)].value[1]; + stokesI := parms[spaste('StokesI.CP',i)].value[1]; + print 'src = ', i, ' ra = ', ra, ' dec = ', dec, + ' I = ', stokesI; + + annotator.change_marker_size(src_mrk[i], stokesI*100); if (is_fail(annotator)) fail; - annotator.add_marker(j, real(ra), real(dec), j==nrpos); + annotator.add_marker(i, real(ra), real(dec), i==nrpos); } - annotator.release(); + annotator.release(); } - sleep_cmd := spaste('sleep ', sleeptime); - if (sleep) shell(sleep_cmd); + + for (i in [1:niter]) { + print "iteration", i; + + srec := mc.solve(); + solverec[spaste("iter",i)] := srec; + + parms := mc.getparms("RA.* DEC.* StokesI.*"); + nrpos := len(parms) / 3; + if (nrpos > 0) { + annotator.hold(); + for (j in [1:nrpos]) { + ra := parms[spaste('RA.CP',j)].value[1]; + dec := parms[spaste('DEC.CP',j)].value[1]; + stokesI := parms[spaste('StokesI.CP',j)].value[1]; + print 'src = ', j, ' ra = ', ra, ' dec = ', dec, + ' I = ', stokesI; + + annotator.change_marker_size(src_mrk[j], stokesI*100); + if (is_fail(annotator)) fail; + annotator.add_marker(j, real(ra), real(dec), j==nrpos); + } + annotator.release(); + } + sleep_cmd := spaste('sleep ', sleeptime); + if (sleep) shell(sleep_cmd); + } + print mc.getstatistics(); + mc.saveresidualdata(); + mc.saveparms(); } - print mc.getstatistics(); - mc.saveresidualdata(); - mc.saveparms(); - } mc.done(); mssel := ''; if (len(ant) > 0) { - ant +:=1; # msselect adds 1 to ANTENNA1,2 - mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); + ant +:=1; # msselect adds 1 to ANTENNA1,2 + mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } mkimg(spaste(fname, '.MS'), spaste(fname, '.img', src[1]), - msselect=mssel, type='corrected'); + msselect=mssel, type='corrected'); return T; # return ref annotator; } peelej := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], - datacol='MODEL_DATA', mapnr='', - modeltype='LOFAR', calcuvw=F, - niter=1, sleep=F, sleeptime=2, wait=F) + datacol='MODEL_DATA', mapnr='', + modeltype='LOFAR', calcuvw=F, + niter=1, sleep=F, sleeptime=2, wait=F) { # # Create calibrater object # mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant, datacolname=datacol, - modeltype=modeltype, calcuvw=calcuvw); + ant=ant, datacolname=datacol, + modeltype=modeltype, calcuvw=calcuvw); if (wait) { - print "Press RETURN to continue."; - shell("read"); + print "Press RETURN to continue."; + shell("read"); } if (is_fail(mc)) { @@ -425,7 +425,7 @@ peelej := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], mc.settimeinterval(3600); # calibrate per 3600 seconds mc.clearsolvableparms(); for (i in src) { - mc.setsolvableparms(spaste('EJ11.*.CP', src, '*')); + mc.setsolvableparms(spaste('EJ11.*.CP', src, '*')); # mc.setsolvableparms(spaste('EJ11.SR{5,9,13,17,21,25,29,33,37}.CP', src)); } solverec := [=]; @@ -434,49 +434,49 @@ peelej := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], if (sleep) shell(sleep_cmd); mc.resetiterator() - while (mc.nextinterval()) - { - mc.peel (src-1, predsrc-1); + while (mc.nextinterval()) + { + mc.peel (src-1, predsrc-1); - d := mc.getsolvedomain(); - print 'solvedomain = ', d; + d := mc.getsolvedomain(); + print 'solvedomain = ', d; - parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); + parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); # parms := mc.getparms("EJ*"); - print parms - - for (i in [1:niter]) { - print "iteration", i; + print parms + + for (i in [1:niter]) { + print "iteration", i; - srec := mc.solve(); - solverec[spaste("iter",i)] := srec; - - for (j in src) { - parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); + srec := mc.solve(); + solverec[spaste("iter",i)] := srec; + + for (j in src) { + parms := mc.getparms("EJ11.{real,imag}.SR{5,9}.*"); # parms := mc.getparms(spaste('EJ*.CP', src, '*')); - nrp := len(parms); - if (nrp > 0) { - for (nm in field_names(parms)) { - print nm, '=', parms[nm].value; - } + nrp := len(parms); + if (nrp > 0) { + for (nm in field_names(parms)) { + print nm, '=', parms[nm].value; + } + } + } + sleep_cmd := spaste('sleep ', sleeptime); + if (sleep) shell(sleep_cmd); } - } - sleep_cmd := spaste('sleep ', sleeptime); - if (sleep) shell(sleep_cmd); + print mc.getstatistics(); + mc.saveresidualdata(); + mc.saveparms(); } - print mc.getstatistics(); - mc.saveresidualdata(); - mc.saveparms(); - } mc.done(); mssel := ''; if (len(ant) > 0) { - ant +:=1; # msselect adds 1 to ANTENNA1,2 - mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); + ant +:=1; # msselect adds 1 to ANTENNA1,2 + mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } mkimg(spaste(fname, '.MS'), spaste(fname, '.img', src[1]), - msselect=mssel, type='corrected'); + msselect=mssel, type='corrected'); return T; # return ref annotator; } @@ -484,9 +484,9 @@ peelej := function(fname='demo', src=1, predsrc=[], ant=4*[0:20], solvepos := function(fname='demo', ant=4*[0:20], niter=1) { annotator := imgannotator(spaste(fname, '.img'), 'raster'); - + mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant); + ant=ant); if (is_fail(mc)) { print "meqcalibratertest(): could not instantiate meqcalibrater"; fail; @@ -501,7 +501,7 @@ solvepos := function(fname='demo', ant=4*[0:20], niter=1) dec := parms[spaste('GSM.',i,'.DEC')].value[1]; print 'src = ', i, ' ra = ', ra, ' dec = ', dec; - if (is_fail(annotator)) fail; + if (is_fail(annotator)) fail; annotator.add_marker(i, real(ra), real(dec), i==nrpos); } } @@ -509,7 +509,7 @@ solvepos := function(fname='demo', ant=4*[0:20], niter=1) mc.settimeinterval(3600); # calibrate per 3600 seconds (1 timeslot = 2 sec) for (i in [1:niter]) { print "iteration", i - mc.clearsolvableparms(); + mc.clearsolvableparms(); print 'Solving for RA ...'; mc.setsolvableparms("GSM.*.RA"); @@ -551,11 +551,11 @@ solvepos := function(fname='demo', ant=4*[0:20], niter=1) solvegain := function(fname='demo', ant=4*[0:20], niter=1) { mc := meqcalibrater(spaste(fname,'.MS'), fname, spaste(fname,'_gsm'), - ant=ant); + ant=ant); if (is_fail(mc)) { - print "meqcalibratertest(): could not instantiate meqcalibrater" - fail - } + print "meqcalibratertest(): could not instantiate meqcalibrater" + fail + } # mc.select('', 0, 0); @@ -570,32 +570,32 @@ solvegain := function(fname='demo', ant=4*[0:20], niter=1) # print 'ALL PARMS = ', parms; mc.resetiterator() - while (mc.nextinterval()) - { - d := mc.getsolvedomain(); - print 'solvedomain = ', d; - - for (i in [1:niter]) - { - print "iteration", i - mc.solve(); - - parms := mc.getparms("gain.11.*"); - print 'GAIN SOLUTION = '; - for (k in [1:5]) - { - print parms[k].value[1]; - } - - parms := mc.getparms("gain.22.*"); - print 'GAIN SOLUTION = '; - for (k in [1:5]) - { - print parms[k].value[1]; - } - } - print mc.getstatistics() - } + while (mc.nextinterval()) + { + d := mc.getsolvedomain(); + print 'solvedomain = ', d; + + for (i in [1:niter]) + { + print "iteration", i + mc.solve(); + + parms := mc.getparms("gain.11.*"); + print 'GAIN SOLUTION = '; + for (k in [1:5]) + { + print parms[k].value[1]; + } + + parms := mc.getparms("gain.22.*"); + print 'GAIN SOLUTION = '; + for (k in [1:5]) + { + print parms[k].value[1]; + } + } + print mc.getstatistics() + } parms := mc.getparms("gain.11.*"); print 'SOLUTION 11: ', parms; @@ -607,90 +607,90 @@ solvegain := function(fname='demo', ant=4*[0:20], niter=1) initparms := function(fname='demo') { - pt := parmtable(spaste(fname,'.MEP'), T); - if (is_fail(pt)) fail; - pt.putinit ('frot', values=0); - pt.putinit ('drot', values=0); - pt.putinit ('dell', values=0); - pt.putinit ('gain.11', values=1); - pt.putinit ('gain.22', values=1); - pt.putinit ('EJ11.real.SR1.CP1', values=1); - pt.putinit ('EJ11.real.SR1.CP2', values=1); - pt.putinit ('EJ11.real.SR1.CP3', values=1); - pt.putinit ('EJ11.real.SR2.CP1', values=1); - pt.putinit ('EJ11.real.SR2.CP2', values=1); - pt.putinit ('EJ11.real.SR2.CP3', values=1); - pt.putinit ('EJ11.real.SR5.CP1', values=1); - pt.putinit ('EJ11.real.SR5.CP2', values=1); - pt.putinit ('EJ11.real.SR5.CP3', values=1); - pt.putinit ('EJ11.real.SR9.CP1', values=1); - pt.putinit ('EJ11.real.SR9.CP2', values=1); - pt.putinit ('EJ11.real.SR9.CP3', values=1); - pt.putinit ('EJ11.real', values=1); - pt.putinit ('EJ12.real', values=0); - pt.putinit ('EJ21.real', values=0); - pt.putinit ('EJ22.real', values=1); - pt.putinit ('EJ11.imag', values=1); - pt.putinit ('EJ12.imag', values=0); - pt.putinit ('EJ21.imag', values=0); - pt.putinit ('EJ22.imag', values=1); - pt.done(); + pt := parmtable(spaste(fname,'.MEP'), T); + if (is_fail(pt)) fail; + pt.putinit ('frot', values=0); + pt.putinit ('drot', values=0); + pt.putinit ('dell', values=0); + pt.putinit ('gain.11', values=1); + pt.putinit ('gain.22', values=1); + pt.putinit ('EJ11.real.SR1.CP1', values=1); + pt.putinit ('EJ11.real.SR1.CP2', values=1); + pt.putinit ('EJ11.real.SR1.CP3', values=1); + pt.putinit ('EJ11.real.SR2.CP1', values=1); + pt.putinit ('EJ11.real.SR2.CP2', values=1); + pt.putinit ('EJ11.real.SR2.CP3', values=1); + pt.putinit ('EJ11.real.SR5.CP1', values=1); + pt.putinit ('EJ11.real.SR5.CP2', values=1); + pt.putinit ('EJ11.real.SR5.CP3', values=1); + pt.putinit ('EJ11.real.SR9.CP1', values=1); + pt.putinit ('EJ11.real.SR9.CP2', values=1); + pt.putinit ('EJ11.real.SR9.CP3', values=1); + pt.putinit ('EJ11.real', values=1); + pt.putinit ('EJ12.real', values=0); + pt.putinit ('EJ21.real', values=0); + pt.putinit ('EJ22.real', values=1); + pt.putinit ('EJ11.imag', values=1); + pt.putinit ('EJ12.imag', values=0); + pt.putinit ('EJ21.imag', values=0); + pt.putinit ('EJ22.imag', values=1); + pt.done(); } setparms := function(fname='demo') { - pt := parmtable(spaste(fname,'.MEP'), T); - if (is_fail(pt)) fail; - pt.putinit ('frot', values=0); - pt.putinit ('drot', values=0); - pt.putinit ('dell', values=0); - pt.putinit ('gain.11', values=1); - pt.putinit ('gain.22', values=0); - pt.putinit ('EJ11.real', values=1); - pt.putinit ('EJ12.real', values=0); - pt.putinit ('EJ21.real', values=0); - pt.putinit ('EJ22.real', values=1); - pt.putinit ('EJ11.imag', values=1); - pt.putinit ('EJ12.imag', values=0); - pt.putinit ('EJ21.imag', values=0); - pt.putinit ('EJ22.imag', values=1); - pt.done(); + pt := parmtable(spaste(fname,'.MEP'), T); + if (is_fail(pt)) fail; + pt.putinit ('frot', values=0); + pt.putinit ('drot', values=0); + pt.putinit ('dell', values=0); + pt.putinit ('gain.11', values=1); + pt.putinit ('gain.22', values=0); + pt.putinit ('EJ11.real', values=1); + pt.putinit ('EJ12.real', values=0); + pt.putinit ('EJ21.real', values=0); + pt.putinit ('EJ22.real', values=1); + pt.putinit ('EJ11.imag', values=1); + pt.putinit ('EJ12.imag', values=0); + pt.putinit ('EJ21.imag', values=0); + pt.putinit ('EJ22.imag', values=1); + pt.done(); } initgsm := function(fname='demo') { - tg := gsm(spaste(fname,'.GSM'), T); - if (is_fail(tg)) fail; - tg.addpointsource ('CP1', [0,1e20], [0,1e20], - 2.734, 0.45379, 1, 0, 0, 0); - tg.addpointsource ('CP2', [0,1e20], [0,1e20], - 2.73402, 0.45369, 0.5, 0, 0, 0); - tg.addpointsource ('CP3', [0,1e20], [0,1e20], - 2.73398, 0.45375, 0.3, 0, 0, 0); - tg.done() - pt := parmtable(spaste(fname,'_gsm.MEP'), T); - pt.loadgsm (spaste(fname,'.GSM')); - pt.done(); + tg := gsm(spaste(fname,'.GSM'), T); + if (is_fail(tg)) fail; + tg.addpointsource ('CP1', [0,1e20], [0,1e20], + 2.734, 0.45379, 1, 0, 0, 0); + tg.addpointsource ('CP2', [0,1e20], [0,1e20], + 2.73402, 0.45369, 0.5, 0, 0, 0); + tg.addpointsource ('CP3', [0,1e20], [0,1e20], + 2.73398, 0.45375, 0.3, 0, 0, 0); + tg.done() + pt := parmtable(spaste(fname,'_gsm.MEP'), T); + pt.loadgsm (spaste(fname,'.GSM')); + pt.done(); } setej := function(fname='demo', mode=1, fact=1) { - pt := parmtable(spaste(fname,'.MEP')); + pt := parmtable(spaste(fname,'.MEP')); # pt.perturb ('NAME==pattern("EJ11.real*")', 0.5, F); - pt.perturb ('NAME==pattern("EJ11.imag*")', -0.5, F); + pt.perturb ('NAME==pattern("EJ11.imag*")', -0.5, F); # pt.perturb ('NAME==pattern("EJ11.SR{1,5,9}.CP{1}")', 0.5*fact, F); - if (mode==2) { + if (mode==2) { # pt.perturb ('NAME==pattern("EJ11.SR2.CP1")', 0.5000015, F); # pt.perturb ('NAME==pattern("EJ11.SR2.CP2")', 0.5000015, F); - pt.perturb ('NAME==pattern("EJ11.SR2.CP3")', 0.5000015, F); - } - pt.done(); + pt.perturb ('NAME==pattern("EJ11.SR2.CP3")', 0.5000015, F); + } + pt.done(); } setgsm := function(fname='demo') { - initgsm(); - + initgsm(fname); + # tg := gsm(spaste(fname,'.GSM'), T); # if (is_fail(tg)) fail; # tg.addpointsource ('CP1', [0,1e20], [0,1e20], @@ -700,65 +700,65 @@ setgsm := function(fname='demo') # tg.addpointsource ('CP3', [0,1e20], [0,1e20], # 2.73399, 0.4537525, 1, 0, 0, 0); # tg.done(); - pt := parmtable(spaste(fname,'_gsm.MEP')); - pt.perturb ('NAME=="RA.CP1"', 0.00003, F); - pt.perturb ('NAME=="RA.CP2"', 0.000005, F); - pt.perturb ('NAME=="RA.CP3"', 0.00001, F); - pt.perturb ('NAME=="DEC.CP1"', -0.000005, F); - pt.perturb ('NAME=="DEC.CP2"', -0.0000025, F); - pt.perturb ('NAME=="DEC.CP3"', 0.0000025, F); - pt.perturb ('NAME=="StokesI.CP1"', 0, F); - pt.perturb ('NAME=="StokesI.CP2"', 0.5, F); - pt.perturb ('NAME=="StokesI.CP3"', 0.7, F); - pt.done(); + pt := parmtable(spaste(fname,'_gsm.MEP')); + pt.perturb ('NAME=="RA.CP1"', 0.00003, F); + pt.perturb ('NAME=="RA.CP2"', 0.000005, F); + pt.perturb ('NAME=="RA.CP3"', 0.00001, F); + pt.perturb ('NAME=="DEC.CP1"', -0.000005, F); + pt.perturb ('NAME=="DEC.CP2"', -0.0000025, F); + pt.perturb ('NAME=="DEC.CP3"', 0.0000025, F); + pt.perturb ('NAME=="StokesI.CP1"', 0, F); + pt.perturb ('NAME=="StokesI.CP2"', 0.5, F); + pt.perturb ('NAME=="StokesI.CP3"', 0.7, F); + pt.done(); } initchan := function(fname='demo') { - t := table(spaste(fname, '.MS/SPECTRAL_WINDOW'),readonly=F) - if (is_fail(t)) fail; - t.putcell ('CHAN_WIDTH', 1, array(1e7,1)); - t.close() -} + t := table(spaste(fname, '.MS/SPECTRAL_WINDOW'),readonly=F) + if (is_fail(t)) fail; + t.putcell ('CHAN_WIDTH', 1, array(1e7,1)); + t.close() + } setchan := function(fname='demo') { - t := table(spaste(fname, '.MS/SPECTRAL_WINDOW'),readonly=F) - if (is_fail(t)) fail; - t.putcell ('CHAN_WIDTH', 1, array(1e6,1)); - t.close() -} + t := table(spaste(fname, '.MS/SPECTRAL_WINDOW'),readonly=F) + if (is_fail(t)) fail; + t.putcell ('CHAN_WIDTH', 1, array(1e6,1)); + t.close() + } init := function(fname='demo') { - initparms(fname=fname); - initgsm(fname=fname); - predict(fname=fname); + initparms(fname=fname); + initgsm(fname=fname); + predict(fname=fname); } demo := function(solver=0,niter=10,fname='demo',sleep=F,sleeptime=2,wait=F) { - global annotator; - - setparms(fname); - setgsm(fname); - if (0 == solver) - { - annotator:=solve(fname=fname, niter=niter, - sleep=sleep, sleeptime=sleeptime, wait=wait); - } - else - { - annotator:=solvepos(fname=fname, niter=niter); - } + global annotator; + + setparms(fname); + setgsm(fname); + if (0 == solver) + { + annotator:=solve(fname=fname, niter=niter, + sleep=sleep, sleeptime=sleeptime, wait=wait); + } + else + { + annotator:=solvepos(fname=fname, niter=niter); + } } demogain := function(niter=10,fname='demo') { - setparms(fname); - initgsm(fname); + setparms(fname); + initgsm(fname); - solvegain(niter=niter,fname=fname); + solvegain(niter=niter,fname=fname); } @@ -855,3 +855,4 @@ testpeelej2a := function() peelej(ant=[], src=3, niter=5, datacol='CORRECTED_DATA', mapnr='2'); return solveej(ant=[], niter=20); } + diff --git a/CEP/CPA/PSS3/CAL/src/mkimg.g b/CEP/CPA/PSS3/CAL/src/mkimg.g index 4c10de2cde8dfd723a3eba89483844ed1e19d6ea..cfd4163eabad960d5e0aa4cc1637d0f67f7f1cbb 100644 --- a/CEP/CPA/PSS3/CAL/src/mkimg.g +++ b/CEP/CPA/PSS3/CAL/src/mkimg.g @@ -24,8 +24,8 @@ todms := function (rad) mkimg := function (msname, imgname, type='model', npix=500, nchan=0, - start=1, step=1, msselect='', mode='mfs', - cellx='0.1arcsec', celly='0.1arcsec') + start=1, step=1, msselect='', mode='mfs', + cellx='0.1arcsec', celly='0.1arcsec') { if (nchan == 0) { t:=table(msname, readonly=F); diff --git a/CEP/CPA/PSS3/CAL/src/pss3-demo.g b/CEP/CPA/PSS3/CAL/src/pss3-demo.g index 59f336a35fad5d59656dd7a859d6420464f1cb2e..e75c2a39eec2782ee22e8af83ebffb9216190c42 100644 --- a/CEP/CPA/PSS3/CAL/src/pss3-demo.g +++ b/CEP/CPA/PSS3/CAL/src/pss3-demo.g @@ -12,7 +12,7 @@ include 'mkimg.g'; # # Demo function showing the predict functionality and creating an image of it. # -predict := function(fname='demo', ant=4*[0:20], +predict := function(fname='demo', ant=5*[0:19], modeltype='LOFAR', calcuvw=F, trace=T) { @@ -45,7 +45,8 @@ predict := function(fname='demo', ant=4*[0:20], mssel := spaste('all([ANTENNA1,ANTENNA2] in ',substitutevar(ant), ')'); } print mssel; - mkimg(spaste(fname, '.MS'), spaste(fname, '.img'), msselect=mssel); + mkimg(spaste(fname, '.MS'), spaste(fname, '.img'), msselect=mssel, + cellx='0.25arcsec', celly='0.25arcsec', npix=1000, type='model'); return T; } @@ -667,18 +668,18 @@ initparms := function(fname='demo') pt.putinit ('dell', values=0); pt.putinit ('gain.11', values=1); pt.putinit ('gain.22', values=1); - pt.putinit ('EJ11.real.SR1.CP1', values=1); - pt.putinit ('EJ11.real.SR1.CP2', values=1); - pt.putinit ('EJ11.real.SR1.CP3', values=1); - pt.putinit ('EJ11.real.SR2.CP1', values=1); - pt.putinit ('EJ11.real.SR2.CP2', values=1); - pt.putinit ('EJ11.real.SR2.CP3', values=1); - pt.putinit ('EJ11.real.SR5.CP1', values=1); - pt.putinit ('EJ11.real.SR5.CP2', values=1); - pt.putinit ('EJ11.real.SR5.CP3', values=1); - pt.putinit ('EJ11.real.SR9.CP1', values=1); - pt.putinit ('EJ11.real.SR9.CP2', values=1); - pt.putinit ('EJ11.real.SR9.CP3', values=1); +# pt.putinit ('EJ11.real.SR1.CP1', values=1); +# pt.putinit ('EJ11.real.SR1.CP2', values=1); +# pt.putinit ('EJ11.real.SR1.CP3', values=1); +# pt.putinit ('EJ11.real.SR2.CP1', values=1); +# pt.putinit ('EJ11.real.SR2.CP2', values=1); +# pt.putinit ('EJ11.real.SR2.CP3', values=1); +# pt.putinit ('EJ11.real.SR5.CP1', values=1); +# pt.putinit ('EJ11.real.SR5.CP2', values=1); +# pt.putinit ('EJ11.real.SR5.CP3', values=1); +# pt.putinit ('EJ11.real.SR9.CP1', values=1); +# pt.putinit ('EJ11.real.SR9.CP2', values=1); +# pt.putinit ('EJ11.real.SR9.CP3', values=1); pt.putinit ('EJ11.real', values=1); pt.putinit ('EJ12.real', values=0); pt.putinit ('EJ21.real', values=0); @@ -763,36 +764,43 @@ initgsm := function(fname='demo') arcsec_per_radian := 180.0*3600.0/pi; alpha0 := (10.0+26.0/60 + 36.3/3600)*15.0*pi/180; delta0 := 26.0*pi/180.0; + print 'alpha0 = ', alpha0; - l := -117.42/arcsec_per_radian; + l := 117.42/arcsec_per_radian; m := 28.38/arcsec_per_radian; - pos := lm_from_alpha_delta(l,m, apha0, delta0); + pos := lm_to_alpha_delta(l,m, alpha0, delta0); tg.addpointsource ('CP1', [0,1e20], [0,1e20], pos[1], pos[2], 1, 0, 0, 0); + print pos; - l := 34.74/arcsec_per_radian; + + l := -34.74/arcsec_per_radian; m := -30.33/arcsec_per_radian; - pos := lm_from_alpha_delta(l,m, apha0, delta0); + pos := lm_to_alpha_delta(l,m, alpha0, delta0); tg.addpointsource ('CP2', [0,1e20], [0,1e20], pos[1], pos[2], 1, 0, 0, 0); + print pos; - l := 65.56/arcsec_per_radian; + l := -65.56/arcsec_per_radian; m := 66.54/arcsec_per_radian; - pos := lm_from_alpha_delta(l,m, apha0, delta0); + pos := lm_to_alpha_delta(l,m, alpha0, delta0); tg.addpointsource ('CP3', [0,1e20], [0,1e20], pos[1], pos[2], 1, 0, 0, 0); + print pos; - l := 93.44/arcsec_per_radian; + l := -93.44/arcsec_per_radian; m := 65.07/arcsec_per_radian; - pos := lm_from_alpha_delta(l,m, apha0, delta0); + pos := lm_to_alpha_delta(l,m, alpha0, delta0); tg.addpointsource ('CP4', [0,1e20], [0,1e20], pos[1], pos[2], 1, 0, 0, 0); + print pos; - l := -8.81/arcsec_per_radian; + l := 8.81/arcsec_per_radian; m := -56.26/arcsec_per_radian; - pos := lm_from_alpha_delta(l,m, apha0, delta0); + pos := lm_to_alpha_delta(l,m, alpha0, delta0); tg.addpointsource ('CP5', [0,1e20], [0,1e20], pos[1], pos[2], 1, 0, 0, 0); + print pos; @@ -830,7 +838,7 @@ setej := function(fname='demo', mode=1, fact=1) setgsm := function(fname='demo') { - initgsm(); + initgsm(fname); # pt := parmtable(spaste(fname,'_gsm.MEP')); # pt.perturb ('NAME=="RA.CP1"', 0.00003, F);