diff --git a/.gitattributes b/.gitattributes
index 9d938e2e0aac871f2e3ceea76ad20f3c83eee967..799f774fb6f6ba189e981d11388316b7fcadc91a 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -954,8 +954,6 @@ CEP/Pipeline/test/test_framework/fixture/pyrap/__init__.py -text
CEP/Pipeline/test/test_framework/fixture/pyrap/tables.py -text
CEP/Pipeline/test/test_framework/fixture/read.me -text
CEP/Pipeline/test/test_framework/unittest_runner.py -text
-CEP/pyparmdb/src/gsmutils.py -text
-CEP/pyparmdb/src/lsm.py -text
CMake/FindAskapSoft.cmake -text
CMake/FindCasarest.cmake -text
CMake/FindJNI.cmake -text
diff --git a/CEP/CMakeLists.txt b/CEP/CMakeLists.txt
index 320f2a7fb22f268fd8f5d1f7f43fdd362f2cf240..7adda5279f50ea72646536fea7b4166dd1612eab 100644
--- a/CEP/CMakeLists.txt
+++ b/CEP/CMakeLists.txt
@@ -2,6 +2,7 @@
lofar_add_package(Calibration)
lofar_add_package(DP3)
+lofar_add_package(GSM)
lofar_add_package(Imager)
lofar_add_package(LMWCommon)
lofar_add_package(MS)
diff --git a/CEP/GSM/CMakeLists.txt b/CEP/GSM/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..669ca137931ca6fbc1806410dde6258f63258cf3
--- /dev/null
+++ b/CEP/GSM/CMakeLists.txt
@@ -0,0 +1,9 @@
+# $Id$
+
+lofar_package(GSM 1.0)
+
+include(LofarFindPackage)
+lofar_find_package(Python REQUIRED)
+
+add_subdirectory(src)
+#add_subdirectory(db)
diff --git a/CEP/GSM/src/CMakeLists.txt b/CEP/GSM/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0090323e0e26c910f4c102f4d8e8d8e92af73896
--- /dev/null
+++ b/CEP/GSM/src/CMakeLists.txt
@@ -0,0 +1,6 @@
+# $Id$
+
+include(PythonInstall)
+
+# Install Python modules
+python_install(gsmutils.py DESTINATION lofar/gsm)
diff --git a/CEP/pyparmdb/src/CMakeLists.txt b/CEP/pyparmdb/src/CMakeLists.txt
index b34a16acc50d3ecbe7380a6ad16de530582546f5..c4f14092f6c1e94ea94fe1ee0ad15c773747f86a 100644
--- a/CEP/pyparmdb/src/CMakeLists.txt
+++ b/CEP/pyparmdb/src/CMakeLists.txt
@@ -26,8 +26,4 @@ lofar_add_bin_program(versionpyparmdb versionpyparmdb.cc)
# Install Python modules
python_install(__init__.py DESTINATION lofar/parmdb)
-python_install(gsmutils.py DESTINATION lofar)
-install(PROGRAMS
- gsm.py
- DESTINATION bin)
diff --git a/CEP/pyparmdb/src/gsm.py b/CEP/pyparmdb/src/gsm.py
deleted file mode 100755
index 560f6a493f31fb95b8799fcc82d8e3eba3ad4eea..0000000000000000000000000000000000000000
--- a/CEP/pyparmdb/src/gsm.py
+++ /dev/null
@@ -1,287 +0,0 @@
-#!/usr/bin/env python
-
-# Select a subset from the GSM and write into a makesourcdb input file.
-# ra and dec both in degrees
-# integrated flux [in Jy] above which the sources will be selected
-# You should know that:
-# cat_id = 3 => NVSS
-# cat_id = 4 => VLSS
-# cat_id = 5 => WENSS
-def doSelect (outfile, dirSelect, fluxi_mins, cat_ids):
-
- import os
- import monetdb.sql as db
- from monetdb.sql import Error as Error
-
- db_type = "MonetDB"
- db_host = "ldb001"
- db_port = 50000
- db_dbase = "gsm"
- db_user = "gsm"
- db_passwd = "msss"
- #db_lang = "sql"
-
- # If a single catalog, use =
- if len(cat_ids) == 1:
- catSelect = 'c1.cat_id = %d AND c1.i_int_avg >= %f' % (cat_ids[0], fluxi_mins[0])
- else:
- # Multiple catalogs.
- # If a single minimum flux, use it for all catalogs.
- # Note that the IN clause needs (), not [].
- if len(fluxi_mins) == 1:
- catSelect = 'c1.cat_id IN (%s) AND c1.i_int_avg >= %f' % (str(cat_ids)[1:-1], fluxi_mins[0])
- else:
- # There is a minimum flux per catalog.
- catSelect = ''
- for i in range(len(cat_ids)):
- if i > 0:
- catSelect += ' OR '
- catSelect += '(c1.cat_id = %d AND c1.i_int_avg >= %f)' % (cat_ids[i], fluxi_mins[i])
- # Form the entire where clause.
- where = 'WHERE (%s) AND (%s)' % (catSelect, dirSelect)
-
- conn = db.connect(hostname=db_host \
- ,port=db_port \
- ,database=db_dbase \
- ,username=db_user \
- ,password=db_passwd \
- )
-
- # Remove output file if already existing.
- if os.path.isfile(outfile):
- os.remove(outfile)
-
- try:
- cursor = conn.cursor()
- # This query concatenates the requested columns per row in a single string in the correct BBS format.
- cursor.execute(#"COPY " + \
- "SELECT t.line " + \
- " FROM (SELECT CAST(CONCAT(t0.catsrcname, CONCAT(',', " + \
- " CONCAT(t0.src_type, CONCAT(',', " + \
- " CONCAT(ra2bbshms(t0.ra), CONCAT(',', " + \
- " CONCAT(decl2bbsdms(t0.decl), CONCAT(',', " + \
- " CONCAT(t0.i, CONCAT(',', " + \
- " CONCAT(t0.q, CONCAT(',', " + \
- " CONCAT(t0.u, CONCAT(',', " + \
- " CONCAT(t0.v, CONCAT(',', " + \
- " CONCAT(t0.MajorAxis, CONCAT(',', " + \
- " CONCAT(t0.MinorAxis, CONCAT(',', " + \
- " CONCAT(t0.Orientation, CONCAT(',', " + \
- " CONCAT(t0.ReferenceFrequency, CONCAT(',', " + \
- " CASE WHEN t0.SpectralIndexDegree = 0" + \
- " THEN '[]'" + \
- " ELSE CASE WHEN t0.SpectralIndexDegree = 1" + \
- " THEN CONCAT('[', CONCAT(t0.SpectralIndex_0, ']'))" + \
- " ELSE CONCAT('[', CONCAT(t0.SpectralIndex_0, CONCAT(',', CONCAT(t0.SPECTRALINDEX_1, ']'))))" + \
- " END" + \
- " END" + \
- ")))))))))))))))))))))))" + \
- " ) AS VARCHAR(300)) AS line " + \
- " FROM (SELECT CAST(TRIM(c1.catsrcname) AS VARCHAR(23)) AS catsrcname " + \
- " ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('POINT' AS VARCHAR(23)) " + \
- " ELSE CAST('GAUSSIAN' AS VARCHAR(23)) " + \
- " END AS src_type " + \
- " ,CAST(c1.ra AS VARCHAR(23)) AS ra " + \
- " ,CAST(c1.decl AS VARCHAR(23)) AS decl " + \
- " ,CAST(c1.i_int_avg AS VARCHAR(23)) AS i " + \
- " ,CAST(0 AS VARCHAR(23)) AS q " + \
- " ,CAST(0 AS VARCHAR(23)) AS u " + \
- " ,CAST(0 AS VARCHAR(23)) AS v " + \
- " ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(23)) " + \
- " ELSE CASE WHEN c1.major IS NULL " + \
- " THEN CAST('' AS VARCHAR(23)) " + \
- " ELSE CAST(c1.major AS varchar(23)) " + \
- " END " + \
- " END AS MajorAxis " + \
- " ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(23)) " + \
- " ELSE CASE WHEN c1.minor IS NULL " + \
- " THEN CAST('' AS VARCHAR(23)) " + \
- " ELSE CAST(c1.minor AS varchar(23)) " + \
- " END " + \
- " END AS MinorAxis " + \
- " ,CASE WHEN c1.pa IS NULL " + \
- " THEN CAST('' AS VARCHAR(23)) " + \
- " ELSE CAST(c1.pa AS varchar(23)) " + \
- " END AS Orientation " + \
- " ,CAST(c1.freq_eff AS VARCHAR(23)) AS ReferenceFrequency " + \
- " ,CASE WHEN si.spindx_degree IS NULL " + \
- " THEN 0 " + \
- " ELSE si.spindx_degree " + \
- " END AS SpectralIndexDegree " + \
- " ,CASE WHEN si.c0 IS NULL " + \
- " THEN 0 " + \
- " ELSE si.c0 " + \
- " END AS SpectralIndex_0 " + \
- " ,CASE WHEN si.c1 IS NULL " + \
- " THEN 0 " + \
- " ELSE si.c1 " + \
- " END AS SpectralIndex_1 " + \
- " FROM catalogedsources c1 " + \
- " LEFT OUTER JOIN spectralindices si ON c1.catsrcid = si.catsrc_id " + \
- where +
- " ) t0 " + \
- " ) t " , ())
- y = cursor.fetchall()
- cursor.close()
- #print "y:", y
-
- except db.Error, e:
- import logging
- logging.warn("gsmSelect to file %s failed " % (outfile))
- logging.warn("Failed on query for reason: %s " % (e))
- logging.debug("Failed select query: %s" % (e))
- return -1
-
- file = open(outfile,'w')
- file.write ("format = Name, Type, Ra, Dec, I, Q, U, V, MajorAxis, MinorAxis, Orientation, ReferenceFrequency='60e6', SpectralIndex='[]'\n")
- for i in range(len(y)):
- file.write(y[i][0] + '\n')
- file.close()
- return len(y)
-
-
-def gsmSelectCone (outfile, ra, dec, radius, fluxi_mins, cat_ids):
- if ra < 0 or ra > 360:
- print 'RA', ra, 'is invalid: <0 or >360 degrees'
- return -1
- if dec < -90 or dec > 90:
- print 'DEC', dec, 'is invalid: <-90 or >90 degrees'
- return -1
- if radius < 0 or radius > 90:
- print 'radius', radius, 'is invalid: <0 or >90 degrees'
- return -1
- # Determine DEC boundaries (for hopefully faster query).
- # Convert degrees to radians.
- import math
- mindec = dec - radius
- maxdec = dec + radius
- d2r = math.pi / 180.
- ra *= d2r
- dec *= d2r
- radius *= d2r
- # Source (srcRA,srcDEC) is inside the cone if:
- # sin(dec)*sin(srcDEC) + cos(dec)*cos(srcDEC)*cos(ra-srcRA) >= cos(radius)
- select = '(c1.decl BETWEEN %17.12f AND %17.12f) AND (%17.14f*sin(c1.decl*%17.14f) + %17.14f*cos(c1.decl*%17.14f)*cos(%17.14f-c1.ra*%17.14f) > %17.14f)' % (mindec, maxdec, math.sin(dec), d2r, math.cos(dec), d2r, ra, d2r, math.cos(radius))
- return doSelect (outfile, select, fluxi_mins, cat_ids)
-
-def gsmSelectBox (outfile, ra_st, ra_end, dec_st, dec_end, fluxi_mins, cat_ids):
- for ra in [ra_st,ra_end]:
- if ra < 0 or ra > 360:
- print 'RA', ra, 'is invalid: <0 or >360 degrees'
- return -1
- for dec in [dec_st,dec_end]:
- if dec < -90 or dec > 90:
- print 'DEC', dec, 'is invalid: <-90 or >90 degrees'
- return -1
- if dec_st > dec_end:
- print 'invalid DEC: start', dec_st, '> end', dec_end
- return -1
- # Form the selection for the RA,DEC.
- # If RA crosses 360 degrees, split in two parts.
- if ra_st <= ra_end:
- select = '(c1.ra BETWEEN %17.12f AND %17.12f)' % (ra_st, ra_end)
- else:
- select = '(c1.ra BETWEEN %17.12f AND 360 OR c1.ra BETWEEN 0 AND %17.12f)' % (ra_st, ra_end)
- select += ' AND (c1.decl BETWEEN %17.12f AND %17.12f)' % (dec_st, dec_end)
- return doSelect (outfile, select, fluxi_mins, cat_ids)
-
-# Interpret the arguments and do the selection.
-def gsmMain (name, argv):
- useBox = False
- starg = 0
- optarg = 4
- if len(argv) > 0 and argv[0] == '-b':
- # select using a box
- useBox = True
- starg = 1
- optarg = 6
- if len(argv) < optarg:
- print ''
- print 'Insufficient arguments given; run as:'
- print ''
- print ' %s outfile RA DEC radius [minFluxI [catalogs]]' % name
- print 'to select using a cone'
- print ''
- print ' %s -b outfile stRA endRA stDEC endDEC [minFluxI [catalogs]]' % name
- print 'to select using a box'
- print ''
- print ' outfile path-name of the output file'
- print ' It will be overwritten if already existing'
- print ' RA cone center Right Ascension (J2000, degrees)'
- print ' DEC cone center Declination (J2000, degrees)'
- print ' radius cone radius (degrees)'
- print ' stRA box start Right Ascension (J2000, degrees)'
- print ' endRA box end Right Ascension (J2000, degrees)'
- print ' stRA can be > endRA indicating crossing 360 deg'
- print ' minDEC box start Declination (J2000, degrees)'
- print ' maxDEC box end Declination (J2000, degrees)'
- print ' minFluxI minimum flux (integrated Stokes I flux in Jy)'
- print ' If a single value is given, it is used for all catalogs.'
- print ' If multiple values given (separated by commas), each applies'
- print ' to the corresponding catalog.'
- print ' Default is 0.5'
- print ' catalogs names of catalogs to search (case-insensitive)'
- print ' (NVSS, VLSS, and/or WENSS)'
- print ' If multiple, separate by commas'
- print ' Default is WENSS'
- print ''
- return False
- # Get the arguments.
- outfile = argv[starg]
- a1 = float(argv[starg+1])
- a2 = float(argv[starg+2])
- a3 = float(argv[starg+3])
- if useBox:
- a4 = float(argv[starg+4])
- minFlux = [0.5]
- if len(argv) > optarg and len(argv[optarg]) > 0:
- minFlux = [float(x) for x in argv[optarg].split(',')]
- # cat_id = 3 => NVSS
- # cat_id = 4 => VLSS
- # cat_id = 5 => WENSS
- cats = ['WENSS']
- if len(argv) > optarg+1 and len(argv[optarg+1]) > 0:
- cats = argv[optarg+1].split (',')
- cat_ids = []
- for cat in cats:
- cat = cat.upper()
- if cat == 'NVSS':
- cat_ids.append (3)
- elif cat == 'VLSS':
- cat_ids.append (4)
- elif cat == 'WENSS':
- cat_ids.append (5)
- else:
- print cat, 'is an invalid catalog name'
- return False
- if len(minFlux) != 1 and len(minFlux) != len(cat_ids):
- print 'Nr of minFlux values (%s) must be 1 or match nr of catalogs' % argv[5]
- return False;
- # Do the selection and create the makesourcedb file.
- if useBox:
- nr = gsmSelectBox (outfile, a1, a2, a3, a4, minFlux, cat_ids)
- else:
- nr = gsmSelectCone (outfile, a1, a2, a3, minFlux, cat_ids)
- if nr < 0:
- return False
- elif nr == 0:
- print 'Warning: no matching sources found in GSM for'
- if useBox:
- print ' stRA=%f endRA=%f stDEC=%f endDEC=%f' % (a1,a2,a3,a4)
- else:
- print ' RA=%f DEC=%f radius=%f' % (a1,a2,a3)
- print ' minFluxI=%s catalogs=%s' % (str(minFlux), cats)
- else:
- print '%d sources written into %s' % (nr, outfile)
- return True
-
-
-# This is the main entry.
-if __name__ == "__main__":
- import sys
- if gsmMain (sys.argv[0], sys.argv[1:]):
- sys.exit(0)
- sys.exit(1)
diff --git a/CEP/pyparmdb/src/gsmutils.py b/CEP/pyparmdb/src/gsmutils.py
deleted file mode 100644
index a2c372ddeb4dcd9826ea3113f869af77ad099326..0000000000000000000000000000000000000000
--- a/CEP/pyparmdb/src/gsmutils.py
+++ /dev/null
@@ -1,537 +0,0 @@
-# LOFAR IMAGING PIPELINE
-#
-# BBS Source Catalogue List
-# Bart Scheers, 2011
-# L.H.A.Scheers@uva.nl
-# ------------------------------------------------------------------------------
-
-import sys, string
-import numpy as np
-import monetdb.sql as db
-import logging
-
-def expected_fluxes_in_fov(conn, ra_central, decl_central, fov_radius, assoc_theta, bbsfile, storespectraplots = False, deruiter_radius = 0.):
- """Search for VLSS, WENSS and NVSS sources that
- are in the given FoV. The FoV is set by its central position
- (ra_central, decl_central) out to a radius of fov_radius.
- The query looks for cross-matches around the sources, out
- to a radius of assoc_theta.
-
- All units are in degrees.
-
- deruiter_radius is a measure for the association uncertainty that takes
- position errors into account (see thesis Bart Scheers). If not given
- as a positive value, it is read from the TKP config file. If not
- available, it defaults to 3.717.
-
- The query returns all vlss sources (id) that are in the FoV.
- If so, the counterparts from other catalogues are returned as well
- (also their ids).
- """
-
- DERUITER_R = deruiter_radius
- if DERUITER_R <= 0:
- try:
- from tkp.config import config
- DERUITER_R = config['source_association']['deruiter_radius']
- print "DERUITER_R =", DERUITER_R
- except:
- DERUITER_R = 3.717
-
- if ra_central + alpha(fov_radius, decl_central) > 360:
- "This will be implemented soon"
- raise BaseException("ra = %s > 360 degrees, not implemented yet" % str(ra_central + alpha(fov_radius, decl_central)))
-
- skymodel = open(bbsfile, 'w')
- header = "# (Name, Type, Ra, Dec, I, Q, U, V, ReferenceFrequency='60e6', SpectralIndex='[0.0]', MajorAxis, MinorAxis, Orientation) = format\n\n"
- skymodel.write(header)
-
- # This is dimensionless search radius that takes into account
- # the ra and decl difference between two sources weighted by
- # their positional errors.
- deRuiter_reduced = DERUITER_R / 3600.
- try:
- cursor = conn.cursor()
- query = """\
-SELECT t0.v_catsrcid
- ,t0.catsrcname
- ,t1.wm_catsrcid
- ,t2.wp_catsrcid
- ,t3.n_catsrcid
- ,t0.v_flux
- ,t1.wm_flux
- ,t2.wp_flux
- ,t3.n_flux
- ,t0.v_flux_err
- ,t1.wm_flux_err
- ,t2.wp_flux_err
- ,t3.n_flux_err
- ,t1.wm_assoc_distance_arcsec
- ,t1.wm_assoc_r
- ,t2.wp_assoc_distance_arcsec
- ,t2.wp_assoc_r
- ,t3.n_assoc_distance_arcsec
- ,t3.n_assoc_r
- ,t0.pa
- ,t0.major
- ,t0.minor
- ,t0.ra
- ,t0.decl
- FROM (SELECT c1.catsrcid AS v_catsrcid
- ,c1.catsrcname
- ,c1.ra
- ,c1.decl
- ,c1.i_int_avg AS v_flux
- ,c1.i_int_avg_err AS v_flux_err
- ,c1.pa
- ,c1.major
- ,c1.minor
- FROM (SELECT catsrcid
- ,catsrcname
- ,ra
- ,decl
- ,pa
- ,major
- ,minor
- ,i_int_avg
- ,i_int_avg_err
- FROM catalogedsources
- WHERE cat_id = 4
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c1
- ) t0
- FULL OUTER JOIN
- (SELECT c1.catsrcid AS v_catsrcid
- ,c2.catsrcid AS wm_catsrcid
- ,c2.i_int_avg AS wm_flux
- ,c2.i_int_avg_err AS wm_flux_err
- ,3600 * DEGREES(2 * ASIN(SQRT( (c1.x - c2.x) * (c1.x - c2.x)
- + (c1.y - c2.y) * (c1.y - c2.y)
- + (c1.z - c2.z) * (c1.z - c2.z)
- ) / 2)
- ) AS wm_assoc_distance_arcsec
- ,SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))
- ) AS wm_assoc_r
- FROM (SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- FROM catalogedsources
- WHERE cat_id = 4
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c1
- ,(SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- ,i_int_avg
- ,i_int_avg_err
- FROM catalogedsources
- WHERE cat_id = 5
- AND (src_type = 'S' OR src_type = 'C')
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c2
- WHERE c1.x * c2.x + c1.y * c2.y + c1.z * c2.z > COS(RADIANS(%s))
- AND SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))) < %s
- ) t1
- ON t0.v_catsrcid = t1.v_catsrcid
- FULL OUTER JOIN
- (SELECT c1.catsrcid AS v_catsrcid
- ,c2.catsrcid AS wp_catsrcid
- ,c2.i_int_avg AS wp_flux
- ,c2.i_int_avg_err AS wp_flux_err
- ,3600 * DEGREES(2 * ASIN(SQRT( (c1.x - c2.x) * (c1.x - c2.x)
- + (c1.y - c2.y) * (c1.y - c2.y)
- + (c1.z - c2.z) * (c1.z - c2.z)
- ) / 2)
- ) AS wp_assoc_distance_arcsec
- ,SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))
- ) AS wp_assoc_r
- FROM (SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- FROM catalogedsources
- WHERE cat_id = 4
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c1
- ,(SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- ,i_int_avg
- ,i_int_avg_err
- FROM catalogedsources
- WHERE cat_id = 6
- AND (src_type = 'S' OR src_type = 'C')
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c2
- WHERE c1.x * c2.x + c1.y * c2.y + c1.z * c2.z > COS(RADIANS(%s))
- AND SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))) < %s
- ) t2
- ON t0.v_catsrcid = t2.v_catsrcid
- FULL OUTER JOIN
- (SELECT c1.catsrcid AS v_catsrcid
- ,c2.catsrcid AS n_catsrcid
- ,c2.i_int_avg AS n_flux
- ,c2.i_int_avg_err AS n_flux_err
- ,3600 * DEGREES(2 * ASIN(SQRT( (c1.x - c2.x) * (c1.x - c2.x)
- + (c1.y - c2.y) * (c1.y - c2.y)
- + (c1.z - c2.z) * (c1.z - c2.z)
- ) / 2)
- ) AS n_assoc_distance_arcsec
- ,SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))
- ) AS n_assoc_r
- FROM (SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- FROM catalogedsources
- WHERE cat_id = 4
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c1
- ,(SELECT catsrcid
- ,ra
- ,decl
- ,ra_err
- ,decl_err
- ,x
- ,y
- ,z
- ,i_int_avg
- ,i_int_avg_err
- FROM catalogedsources
- WHERE cat_id = 3
- AND zone BETWEEN CAST(FLOOR(CAST(%s AS DOUBLE) - %s) AS INTEGER)
- AND CAST(FLOOR(CAST(%s AS DOUBLE) + %s) AS INTEGER)
- AND decl BETWEEN CAST(%s AS DOUBLE) - %s
- AND CAST(%s AS DOUBLE) + %s
- AND ra BETWEEN CAST(%s AS DOUBLE) - alpha(%s, %s)
- AND CAST(%s AS DOUBLE) + alpha(%s, %s)
- AND x * COS(RADIANS(%s)) * COS(RADIANS(%s))
- + y * COS(RADIANS(%s)) * SIN(RADIANS(%s))
- + z * SIN(RADIANS(%s)) > COS(RADIANS(%s))
- ) c2
- WHERE c1.x * c2.x + c1.y * c2.y + c1.z * c2.z > COS(RADIANS(%s))
- AND SQRT(((c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- * (c1.ra * COS(RADIANS(c1.decl)) - c2.ra * COS(RADIANS(c2.decl)))
- / (c1.ra_err * c1.ra_err + c2.ra_err * c2.ra_err))
- + ((c1.decl - c2.decl) * (c1.decl - c2.decl)
- / (c1.decl_err * c1.decl_err + c2.decl_err * c2.decl_err))) < %s
- ) t3
- ON t0.v_catsrcid = t3.v_catsrcid
- """
- cursor.execute(query, (
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- assoc_theta, deRuiter_reduced,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- assoc_theta, deRuiter_reduced,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius, decl_central, fov_radius,
- ra_central, fov_radius, decl_central, ra_central, fov_radius, decl_central,
- decl_central, ra_central, decl_central, ra_central, decl_central, fov_radius,
- assoc_theta, deRuiter_reduced
- ))
- results = zip(*cursor.fetchall())
- if len(results) != 0:
- vlss_catsrcid = results[0]
- vlss_name = results[1]
- wenssm_catsrcid = results[2]
- wenssp_catsrcid = results[3]
- nvss_catsrcid = results[4]
- v_flux = results[5]
- wm_flux = results[6]
- wp_flux = results[7]
- n_flux = results[8]
- v_flux_err = results[9]
- wm_flux_err = results[10]
- wp_flux_err = results[11]
- n_flux_err = results[12]
- wm_assoc_distance_arcsec = results[13]
- wm_assoc_r = results[14]
- wp_assoc_distance_arcsec = results[15]
- wp_assoc_r = results[16]
- n_assoc_distance_arcsec = results[17]
- n_assoc_r = results[18]
- pa = results[19]
- major = results[20]
- minor = results[21]
- ra = results[22]
- decl = results[23]
- spectrumfiles = []
- for i in range(len(vlss_catsrcid)):
- print "\ni = ", i
- bbsrow = ""
- # Here we check the cases for the degree of the polynomial spectral index fit
- print vlss_catsrcid[i], wenssm_catsrcid[i], wenssp_catsrcid[i], nvss_catsrcid[i]
- #print "VLSS",vlss_name[i]
- bbsrow += vlss_name[i] + ", "
- # According to Jess, only sources that have values for all
- # three are considered as GAUSSIAN
- if pa[i] is not None and major[i] is not None and minor[i] is not None:
- #print "Gaussian:", pa[i], major[i], minor[i]
- bbsrow += "GAUSSIAN, "
- else:
- #print "POINT"
- bbsrow += "POINT, "
- #print "ra = ", ra[i], "; decl = ", decl[i]
- #print "BBS ra = ", ra2bbshms(ra[i]), "; BBS decl = ", decl2bbsdms(decl[i])
- bbsrow += ra2bbshms(ra[i]) + ", " + decl2bbsdms(decl[i]) + ", "
- # Stokes I id default, so filed is empty
- #bbsrow += ", "
- lognu = []
- logflux = []
- lognu.append(np.log10(74.0 / 60.0))
- logflux.append(np.log10(v_flux[i]))
- if wenssm_catsrcid[i] is not None:
- lognu.append(np.log10(325.0 / 60.0))
- logflux.append(np.log10(wm_flux[i]))
- if wenssp_catsrcid[i] is not None:
- lognu.append(np.log10(352.0 / 60.0))
- logflux.append(np.log10(wp_flux[i]))
- if nvss_catsrcid[i] is not None:
- lognu.append(np.log10(1400.0 / 60.0))
- logflux.append(np.log10(n_flux[i]))
- f = ""
- for j in range(len(logflux)):
- f += str(10 ** logflux[j]) + "; "
- print f
- #print "len(lognu) = ",len(lognu), "nvss_catsrcid[",i,"] =", nvss_catsrcid[i]
- # Here we write the expected flux values at 60 MHz, and the fitted spectral index and
- # and curvature term
- if len(lognu) == 1:
- #print "Exp. flux:", 10**(np.log10(v_flux[i]) + 0.7 * np.log10(74.0/60.0))
- #print "Default -0.7"
- bbsrow += str(round(10 ** (np.log10(v_flux[i]) + 0.7 * np.log10(74.0 / 60.0)), 2)) + ", , , , , "
- bbsrow += "[-0.7]"
- elif len(lognu) == 2 or (len(lognu) == 3 and nvss_catsrcid[i] is None):
- #print "Do a 1-degree polynomial fit"
- # p has form : p(x) = p[0] + p[1]*x
- p = np.poly1d(np.polyfit(np.array(lognu), np.array(logflux), 1))
- #print p
- if storespectraplots:
- spectrumfile = plotSpectrum(np.array(lognu), np.array(logflux), p, "spectrum_%s.eps" % vlss_name[i])
- spectrumfiles.append(spectrumfile)
- # Default reference frequency is reported, so we leave it empty here;
- # Catalogues just report on Stokes I, so others are empty.
- bbsrow += str(round(10 ** p[0], 4)) + ", , , , , "
- bbsrow += "[" + str(round(p[1], 4)) + "]"
- elif (len(lognu) == 3 and nvss_catsrcid[i] is not None) or len(lognu) == 4:
- #print "Do a 2-degree polynomial fit"
- # p has form : p(x) = p[0] + p[1]*x + p[2]*x**2
- p = np.poly1d(np.polyfit(np.array(lognu), np.array(logflux), 2))
- #print p
- if storespectraplots:
- spectrumfile = plotSpectrum(np.array(lognu), np.array(logflux), p, "spectrum_%s.eps" % vlss_name[i])
- spectrumfiles.append(spectrumfile)
- # Default reference frequency is reported, so we leave it empty here
- bbsrow += str(round(10 ** p[0], 4)) + ", , , , , "
- bbsrow += "[" + str(round(p[1], 4)) + ", " + str(round(p[2], 4)) + "]"
- if pa[i] is not None and major[i] is not None and minor[i] is not None:
- # Gaussian source:
- bbsrow += ", " + str(round(major[i], 2)) + ", " + str(round(minor[i], 2)) + ", " + str(round(pa[i], 2))
- #print bbsrow
- skymodel.write(bbsrow + '\n')
-
- if storespectraplots:
- print "Spectra available in:", spectrumfiles
-
- skymodel.close()
- print "Sky model stored in source table:", bbsfile
-
- except db.Error, e:
- logging.warn("Failed on query nr %s; for reason %s" % (query, e))
- raise
- finally:
- cursor.close()
-
-def plotSpectrum(x, y, p, f):
- import pylab
- expflux = "Exp. flux: " + str(round(10 ** p(0), 3)) + " Jy"
- fig = pylab.figure()
- ax = fig.add_subplot(111)
- for i in range(len(ax.get_xticklabels())):
- ax.get_xticklabels()[i].set_size('x-large')
- for i in range(len(ax.get_yticklabels())):
- ax.get_yticklabels()[i].set_size('x-large')
- ax.set_xlabel(r'$\log \nu/\nu_0$', size = 'x-large')
- ax.set_ylabel('$\log S$', size = 'x-large')
- # Roughly between log10(30/60) and log10(1500/60)
- xp = np.linspace(-0.3, 1.5, 100)
- ax.plot(x, y, 'o', label = 'cat fluxes')
- ax.plot(0.0, p(0), 'o', color = 'k', label = expflux)
- ax.plot(xp, p(xp), linestyle = '--', linewidth = 2, label = 'fit')
- pylab.legend(numpoints = 1, loc = 'best')
- pylab.grid(True)
- pylab.savefig(f, dpi = 600)
- return f
-
-
-def decl2bbsdms(d):
- """Based on function deg2dec Written by Enno Middelberg 2001
- http://www.atnf.csiro.au/people/Enno.Middelberg/python/python.html
- """
- deg = float(d)
- sign = "+"
-
- # test whether the input numbers are sane:
-
- # if negative, store "-" in sign and continue calulation
- # with positive value
-
- if deg < 0:
- sign = "-"
- deg = deg * (-1)
-
- #if deg > 180:
- # logging.warn("%s: inputs may not exceed 180!" % deg)
- # raise
-
- #if deg > 90:
- # print `deg`+" exceeds 90, will convert it to negative dec\n"
- # deg=deg-90
- # sign="-"
-
- if deg < -90 or deg > 90:
- logging.warn("%s: inputs may not exceed 90 degrees!" % deg)
-
- hh = int(deg)
- mm = int((deg - int(deg)) * 60)
- ss = '%10.8f' % (((deg - int(deg)) * 60 - mm) * 60)
-
- #print '\t'+sign+string.zfill(`hh`,2)+':'+string.zfill(`mm`,2)+':'+'%10.8f' % ss
- #print '\t'+sign+string.zfill(`hh`,2)+' '+string.zfill(`mm`,2)+' '+'%10.8f' % ss
- #print '\t'+sign+string.zfill(`hh`,2)+'h'+string.zfill(`mm`,2)+'m'+'%10.8fs\n' % ss
- return sign + string.zfill(`hh`, 2) + '.' + string.zfill(`mm`, 2) + '.' + string.zfill(ss, 11)
-
-def ra2bbshms(a):
- deg = float(a)
-
- # test whether the input numbers are sane:
-
- if deg < 0 or deg > 360:
- logging.warn("%s: inputs may not exceed 90 degrees!" % deg)
-
- hh = int(deg / 15)
- mm = int((deg - 15 * hh) * 4)
- ss = '%10.8f' % ((4 * deg - 60 * hh - mm) * 60)
-
- #print '\t'+string.zfill(`hh`,2)+':'+string.zfill(`mm`,2)+':'+'%10.8f' % ss
- #print '\t'+string.zfill(`hh`,2)+' '+string.zfill(`mm`,2)+' '+'%10.8f' % ss
- #print '\t'+string.zfill(`hh`,2)+'h'+string.zfill(`mm`,2)+'m'+'%10.8fs\n' % ss
- return string.zfill(`hh`, 2) + ':' + string.zfill(`mm`, 2) + ':' + string.zfill(ss, 11)
-
-def alpha(theta, decl):
- if abs(decl) + theta > 89.9:
- return 180.0
- else:
- return degrees(abs(np.arctan(np.sin(radians(theta)) / np.sqrt(abs(np.cos(radians(decl - theta)) * np.cos(radians(decl + theta)))))))
-
-def degrees(r):
- return r * 180 / np.pi
-
-def radians(d):
- return d * np.pi / 180
diff --git a/CEP/pyparmdb/src/lsm.py b/CEP/pyparmdb/src/lsm.py
deleted file mode 100644
index f943b2c8b4fb8d8b35a836d8c2f323d864028f78..0000000000000000000000000000000000000000
--- a/CEP/pyparmdb/src/lsm.py
+++ /dev/null
@@ -1,43 +0,0 @@
-#!/usr/bin/python
-
-#
-# Script used for testing the bbs skymodel files
-#
-import sys, os, time
-
-import monetdb
-import monetdb.sql as db
-import gsmutils as gsm
-
-db_host = "ldb001"
-db_dbase = "gsm"
-db_user = "gsm"
-db_passwd = "msss"
-db_port = 50000
-
-try:
- conn = monetdb.sql.connect(hostname=db_host, database=db_dbase, username=db_user, password=db_passwd, port=db_port)
-except db.Error, e:
- raise
-
-#ra_c = 289.89258333333333
-#decl_c = 0.0017444444444444445
-
-ra_c = 287.80216666666666
-decl_c = 9.096861111111112
-
-#ra_c = 362
-#decl_c = 10
-
-#ra_c = 286.87425
-#decl_c = 7.1466111111111115
-
-#ra_c = 70.0
-#decl_c = 33.0
-fov_radius = 5.0
-assoc_theta = 0.025
-
-gsm.expected_fluxes_in_fov(conn, ra_c, decl_c, fov_radius, assoc_theta, 'bbs.skymodel.test', storespectraplots=False)
-
-conn.close()
-
diff --git a/CMake/LofarPackageList.cmake b/CMake/LofarPackageList.cmake
index c3275e3d7a03ecac9ea5be31ea5ae4807b01fe3e..613bb07f07cc3472c1cf10a7da1ea3ae0ceecaa1 100644
--- a/CMake/LofarPackageList.cmake
+++ b/CMake/LofarPackageList.cmake
@@ -16,6 +16,7 @@ if(NOT DEFINED LOFAR_PACKAGE_LIST_INCLUDED)
set(LOFAR_PACKAGE_LIST_INCLUDED TRUE)
set(Calibration_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/Calibration)
set(DP3_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3)
+ set(GSM_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/GSM)
set(Imager_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/Imager)
set(LMWCommon_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/LMWCommon)
set(MS_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/MS)