Initial commit

7 files changed, 355 insertions, 0 deletions

diff --git a/noao/obsutil/src/pairmass/airmass.x b/noao/obsutil/src/pairmass/airmass.x
new file mode 100644
index 00000000..8e9c585d
--- /dev/null
+++ b/noao/obsutil/src/pairmass/airmass.x
@@ -0,0 +1,23 @@
+include	<math.h>
+
+# AIRMASS -- Compute airmass from DEC, LATITUDE and HA
+
+# Airmass formulation from Allen "Astrophysical Quantities" 1973 p.125,133.
+# and John Ball's book on Algorithms for the HP-45
+
+double procedure airmass (ha, dec, lat)
+
+double  ha, dec, lat, cos_zd, x
+
+define  SCALE   750.0d0                 # Atmospheric scale height
+
+begin
+        if (IS_INDEFD (ha) || IS_INDEFD (dec) || IS_INDEFD (lat))
+            call error (1, "Can't determine airmass")
+
+        cos_zd = sin(DEGTORAD(lat)) * sin(DEGTORAD(dec)) +
+                 cos(DEGTORAD(lat)) * cos(DEGTORAD(dec)) * cos(DEGTORAD(ha*15.))
+        x  = SCALE * cos_zd
+
+        return (sqrt (x**2 + 2*SCALE + 1) - x)
+end
diff --git a/noao/obsutil/src/pairmass/drawvector.x b/noao/obsutil/src/pairmass/drawvector.x
new file mode 100644
index 00000000..770689d3
--- /dev/null
+++ b/noao/obsutil/src/pairmass/drawvector.x
@@ -0,0 +1,121 @@
+# DRAW_VECTOR -- Draw the projected vector to the screen.
+
+include <gset.h>
+include	<mach.h>
+
+procedure draw_vector (def_title, timesys, xvec, yvec, n,
+	xmin, xmax, ymin, ymax)
+
+char	def_title[ARB]				#I default plot title
+char	timesys[ARB]				#I time system
+real	xvec[n], yvec[n]			#I vectors to plot
+int	n					#I npts in vectors
+real	xmin, xmax				#I x vector min & max
+real	ymin, ymax				#I y vector min & max
+
+pointer	sp, gp
+pointer	device, marker, xlabel, ylabel, title, suffix
+real	wx1, wx2, wy1, wy2, vx1, vx2, vy1, vy2, szm, tol
+int	mode, imark
+bool	pointmode
+
+pointer	gopen()
+real	clgetr()
+bool	clgetb(), streq()
+int	btoi(), clgeti()
+
+begin
+	call smark (sp)
+	call salloc (device, SZ_FNAME, TY_CHAR)
+	call salloc (marker, SZ_FNAME, TY_CHAR)
+	call salloc (xlabel, SZ_LINE,  TY_CHAR)
+	call salloc (ylabel, SZ_LINE,  TY_CHAR)
+	call salloc (title,  SZ_LINE,  TY_CHAR)
+	call salloc (suffix, SZ_FNAME, TY_CHAR)
+
+	call clgstr ("device", Memc[device], SZ_FNAME)
+	mode = NEW_FILE
+	if (clgetb ("append"))
+	    mode = APPEND
+
+	gp = gopen (Memc[device], mode, STDGRAPH)
+	tol = 10. * EPSILONR
+
+	if (mode != APPEND) {
+	    # Establish window.
+	    wx1 = clgetr ("wx1")
+	    wx2 = clgetr ("wx2")
+	    wy1 = clgetr ("wy1")
+	    wy2 = clgetr ("wy2")
+
+	    # Set window limits to defaults if not specified by user.
+	    if ((wx2 - wx1) < tol) {
+	        wx1 = xmin
+	        wx2 = xmax
+	    }
+
+	    if ((wy2 - wy1) < tol) {
+	        wy1 = ymin
+	        wy2 = ymax
+	    }
+
+	    call gswind (gp, wx1, wx2, wy1, wy2)
+    
+	    # Establish viewport.
+	    vx1 = clgetr ("vx1")
+	    vx2 = clgetr ("vx2")
+	    vy1 = clgetr ("vy1")
+	    vy2 = clgetr ("vy2")
+
+	    # Set viewport only if specified by user.
+	    if ((vx2 - vx1) > tol && (vy2 - vy1) > tol)
+	        call gsview (gp, vx1, vx2, vy1, vy2)
+	    else {
+		if (!clgetb ("fill"))
+		    call gseti (gp, G_ASPECT, 1)
+	    }
+   
+            call clgstr ("xlabel", Memc[xlabel], SZ_LINE)
+            call clgstr ("ylabel", Memc[ylabel], SZ_LINE)
+            call clgstr ("title",  Memc[title],  SZ_LINE)
+
+            if (streq (Memc[title], "default"))
+                call strcpy (def_title, Memc[title], SZ_LINE)
+            if (streq (Memc[xlabel], "default")) {
+                call sprintf (Memc[xlabel], SZ_LINE, "%s Time")
+		    call pargstr (timesys)
+	    }
+    
+	    call gseti (gp, G_XNMAJOR, clgeti ("majrx"))
+	    call gseti (gp, G_XNMINOR, clgeti ("minrx"))
+	    call gseti (gp, G_YNMAJOR, clgeti ("majry"))
+	    call gseti (gp, G_YNMINOR, clgeti ("minry"))
+
+	    call gseti (gp, G_ROUND, btoi (clgetb ("round")))
+
+	    if (clgetb ("logx"))
+	        call gseti (gp, G_XTRAN, GW_LOG)
+	    if (clgetb ("logy"))
+	        call gseti (gp, G_YTRAN, GW_LOG)
+
+	    # Draw axes using all this information.
+	    call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel])
+	}
+    
+	pointmode = clgetb ("pointmode")
+        if (pointmode) {
+            call clgstr ("marker", Memc[marker], SZ_FNAME)
+            szm = clgetr ("szmarker")
+            call init_marker (Memc[marker], imark)
+        }
+
+        # Now to actually draw the plot.
+        if (pointmode)
+            call gpmark (gp, xvec, yvec, n, imark, szm, szm)
+        else
+            call gpline (gp, xvec, yvec, n)
+       
+	call gflush (gp)
+        call gclose (gp)
+	call sfree (sp)
+end
diff --git a/noao/obsutil/src/pairmass/initmarker.x b/noao/obsutil/src/pairmass/initmarker.x
new file mode 100644
index 00000000..c506ecbb
--- /dev/null
+++ b/noao/obsutil/src/pairmass/initmarker.x
@@ -0,0 +1,39 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<gset.h>
+
+# INIT_MARKER -- Returns integers code for marker type string.
+
+procedure init_marker (marker, imark)
+
+char	marker[SZ_FNAME]	# Marker type as a string
+int	imark			# Integer code for marker - returned
+
+bool	streq()
+
+begin
+	if (streq (marker, "point"))
+	    imark = GM_POINT
+	else if (streq (marker,   "box"))
+	    imark = GM_BOX
+	else if (streq (marker,  "plus"))
+	    imark = GM_PLUS
+	else if (streq (marker, "cross"))
+	    imark = GM_CROSS
+	else if (streq (marker, "circle"))
+	    imark = GM_CIRCLE
+	else if (streq (marker, "hebar"))
+	    imark = GM_HEBAR
+	else if (streq (marker, "vebar"))
+	    imark = GM_VEBAR
+	else if (streq (marker, "hline"))
+	    imark = GM_HLINE
+	else if (streq (marker, "vline"))
+	    imark = GM_VLINE
+	else if (streq (marker, "diamond"))
+	    imark = GM_DIAMOND
+	else {
+	    call eprintf ("Unrecognized marker type, using 'box'\n")
+	    imark = GM_BOX
+	}
+end
diff --git a/noao/obsutil/src/pairmass/mkpkg b/noao/obsutil/src/pairmass/mkpkg
new file mode 100644
index 00000000..bfc564e5
--- /dev/null
+++ b/noao/obsutil/src/pairmass/mkpkg
@@ -0,0 +1,19 @@
+# Make the PAIRMASS task.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+standalone:
+	$update	libpkg.a
+	$omake	x_pairmass.x
+	$link	x_pairmass.o libpkg.a -lxtools -lasttools -o xx_pairmass.e
+	;
+
+libpkg.a:
+	airmass.x	<math.h>
+	drawvector.x	<gset.h> <mach.h>
+	initmarker.x	<gset.h>
+	t_pairmass.x
+	;
diff --git a/noao/obsutil/src/pairmass/pairmass.par b/noao/obsutil/src/pairmass/pairmass.par
new file mode 100644
index 00000000..1b0f28d6
--- /dev/null
+++ b/noao/obsutil/src/pairmass/pairmass.par
@@ -0,0 +1,40 @@
+ra,r,h,,0.,24.,Right Ascension of the object
+dec,r,h,,-90.,90.,Declination of the object
+epoch,r,h,INDEF,,,Epoch of the coordinates
+
+year,i,h,,,,Year of the observation
+month,i,h,,1,12,Numerical month specification
+day,i,h,,1,31,Day of the month
+
+observatory,s,h,"observatory",,,Observatory
+
+timesys,s,h,"Standard","|Universal|Standard|Siderial|",,Time system
+resolution,r,h,4,0.25,,Number of UT points per hour
+listout,b,h,no,,,"List, rather than plot, the airmass vs. UT"
+
+wx1,r,h,-7.,,,left user x-coord if not autoscaling
+wx2,r,h,7.,,,right user x-coord if not autoscaling
+wy1,r,h,0.,,,lower user y-coord if not autoscaling
+wy2,r,h,5.,,,upper user y-coord if not autoscaling
+pointmode,b,h,no,,,plot points instead of lines
+marker,s,h,"box",\
+  "point|box|plus|cross|circle|hebar|vebar|hline|vline|diamond",\
+  ,point marker character 
+szmarker,r,h,5E-3,,,marker size (0 for list input)
+logx,b,h,no,,,log scale x-axis
+logy,b,h,no,,,log scale y-axis
+xlabel,s,h,"default",,,x-axis label
+ylabel,s,h,"Airmass",,,y-axis label
+title,s,h,"default",,,title for plot
+vx1,r,h,0.,,,left limit of device window (ndc coords)
+vx2,r,h,0.,,,right limit of device window (ndc coords)
+vy1,r,h,0.,,,bottom limit of device window (ndc coords)
+vy2,r,h,0.,,,upper limit of device window (ndc coords)
+majrx,i,h,5,,,number of major divisions along x grid
+minrx,i,h,5,,,number of minor divisions along x grid
+majry,i,h,5,,,number of major divisions along y grid
+minry,i,h,5,,,number of minor divisions along y grid
+round,b,h,no,,,round axes to nice values
+fill,b,h,yes,,,fill device viewport regardless of aspect ratio?
+append,b,h,no,,,append to existing plot
+device,s,h,"stdgraph",,,output device
diff --git a/noao/obsutil/src/pairmass/t_pairmass.x b/noao/obsutil/src/pairmass/t_pairmass.x
new file mode 100644
index 00000000..7b9294a1
--- /dev/null
+++ b/noao/obsutil/src/pairmass/t_pairmass.x
@@ -0,0 +1,112 @@
+# PAIRMASS -- plot the airmass for a given RA & Dec on a given date.
+
+procedure t_pairmass()
+
+pointer	sp, observat, timesys, ut, air, title
+pointer	obs
+double	ra, dec, epoch, ra0, dec0, epoch0
+double	longitude, latitude, zone, st, ha, utd, resolution
+int	day, month, year, nsteps, i, tsys
+real	amin, amax
+
+pointer	obsopen()
+double	clgetd(), obsgetd(), ast_mst(), airmass()
+bool	clgetb()
+int	clgeti(), clgwrd()
+
+begin
+	call smark (sp)
+	call salloc (observat, SZ_FNAME, TY_CHAR)
+	call salloc (timesys, SZ_FNAME, TY_CHAR)
+
+	ra0 = clgetd ("ra")
+	dec0 = clgetd ("dec")
+	epoch0 = clgetd ("epoch")
+
+	year = clgeti ("year")
+	month = clgeti ("month")
+	day = clgeti ("day")
+
+	# Set time to plot.
+	tsys = clgwrd ("timesys", Memc[timesys], SZ_FNAME,
+	    "|Universal|Standard|Siderial|")
+
+	# Get observatory information.
+	call clgstr ("observatory", Memc[observat], SZ_FNAME)
+	obs = obsopen (Memc[observat])
+	#call obslog (obs, "PAIRMASS", "latitude longitude timezone", STDOUT)
+	call obsgstr (obs, "name", Memc[observat], SZ_FNAME)
+	longitude = obsgetd (obs, "longitude")
+	latitude = obsgetd (obs, "latitude")
+	zone = obsgetd (obs, "timezone")
+	call obsclose (obs)
+
+	resolution = clgetd ("resolution")
+	nsteps = nint (24 * resolution) + 1
+
+	call salloc (ut, 3*nsteps, TY_REAL)
+	call salloc (air, 3*nsteps, TY_REAL)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	call ast_date_to_epoch (year, month, day, 12.d0, epoch)
+	call ast_precess (ra0, dec0, epoch0, ra, dec, epoch)
+
+	do i = 1, nsteps {
+	    utd = (i-1) / resolution
+	    call ast_date_to_epoch (year, month, day, utd, epoch)
+	    st = ast_mst (epoch, longitude)
+	    ha = st - ra
+
+	    switch (tsys) {
+	    case 1:
+		Memr[ut+i-1] = utd
+	    case 2:
+		if (utd < zone)
+		    Memr[ut+i-1] = utd - zone + 24
+		else
+		    Memr[ut+i-1] = utd - zone
+	    case 3:
+		Memr[ut+i-1] = st
+	    }
+	    Memr[air+i-1] = real (airmass (ha, dec, latitude))
+
+	    Memr[ut+i-1+nsteps] = Memr[ut+i-1] - 24
+	    Memr[ut+i-1+2*nsteps] = Memr[ut+i-1] + 24
+	    Memr[air+i-1+nsteps] = Memr[air+i-1]
+	    Memr[air+i-1+2*nsteps] = Memr[air+i-1]
+	}
+	call xt_sort2 (Memr[ut], Memr[air], 3*nsteps)
+	call alimr (Memr[air], 3*nsteps, amin, amax)
+
+	call sprintf (Memc[title], SZ_LINE,
+	    "Airmass for %d/%d/%d\n%s\nRA=%h, Dec=%h (%g)")
+	    call pargi (month)
+	    call pargi (day)
+	    call pargi (year)
+	    call pargstr (Memc[observat])
+	    call pargd (ra0)
+	    call pargd (dec0)
+	    call pargd (epoch0)
+
+	if (clgetb ("listout")) {
+	    call printf ("%s\nTime System=%s\n\n")
+		call pargstr (Memc[title])
+		call pargstr (Memc[timesys])
+	    amax = clgetd ("wy2")
+	    if (amax < 1)
+		amax = 5
+	    do i = 1, 3*nsteps {
+		if (Memr[ut+i-1] < 0. || Memr[ut+i-1] >= 24.)
+		    next
+		if (Memr[air+i-1] > amax)
+		    next
+		call printf ("%6.0m\t%8.4f\n")
+		    call pargr (Memr[ut+i-1])
+		    call pargr (Memr[air+i-1])
+	    }
+	} else
+	    call draw_vector (Memc[title], Memc[timesys], Memr[ut], Memr[air],
+		3*nsteps, 0., 24., amin, amax)
+
+	call sfree (sp)
+end
diff --git a/noao/obsutil/src/pairmass/x_pairmass.x b/noao/obsutil/src/pairmass/x_pairmass.x
new file mode 100644
index 00000000..82f80509
--- /dev/null
+++ b/noao/obsutil/src/pairmass/x_pairmass.x
@@ -0,0 +1 @@
+task	pairmass = t_pairmass