% 2007.09.18 NOTE: putting under svn control. % % $Date: 2008-05-02 11:48:39 -0400 (Fri, 02 May 2008) $ % $Author: dph $ % $Rev: 188 $ % % % This file is part of findzo.sl (finding zeroth order). % Copyright (C) 2006-2008 Massachusetts Institute of Technology % % This software was developed by the MIT Center for Space Research under % contract SV3-73016 from the Smithsonian Institution. % % 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 Mass Ave, Cambridge, MA 02139, USA % %----------------------------------------------------------------------- % % File: findzo.sl % % Version: % % Authors: Bish Ishibashi % Date Created: 21 APR 2005 % % % HISTORY: % This script has been modified numerous times. % % September 21, 2005: Repair for LETG+HRC-S. % Hidden / unnecessary parameters removed. % Then some unnecessary parameters added as a placeholder % for future support of HRC-I and ACIS-I. % % September 23, 2005: Add a support for cfitsio subspace filtering. % September 30, 2005: Remove the namespace setting since sherpa % imports isis by default and is no % longer necessary. % % March 21, 2006: Changes in parameter names and correction % on cfitsio subspace filtering. % % August 24, 2006: Add a check on subarray condition. Scale rad_strc if necessary. % August 31, 2006: Replace "static" with "private". Add "provide()" at the end. % October 3, 2006: Add a better bowl_strc condition in the subarray situation. % Allow NROWS to be set to 1024 when no subarray is used. % October 17, 2006: Correct an error on deriving the initial guess [gx, gy] for % ftype==1. I sloppily omitted cos(decnom) in evaluating an % angular separation in RA. % October 18, 2006: Per JED's advice, implement the conditions to control verbosity % (i.e., to enable users to make it less chatty). The executable % script testzo is also adjusted for this. A library called % set_findzo_verbose is added to utilize this feature. % October 23, 2006: Change in the grating angle parameters (consistent with the % geometry update reported by me in July 2006. % December 27, 2006: Adjust verbose mode settings and messages. % % March 22, 2007: [dph] fix plot of sub-sampled % data. Picking first N events is not % reliable because sometimes there are % bad-aspect points. Better to simply % sample every nth via [[::nth]]. % % August 21 2007 [dph] changed scaling of zoomed plot to % show more zo detail. (zoomed in) % %========================= % % Purpose: This program finds the centroid position (SKY_X,SKY_Y) of zeroth order % image by fitting lines to the MEG and readout streaks, and then determining % the intersection of the lines. % % Syntax: (x,y)=findzo(evtfile_name, grating_type [,guessing_mode, gx, gy]) % % Input: % evtfile_name == Event L1.5 or 2 file. % grating_type == "h" (HEG), "m", or "l" (LEG). % guessing_mode == 1 = astrometry based % 2 = median filter based % 3 = dumb guess ([X,Y]=[4096.5,4096.5]) % 4 = user specified guesses (sky x and y) % gx, gy == initial guess on the source position (in sky) % % Example % by default, flag == 1: % (x,y)=findzo("acisf0001N0001_evt1.fits","m") % % Or try something fancier. When using CFITSIO filtering, % do not specify the extention block name since it is set to % [EVENTS] by default (always true for Chandra Event files). % % (x,y)=findzo("acisf0001N0001_evt1.fits[energy=2000:8000]","h",4,4090.2,4061.72) % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Version: commented out for now. % private variable _version = [0, 9, 3 ] ; % major, minor, patch % % % Verbosity default value % private variable FINDZO_Verbose = 5; private define vmsg () { variable args = __pop_args (_NARGS-1); variable vb = (); variable fp = stdout; if (FINDZO_Verbose >= vb) () = fprintf(fp, __push_args (args)); } define set_findzo_verbose (vb) { FINDZO_Verbose = vb; } % % Control Display Setting % private variable FINDZO_Plot_Dev = "/xwin"; private define pdev () { variable pb = (); variable pw; if (FINDZO_Plot_Dev == "/xwin") { if (pb == 3) return; pw = open_plot("/xwin");resize(20,0.9);erase; window(pw); } else if (FINDZO_Plot_Dev == "/cps") { if (pb == 1) pw = open_plot("findzo.ps/cps"); if (pb == 3) close_plot; else return; } else pw = open_plot("/xwin"); } define set_findzo_plot_dev (pb) { FINDZO_Plot_Dev = pb; } % % Subroutines: % %------------------------------------------------------------------------ % bfun: % % S-Lang function: y = f(x,pars) = grang*(x + pars[0]) + b. % % The function will be fit to derive the best-fit parameter value "b". % The slope of the function f(x,pars) is fixed in this function. % private define b_fun (x, pars) { return pars[0] + x ; } %------------------------------------------------------------------------ % garm_fit: % % Fit a line with a specificied slope value. % % This function depends on a S-lang function b_fun. % private define garm_fit (x,y,angle) { variable pars, pars_min,pars_max; variable out,stat; pars=4000; % This is rather ACIS-S specific. pars_min=-1e5; pars_max=1e5; (out,stat) = array_fit (x*angle,y, NULL, pars, pars_min, pars_max, &b_fun); return [out,angle], NULL; } %------------------------------------------------------------------------ % clip_array: % % A customized S-lang function to perform clipping of array % elements found NSIG x SDEV times above or below the AVERAGE % value of the array Y. This function will return the indices % of the array Y that meets the given criterion NSIG. % private define clip_array( y, nsig ) { variable l, m, n, lprev; l = where( y == y ); lprev = 0 ; while ( length(l) - lprev ) { lprev = length(l); m = moment( y[l] ); l = where ( abs( y - m.ave ) < (nsig * m.sdev) ) ; } return l ; } %------------------------------------------------------------------------ % rotxy: % % A simple S-Lang function for 2-dimensional array rotation. % % The function rotates the input array (x,y) by the value ANGLE around % the pivot point (x0, y0). For a given positive ANGLE value, the array % is rotated in a CLOCK-WISE direction. % private define rotxy( x, y, x0, y0, angle ) { if (_NARGS !=5 ) { message("USAGE: (newx, newy) = rotxy( x, y, x0, y0, angle );"); message(" Rotate axes clockwise by angle (in radians)"); message(" with center, x0, y0"); return -1; } variable ca, sa, xca, xsa, yca, ysa ; ca = cos(angle); sa = sin(angle); xca = (x-x0)*ca; xsa = (x-x0)*sa; yca = (y-y0)*ca; ysa = (y-y0)*sa; return xca+ysa + x0, -xsa+yca + y0 ; } %------------------------------------------------------------------------ % median_peak: % % A special S-Lang function to trace CHANDRA grating arms detected % by ACIS and HRC detectors. % % This function performs the following steps: % % (1) read in two data arrays (X,Y); the elements in Y is binned up % along X with a specified bin width BIN. % % (2) Use the clip_array function to clip data points that are considered % as insignificant (its significance controlled via DEL). % % (3) then derive the median value of data points in each data bin along X. % % (4) return the median value in each bin. % % This function is actually NOT a mission-specific one. It could work for % any 2-D data array. % private define median_peak(x,y,bin,del) { variable xp, yp, nx, j, k, l,m; nx=int((max(x)-min(x))/bin); xp=[0:nx-1:1]*bin +min(x)+bin*0.5; yp=[0:nx-1:1]*0.0-32768.; % Silly little hack. % Sometimes the median values % are not found in the routine; % when that happens, the value YP % stays as signed integer -32768. % The array elements with this value % are rejected and not returned. for (k=0; k < length(xp); k++) { j=where((x >= xp[k]-bin*0.5) and (x < xp[k]+bin*0.5)); (m,,)=array_info(j); if (m[0] >= 2) { l=clip_array(y[j],del); (m,,)=array_info(l); if (m[0] >= 2) { yp[k]=median(y[j[l]]); } } } m=where(yp != -32768.); return xp[m],yp[m]; } %------------------------------------------------------------------------ % gparam_io: % % A S-Lang function to set up both required and hidden parameters % to make the function "findzo". % % This function is essentially Bish's version of "param_io" that controls % the parameter setting to run the function FINDZO. This one is CHANDRA % mission specific. Almost all the hard-coded parameter values are found % in this function (hence, this is where the users would want to mess with % if he/she wants to customize parameters...I don't see why, though). % private define gparam_io (fitsfile, fitssub, gtype, ftype, x0 ,y0) { % % Initialize Variables variable fp, nx, ny; variable rotkey="ROLL_NOM"; variable grelem="GRATING"; variable detnam="INSTRUME"; variable ranom, denom, ratg, detg, rarf, derf; variable deg_rad=PI/180.; variable gx, gy, gflag; variable x, y, en, cy; variable rotang, grang; variable subarr; % % note on angles: positive angle == clockwise rotation angle in space craft coordinate. % variable strc = 0.0; % Initialization variable strc1 =-0.011; % Streak offset position for MEG variable strc2 = 0.05; % Streak offset position for HEG variable strc3 = 0.00; % Streak offset position for LEG (zero) variable acisscor = 0.03; % empirical acis-s rotation correction variable hrcscor = 0.2744; % empirical hrc-s rotation correction variable acisicor = 0.0; % empirical acis-i rotation correction (dummy placeholder) variable hrcicor = 0.0; % empirical hrc-i rotation correction (dummy placeholder) variable alpha1 = 4.755; % dispersion angle for MEG wrt CHIP_Y (from CALDB) variable alpha2 =-5.205; % dispersion angle for HEG wrt CHIP_Y (from CALDB) variable alpha3 = 0.07; % dispersion angle for LEG wrt CHIP_Y (from CALDB) variable corr_a1 =-0.006; % supplemental correction factor for alpha1 variable corr_a2 = 0.001; % supplemental correction factor for alpha2 variable corr_a3 = 0.00; % supplemental correction factor for alpha3 (dummy placeholder) variable tmp_a1 = 0.00; % ad-hoc correction factor allowed for alpha1 (dummy placeholder) variable tmp_a2 = 0.00; % ad-hoc correction factor allowed for alpha2 (dummy placeholder) variable tmp_a3 =-0.089; % ad-hoc correction factor allowed for alpha3 % ACIS-S or LETG arm is wicked and need % an additional bogus correction factor. variable tmp_a4 = 0.00; % ad-hoc correction factor allowed for alpha4 % (LETG+HRC-S: dummy placeholder) % % Plate scale (in SKY) % variable acis_plate = (0.492/3600.)*deg_rad; % in radian variable hrc_plate = (0.1318/3600.)*deg_rad; % % % First, check its grating configulation. And then check the nominal roll angle. % If these information are not found (e.g., it is not a grating observation), % then the NULL values will be returned. % if (fits_key_exists(fitsfile,grelem)) { grelem=fits_read_key(fitsfile+fitssub,grelem); } else return NULL,NULL; % no rotkey then return null if (fits_key_exists(fitsfile,rotkey)) { fp=fits_open_file(fitsfile+fitssub,"r"); ()=_fits_get_colnum(fp,"x",&nx); ()=_fits_get_colnum(fp,"y",&ny); ()=_fits_close_file(fp); (rotang,ranom,denom,ratg,detg,rarf,derf,detnam)= fits_read_key(fitsfile+fitssub,rotkey,"TCRVL"+string(nx),"TCRVL"+string(ny), "RA_TARG","DEC_TARG","TCRPX"+string(nx),"TCRPX"+string(ny),detnam); } else return NULL,NULL; % no grating then return null % % Case 1: HETG % % % The lines below seem too complicated? Here is a synopsis. % % First for ACIS-S. The ACIS-S array happens to be tilted by -0.03deg, i.e., % transfer streak of zeroth order is tilted by +0.03 deg in Sky coordinate. % And yet the alpha clock angles of HEG and MEG arms are -5.205 and 4.755 degs % with respect to CHIP_Y axis across the detector array (reminder, positive % angle goes in clockwise direction; i.e., -5.205 deg == 5.205 deg in ccw % from CHIP_Y axis). So if the whole array is rotated by -0.03 such that % the streak is now vertical in the mapped frame, then the alpha clocking % angles must be rotated by -0.03 degrees. And that's what is happening. % % The same goes with LETG and ACIS-S, though there is an ad-hoc parameter % to correct inconsistency (that rotation problem that plague LETG). % % Also a similar issue emerges in LETG+HRC-S. But things get much muddier % in that configuration. Really really nasty and I have no idea how accurate % the existing geometrical calibration on this particular configuration is. % % % rotang, strc, and grang are determined here. % if (grelem == "HETG") { alpha1 = alpha1 + corr_a1; % Ad-hoc correction (empirical value) alpha2 = alpha2 + corr_a2; % Ad-hoc correction (empirical value) % maybe set to zero (Sep 20, 2005) rotang = rotang + acisscor; % compensating a small ACIS-S detector % rotation so that the streak is vertical. if (gtype == "h") { strc = strc2; grang = alpha2 - acisscor + tmp_a2; % this is the corrected alpha clocking angle % in the final rotated frame (streak vertical). } if (gtype == "m") { strc = strc1; grang = alpha1 - acisscor + tmp_a1; % Ditto } gflag=12; } % the end of if (grelem == "HETG") % % Case 2: LETG + HRC-S or ACIS-S. This section is a mess. % if (grelem == "LETG") { gflag=3; alpha3 = alpha3 + corr_a3; % Ad-hoc alpha correction (empirical value) if (detnam == "HRC") { strc = strc3; rotang = rotang + hrcscor; grang = alpha3 + tmp_a4; } if (detnam == "ACIS") { strc = strc3; rotang = rotang + acisscor; % compensating a small ACIS-S detector rotation wrt Roll. grang = alpha3 + tmp_a3; % ad-hoc correction for detector mismatch (acis vs. hrcs w/ letg) } if (detnam == "ACIS-I") % This is a dummy placeholder { message("This Mode is not Supported.\n"); return -1; } if (detnam == "HRC-I") % This is a dummy placeholder { message("This Mode is not Supported.\n"); return -1; } } if (grelem == "NONE") { message("This is not a grating dataset.\n"); % This should not happen. return -1; } % % Do a quick sanity check to see if the right grating type is selected. % % gflag == 12 if the grating modes are 1 (HEG) or 2 (MEG), % or == 3 if the grating mode is 3 (LEG). % if ((gflag == 12 and gtype == "l") or (gflag == 3 and (gtype == "h" or gtype == "m"))) { message("Grating Type Mismatch. Please check your brain.\n"); return -1; } % % Initial Guess (x0,y0) from _TARG coordinate % % Now read in the raw data points and re-rotate back such that % the transfer streak is always vertical across the field in the resulted map. % if (detnam == "ACIS") { (x,y,en,cy)=fits_read_col(fitsfile + fitssub,"x","y","energy","CHIPY"); } if (detnam == "HRC") { (x,y,cy)=fits_read_col(fitsfile + fitssub,"x","y","CHIPY"); en = x * 0.0 + 2000.00; % Fudge the energy column since HRC doesn't have one. } if (ftype == 1) { vmsg(2, "Assuming that the RA_TARG and DEC_TARG values are correct\n"); vmsg(2, " for this particular target. \n"); if (detnam == "ACIS") { gx = -(atan((ratg - ranom)*deg_rad)/(acis_plate))*cos(denom*deg_rad)+rarf; gy = (atan((detg - denom)*deg_rad)/(acis_plate))+derf; } if (detnam == "HRC") { gx = -(atan((ratg - ranom)*deg_rad)/(hrc_plate ))*cos(denom*deg_rad)+rarf; gy = (atan((detg - denom)*deg_rad)/(hrc_plate ))+derf; } } if (ftype == 2) { vmsg(2, "Making the guess based on Median Sky X and Y values\n"); gx = median(x); gy= median(y); } if (ftype == 3) { if (detnam == "ACIS") { vmsg(2, "Assume dumb coordinate values [4096.5, 4096.5]\n"); gx = 4096.5; gy = 4096.5; % default dummy x,y values } if (detnam == "HRC") { vmsg(1, "Assume dumb coordinate values [32768.5,32768.5]\n"); gx = 32768.5; gy = 32768.5; % default dummy x,y values } } if (ftype == 4) { vmsg(2, "User specified values will be used\n"); gx = x0; % user specified values. gy = y0; % } % % Define the REGION BOX for filtering events. % variable boxw_strc = 80.0; % Box width for streak region variable boxl_strc = 2000.0; % Box length for streak region variable boxw_garm = 50.0; % Box width for grating arm region variable boxl_garm; % Box length for grating arm region variable bin_garm; % Bin size along X variable rad_garm = 50.; % Circle Radius for grating arm region variable rad_strc = 50.0; % Circle Radius for streak region % Check the array usage subarr = max(cy) - min(cy); if (subarr <= 768) % if 3/4 or less CCD area is used (subarray) { % then use a smaller radius of circle for streak. rad_strc *= 0.4; boxw_strc *= 0.5; } % % Set grating arm region % if (detnam == "ACIS") { boxl_garm = 5000.0; bin_garm = 20; } if (detnam == "HRC") { boxl_garm = 12000.0; bin_garm = 500; } % % Control the display box size, etc. % variable disp_n_elem = 50000; % Number of elements displayed. variable disp_bx_size = 1000; % Define the range of (X,Y) box size % in the plot display % % Misc. filtering (energy cut, significance setting for data clipping) % variable n_del = 1.3; % significance setting for array clipping variable n_eng_hi = 4000; % energy filtering (high threshold) variable n_eng_lo = 500; % energy filtering (low threshold) % % and that's it. Return the parameter values required for findzo. % return, gx, gy, x, y, en, grang, rotang, detnam, strc, boxw_strc, boxl_strc, rad_strc, boxw_garm, boxl_garm, rad_garm, bin_garm, disp_bx_size, disp_n_elem, n_del, n_eng_hi, n_eng_lo; } % % Main Routine: % %------------------------------------------------------------------------ % findzo: % % A S-Lang function to determine the position of zeroth-order % centroid (CHANDRA/HETGS and LETGS specific) by determining % the intersection of grating arms and data transfer streak % (or fine support structure diffraction pattern in LETG/HRC). % % % Usage: (x,y) = findzo(fits_event_file,gtype [,ftype, x0, y0]); % % where: % % gtype == h, m or l. % ftype == 1, 2, 3, or 4. % % 1: use RA_ and DEC_TARG to derive the guesstimated source position [gx, gy] % 2: use the median values as [gx, gy] % 3: use the dumb guess [x0,y0] as the initial guesses for [gx, gy] % 4: use the user specified coordinate [x0, y0] as [gx, gy] % define findzo( ) { vmsg(1, "Running FINDZO version "+ string(_version[0])+ "."+string(_version[1])+"."+string(_version[2])+"\n"); % % here is one hack for isis 1.2.6 or later. % set_fit_method("marquardt;delta=1.0e-12"); % % Initialize Variables % variable fitsfile, gtype, ftype; variable gx, gy; variable x0, y0; variable boxl_strc, boxw_strc, rad_strc; variable boxl_garm, boxw_garm, rad_garm, bin_garm; variable detnam="INSTRUME"; variable deg_rad=3.1416/180.; variable en, x, y, xp, yp, xpp, ypp; variable sxp, syp, sxp2, syp2, cnt, cnt2; variable cltmp, l, ll, mxp, myp, strx, slp, islp, zx, zy; variable nx, ny; variable disp_bx_size, disp_n_elem; variable n_del, n_eng_hi, n_eng_lo; variable p1, p2; variable strc; % % note on angles: positive angle == clockwise rotation angle. % variable grang, rotang; % % Read input from command line % if ((_NARGS < 2) or (_NARGS == 4) or (_NARGS > 5)) { vmsg(1, " \n"); vmsg(1, "USAGE in ISIS: (x, y) = findzo(fits_file,\"l\",[guessing_mode,gX,gY]);\n"); vmsg(1, " where \"l\" (LEG) can be replaced with \"h\" (HEG)\n" ); vmsg(1, " or \"m\" (MEG). guessing_mode, gx and gy are optional.\n" ); vmsg(1, " \n"); vmsg(1, "USAAGE in CIAO: tg_findzo infile l guessing_mode gX gY\n"); return -1; } if (_NARGS ==2 ) { (fitsfile,gtype)=(); ftype = 1; x0=4096.5; y0=4096.5; % These are dummy values. Not actually used. } if (_NARGS ==3 ) { (fitsfile,gtype,ftype)=(); % here, ftype can be any integer [1,2,3, or 4]. if (ftype == NULL) { vmsg(1,"ftype = NULL invalid\n"); return -1; } x0=4096.5; y0=4096.5; % Initial guesses. Will not be used. } if (_NARGS ==5 ) { (fitsfile,gtype,ftype,x0,y0)=(); % here ftype should be 4. [x0, y0] needs % to be specified. } % % Now allowing to accept the use of cfitsio extended name syntax. % % % The extention block name [EVENTS] should always be used for Chandra Event files. % variable tmp_name = strchop(fitsfile,'[',0); variable fitsname = tmp_name[0]; variable fitssub = "[EVENTS]"; if ( length(tmp_name) > 1 ) { tmp_name = tmp_name[[1:length(tmp_name)-1]]; % remove the first element since % they should be a fits filename % % Dummy proof % if ((strup(tmp_name[0]) == "EVENTS]") or (strup(tmp_name[0]) == "1]")) { vmsg(2," \n"); vmsg(2,"#############################################\n"); vmsg(2,"# Please do not specify the extention name. #\n"); vmsg(2,"#############################################\n"); vmsg(2," \n"); tmp_name = tmp_name[[1:length(tmp_name)-1]]; } % % vmsg(2," \n"); vmsg(2,"The extention block name is always set to be [EVENTS] by default. \n"); vmsg(2," \n"); if ( length(tmp_name) > 0 ) { fitssub = fitssub + "[" + strjoin(tmp_name,"["); } } % % Initialize the required parameters specific to this configuration. % % gx, gy == guesstimated coordinate position of the source (in SKY) % x, y == event data arrays % en == event energy % grang == grating arm angle % rotang == ajusted ROLL_NOM angle % detnam == DETECTOR name. % strc == Streak Correction Factor % box[w,l]_strc == Region Box Size for transfer streak (for fitting) % rad_strc == Circle of radius exclusion zone around the zeroth order image. % box[w,l]_garm == Region Box Size for grating arms (for fitting) % rad_garm == Circle of radius exclusion zone around the zeroth order image. % bin_garm == binning width for histogramming. % disp_bx_size == display box size to show the fidelity of fit % disp_n_elem == a total number of events shown in the pgplot screen. % n_del = a level of statistical significance used in filtering (e.g., clip_array). % n_eng_[hi,low] == Upper and lower threshold levels for energy filtering. % (gx, gy, x, y, en, grang, rotang, detnam, strc, boxw_strc, boxl_strc, rad_strc, boxw_garm, boxl_garm, rad_garm, bin_garm, disp_bx_size, disp_n_elem, n_del, n_eng_hi, n_eng_lo) = gparam_io(fitsname, fitssub, gtype, ftype, x0 ,y0); % % Rotate the SKY event array (X,Y) by (-rotang) degree around (gx,gy). % This will make the transfer streak go vertical in the new frame. % (xp,yp)=rotxy(x,y,gx,gy,-rotang*deg_rad); % % And draw a first region box to isolate transfer streak events % without zero order and with energy filtering. % l=where(xp > (gx-boxw_strc/2) and xp < (gx+boxw_strc/2) and yp > (gy-boxl_strc/2) and yp < (gy+boxl_strc/2) and x > 0.0 and y > 0.0 and hypot(xp-gx,yp-gy)>rad_strc and (en > n_eng_lo and en < n_eng_hi)); % % Fit the streak profile % (myp,mxp) = median_peak(yp[l], xp[l],20,n_del); (slp,) = garm_fit(myp,mxp,0.0); islp = 0; while (islp != slp[0]) { islp=slp[0]; cltmp=clip_array(mxp-myp*slp[1]+slp[0],2); % loosely filtered here. hardcode OK. (slp,)=garm_fit(myp[cltmp],mxp[cltmp],0.0); } strx = slp[0]+strc; % % Save some parameter values for later use % sxp=strx;syp=gy;sxp2=0.;syp2=0.;cnt=0; cnt2=0; % % Draw the events in the selected region box % pdev(1); connect_points(0); ll=where(x > 0.0 and y > 0.0); xrange(sxp-disp_bx_size,sxp+disp_bx_size); yrange(syp-disp_bx_size,syp+disp_bx_size); %%% dph 2007.03.22 % plot(xp[ll[[0:disp_n_elem]]],yp[ll[[0:disp_n_elem]]]); variable nth = max( [ length( ll ) / disp_n_elem, 1 ] ) ; plot( xp[ ll[ [::nth] ] ], yp[ ll[ [::nth] ] ] ) ; % % here, the while loop, is to check whether the position % solution points to the same place or not. Until it does, % it'll keep runing (til 20th cycle). % while ((sxp2 != sxp or syp2 != syp) and (cnt < 20)) { sxp2=sxp; syp2=syp; % % Rotate the region box again to rectify the selected grating arm. % (xpp,ypp)=rotxy(xp,yp,sxp,syp,-grang*deg_rad); % % Then index the events in the newly selected region box. % Here, the energy filtering is hardcoded (and the user % does not need to control this). % l=where(xpp>(sxp-boxl_garm/2) and xpp<(sxp+boxl_garm/2) and ypp>(syp-boxw_garm/2) and ypp<(syp+boxw_garm/2) and x > 0.0 and y > 0.0 and hypot(xpp-sxp,ypp-syp)>rad_garm and (en < 8000.)); % % Call median_peak function. % % The output are the point array (mxp,myp) (median point along the arm). % (mxp,myp)=median_peak(xp[l],yp[l],bin_garm,n_del); % % Fit the vector (mxp, myp). % (slp,)=garm_fit(mxp,myp,-grang*deg_rad); % % Do interative filtering; this may be redundant. % cltmp=clip_array(myp-mxp*slp[1]+slp[0],2); % again, hardcode OK (slp,)=garm_fit(mxp[cltmp],myp[cltmp],-grang*deg_rad); islp=slp[1]; variable flag=1; while (flag and cnt2 < 20) { cltmp=clip_array(myp-mxp*slp[1]+slp[0],2); % hardcode ok. (slp,)=garm_fit(mxp[cltmp],myp[cltmp],-grang*deg_rad); if (slp[1] == islp) { flag = 0; } islp=slp[1]; cnt2++; } if (cnt2 == 20) { vmsg(5,"the slope solution did not converge after "+string(cnt2)+" trials.\n"); } connect_points(1); sxp = strx; syp = slp[1]*sxp+slp[0]; cnt++; } % % Plot out the resulting line fits. % oplot([strx,strx],[min(yp),max(yp)]); oplot(mxp,mxp*slp[1]+slp[0]); % % Zoom in onto the zeroth order region. % % dph mod: 2007.08.16, zoom factor, make larger (smaller factor) variable zfact = 0.05 ; pdev(2); connect_points(0); % xrange(sxp-disp_bx_size*0.15,sxp+disp_bx_size*0.15); % yrange(syp-disp_bx_size*0.15,syp+disp_bx_size*0.15); xrange( sxp - disp_bx_size*zfact, sxp + disp_bx_size*zfact); yrange( syp - disp_bx_size*zfact, syp + disp_bx_size*zfact); ll=where(x > 0.0 and y > 0.0); % 2007.03.22 dph % plot(xp[ll[[0:disp_n_elem]]],yp[ll[[0:disp_n_elem]]]); nth = max( [ length( ll ) / disp_n_elem, 1 ] ) ; plot( xp[ ll[ [ ::nth ] ] ], yp[ ll[ [::nth] ] ] ) ; connect_points(1); oplot( [strx,strx], [ min(yp), max(yp)] ); oplot( mxp, mxp*slp[1]+slp[0] ); pdev(3); % % Return error message if the fit did not converge. % if (cnt == 20) { vmsg(5,"the position solution did not converge after "+string(cnt)+" trials.\n"); } % % Rotate back to the original [X,Y] coordinate frame. % (zx,zy)=rotxy([sxp],[syp],gx,gy,rotang*deg_rad); % % return the values. DONE. % vmsg(5, "Derived zeroth order position: X0 = "+string(zx[0])+" Y0 = "+string(zy[0])+" \n"); return zx[0], zy[0]; } provide("findzo"); provide("set_findzo_verbose"); provide("set_findzo_plot_dev"); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%