aboutsummaryrefslogtreecommitdiff
path: root/noao/nproto/t_linpol.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/nproto/t_linpol.x')
-rw-r--r--noao/nproto/t_linpol.x547
1 files changed, 547 insertions, 0 deletions
diff --git a/noao/nproto/t_linpol.x b/noao/nproto/t_linpol.x
new file mode 100644
index 00000000..d41270ec
--- /dev/null
+++ b/noao/nproto/t_linpol.x
@@ -0,0 +1,547 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <error.h>
+
+define MAX_IMAGES 4
+
+define LP_TITLE "Linear polarization image"
+define LP_IMKEY "POL" # keyword prefix for input images
+
+define LP_PBAND 1
+define LP_PKEY "POLAR"
+define LP_PSTR "Band 1 is the percent polarization"
+
+define LP_ABAND 2
+define LP_AKEY "ANGLE"
+define LP_ASTR "Band 2 is the polarization angle"
+
+define LP_IBAND 3
+define LP_IKEY "I-STOKES"
+define LP_ISTR "Band 3 is the Stokes I parameter"
+
+define LP_QBAND 4
+define LP_QKEY "Q-STOKES"
+define LP_QSTR "Band 4 is the Stokes Q parameter"
+define LP_QSTRN "Band 4 is the normalized Stokes Q parameter"
+
+define LP_UBAND 5
+define LP_UKEY "U-STOKES"
+define LP_USTR "Band 5 is the Stokes U parameter"
+define LP_USTRN "Band 5 is the normalized Stokes U parameter"
+
+
+# LINPOL -- Calculate the percent polarization and the polarization
+# angle images for the simplest linear polarization cases, 0-45-90 or
+# 0-45-90-135 polarizer positions.
+
+procedure t_linpol ()
+
+pointer inlist, output, in[MAX_IMAGES], out, key, sp
+bool dflag, sflag, nflag
+int len
+
+int imtopenp(), imtlen()
+bool clgetb()
+
+errchk lp_map, lp_polarize
+
+begin
+ call smark (sp)
+ call salloc (output, SZ_FNAME, TY_CHAR)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+
+ # get the input image template
+ inlist = imtopenp ("input")
+ len = imtlen (inlist)
+ if (len != 3 && len != 4) {
+ call imtclose (inlist)
+ call sfree (sp)
+ call error (1, "Must supply either three or four images.")
+ }
+
+ # get the output image stack name
+ call clgstr ("output", Memc[output], SZ_FNAME)
+
+ sflag = clgetb ("stokes")
+ dflag = clgetb ("degrees")
+ nflag = clgetb ("normalize")
+
+ # keyword for polarizer angle - UPPERcase for neatness
+ call clgstr ("keyword", Memc[key], SZ_FNAME)
+ call strupr (Memc[key])
+
+ iferr {
+ # pass the number of possible frames (4) explicitly in
+ # hopes of later relaxing the 45 degree restriction
+ call lp_map (inlist, in, MAX_IMAGES,
+ Memc[key], Memc[output], out, sflag, nflag)
+ call lp_polarize (in, MAX_IMAGES, out, dflag, sflag, nflag)
+ } then {
+ call lp_unmap (in, MAX_IMAGES, out)
+ call imtclose (inlist)
+ call sfree (sp)
+ call erract (EA_ERROR)
+ }
+
+ call lp_unmap (in, MAX_IMAGES, out)
+ call imtclose (inlist)
+ call sfree (sp)
+end
+
+
+# LP_MAP -- map the set of input images.
+
+procedure lp_map (inlist, in, nin, key, output, out, sflag, nflag)
+
+pointer inlist #I input image template
+pointer in[nin] #O input image descriptor array
+int nin #I size of the array (4)
+char key[ARB] #I keyword for polarizer angle
+char output[ARB] #I output image name
+pointer out #O output image descriptor
+bool sflag #I include stokes frames in output?
+bool nflag #I normalize the stokes frames?
+
+pointer input, im_tmp, sp
+real pol
+int i, j, ipol, ndim
+long axis[IM_MAXDIM]
+bool firsttime
+
+int imtgetim()
+real imgetr()
+pointer immap()
+bool fp_equalr()
+
+errchk immap, imgetr, imdelf
+
+begin
+ # for graceful error recovery
+ im_tmp = NULL
+ out = NULL
+ do i = 1, nin
+ in[i] = NULL
+
+
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+
+ iferr {
+ while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+ im_tmp = immap (Memc[input], READ_ONLY, 0)
+
+ if (IM_NDIM(im_tmp) > 2)
+ call error (1, "only 1 or 2 dimensional images allowed")
+
+ pol = imgetr (im_tmp, key)
+
+ if (pol < 0 || pol > 135 || mod (nint(pol), 45) != 0 ||
+ ! fp_equalr (pol, real(nint(pol)))) {
+ call eprintf ("image %s, %s must be 0,45,90,135 degrees\n")
+ call pargstr (Memc[input])
+ call pargstr (key)
+ call flush (STDERR)
+ call error (1, "task LINPOL")
+ }
+
+ # index into in pointer array
+ ipol = max (1, min (nin, 1 + int(pol) / 45))
+
+ if (in[ipol] == NULL) {
+ in[ipol] = im_tmp
+ im_tmp = NULL
+ } else {
+ call eprintf ("multiple images specified at %d degrees\n")
+ call pargi ((ipol-1) * 45)
+ call flush (STDERR)
+ call error (1, "task JOIN")
+ }
+ }
+
+ # check dimensionality
+ firsttime = true
+ do i = 1, nin {
+ if (in[i] == NULL)
+ next
+
+ if (firsttime) {
+ ndim = IM_NDIM(in[i])
+ do j = 1, IM_MAXDIM
+ axis[j] = IM_LEN(in[i],j)
+ firsttime = false
+ next
+ }
+
+ if (IM_NDIM(in[i]) != ndim)
+ call error (1, "images are different sizes")
+
+ do j = 1, ndim
+ if (IM_LEN(in[i],j) != axis[j])
+ call error (1, "images are different sizes")
+ }
+
+ # create the output polarization (hyper) cube
+ # just copy header from first image available
+ do i = 1, nin
+ if (in[i] != NULL) {
+ out = immap (output, NEW_COPY, in[i])
+ break
+ }
+
+ # increase the image's girth
+ IM_NDIM(out) = ndim + 1
+ for (i=1; i <= ndim; i=i+1)
+ IM_LEN(out,i) = axis[i]
+
+ if (sflag)
+ IM_LEN(out,i) = 5
+ else
+ IM_LEN(out,i) = 2
+
+ call strcpy (LP_TITLE, IM_TITLE(out), SZ_IMTITLE)
+
+ # delete the polarizer angle keyword
+ call imdelf (out, key)
+
+ # add keywords naming the input images
+ do i = 1, nin {
+ if (in[i] == NULL)
+ next
+
+ call sprintf (Memc[input], SZ_FNAME, "%s%d")
+ call pargstr (LP_IMKEY)
+ call pargi (45*(i-1))
+
+ call imastr (out, Memc[input], IM_HDRFILE(in[i]))
+ }
+
+ # add keywords to index output frames
+ call imastr (out, LP_PKEY, LP_PSTR)
+ call imastr (out, LP_AKEY, LP_ASTR)
+
+ if (sflag) {
+ call imastr (out, LP_IKEY, LP_ISTR)
+ if (nflag) {
+ call imastr (out, LP_QKEY, LP_QSTRN)
+ call imastr (out, LP_UKEY, LP_USTRN)
+ } else {
+ call imastr (out, LP_QKEY, LP_QSTR)
+ call imastr (out, LP_UKEY, LP_USTR)
+ }
+ }
+
+ } then {
+ if (im_tmp != NULL)
+ call imunmap (im_tmp)
+ call sfree (sp)
+ call erract (EA_ERROR)
+ }
+
+ # start off with a clean slate
+ call imflush (out)
+ call sfree (sp)
+end
+
+
+# LP_UNMAP -- unmap the set of input images.
+
+procedure lp_unmap (in, nin, out)
+
+pointer in[nin] #U input image pointer array
+int nin #I size of the array (4)
+pointer out #U output image pointer
+
+int i
+
+begin
+ do i = 1, nin
+ if (in[i] != NULL)
+ call imunmap (in[i])
+
+ if (out != NULL)
+ call imunmap (out)
+end
+
+
+# LP_POLARIZE -- calculate the polarization given at least 3 of the 4
+# possible frames taken with the polarizer at 45 degree increments.
+
+procedure lp_polarize (in, nin, out, dflag, sflag, nflag)
+
+pointer in[nin] #I input image pointer array
+int nin #I size of the array (4)
+pointer out #I output image pointer
+bool dflag #I report the angle in degrees?
+bool sflag #I include stokes frames in output?
+bool nflag #I normalize the stokes frames?
+
+pointer ibuf, qbuf, ubuf, sp
+pointer buf1, buf2, buf3, buf4
+long v1[IM_MAXDIM], v2[IM_MAXDIM], v3[IM_MAXDIM], v4[IM_MAXDIM]
+int line, npix, skip, i
+
+int imgnlr()
+
+begin
+ npix = IM_LEN(out,1)
+
+ call smark (sp)
+ call salloc (ibuf, npix, TY_REAL)
+ call salloc (qbuf, npix, TY_REAL)
+ call salloc (ubuf, npix, TY_REAL)
+
+ call amovkl (long(1), v1, IM_MAXDIM)
+ call amovkl (long(1), v2, IM_MAXDIM)
+ call amovkl (long(1), v3, IM_MAXDIM)
+ call amovkl (long(1), v4, IM_MAXDIM)
+
+ # choose the combining scheme
+ skip = 0
+ do i = 1, nin
+ if (in[i] == NULL) {
+ skip = i
+ break
+ }
+
+ # not worth generalizing the method, just duplicate the code...
+ switch (skip) {
+
+ case 0:
+ # I = (im0 + im45 + im90 + im135) / 4
+ # Q = (im0 - im90) / 2
+ # U = (im45 - im135) / 2
+ while (imgnlr (in[1], buf1, v1) != EOF &&
+ imgnlr (in[2], buf2, v2) != EOF &&
+ imgnlr (in[3], buf3, v3) != EOF &&
+ imgnlr (in[4], buf4, v4) != EOF) {
+
+ call aaddr (Memr[buf1], Memr[buf2], Memr[ibuf], npix)
+ call aaddr (Memr[buf3], Memr[ibuf], Memr[ibuf], npix)
+ call aaddr (Memr[buf4], Memr[ibuf], Memr[ibuf], npix)
+ call adivkr (Memr[ibuf], 4., Memr[ibuf], npix)
+
+ call asubr (Memr[buf1], Memr[buf3], Memr[qbuf], npix)
+ call adivkr (Memr[qbuf], 2., Memr[qbuf], npix)
+
+ call asubr (Memr[buf2], Memr[buf4], Memr[ubuf], npix)
+ call adivkr (Memr[ubuf], 2., Memr[ubuf], npix)
+
+ line = int(v1[2]) - 1
+ call lp_stokes (Memr[ibuf], Memr[qbuf], Memr[ubuf],
+ npix, out, line, dflag, sflag, nflag)
+ }
+
+ case 1:
+ # I = (im45 + im135) / 2
+ # Q = I - im90
+ # U = (im45 - im135) / 2
+ while (imgnlr (in[2], buf2, v2) != EOF &&
+ imgnlr (in[3], buf3, v3) != EOF &&
+ imgnlr (in[4], buf4, v4) != EOF) {
+
+ call aaddr (Memr[buf2], Memr[buf4], Memr[ibuf], npix)
+ call adivkr (Memr[ibuf], 2., Memr[ibuf], npix)
+
+ call asubr (Memr[ibuf], Memr[buf3], Memr[qbuf], npix)
+
+ call asubr (Memr[buf2], Memr[buf4], Memr[ubuf], npix)
+ call adivkr (Memr[ubuf], 2., Memr[ubuf], npix)
+
+ line = int(v2[2]) - 1
+ call lp_stokes (Memr[ibuf], Memr[qbuf], Memr[ubuf],
+ npix, out, line, dflag, sflag, nflag)
+ }
+
+ case 2:
+ # I = (im0 + im90) / 2
+ # Q = (im0 - im90) / 2
+ # U = I - im135
+ while (imgnlr (in[1], buf1, v1) != EOF &&
+ imgnlr (in[3], buf3, v3) != EOF &&
+ imgnlr (in[4], buf4, v4) != EOF) {
+
+ call aaddr (Memr[buf1], Memr[buf3], Memr[ibuf], npix)
+ call adivkr (Memr[ibuf], 2., Memr[ibuf], npix)
+
+ call asubr (Memr[buf1], Memr[buf3], Memr[qbuf], npix)
+ call adivkr (Memr[qbuf], 2., Memr[qbuf], npix)
+
+ call asubr (Memr[ibuf], Memr[buf4], Memr[ubuf], npix)
+
+ line = int(v1[2]) - 1
+ call lp_stokes (Memr[ibuf], Memr[qbuf], Memr[ubuf],
+ npix, out, line, dflag, sflag, nflag)
+ }
+
+ case 3:
+ # I = (im45 + im135) / 2
+ # Q = im0 - I
+ # U = (im45 - im135) / 2
+ while (imgnlr (in[1], buf1, v1) != EOF &&
+ imgnlr (in[2], buf2, v2) != EOF &&
+ imgnlr (in[4], buf4, v4) != EOF) {
+
+ call aaddr (Memr[buf2], Memr[buf4], Memr[ibuf], npix)
+ call adivkr (Memr[ibuf], 2., Memr[ibuf], npix)
+
+ call asubr (Memr[buf1], Memr[ibuf], Memr[qbuf], npix)
+
+ call asubr (Memr[buf2], Memr[buf4], Memr[ubuf], npix)
+ call adivkr (Memr[ubuf], 2., Memr[ubuf], npix)
+
+ line = int(v1[2]) - 1
+ call lp_stokes (Memr[ibuf], Memr[qbuf], Memr[ubuf],
+ npix, out, line, dflag, sflag, nflag)
+ }
+
+ case 4:
+ # I = (im0 + im90) / 2
+ # Q = (im0 - im90) / 2
+ # U = im45 - I
+ while (imgnlr (in[1], buf1, v1) != EOF &&
+ imgnlr (in[2], buf2, v2) != EOF &&
+ imgnlr (in[3], buf3, v3) != EOF) {
+
+ call aaddr (Memr[buf1], Memr[buf3], Memr[ibuf], npix)
+ call adivkr (Memr[ibuf], 2., Memr[ibuf], npix)
+
+ call asubr (Memr[buf1], Memr[buf3], Memr[qbuf], npix)
+ call adivkr (Memr[qbuf], 2., Memr[qbuf], npix)
+
+ call asubr (Memr[buf2], Memr[ibuf], Memr[ubuf], npix)
+
+ line = int(v1[2]) - 1
+ call lp_stokes (Memr[ibuf], Memr[qbuf], Memr[ubuf],
+ npix, out, line, dflag, sflag, nflag)
+ }
+
+ }
+
+ call sfree(sp)
+end
+
+
+# LP_STOKES -- calculate the fractional polarization and angle for a
+# specific line (from the stokes parameters) and output the results.
+
+procedure lp_stokes (i, q, u, npix, out, line, dflag, sflag, nflag)
+
+real i[ARB] #I Stokes I vector
+real q[ARB] #I Stokes Q vector
+real u[ARB] #I Stokes U vector
+int npix #I length of the vectors
+pointer out #I output image descriptor
+int line #I line number
+bool dflag #I convert to degrees?
+bool sflag #I include stokes frames in output?
+bool nflag #I normalize the stokes frames?
+
+pointer pbuf, abuf, sp
+
+pointer impl3r()
+real lp_errfcn()
+extern lp_errfcn
+
+begin
+ call smark (sp)
+ call salloc (pbuf, npix, TY_REAL)
+ call salloc (abuf, npix, TY_REAL)
+
+ call lp_pol (i, q, u, Memr[pbuf], npix)
+ call lp_ang (q, u, Memr[abuf], npix, dflag)
+
+ call amovr (Memr[pbuf], Memr[impl3r (out, line, LP_PBAND)], npix)
+ call amovr (Memr[abuf], Memr[impl3r (out, line, LP_ABAND)], npix)
+
+ if (sflag) {
+ call amovr (i, Memr[impl3r (out, line, LP_IBAND)], npix)
+ if (nflag) {
+ call advzr (q, i, q, npix, lp_errfcn)
+ call advzr (u, i, u, npix, lp_errfcn)
+ }
+ call amovr (q, Memr[impl3r (out, line, LP_QBAND)], npix)
+ call amovr (u, Memr[impl3r (out, line, LP_UBAND)], npix)
+ }
+
+ call sfree (sp)
+end
+
+
+# LP_POL -- calculate the fractional linear polarization for a vector,
+# given the stokes I, Q, and U vectors.
+
+procedure lp_pol (i, q, u, p, npix)
+
+real i[ARB] #I Stokes I vector
+real q[ARB] #I Stokes Q vector
+real u[ARB] #I Stokes U vector
+real p[ARB] #O fractional polarization vector
+int npix #I length of the vectors
+
+pointer tmp, sp
+
+real lp_errfcn()
+extern lp_errfcn
+
+begin
+ call smark (sp)
+ call salloc (tmp, npix, TY_REAL)
+
+ call amulr (q, q, p, npix)
+ call amulr (u, u, Memr[tmp], npix)
+ call aaddr (p, Memr[tmp], p, npix)
+ call asqrr (p, p, npix, lp_errfcn)
+ call advzr (p, i, p, npix, lp_errfcn)
+
+ call sfree (sp)
+end
+
+
+# LP_ERRFCN -- error function for the square root of negative numbers.
+
+real procedure lp_errfcn (x)
+
+real x
+
+begin
+ return (0.)
+end
+
+
+# LP_ANG -- calculate the polarization angle, given the Stokes params.
+
+procedure lp_ang (q, u, a, npix, dflag)
+
+real q[ARB] #I Stokes Q vector
+real u[ARB] #I Stokes U vector
+real a[ARB] #O polarization angle vector
+int npix #I length of the vectors
+bool dflag #I convert to degrees?
+
+define PI 3.14159265358979
+define RAD2DEG (180./PI)
+
+begin
+ call lp_aatn2r (u, q, a, npix)
+ call adivkr (a, 2., a, npix)
+
+ if (dflag)
+ call amulkr (a, RAD2DEG, a, npix)
+end
+
+
+# LP_AATN2R -- calculate the arctangent in the proper quadrant.
+
+procedure lp_aatn2r (y, x, a, npix)
+
+real y[ARB], x[ARB] #I numerator and denominator, respectively
+real a[ARB] #O arctangent vector (radians)
+
+int npix, i
+
+begin
+ do i = 1, npix {
+ a[i] = atan2 (y[i], x[i])
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 16:21:00 -0400