path: root/vo/src/OLD/sloanspec.cl

#{  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