From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/digiphot/apphot/center/apcplot.x | 304 ++++++++++++++++++++++++++++++++++ 1 file changed, 304 insertions(+) create mode 100644 noao/digiphot/apphot/center/apcplot.x (limited to 'noao/digiphot/apphot/center/apcplot.x') diff --git a/noao/digiphot/apphot/center/apcplot.x b/noao/digiphot/apphot/center/apcplot.x new file mode 100644 index 00000000..6e2845c3 --- /dev/null +++ b/noao/digiphot/apphot/center/apcplot.x @@ -0,0 +1,304 @@ +include +include +include +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +# AP_CPLOT -- Procedure to compute radial profile plots for the centering +# routine. + +procedure ap_cplot (ap, sid, gd, makeplot) + +pointer ap # the pointer to the apphot structure +int sid # the output sequence number of the star +pointer gd # the graphics stream +int makeplot # the make a plot ? + +int nx, ny +pointer ctr, sp, str, r, gt +real xcenter, ycenter, xc, yc, rmin, rmax, imin, imax +real u1, u2, v1, v2, x1, x2, y1, y2 +int apstati() +pointer ap_gtinit() +real apstatr() + +begin + # Check for enabled graphics stream. + if (gd == NULL || makeplot == NO) + return + + # Check for defined center. + xcenter = apstatr (ap, XCENTER) + ycenter = apstatr (ap, YCENTER) + if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter)) + return + + # Check for a centering algorithm. + if (apstati (ap, CENTERFUNCTION) == AP_NONE) + return + + # Get the pixel buffer parameters. + ctr = AP_PCENTER(ap) + nx = AP_CNX(ctr) + ny = AP_CNY(ctr) + if (AP_CTRPIX(ctr) == NULL || nx <= 0 || ny <= 0) + return + xc = AP_CXC(ctr) + (xcenter - apstatr (ap, CXCUR)) + yc = AP_CYC(ctr) + (ycenter - apstatr (ap, CYCUR)) + + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (r, nx * ny, TY_REAL) + + # Compute the radii and the plot limits. + call ap_ijtor2 (Memr[r], nx, ny, xc, yc) + call alimr (Memr[r], nx * ny, rmin, rmax) + call alimr (Memr[AP_CTRPIX(ctr)], nx * ny, imin, imax) + + # Reactivate the work station. + call greactivate (gd, 0) + + # Save old viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Initialize the plot. + call apstats (ap, IMROOT, Memc[str], SZ_LINE) + call sprintf (Memc[str], SZ_LINE, "%s Star %d") + call pargstr (Memc[str]) + call pargi (sid) + gt = ap_gtinit (Memc[str], apstatr (ap,OXINIT), apstatr (ap, OYINIT)) + + # Make the plot. + call gclear (gd) + call ap_cpset (gd, gt, ap, rmin, rmax, imin, imax) + if (apstati (ap, POSITIVE) == YES) + call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny, + rmin, rmax, imin + EPSILONR, imax, "plus") + else + call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny, + rmin, rmax, imin, imax - EPSILONR, "plus") + call ap_cpreset (gd, gt, ap, rmin, rmax, imin, imax) + call ap_cpannotate (gd, ap) + + # Restore the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + # Free the space. + call ap_gtfree (gt) + call gdeactivate (gd, 0) + call sfree (sp) +end + + +# AP_CPSET -- Procedure to set up the parameters for the center radial profile +# plot. + +procedure ap_cpset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # the graphics stream +pointer gt # the gtools pointer +pointer ap # the apphot pointer +real xmin, xmax # the minimum and maximum radial distance +real ymin, ymax # the minimum and maximum of y axis (ymin not used) + +int fd +pointer sp, str, title +real scale, aspect, datalimit, skysigma, threshold, vx1, vx2, vy1, vy2 +real cthreshold +int stropen(), apstati() +real apstatr(), gstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + + # Encode the parameter string. + fd = stropen (Memc[str], SZ_LINE, WRITE_ONLY) + call sysid (Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + call fprintf (fd, + "Center: xc=%0.2f yc=%0.2f xerr=%0.2f yerr=%0.2f\n") + call pargr (apstatr (ap, OXCENTER)) + call pargr (apstatr (ap, OYCENTER)) + call pargr (apstatr (ap, XERR)) + call pargr (apstatr (ap, YERR)) + call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + call strclose (fd) + + # Store some default plotting parameters. + scale = apstatr (ap, SCALE) + aspect = gstatr (gd, G_ASPECT) + call gsetr (gd, G_ASPECT, 0.75) + + # Set the labels and window. + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + if (apstati (ap, POSITIVE) == YES) { + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, datalimit + threshold, + datalimit + threshold + ymax) + call glabax (gd, Memc[str], "", "Intensity") + call gseti (gd, G_YDRAWAXES, 0) + call gseti (gd, G_XDRAWAXES, 1) + call ggview (gd, vx1, vx2, vy1, vy2) + call gswind (gd, vx1, vx2, vy1, vy2) + call gswind (gd, xmin, xmax, datalimit + threshold, + datalimit + threshold + ymax) + call glabax (gd, "", + "Radial Distance (lower-pixels, upper-scale units)", "") + } else { + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, datalimit - + threshold - ymax, datalimit - threshold) + call glabax (gd, Memc[str], "", "Intensity") + call gseti (gd, G_YDRAWAXES, 0) + call gseti (gd, G_XDRAWAXES, 1) + call ggview (gd, vx1, vx2, vy1, vy2) + call gswind (gd, vx1, vx2, vy1, vy2) + call gswind (gd, xmin, xmax, datalimit - threshold - ymax, + datalimit - threshold) + call glabax (gd, "", + "Radial Distance (lower-pixels, upper-scale units)", "") + } + + # Set the window up for plotting the data. + call gseti (gd, G_YDRAWAXES, 3) + call gseti (gd, G_XDRAWAXES, 3) + call gsetr (gd, G_ASPECT, aspect) + call gt_sets (gt, GTTYPE, "mark") + call gt_setr (gt, GTXMIN, xmin) + call gt_setr (gt, GTXMAX, xmax) + if (apstati (ap, POSITIVE) == YES) { + call gt_setr (gt, GTYMIN, ymin) + call gt_setr (gt, GTYMAX, ymax) + } else { + call gt_setr (gt, GTYMIN, ymin) + call gt_setr (gt, GTYMAX, ymax) + } + call gt_swind (gd, gt) + + call sfree (sp) +end + + +# AP_CPANNOTATE -- Procedure to annotate the radial plot in center. + +procedure ap_cpannotate (gd, ap) + +pointer gd # the graphics stream +pointer ap # the apphot structure + +pointer sp, str +real fwhmpsf, capert, datalimit, threshold, skysigma, cthreshold +real xmin, xmax, ymin, ymax +int apstati() +real apstatr() + +begin + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call gseti (gd, G_PLTYPE, GL_DASHED) + call ggwind (gd, xmin, xmax, ymin, ymax) + + fwhmpsf = 0.5 * apstatr (ap, FWHMPSF) * apstatr (ap, SCALE) + capert = 2.0 * fwhmpsf * apstatr (ap, CAPERT) + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + if (apstati (ap, POSITIVE) == YES) + threshold = datalimit + threshold + else + threshold = datalimit - threshold + + # Plot the full width half maximum of the radial profile. + if (fwhmpsf >= xmin && fwhmpsf <= xmax) { + call gamove (gd, fwhmpsf, ymin) + call gadraw (gd, fwhmpsf, ymax) + call sprintf (Memc[str], SZ_LINE, "hwhm = %0.2f") + call pargr (fwhmpsf) + call gtext (gd, fwhmpsf, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the centering aperture. + if (capert >= xmin && capert <= xmax) { + call gamove (gd, capert, ymin) + call gadraw (gd, capert, ymax) + call sprintf (Memc[str], SZ_LINE, "cbox half-width = %0.2f") + call pargr (capert) + call gtext (gd, capert, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the threshold level for centering. + call sprintf (Memc[str], SZ_LINE, "threshold = %g") + call pargr (threshold) + call gtext (gd, xmin, ymin, Memc[str], "q=h") + + # Mark the sky sigma if defined. + if (! IS_INDEFR(skysigma) && skysigma >= ymin && skysigma <= ymax) { + call gmark (gd, (xmin + xmax) / 2.0, (ymin + ymax) / 2.0, + GM_VEBAR, -0.25, -skysigma) + call sprintf (Memc[str], SZ_LINE, "sigma = %g") + call pargr (skysigma) + call gtext (gd, (xmin + xmax) / 2.0, (ymin + ymax + skysigma) / 2.0, + Memc[str], "q=h;h=c") + } + + call sfree (sp) +end + + +# AP_CPRSET -- Procedure to reset the plot window after the data points +# have been plotted. + +procedure ap_cpreset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # the graphics stream +pointer gt # the gtools pointer +pointer ap # the apphot pointer +real xmin, xmax # the minimum and maximum radial distance +real ymin, ymax # the minimum and maximum of the y axis (ymin not used) + +real datalimit, skysigma, threshold, cthreshold +int apstati() +real apstatr() + +begin + # Set the data window. + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + call gt_setr (gt, GTXMIN, xmin) + call gt_setr (gt, GTXMAX, xmax) + if (apstati (ap, POSITIVE) == YES) { + call gt_setr (gt, GTYMIN, datalimit + threshold) + call gt_setr (gt, GTYMAX, ymax + datalimit + threshold) + } else { + call gt_setr (gt, GTYMIN, datalimit - ymax - threshold) + call gt_setr (gt, GTYMAX, datalimit - threshold) + } + call gt_swind (gd, gt) +end -- cgit