; RMFIT ;--------------------------------------------------------------- ;! Fits 1-dimensional polarization spectrum to Q/U cube ;# TASK ANALYSIS PLOT SPECTRAL ONED MODELING TV-APPL POLARIZATION ;----------------------------------------------------------------------- ;; Copyright (C) 2013 ;; Associated Universities, Inc. Washington DC, USA. ;; ;; This program is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation; either version 2 of ;; the License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public ;; License along with this program; if not, write to the Free ;; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, ;; MA 02139, USA. ;; ;; Correspondence concerning AIPS should be addressed as follows: ;; Internet email: aipsmail@nrao.edu. ;; Postal address: AIPS Project Office ;; National Radio Astronomy Observatory ;; 520 Edgemont Road ;; Charlottesville, VA 22903-2475 USA ;----------------------------------------------------------------------- RMFIT LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC RMFIT Fits polarization spectrum to Q/U cubes. INNAME Input Q image name (name) INCLASS Input Q image name (class) INSEQ 0.0 9999.0 Input Q image name (seq. #) INDISK 0.0 9.0 Input Q image disk unit # INVERS 0.0 RM table version number IN2NAME Input U image name (name) IN2CLASS Input U image name (class) IN2SEQ 0.0 9999.0 Input U image name (seq. #) IN2DISK 0.0 9.0 Input U image disk unit # IN3NAME Input I image name (name) IN3CLASS Input I image name (class) IN3SEQ 0.0 9999.0 Input I image name (seq. #) IN3DISK 0.0 9.0 Input I image disk unit # IN4NAME FARS real image name (name) IN4CLASS FARS real image name (class) IN4SEQ 0.0 9999.0 FARS real image name (seq. #) IN4DISK 0.0 9.0 FARS real image disk unit # IN5NAME FARS imag image name (name) IN5CLASS FARS imag image name (class) IN5SEQ 0.0 9999.0 FARS imag image name (seq. #) IN5DISK 0.0 9.0 FARS imag image disk unit # OUTNAME Output image name (name) OUTSEQ -1.0 9999.0 Output image name (seq. #) OUTDISK 0.0 9.0 Output image disk unit #. BLC Bottom left corner of input TRC Top right corner of input YINC 0.0 128.0 Do every YINCth row ZINC 0.0 128.0 Do every ZINCth plane PCUT 0.0 Flux cutoff in total polarization flux ICUT 0.0 Flux cutoff in unpolarized I DOOUTPUT -1.0 3.0 Write residual image: 1,3 write parameter img: 2,3 NITER 0.0 4000.0 Max # fit iterations 0 -> 100. NGAUSS 1.0 4.0 Max number components DOSPIX -1.0 1.0 > 0 => fit spectral index else keep at zero. DOWEIGHT -1.0 2.0 > 0 weighted solution > 1.5 RMFIT finds weights INFILE Text file with weights BADDISK Disk to avoid for scratch ---------------------------------------------------------------- RMFIT Task: Fits polarization spectra for one or more components of total polarization, polarization angle at 0 wavelength, and rotation measure. RMFIT will be used on Q and U image cubes with frequency or frequency ID as the first axis and latitude and longitude as the second and third axes. It requires the real and imaginary or amplitude and phase images output by FARS, with rotation measure as the first axis and latitude and longitude as axes 2 and 3. It also requires an image of I polarization which can be a corresponding transposed cube or a single latitude/longitude plane. All images must have the same coordinates in all corresponding axes or the results will be rubish. The task fits up to 4 polarization components to each spatial pixel. Optionally, it writes up to 10*NGAUSS n-1 dimensional images containing the fit parameters and their uncertainties and/or the residual image. An interactive mode (using the TV graphics planes) is required. With NGAUS = 1, it is often possible to turn off the interactivity for large parts of the fitting if it is going well. The interactivity will be turned back on automatically if the answers are obviously wrong. Note that, since the TV is not allowed in BATCH jobs, RMFIT may not be run in BATCH. Color differentiates what is plotted. Graphics channel 1 (usually yellow) is used for axis labels and the data. The FARS output is plotted in amplitude/phase initially. The initial guess is shown as a large X in graphics channel 2 (usually green). This screen may be used to reset the initial guess. Then the Q and U data are plotted. Since the Q and U data may be on a FQID axis and so may not actually be in frequency order, the Q and U data are plotted as small X's. The initial guess is shown as a smooth curve in graphics channel 2 (usually green). The fit is then shown as somewhat large X's in graphics channel 4 (usually cyan). RMFIT started with the transformed XGAUS to be the task it should be. When run in interactive mode, no output files are created except for a large table file used to hold the results of the fitting (an RM file). The table is initialized with pixel positions and fluxes before any fitting is done. Then the fitting is done sequentially through the input cube. After that, in interactive mode, a menu-driven routine is brought up to allow you to review the fits including re-doing poor ones, rearranging the numbering of the components, and other editing of the results. Finally output files are written if requested. RMFIT may be re-started using the RM file either resuming the sequential fitting (if it was not completed) or dropping into the edit stage directly. You may create the RM file with expansive ranges for BLC(2) - TRC(2), BLC(3) - TRC(3) and then quit as soon as the first spectrum to be fit is displayed. You may then re-start using rather smaller ranges in Y and Z covering particular regions of similar polarization. The output images will cover the entire range initially specified. When doing this, it is a good idea not to change BLC and TRC (1,4,5,6,7). Details of the interactive operation and options are described in the EXPLAIN file. The model fit is: Q(i) = SUM_j [ P_j (nu(i) / nu_1)^alpha_j cos ( 2 Theta_j + 2 RM_j lambda(i)^2 ] U(i) = SUM_j [ P_j (nu(i) / nu_1)^alpha_j sin ( 2 Theta_j + 2 RM_j lambda(i)^2 ] where P_j is the polarizaed flux at 1 GHz, nu(i) is frequency, nu_1 is 1 GHz, lambda(i) is wavelength, alpha_j is spectral index. Output images are written of the P_j, Theta_j, RM_j, alpha_j, Q0_j, and U0_j and their uncertainties, where Q0_j = P_j cos ( 2 Theta_j ) U0_j = P_j sin ( 2 Theta_j ) are Q and U at 1 GHz without rotation measure. Adverbs: INNAME.....Input Q image name (name). Standard defaults. INCLASS....Input Q image name (class). Standard defaults. INSEQ......Input Q image name (seq. #). 0 => highest. INDISK.....Disk drive # of input Q image. 0 => any. INVERS.....Input RM table version. 0 => a new one. Set this adverb if you wish to re-start a fitting session. IN2NAME....Input U image name (name). Standard defaults. IN2CLASS...Input U image name (class). Standard defaults. IN2SEQ.....Input U image name (seq. #). 0 => highest. IN2DISK....Disk drive # of input U image. 0 => any. The I image may either be a transposed cube matching the Q and U cubes, or it may be a single plane with the latitude/longitude axes 1 and 2 matching the latitude/longitude axes 2/3 of the Q, U, real, and imaginary cubes. IN3NAME and IN3CLASS both blank means to ignore the test on I entirely. IN3NAME....Input I image name (name). Standard defaults. IN3CLASS...Input I image name (class). Standard defaults. IN3SEQ.....Input I image name (seq. #). 0 => highest. IN3DISK....Disk drive # of input I image. 0 => any. IN4NAME....Input FARS real/amplitude image name (name). Standard defaults. IN4CLASS...Input FARS real/amplitude image name (class). Standard defaults. IN4SEQ.....Input FARS real/amplitude image name (seq. #). 0 => highest. IN4DISK....Disk drive # of input FARS real/amplitude image. 0 => any. IN5NAME....Input FARS imaginary/phase image name (name). Standard defaults. IN5CLASS...Input FARS imaginary/phase image name (class). Standard defaults. IN5SEQ.....Input FARS imaginary/phase image name (seq. #). 0 => highest. IN5DISK....Disk drive # of input FARS imaginary/phase image. 0 => any. OUTNAME....Output image name (name). Standard defaults. OUTCLASS is set for each image by the task: QRESID, URESID, PPOLn, THETAn, ROTMEn, Q0_n, U_n, DPPOLn, DTHETn, DROTMn, DQ0_n, DU0_n for n = 1 to NGAUS. OUTSEQ.....Output image name (seq. #). 0 => highest unique. OUTDISK....Disk drive # of output image. 0 => highest number with sufficient space. BLC........Bottom right corner in input image of desired subimage. Default is entire image. If INVERS is specified, this BLC - TRC region in axes 2 and 3 must match or be interior to the region specified when the RM table was created. Processing is restricted to the current BLC to TRC. TRC........Top right corner in input image of desired subimage. Default is entire image. YINC.......Do initially every YINC'th row (beginning at BLC(2)). All rows will be done in the second stage. ZINC.......Do initially every ZINC'th plane (beginning at BLC(3)). All planes will be done in the second stage. PCUT.......A polarized flux cutoff in the same units as the input image (i.e. Jy/beam). If a row does not have an average polarized flux above this level, no component is fit to the row. It is also used to limit the points which determine the initial guess. If you do a re-start with a new lower value of PCUT, RMFIT will go into a special sequential mode to fit all rows with values between the previous and new PCUT levels. ICUT.......A total intensity flux cutoff in the same units as the input image (i.e. Jy/beam). If a celestial-coordinate pixel does not have an average total intensity flux above this level, no component is fit to the pixel. If you do a re-start with a new lower value of ICUT, RMFIT will go into a special sequential mode to fit all rows with values between the previous and new ICUT levels. DOOUTPUT...= 1 or 3 requests that the residual (data-model) images be written out as cataloged images = 2 or 3 requests that the model images be written out as cataloged images Note that DOOUTPUT can be changed interactively during the imaging portion of the task and the value of the adverb at exit is all that matters. NITER......Maximum function evaluations during the fit of each row. (< 10 -> 100 for NGAUSS > 1, < 100 -> 150 for 1 Component) NGAUSS.....Number of Components to fit (between 1 and 4). DOSPIX.....If > 0, solve for spectral index in the corresponding component, otherwise take the spectral index to be 0. DOWEIGHT...> 0 => Do a weighted solution. If INFILE is not blank, and DOWEIGHT < 1.5, read the weights from INFILE. Otherwise compute the weights from the rms of each image plane found through robust methods. <= 0 => Use all weights = 1. INFILE.....Text file containing the weights, like the file in FARS but allowing weights of any positive value (not just integers from 1 to 99). Note that RMFIT expects a weight for every pixel on the first axis of INNAME but uses only those numbered BLC(1) through TRC(1). If there are more weights than the total number of pixels on the first axis, then the second half of the file is assumed to be for U weights. Any Q weight <= 0 is replaced with 1.0. Any U weight <= 0 is replaced with the corresponding Q weight. Weights may have any positive value, not just integers and more than one value may appear on each line of INFILE. If a line contains a semi-colon, then the content of the line from the semi-colon to the end is treated as a comment. Note that the weights actually used are scaled so that their sum is 2N where N is the number of spectral channels used in the fit. They are recorded in the history file of the output images. BADDISK....Disk drives to avoid for scratch files. ---------------------------------------------------------------- IMAGR: One-dimensional polarization fitting of image cubes Documenter: E. W. Greisen NRAO Related Programs: IMFIT, JMFIT, SAD, SLFIT, ZEMAN, XGAUS RMFIT is a new task in 31DEC13 based on the overhauled Gaussian task XGAUS. At present, RMFIT assumes that there are at most 3 axes with more than one pixel in the input images. RMFIT begins by making an RM table containing one row for every pixel in the output plane from BLC(2,3) to TRC(2,3). Each row of the table is initialized with the average value of total polarization and total intensity for each output (spatial) pixel. Then RMFIT begins fitting those rows having average values above PCUT and ICUT, doing every YINC pixel in the Y axis and ZINC pixel in the Z axis. Initial guesses for each fit are found from the FARS output data. When NGAUSS = 1 this can happen automatically from the row data or it can come from a previous fit for NGAUSS > 1. After each input image row (output pixel) is fit, the results are added to the RM table immediately. This allows the interactive user to quit at any time and then restart the process later (set INVERS to point at the pre-existing RM table which you intend to modify). After the YINC by ZINC pass through the data, a second pass through the data is done on every output pixel. Those pixels which have already been fit are skipped, but the results from them contribute to the next initial guess. Finally, when all pixels have been fit, RMFIT enters an interactive routine designed to improve the results before any output images are written. However, before describing that function, we need to provide details of the interactive fitting of each row of the input image. This process is followed whenever a row is fit in the first YINC/ZINC function, in the second every pixel function, and in the third result editing function. For each spectrum, the first step in the fitting process is to determine an initial guess non-interactively. If RMFIT is currently in interactive mode, the next steps are: 1. The input FARS output data are plotted on the TV as amplitude and phase. The initial guess is added to the plot of the data as a large X in both the amplitude and phase portions of the plot. 2. If either NGAUSS > 1 or a fitting has been forced, RMFIT asks in the AIPS window for instructions. If your answer begins with E or e, RMFIT will ask you to enter a new initial guess (step 3). If the answer begins with B or b, RMFIT will mark this pixel as bad and, if the answer begins with Q or q, RMFIT will simply exit, allowing you to restart it at a later time. If the answer is anything else, e.g. a simple carriage return ("Enter"), RMFIT goes on to step 4 below. 3. To enter a new initial guess, you will be prompted to position the TV cursor at the rotation measure of each Component (note that only the X position matters). To mark each point, hit any button (A, B, C, or D). RMFIT then returns to step 1 above to plot the new guess and ask again. 4. Once an acceptable initial guess has been found, RMFIT proceeds to call a non-linear least squares fitting routine to determine the polarization parameters that appear to fit the data best. The answers are then checked to see if they are "reasonable" - negative and very large components and components with rotation measures outside the input data are unreasonable. 5. If the result is unreasonable and the current mode is not interactive, the interactive mode is turned back on (with a message and a display of the answers in the message terminal window) and RMFIT returns to step 1. 6. If the current mode is interactive, the final model is added to the plot. In the input terminal, you are then told if the answers are "unreasonable" and are shown the answers, reasonable or not. You are then offered a variety of choices. If your response begins with: a) B or b - the solution is marked as bad and RMFIT goes on to the next input row. b) Q or q - RMFIT exits cleanly, leaving the RM table to be worked on again at a later time. c) R or r - RMFIT returns to step 1 to try again with the current number of Gaussians (<= NGAUSS) c) E or e - RMFIT returns to step 1 to try again with the current number of Gaussians (<= NGAUSS) e) 1, 2, 3, or 4 - RMFIT returns to step 1 to try again with the specified number of Components (note that only numbers <= NGAUSS are respected) f) H or h - RMFIT prompts you to enter the polarization parameters for each component. You enter the total polarization in input image units, the polarization position angle at zero wavelength in degrees (-90 to +90), and the rotation measure in radians per meter squared. All 3 numbers must appear in one line. It will prompt for all NGAUSS components, but will change the prompt if the current number of Components is < NGAUSS. If you enter something besides 0 in that case, the the current number of Components is increased appropriately. RMFIT then returns to the start of this step 6 to allow you to see if you made a good guess. g) T or t - RMFIT turns off the interactive mode and does solutions in a batch-like fashion until an unreasonable solution is found or it finishes the current function. h) P or p - RMFIT prompts you for the Q and U display ranges to re-display the data. Enter Qmin, Qmax, Umin, Umax in that order in one line. If you enter only two values, then Umin and Umax are taken to equal Qmin and Qmax. If you enter no values, then the plot reverts to self-scaling. The Q plot self-scales if Qmin >= Qmax; the U plot self-scales if Umin >= Umax. Note that on the re-plot, the smooth "initial guess" line is now actually the result of the fit. i) Anything else - RMFIT saves the answer (with uncertainties) and goes on to the next input spectrum. Note, on retries (R, E, 1, 2, 3, and 4), step 2 will skip the questions and go on to step 3 directly (as if you had answered E). After all pixels that are strong enough have been given a solution, RMFIT enters a menu-driven function. The menu has in the left column: --------------------- | EXIT | Exit RMFIT, writing output images if DOOUTPUT is now > 0. | SET MIN S/N | Set minimum total polarization S/N(s) for okay solutions | SET MAX RESID | Set maximum residual for okay solutions | SET RM RANGE | Set rotation measure range(s) for okay solutions | SET MAX THETA ERR | Set maximum error(s) in polarization angle for okay solutions | REDO ALL | Re-do all solutions which are not okay | FLAG ALL | Mark bad all solutions which are not okay | OFF ZOOM | Turn of TV zoom | OFF TRANSFER | Turn off black & white and color TV enhancements | SET DOOUTPUT | Increment DOOUTPUT in loop 0-3 - with 1 and 3 causing residual images and 2 and 3 causing parameter images to be written on EXIT | ADD TO LIST | Type in output pixel coordinates to add to list | SHOW LIST | Display coordinates in list | REDO LIST | Re-do solutions for all pixels in list | FLAG LIST | Flag solutions for all pixels in list | SWAP LIST 1-2 | Swap solutions for components 1 and 2 for all pixels in list | SWAP LIST 1-3 | Swap solutions for components 1 and 3 for all pixels in list | SWAP LIST 2-3 | Swap solutions for components 2 and 3 for all pixels in list | SWAP LIST 1-4 | Swap solutions for components 1 and 4 for all pixels in list | SWAP LIST 2-4 | Swap solutions for components 2 and 4 for all pixels in list | SWAP LIST 3-4 | Swap solutions for components 3 and 4 for all pixels in list --------------------- Only appropriate SWAP LIST options are actually displayed. Three editing concepts are present in this menu. The first is to establish what parameter values (S/N in peak, peak residual, rotation measures, and error in polarization angle) for each component constitute an "okay" solution. The selected limits are displayed above the menu. Then you can choose to try to REDO ALL solutions which are not okay or simply FLAG ALL, marking them as bad. The second concept is to create a list of up to 1000 output pixels which are thought to have some problem. You can enter these pixels by hand (or remove pixels from the list by hand) with ADD TO LIST. A faster way to add to the list is to show the image of some parameter and then select interesting pixels while doing CURVALUE (see below). The list of pixels can be re-done (REDO LIST) or marked as bad (FLAG LIST). The list is cleared by these operations so you can make a new list as needed. The third editing method available here is to take the pixel list and swap the solutions at those pixels between component N with component M. The right hand column of options include: ------------------ | SHOW IMAGE P0_1 | Enter image interaction with total polarization value of component 1 | SHOW IMAGE TH_1 | Enter image interaction with polarization angle of component 1 | SHOW IMAGE RM_1 | Enter image interaction with rotation measure of component 1 | SHOW IMAGE Q0_1 | Enter image interaction with Q polarization at wavelength 0 of component 1 | SHOW IMAGE U0_1 | Enter image interaction with U polarization at wavelength 0 of component 1 | SHOW IMAGE EP0_1 | Enter image interaction with uncertainty in total polarization of component 1 | SHOW IMAGE ETH_1 | Enter image interaction with uncertainty in polarization angle of component 1 | SHOW IMAGE ERM_1 | Enter image interaction with uncertainty in rotation measure of component 1 | SHOW IMAGE EQ0_1 | Enter image interaction with uncertainty in Q polarization at wavelength 0 of component 1 | SHOW IMAGE EU0_1 | Enter image interaction with uncertainty in U polarization at wavelength 0 of component 1 ------------------ For NGAUSS > 1, appropriate additional choices are offered. When you select one of the above options, the above menus are turned off, the image plane maintained in memory is displayed on the TV screen, and a new menu is displayed. If offers: --------------- | RETURN | Return to the above menus, image stays displayed | LOAD AS SQ | Re-load image with square root transfer function | LOAD AS LG | Re-load image with log transfer function | LOAD AS L2 | Re-load image with extreme log transfer function | LOAD AS LN | Re-load image with linear transfer function | OFF TRANSF | Turn off enhancement done with TVTRANSF | OFF COLOR | Turn off color enhancements | TVTRANSF | Black & white image enhancement | TVPSEUDO | Color enhancement of numerous sorts | TVPHLAME | Color enhancement of flame type, multiple colors | TVZOOM | Interactive zooming and centering of image | CURVALUE | Display value under cursor, mark pixels for list | SWAP 1-2 | Swap solutions for components 1 and 2 interactively | SWAP 1-3 | Swap solutions for components 1 and 3 interactively | SWAP 2-3 | Swap solutions for components 2 and 3 interactively | SWAP 1-4 | Swap solutions for components 1 and 4 interactively | SWAP 2-4 | Swap solutions for components 2 and 4 interactively | SWAP 3-4 | Swap solutions for components 3 and 4 interactively | NEXT WINDOW | Move to next window into large images --------------- For NGAUSS < 4, appropriate SWAP options are suppressed. Only one of the LOAD AS xx options is offered - SQ when the current function type is LN, LG when the current function type is SQ, L2 when the current function type is LG, and LN when the current type is L2. Selecting this option, reloads the image with the newly selection function type and changes the menu option accordingly. The options OFF TRANSF through TVZOOM are essentially the same as the corresponding verbs in AIPS. CURVALUE is like the verb of that name, but, when you press buttons A or B, the pixel under the cursor is added to the editing list. The SWAP n-m options invoke the interactive polygon setting operation used by the AIPS verb TVSTAT. This starts in the "set polygon n" mode where button A selects a vertex in the polygon, B sets the last vertex in the polygon and gets ready to set polygon n+1, C sets the last vertex and enters a vertex editing mode, and D sets the last vertex and exits the polygon setting. In the vertex editing mode, move the cursor to a vertex to be moved and press button A or B to reset its position. After moving it to the desired point, press button A or B to fix that point and restart the vertex editing mode. Press button C to fix that point and then start creating polygon n+1. Press button D to fix that point and exit the polygon setting. Once button D has been set, the selected areas of pixels have their solutions swapped. The images are updated and reloaded automatically. Very large images may not be able to fit on your TV. When the image exceeds the size of the TV, the top line will show the component number followed by a subimage number in parentheses and the NEXT WINDOW option will appear. The first sub-images displayed is called number 0 and shows the full image every n'th pixel in X and Y. Sub-image 1 begins at the lower left, moves right, then back to the left and up, and so forth until the top right is reached. Every pixel is displayed in these sub-images. Use the NEXT WINDOW option to step through the sub-images in a circular fashion.