Initial commit

1 files changed, 413 insertions, 0 deletions

diff --git a/vo/src/OLD/sloanspec.cl b/vo/src/OLD/sloanspec.cl
new file mode 100644
index 00000000..c6f24aed
--- /dev/null
+++ b/vo/src/OLD/sloanspec.cl
@@ -0,0 +1,413 @@
+#{  SLOANSPEC --  Query the SDSS Spectral Data (DR4) given an image footprint,
+#   object name, or equatorial position and search size.
+
+procedure sloanspec (what)
+
+string  what 			{ prompt = "Image or Object name"	  }
+
+real	ra = INDEF		{ prompt = "RA of Field (hours)"	  }
+real	dec = INDEF		{ prompt = "Dec of Field (degrees)"	  }
+real	sr = 0.25		{ prompt = "Search size (degrees)"	  }
+
+int     maxrecords  = 100       { prompt = "Max records to return"        }
+bool	display     = no	{ prompt = "Display result?"		  }
+bool	interactive = yes	{ prompt = "Interactive mode?"		  }
+bool	autosave    = yes	{ prompt = "Auto-save downloaded data?"	  }
+bool	verbose     = yes	{ prompt = "Verbose output?"		  }
+bool	ned         = no	{ prompt = "Overlay NED objects?"	  }
+bool	grid        = no	{ prompt = "Overlay coord grid?"	  }
+bool	download    = no	{ prompt = "Automatically download data?" }
+
+string  spectype    = "all"	{ prompt = "Spectra type (qso|galaxy|all)",
+					enum="qso|galaxy|all"		  }
+real	zlow        = INDEF	{ prompt = "Low redshift filter"	  }
+real	zhigh       = INDEF	{ prompt = "High redshift filter"	  }
+
+real	size = 15.0		{ prompt = "Size of marker (pixels)"	  }
+int	color = 207		{ prompt = "Marker color"		  }
+int	status = 0		{ prompt = "Service status code"	  }
+
+string  *llist
+
+begin
+    string this, img, obj, stab, ttab, mtab, coords, tvcoords, filter
+    string name, class, oclass, ch, cmd, tnames, imlist, sv_img, stype
+    bool   verb, disp, do_grid, interact, saving, do_ned, found, getdata
+    int    maxrec, mkcolor, n, wcs, siap, count
+    real   szmark, d0, d1, xp, yp, z, x1, x2, y1, y2
+    real   lra, ldec, lsr, sra, sdec, sz, sx, sy
+    real   llx, lly, urx, ury, cx, cy, bsz, z1, z2
+    string sname, sid, sclass, mission
+
+
+    lra      = ra
+    ldec     = dec
+    this     = ""
+    if (isindef(lra) && isindef(ldec))
+        this     = what
+
+    lsr      = sr
+    disp     = display
+    maxrec   = maxrecords
+    szmark   = size
+    mkcolor  = color
+    verb     = verbose
+    interact = interactive
+    do_grid  = grid
+    do_ned   = ned
+    saving   = autosave
+    stype    = spectype
+    z1       = zlow
+    z2       = zhigh
+    getdata  = download
+    sv_img   = ""
+
+    # Set the environment.
+    set clobber = yes
+    set imclobber = yes
+
+
+    if (imaccess (this)) {
+        img = this
+	obj = ""
+        iferr { wcsinfo (this) } then {
+            error (0, "Cannot determine image coords for `"//this//"'")
+        } else {
+            ra  = wcsinfo.midx
+            dec = wcsinfo.midy
+            size = max (wcsinfo.width, wcsinfo.height) / 60.0
+        }
+
+    } else { 
+        # Assume it's an object name, resolve to coords.
+        sesame (this)
+        if (sesame.status < 0)
+            error (0, "Cannot resolve object `" // this // "'")
+        lra  = sesame.ra
+        ldec = sesame.dec
+        lsr  = sr
+        img  = ""
+        obj  = this
+    }
+
+    filter = "" 
+    if (stype == "qso") {
+	filter = "Qso" 
+    } else if (stype == "galaxy") {
+	filter = "Galaxy" 
+    }
+    ;
+
+    # Create temp files.
+    stab = mktemp ("tmp$stab")
+    ttab = mktemp ("tmp$ttab")
+
+    # Display an image early on so we have something to look at.
+    if (display) {
+	if (img == "" || imaccess(img) == no) {
+	    # Display an image from the DSS.
+	    img = mktemp ("tmp$sspec") // ".fits"
+	    sv_img = img
+
+    	    # Get an image of the field if the user didn't supply one.
+    	    siap  = dalSiapSvc(regResolver("ivo://nasa.heasarc/skyview/dss1","siap"), lra,ldec,lsr,lsr)
+    
+    	    count = dalRecordCount(siap)
+    	    if (count <= 0)
+   	         error (0, "No optical images found")
+    	    s1 = dalGetData(siap, 0, osfn(img))
+
+            display (img//"[0]", 1, fill+, bpmask="", >& "dev$null")
+	    img = ""
+	} else {
+            display (img, 1, fill+, bpmask="", >& "dev$null")
+	}
+    }
+
+    
+    # Query for the spectra.
+    if (verb) {
+	printf ("#\n")
+	if (img != "")
+            printf ("# Search Params: Img='%s' RA=%.5f Dec=%.5f SR=%.2f\n",
+		img, lra, ldec, lsr)
+	else if (obj != "")
+            printf ("# Search Params: Obj='%s' RA=%.5f Dec=%.5f SR=%.2f\n",
+		obj, lra, ldec, lsr)
+	else
+            printf ("# Search Position: RA=%.5f Dec=%.5f SR=%.2f\n",
+		lra, ldec, lsr)
+	print ("#\n")
+    }
+    ssquery (image=img, object=obj, ra=lra, dec=ldec, sr=(lsr * 60.0), 
+	    zlow=z1, zhigh=z2, output="STDOUT", maxrec=maxrec, verb-) | \
+	match (filter, "STDIN", >& stab)
+
+    status = ssquery.status
+    if (ssquery.status < 0)  {
+	error (0, "No results found.")
+
+    } else if (img != "") {
+	wcsinfo (img)
+	lra = wcsinfo.midx		# decimal degrees
+	ldec = wcsinfo.midy
+
+    } else if (obj != "") {
+	sesame (obj)
+	lra = sesame.ra			# decimal degrees
+	ldec = sesame.dec
+    }
+    #printf("img='%s' obj='%s' ra=%.4f dec=%.4f sr=%.2f\n",img,obj,lra,ldec,lsr)
+
+
+
+    # See whether we display the image or simply download a cutout of the
+    # field we want.
+    if (display) {
+	coords = mktemp ("tmp$coords")
+	tvcoords = mktemp ("tmp$coords")
+
+	if (img == "" || imaccess(img) == no) {
+#	    # Display an image from the DSS.
+#	    img = mktemp ("tmp$sspec") // ".fits"
+	    img = sv_img
+#
+#    	    # Get an image of the field if the user didn't supply one.
+#    	    siap  = dalSiapSvc(regResolver("DSS2R","siap"), lra,ldec,lsr,lsr)
+#    
+#    	    count = dalRecordCount(siap)
+#    	    if (count <= 0)
+#   	         error (0, "No optical images found")
+#    	    s1 = dalGetData(siap, 0, osfn(img))
+#
+#            display (img//"[0]", 1, fill+, bpmask="", >& "dev$null")
+            fields (stab, "2,3", >& coords)
+            wcsctran (coords, tvcoords, img//"[0]", "world", "logical", 
+	        units="h n", columns="1 2", verb-)
+
+	    wcsinfo (img//"[0]")
+
+	    sv_img = img
+	    img    = img // "[0]"
+
+    	    if (!interact && imaccess(sv_img//".fits"))
+	        delete (sv_img//".fits", verify-, >& "dev$null")
+
+	} else {
+	    # Display the image.
+#            display (img, 1, fill+, bpmask="", >& "dev$null")
+            fields (stab, "2,3", >& coords)
+            wcsctran (coords, tvcoords, img, "world", "logical", 
+	        units="h n", verb-)
+	    sv_img = img
+	    wcsinfo (img)
+	}
+
+        if (verb) {
+     	    printf ("#\n")
+            if (img != "")
+     	        printf ("# Image: %s\n", img)
+            else if (obj != "")
+     	        printf ("# Object %s\n", obj)
+
+            if (img == "" && obj == "")
+       	        printf (   "# Search Pos: (%H,%h)  Size: %.2fx%.2f arcmin\n",
+	            ra, dec, lsr, lsr)
+            else
+       	        printf (   "# Search Pos: (%H,%h)  Size: %.2fx%.2f arcmin\n",
+	            wcsinfo.midx, wcsinfo.midy, wcsinfo.width, wcsinfo.height)
+       	    printf ("#\n")
+	    fields (stab, "1-5")
+        }
+
+	# Create local catalog of spectra.
+        joinlines (tvcoords//","//stab) | match ("-", "STDIN", stop+, >& ttab)
+
+    	if (access(coords))   delete (coords, verify-, >& "dev$null")
+    	if (access(tvcoords)) delete (tvcoords, verify-, >& "dev$null")
+
+
+        mtab = mktemp ("tmp$mtab")
+	fields (ttab, "1,2,3", >& mtab)
+type (mtab)
+        tvmark (1, mtab, mark="circle", radii=szmark, col=mkcolor,
+	    lab+, txsize=1)
+
+	if (access(mtab)) delete (mtab, verify-, >& "dev$null")
+
+        # Overlay the NED objects?
+        if (do_ned) 
+            nedoverlay (img)
+
+        # Draw the coordinate overlay grid?
+        if (do_grid) 
+            wcslab (img, 1, use-, fill+, overplot+, append+, 
+                labout-, dev="imdy")
+
+        if (interact) {
+	    # When verbose we printed the table above, otherwise we want some
+	    # kind of summary visible here.
+	    if (!verbose) 
+		type (ttab)
+
+	    print  ("\n\nIMSSPEC Command Summary:")
+	    print  ("b    Batch display w/in a box")
+	    print  ("d    Download spectrum")
+	    print  ("i    Print Info about spectrum")
+	    print  ("n    NED Overlay")
+	    print  ("r    Radio contour overlay")
+	    print  ("s    Show spectrum")
+	    print  ("w    WCS Overlay")
+	    print  ("q    Quit")
+	    printf ("\nCommand? ")
+
+	    while (fscan (imcur, x, y, wcs, cmd) != EOF) {
+		# Find nearest point.  All other commands will depend on the
+		# chosen spectrum.
+		d0 = 999999999.0
+		d1 = 0.0
+		found = no
+		llist = ttab
+	    	while (fscan(llist, xp,yp, name,ra,dec,z,class,line) != EOF) {
+		    d1 = sqrt ((xp-x)*(xp-x) + (yp-y)*(yp-y))
+		    if (d1 < d0) {
+			d0 = d1
+		 	sx = xp    ; sy = yp
+			sra = lra  ; sdec = ldec
+			sz = z
+			sname = name
+			sclass = class
+			found = yes
+		    }
+	    	}
+		llist = ""
+		if (!found)
+		    print ("Error: object not found")
+
+	  	print ("")
+	        switch (cmd) {
+	        case "i": 			# Show info on spectrum 
+	            {   printf ("%s  %H %h  %.5f  %s\n", 
+			    sname, (sra*15.), sdec, sz, sclass)
+	            }   
+
+	        case "b": 			# Get spectrum in the box
+	            {   printf ("Again...")
+			x1 = x
+			y1 = y
+			n = fscan (imcur, x2, y2, wcs, ch)
+			tnames = mktemp ("tmp$tname")
+			imlist = mktemp ("tmp$imlist")
+			tabclip (ttab, "STDOUT", "c1","c2", x1,y1, x2,y2) | \
+			    fields ("STDIN", "1,2,3,9", >& tnames)
+
+			# Draw a box around the selected region
+			printf ("%f %f 101 b\n%f %f 101 b\n", x1,y1,x2,y2) | \
+			tvmark (1, "", commands="STDIN",  lab-, col=200)
+
+			print ("\nDownloading images ....")
+			llist = tnames
+			while (fscan (llist, x, y, sname, sid) != EOF) {
+			    name = substr (sname, 1, 10) // ".fits"
+			    print (name, >>& imlist)
+			    printf ("%f %f %s\n", x, y, sname) | \
+			    tvmark (1, "STDIN", mark="circle", 
+				radii=szmark, col=205, lab+, txsize=1)
+			    if (imaccess(name) == no)
+			        ssget (sid, name, verbose+)
+			    else
+				printf ("Image '%s' exists already...\n", name)
+			}
+			llist = ""
+
+			printf ("Running SPECLOT ....")
+			specplot ("@"//imlist, label="imname",cursor="dev$null")
+
+			if (!saving)
+			    imdelete ("@"//imlist, >& "dev$null")
+			delete (imlist, verify-, >& "dev$null")
+			delete (tnames, verify-, >& "dev$null")
+			gflush
+
+	            }   
+	        case "d": 			# Get spectrum 
+	            {   
+			match (sname, stab) | scan (line)
+			print (line) | fields ("STDIN","1,7") | scan(sname,sid)
+			name = substr (sname, 1, 10) // ".fits"
+			printf ("%f %f %s\n", sx, sy, sname) | \
+			    tvmark (1, "STDIN", mark="circle", radii=szmark,
+				col=205, lab+, txsize=1)
+			if (imaccess(name) == no)
+			    ssget (sid, name, verbose+)
+	            }   
+	        case "o":			# Observation Overlay
+	            {   while (yes) {
+			    printf ("Mission (HST|XMM|Spitzer|Chandra|Rosat)? "
+			    i = scan (mission)
+			    if (i > 0)
+				break
+			}
+			obslogoverlay (img, mission)
+	            }
+	        case "r":			# NVSS Radio contour overlay
+	            {   radiooverlay (img)
+	            }
+	        case "n":			# NED Overlay
+	            {   nedoverlay (img)
+	            }
+	        case "w":			# WCS Overlay
+	            {   wcslab (img, 1, use-, fill+, overplot+, append+, 
+                	    labout-, dev="imdy")
+	            }
+	        case "s": 			# SPLOT spectrum 
+	            {   
+			match (sname, stab) | scan (line)
+			print (line) | fields ("STDIN","1,7") | scan(sname,sid)
+			name = substr (sname, 1, 10) // ".fits"
+			printf ("%f %f %s\n", sx, sy, sname) | \
+			    tvmark (1, "STDIN", mark="circle", radii=szmark,
+				col=205, lab+, txsize=1)
+			if (imaccess(name) == no)
+			    ssget (sid, name, verbose+)
+
+			if (imaccess(name) == yes) {
+			    printf ("Running SPLOT on image '%s'....", name)
+			    splot (name)
+			} else
+			    printf ("Error loading image '%s'\n", name)
+
+			if (!saving)
+			    imdelete (name, >& "dev$null")
+			gflush
+	            }
+	        case "q":			# Quit
+	            {   break
+	            }
+	        }
+
+	        printf ("\nCommand? ")
+	    }
+        }
+
+    } else {
+	type (stab)
+
+	if (getdata) {
+	    spectab (stab, imroot="", verbose+)
+	}
+    }
+
+
+    # Clean up
+    if (access(stab))	 delete (stab, verify-, >& "dev$null")
+    if (access(ttab))	 delete (ttab, verify-, >& "dev$null")
+
+    ra = INDEF		# reset params
+    dec = INDEF
+
+    if (interact && imaccess(sv_img//".fits"))
+	delete (sv_img//".fits", verify-, >& "dev$null")
+end
+
+