diff options
Diffstat (limited to 'noao/obsutil/src/pairmass')
| -rw-r--r-- | noao/obsutil/src/pairmass/airmass.x | 23 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/drawvector.x | 121 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/initmarker.x | 39 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/mkpkg | 19 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/pairmass.par | 40 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/t_pairmass.x | 112 | ||||
| -rw-r--r-- | noao/obsutil/src/pairmass/x_pairmass.x | 1 |
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 |