1 files changed, 112 insertions, 0 deletions
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