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

%------------------------------------------------------------------------

% sum() computes the sum of elements of an array.

define sum (x)
{
variable num = length(x);
variable total, i;

total = 0;
_for (0, num-1, 1)
{
i = ();
total += x[i];
}

return total;
}



% integrate_binned() assumes the following:
% f = struct {val, elo, ehi}
% f.val[k] = \int_f.elo[k]^f.ehi[k] g(e) de
% f.elo[k] < f.ehi[k]
% f.elo[k] < f.elo[k+1]
%
% Integration limits:
% a < b
%
% Usage: result = integrate_binned (f, a, b);
%
% result is the integral of 'f' from a -> b.
%
private define integrate_binned (f, a, b)
{
variable i, result;

i = where ((a < f.elo) and (f.ehi < b));
variable l = length(i);

if (l == 0)
{
result = 0;
return result;
}

result = sum(f.val[i]);

% Include the endpoint corrections:

variable k, ni, nf, frctn;

k = i[0];
if (k > 0)
{
frctn = (f.ehi[k-1] - a) / (f.ehi[k-1] - f.elo[k-1]);
result += f.val[k-1] * frctn;
}

nf = length(f.val);
ni = length(i);

k = i[ni-1];
if (k < nf - 1)
{
frctn = (b - f.elo[k+1]) / (f.ehi[k+1] - f.elo[k+1]);
result += f.val[k+1] * frctn;
}

return result;
}

typedef struct
{
elo, ehi
}
Energy_Band_Type;

private define init_bands (elo, ehi)
{
variable dest = @Energy_Band_Type;
dest.elo = @elo;
dest.ehi = @ehi;

return dest;
}

private variable Standard_Bands = Assoc_Type [Energy_Band_Type];

% The band positions listed here are for illustration only --
% the real values should be chosen using the effective area curves.
%
Standard_Bands["ACIS"] = init_bands ([0.5, 1.0, 1.8, 2.0, 4.0, 6.0, 8.0],
[1.0, 1.8, 2.0, 4.0, 6.0, 8.0, 9.9]);
Standard_Bands["HRC"] = init_bands ([0.5, 2.0, 4.0, 8.0],
[2.0, 4.0, 8.0, 10.0]);
Standard_Bands["SNR"] = init_bands ([0.3, 0.42, 0.54],
[0.42, 0.54, 0.6]);

public define get_bands_with_energy_limit(lo_e, hi_e)
{

variable Boundaries = [0.31, 0.42, 0.54, 0.8, 1.58, 2.08, 4., 6., 9.9];
variable Band, Band1 ;

variable low, high;

variable l, i;

% hi_e = ();
% lo_e = ();

Band = Boundaries[where( Boundaries > lo_e and Boundaries < hi_e)];
% Band = where( Band1 < hi_e);

l = length(Band);

if ( l == 0)
{
Band = [(lo_e+hi_e)/2.];
l = 1;
}

% vmessage ("-- %e %e\n", lo_e, hi_e);

low = Double_Type [l + 1];
high = Double_Type [l + 1];

% for (i = 0; i <= l-1; i++) vmessage (" %e ", Band[i]);;
% vmessage ("-- %d\n", l);

low[0] = lo_e;
low[[1:l]] = Band[[0:l-1]];
high[[0:l-1]] = Band[[0:l-1]];
high[l] = hi_e;

for (i = 0; i <= l; i++) vmessage (" %e ", low[i]);;
vmessage ("-- %d\n", l);

for (i = 0; i <= l; i++) vmessage (" %e ", high[i]);;
vmessage ("-- %d\n", l);

return init_bands(low, high);

}

public define get_standard_bands ()
{
variable detector;

if (_NARGS != 1)
{
message ("Usage: bands = get_standard_bands (\"detector\")");
return;
}

detector = ();

if (0 != assoc_key_exists (Standard_Bands, detector))
return Standard_Bands[detector];

vmessage ("Unrecognized detector name: %s", detector);
}

% calc_instmap_weights assumes that
% f = struct {val, elo, ehi}
% bands = struct {elo, ehi}
%
% Usage: weights = calc_instmap_weights (f, bands);
%
public define calc_instmap_weights ()
{
if (_NARGS != 2)
{
message ("Usage: weights = calc_instmap_weights (spectrum, bands)");
return;
}

variable f, bands;
(f, bands) = ();

variable i, n, weights, norm;

n = length (bands.elo);
weights = Float_Type[n];

_for (0, n-1, 1)
{
i = ();
weights[i] = integrate_binned (f, bands.elo[i], bands.ehi[i]);
}

norm = sum(weights);

if (norm > 0.0)
return weights / norm;

vmessage ("Invalid weights: sum = %f", norm);
return weights;
}

public define write_weights_file ()
{
variable filename, weights, bands;

if (_NARGS != 3)
{
message ("Usage: write_weights_file (\"filename\", weights, bands)");
return;
}

(filename, weights, bands) = ();
variable fp = fopen (filename, "w");

if (NULL == fp)
{
vmessage ("Failed opening file %s for writing", filename);
return;
}

variable i,n = length(weights);

_for (0, n-1, 1)
{
i = ();
() = fprintf (fp, "%10.3e %10.3e\n",
0.5 * (bands.elo[i] + bands.ehi[i]), weights[i]);
}

() = fclose (fp);
}

private define load_model_spectrum (file)
{
variable s = readascii(file); % read the sherpa model file
variable n = s._nrows;

variable f = struct {val, elo, ehi};

% assuming col2 is the bin-integral:

f.elo = s.col1[[0:n-2]]; % bin low edge
f.ehi = s.col1[[1:n-1]]; % bin high edge
f.val = s.col2[[1:n-1]]; % bin integral

return f;
}

public define runtest2 (file, outfile, e1, e2)
{
variable f, b, w;
f = load_model_spectrum (file);
b = get_bands_with_energy_limit (e1, e2);
w = calc_instmap_weights (f, b);
write_weights_file (outfile, w, b);
}