Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/ChandraSNR/G010.0-00.3/746/work/analyze_ltcrv.sl Дата изменения: Mon Aug 15 22:12:02 2005 Дата индексирования: Tue Oct 2 11:55:00 2012 Кодировка: Поисковые слова: annular solar eclipse

%
% $Id: analyze_ltcrv.sl,v 1.1 2003/06/08 22:40:52 jhagler Exp $
%
% Usage:
% analyze_ltcrv( filename )
%
% where filename is the output from the CIAO lightcurve program
% (or any file with TIME and RATE columns)
%
% Example:
% chips> () = evalfile( "analyze_ltcrv.sl" );
% chips> analyze_ltcrv( "lc.fits" );
%
% Thread:
% http://asc.harvard.edu/ciao/threads/filter_ltcrv.thread.html
%
% Aim:
% Calculate the time periods when a lightcurve lies within
% +- 3 sigma of the mean level. Displays these periods both
% in graphical form (in a ChIPS window with green/red points
% showing allowed/rejected times), and as a list of times that
% can be used in a DM filter, or passed to dmgti.
%
% The rejection of points uses a sigma-clipping algorithm,
% iterating through the dataset until no more points are
% removed.
%
% If you wish to clip using a different number of standard
% deviations, change the stddev_limit variable in the
% analyze_ltcrv() function.
%
%

define analyze_ltcrv() { % filename

%
% Takes a lightcurve, reads it in, finds the average and standard
% deviation, then removes all points more than X standard deviations
% from the average and refits. Finally, gives a list of all excluded
% times suitable for dmgti.
%

variable filename, done=0;
variable stddev_limit = 3;
variable ltcrv_data=NULL;
variable lt_avg, lt_std, sigmas;
variable minLength=3; % At least three points in a row must be good.
variable bp, gp, ap, gp_i, i;
variable fieldname, hasrate=0, hascounts=0;

if (_NARGS != 1) {
message("Call as analyze_ltcrv(filename)");
message("where filename is the output from the CIAO lightcurve routine.");
if (_NARGS != 0) _pop_n(_NARGS);
return NULL;
}
filename = ();

variable ltcrv = readfile(filename);
if ( ltcrv == NULL ) {
error( "analyze_ltcrv() stopped as unable to read data from file" );
}

%
% Search for "rate" field, or if not available "counts"
%
foreach (get_struct_field_names (ltcrv)) {
fieldname = ();
if (strcmp (fieldname, "RATE") == 0) hasrate = 1;
if (strcmp (fieldname, "COUNTS") == 0) hascounts = 1;
}

if (hasrate == 1) {
ltcrv_data = @ltcrv.RATE;
} else {
if (hascounts == 1) {
ltcrv_data = @ltcrv.COUNTS;
} else {
vmessage("Datafile %s does not have RATE or COUNTS column; is it a lightcurve file?",filename);
return NULL;
}
}

ap = Integer_Type [length(ltcrv_data)] + 1; % 1's and 0's indicating OK pts

while (done == 0) {
done = 1;
gp = where(ap);

lt_avg = sum( ltcrv_data[gp] ) / length(gp);
lt_std = sqrt( sum(ltcrv_data[gp]*ltcrv_data[gp])/length(gp) - lt_avg*lt_avg );

% vmessage("Average : %f, Std Dev : %f",lt_avg,lt_std);
sigmas = (ltcrv_data - lt_avg)/lt_std;
for (i=0;i if ((abs(sigmas[i]) > stddev_limit)and(ap[i] == 1)) {
ap[i] = 0;
done = 0;
}
}
%
% At this point, gp should be all points that survived.
%

bp = where(1-ap);

() = chips_eval( "clear" );
() = chips_eval( "redraw off" ); % stop flickering

() = curve(ltcrv.TIME[gp], ltcrv_data[gp]);
() = chips_eval("c 1 noline");
() = chips_eval("c 1 symbol bigpoint");
() = chips_eval("c 1 symbol size 1");
() = chips_eval("c 1 symbol green");

if (length(bp) != 0) {
() = curve(ltcrv.TIME[bp], ltcrv_data[bp]);
() = chips_eval("c 2 noline");
() = chips_eval("c 2 symbol bigpoint");
() = chips_eval("c 2 symbol size 1");
() = chips_eval("c 2 symbol red");
}

() = chips_eval( "xlabel 'Time (s)'" );
() = chips_eval( "ylabel 'Rate (/s)'" );
% protect any underscores in the filename
variable title;
( title, ) = strreplace( filename, "_", "\\_", strlen(filename) );
() = chips_eval( "title 'Lightcurve for: " + title + "'" );
() = chips_eval( "title size 1.2" );

() = chips_eval( "redraw on" ); % stop flickering

%
% Now, need to print out the good time intervals
%
variable lowtime=0, hightime=0, findlow=1, iLow;

%%%% LUCA
variable fp = fopen ("lc_bg.time_for_dmcopy", "w");

%%%% LUCA
fprintf (fp, "0:0 ");


for (i=0;i if (ap[i] == 1) { % A good point
if (findlow == 1) {
lowtime = ltcrv.TIME[i];
findlow = 0;
iLow = i;
}
} else {
if (findlow == 0) { % End of good points
hightime = ltcrv.TIME[i];
findlow = 1;
if (i - iLow >= minLength) {
vmessage("((time > %f) && (time < %f)) ; %5.2f ksec",
lowtime, hightime, (hightime-lowtime)/1.e3);

%%%% LUCA
() = fprintf (fp, ",%f:%f", lowtime, hightime);

}
}
}
}
% only worry about end of array if we've found some good points
if ( findlow == 0 ) {
i = i-1;
if (i - iLow >= minLength) {
hightime = ltcrv.TIME[i];
vmessage("((time > %f) && (time < %f)) ; %5.2f ksec",
lowtime, hightime, (hightime-lowtime)/1.e3);

%%%% LUCA
() = fprintf (fp, ",%f:%f", lowtime, hightime);

}
}
}