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: unsubscribeReceived 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