Re: fits in channel space

From: Kazunori Ishibashi <bish_at_email.domain.hidden>
Date: Tue, 04 Aug 2009 04:02:25 +0900
Hi there, 

I don't have much time to elaborate, but isn't it simpler 
-- if it's a Gaussian function you want to fit -- to write
S-lang functions using array_fit? 

Like, 

define g_fun (x,pars)
{
	return pars[3] * exp(-(x-pars[0])^2/pars[1]);
}


and

define g_fit (x,y,pars)
{
variable ....;
variable pars_max, pars_min;
pars_max=...;
pars_min=...;

(out,stat)=array_fit(x,y, NULL, pars, pars_min,pars_max, &g_fun);
return [out,stat],...;
}


Or something like that (I am pretty sure that the examples 
above make John Davis pull his hair...)? 


....


BTW, I just remembered that the very first prototype of tg_findzo
(aka findzero.sl) did what Maurice was trying. I attached the ISIS
script in this mail (is an attachment allowed? I forgot). Take a look 
at functions "load_bindata" and "zerostr" in the script, which may 
give you some ideas how it was done about 3 years ago. 


Cheers, 

Bish 



(09.8.4 3:39 AM), Leutenegger, Maurice A. (GSFC-662.0)[OAK RIDGE ASSOCIATED UNIVERSITIES (ORAU)] wrote:
> Hi Dan,
> 
> actually, loading the data is not my problem. (I think.)
> 
> I have my data in a fits file, and I can use fits_read_table to load it, followed by define_counts.
> 
> As far as I can tell, the problem is that I'm not in wavelength units, so no matter how I load the data, I can't use it.
> 
> Maurice
> ----
> You received this message because you are
> subscribed to the isis-users list.
> To unsubscribe, send a message to
> isis-users-request_at_email.domain.hidden> with the first line of the message as:
> unsubscribe

% File: findzero.sl
%
% Version: 2.0.0
%
% Authors: K. Tibbetts (originally linord.sl)
%	   Bish Ishibashi (modified as findzero.sl)	
% Date Created: 02 JUN 2003
%
% HISTORY:
%	   K. Tibbetts  - Written linord.sl as original.
%	   B. Ishibashi - Modified linord.sl to converge into 
%			  the better solution by ignoring 
%			  the scatter wing component. 
%	   K. Tibbetts  - replace the dmcopy call with 
%			  rotxy.sl.
%	   B. Ishibashi - Embedded rotxy.sl code into this one. 
%			  plus removing unused variables. 
%	   K. Tibbetts  - Added wrapper to linfit to centroid
%			  on x if angle is near 90 degrees.
%	   B. Ishibashi - Fix bugs introduced by KT. 
%	   B. Ishibashi - Add LETG setting. 
%	   B. Ishibashi - Overhaul the entire code and fitting schemes
%			  so that the fit is more accurate.  
%
%=========================
%
% 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)=findzero(evtfile_name,grating_type,guess_x,guess_y)
%
% Input:
%	   evtfile_name == Event L1.5 or 2 file. 
%	   grating_type == "h" (HEG), "m", or "l" (LEG). 
%	   guess_[x,y]  == initial guess for the zeroth order position in SKY coordinate. 
%			   the best if you can guess as accurately as possible. 
% Example
%	   (x,y)=findzero("acisf0001N0001_evt1.fits","h", 4096.5,4096.5)
%
%	 Then try to throw in the value (X,Y) as the new guess. Iterate this
%	 step until the value (X,Y) converges to more or less the same value. 
%
% Caveat:
%	Users need to give it a reasonable initial guess on zeroth order position. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


static define load_bindata(lo, hi, y, err)
%
% load_bindata: load the data array x and y into ISIS 
% so the loaded histogram data can be fitted with ISIS 
% functions. 
%
        {
                variable ff;
                ff=define_flux(lo,hi,y,err);
                plot_unit("A");
        }

static define zerostr(x,xmin, xmax, xn)
%
% Find a positon of zeroth order streak/fine-structure streak via Gaussian fitting. 
% Return the centroid position of Gaussian only.
%
	{
	variable lo, hi, hx, err, hm, lx, l,h, ids, cx;

	(lo,hi)=linear_grid(xmin,xmax,xn);
	hx  = histogram(x,lo,hi); err = (sqrt(hx+0.75));
	
	hm  = where(hx == max(hx)); lx  = lo[hm];

	load_bindata(lo,hi,hx,err); 
	fit_fun("poly(1)+gauss(2)");
	set_par(1,1,0,0,1e8); set_par(2,0,1); set_par(3,0,1);
        set_par(4,1,0,0,1e8);
        set_par(5,lx,0,lx-20,lx+20);
        set_par(6,1,0,0,1e8);

	lx = get_par("gauss(2).center");
        xnotice(1, lx - 20., lx + 20.); 
	()=renorm_flux;
	()=fit_flux;

        (l,h) = conf(6);
        while (l == h) 
                { 
                ()=conf(1,2);
                ()=conf(4,2);
                ()=conf(5,2);
                (l,h)=conf(6,2);
        }

        ()=fit_flux;
	ids = all_data;
        delete_data(ids);

	cx  = get_par("gauss(2).center");
	return cx; 
	}

static define findmed(x,y,xinc)
%
% This is to bin up the samples to locate a median value in data 
% distribution in each bin. 
%
 	{
	variable xmin,xmax, xn, xbin, ybin, i, h; 

	xmin=min(x); xmax=max(x);
	xn=typecast((xmax-xmin)/xinc,Integer_Type);
	xbin=xmin+xinc*0.5+xinc*[0:xn-1:1.0];
	ybin=[0:xn-1:1.0]*0.0;

	for (i=0; i<xn-1; i++)
		{
		h=where(x > xmin+i*xinc and x < xmin+(i+1)*xinc);
		ybin[i]=median(y[h]);
		}

	return xbin,ybin;
	}



static define linfit (x,y)
%
% A good ol' linear least square fit. 
%
  {
    variable ndata,ndatax,ndatay,s,sx,sy,sxx,sxy,delta,a,b,i,sol;

    (ndatax,,)=array_info(x); (ndatay,,)=array_info(y);

    if (ndatax[0] != ndatay[0]) 
	{
	message("x and y must have the same number of elements.\n");
	return 0;
	}
    ndata=ndatax;
    s=ndata; sx=0; sy=0; sxx=0; sxy=0;
    for (i=0; i<ndata[0]; i++)
	{
	sx=sx+x[i];  sy=sy+y[i]; sxx=sxx+x[i]*x[i]; sxy=sxy+x[i]*y[i];
	}
    delta=s*sxx-sx*sx; a=(sxx*sy-sx*sxy)/delta; b=(s*sxy-sx*sy)/delta;

    sol=[a,b];
  return sol;
  }


static define clip_array( y, nsig )
%
% data clipping routine. 
%
{
    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 ;
}


static define rotxy( x, y, x0, y0, angle )
%
% purpose:  rotate coords in CLOCK-WISE direction (and Bish means it) 
%
	{
   	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 ; 
	}


define findzero( )
  {
 	variable fitsfile, gtype, ftype;
  	if (_NARGS < 2 )
   		{
   		   message("USAGE: (x, y) = findzero(fits_file,\"l\",flag); ");
   		   message("    where \"l\" (LEG) can be replaced with \"h\" (HEG)");
   		   message("    or \"m\" (MEG). Flag is optional.");
   		   return -1;
   		}
   	if (_NARGS !=3 )
   		{
		  (fitsfile,gtype)=(); 
   		   message("Assuming that the RA_TARG and DEC_TARG values are correct");
		   message("  for this particular target.                            "); 
		   ftype = 1; 
   		}
	if (_NARGS ==3 )
		{
		(fitsfile,gtype,ftype)=();
		}

    variable gx, gy;
    variable boxl, boxw, rad; 
    variable rotkey="ROLL_NOM";
    variable grelem="GRATING";
    variable detnam="INSTRUME";
    variable ranom, denom, ratg,detg, rarf,derf;
    variable deg_rad=3.1416/180.;
%
% notation on angles: positive angle == clockwise rotation angle. 
%
    variable rotang, alpha1, alpha2, alpha3, acisscor, hrcscor, legcor, grang, gflag; 
    variable x,y,xrot,yrot,xrot2,yrot2;
    variable sx,sy, sx2, sy2, cnt, cnt2; 
    variable xrottmp,cltmp,l, xm, ym, strx, hegx, hegsl, solx, soly;

%
% First, check its grating configulation. 
%
    if (fits_key_exists(fitsfile,grelem))
	{
	grelem=fits_read_key(fitsfile+"[EVENTS]",grelem);

% 
% Then figure out the nominal roll angle. 
%
    	if (fits_key_exists(fitsfile,rotkey))
		{	
		(rotang,ranom,denom,ratg,detg,rarf,derf,detnam)=
			fits_read_key(fitsfile+"[EVENTS]",rotkey,"TCRVL11","TCRVL12",
				"RA_TARG","DEC_TARG","TCRPX11","TCRPX12",detnam);

%
% 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. 
%
	if (grelem == "HETG")
		{ 
			acisscor=0.03;  	% correction for ACIS-S rotation 
			alpha1= 4.755;   	% dispersion angle for MEG wrt CHIP_Y 
			alpha2=-5.205;  	% dispersion angle for HEG wrt CHIP_Y
			rotang=(rotang+acisscor);% compensating a small ACIS-S detector 
					        % rotation so that the streak is vertical. 
		if (gtype == "h")
			{
			grang=alpha2-acisscor;	% this is the corrected alpha clocking angle
			}			% in the final rotated frame (streak vertical). 
		if (gtype == "m")
			{
			grang=alpha1-acisscor; 
			}
			gflag=12; 
		} % grelem
%
% Case 2: LETG + HRC-S or ACIS-S. 
%
	if (grelem == "LETG")
		{
			gflag=3; 
		if (detnam == "HRC") 
			{
			(ranom,denom,ratg,detg,rarf,derf)=
			fits_read_key(fitsfile+"[EVENTS]","TCRVL19","TCRVL20",
				"RA_TARG","DEC_TARG","TCRPX19","TCRPX20");

			hrcscor=0.2744; 	% empirical hrc-s rotation correction
			alpha3=0.07; 		% dispersion angle for LEG wrt CHIP_Y 
			rotang=rotang+hrcscor; 
			legcor=-0.2744;		% again, empirical letg rotation correction
			grang=alpha3-hrcscor-legcor;
			}
		if (detnam == "ACIS")
			{
			acisscor=0.03;
			alpha3=0.07;
			rotang=rotang+acisscor;% compensating a small ACIS-S detector 
					       % rotation wrt Roll.
			legcor=0.089;         % ACIS-S or LETG arm is wicked and need
					       % an additional bogus correction factor. 
			grang=alpha3-legcor; 
			}
		if (detnam == "ACIS-I") % This is a dummy entry 
			{
			alpha3=0.07;
			legcor=0.0;
			}
		if (detnam == "HRC-I")  % This is a dummy entry 
			{
			alpha3=0.07;
			legcor=0.0;
			}
		}
	if (grelem == "NONE")
		{
		message("This is not a grating dataset.\n");
		return NULL, NULL;
		}

% First, a sanity check. 
%

	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 0,0;
		}


%
% 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. 
% 

	(x,y)=fits_read_col(fitsfile,"x","y");

	if (ftype == 1) 
		{
		if (detnam == "ACIS")
			{
			gx = -(atan((ratg - ranom)*deg_rad)/((0.492/3600.)*deg_rad))+rarf;
			gy =  (atan((detg - denom)*deg_rad)/((0.492/3600.)*deg_rad))+derf;
			}
		if (detnam == "HRC")
			{
			gx = -(atan((ratg - ranom)*deg_rad)/((0.1318/3600.)*deg_rad))+rarf;
			gy =  (atan((detg - denom)*deg_rad)/((0.1318/3600.)*deg_rad))+derf;
			}
		}
	if (ftype == 2)
		{
		%
		% just take a median in x & y. 
		%
		gx = median(x);gy= median(y);
		}
	if (ftype == 3) 
		{
		%
		% make a dumb guess
		%
		if (detnam == "ACIS") 
			{
			gx = 4096.5; gy = 4096.5;
			}
		if (detnam == "HRC")
			{
			gx = 32768.5;gy=32768.5;
			}
		}

	(xrot,yrot)=rotxy(x,y,gx,gy,-rotang*deg_rad);

    	boxl=1000.0;  boxw=500.0; rad=30.0;

	l=where(xrot>(gx-boxw/2) and xrot<(gx+boxw/2) and yrot>(gy-boxl/2) and 
		   yrot<(gy+boxl/2) and hypot(xrot-gx,yrot-gy)>rad);

	strx=zerostr(xrot[l],gx-100, gx+100,400);

%
% Now fit the grating arm. 
%

	boxl=5000.0;  boxw=400.0; rad=100.0;

	sx=strx;sy=gy;sx2=0.;sy2=0.;cnt=0; cnt2=0;

	()=open_plot("/xwin");resize(20,0.9);erase;

	connect_points(0);
	l=[0:30000:1];
	plot(xrot[l],yrot[l]);

	while ((sx2 != sx or sy2 != sy) and (cnt < 20)) 
		{
		sx2=sx; sy2=sy; 
		(xrot2,yrot2)=rotxy(xrot,yrot,sx,sy,-grang*deg_rad);
		l=where(xrot2>(sx-boxl/2) and xrot2<(sx+boxl/2) and yrot2>(sy-boxw/2) and 
		   yrot2<(sy+boxw/2) and hypot(xrot2-sx,yrot2-sy)>rad);

		(xm,ym)=findmed(xrot[l],yrot[l],100);
%		oplot(xm,ym);

		hegx=linfit(xm,ym);
		cltmp=clip_array(ym-xm*hegx[1]+hegx[0],2);
		hegx=linfit(xm[cltmp],ym[cltmp]);
		hegsl=hegx[1];

		variable flag=1;
		while (flag and cnt2 < 20) 
			{
			cltmp=clip_array(ym-xm*hegx[1]+hegx[0],2);
			hegx=linfit(xm[cltmp],ym[cltmp]);
			if (hegx[1] == hegsl) 
				{ 
				flag = 0;
				}
				hegsl=hegx[1];
				cnt2++;
			}
		if (cnt2 == 20)
			{
		message("the slope solution did not converge after "+string(cnt2)+" trials.\n");
			}
		connect_points(1);
		sx = strx; sy = hegx[1]*sx+hegx[0];
		cnt++; 
		}
		oplot([strx,strx],[min(yrot),max(yrot)]);
		oplot(xm,xm*hegx[1]+hegx[0]);

	if (cnt == 20)
		{
		message("the position solution did not converge after "+string(cnt)+" trials.\n");
		}
	
	(solx,soly)=rotxy([sx],[sy],gx,gy,rotang*deg_rad);

  return solx[0], soly[0];
 		} 
else return NULL,NULL; % no rotkey then return null
	}	  
else return NULL,NULL; % no grating then return null 
 }








----
You received this message because you are
subscribed to the isis-users list.
To unsubscribe, send a message to
isis-users-request_at_email.domain.hiddenwith the first line of the message as:
unsubscribe
Received on Mon Aug 03 2009 - 15:05:44 EDT

This archive was generated by hypermail 2.3.0 : Fri May 02 2014 - 08:35:46 EDT