ISIS Equivalent Width Function

In response to Juan Luna's request for an XSPEC-like equivalent width 
function for ISIS, the below attached code will do it.  It tries to 
replicate the XSPEC definition as closely as possible, which may or may 
not be a good thing, as that a) will calculate an equivalent width for 
*anything*, whether it's sensible or not, by merely subtracting off 
that component, differencing with the original model, and doing a 
ratio, and b) can depend upon the energy/wavelength range covered by 
the data (or by the energy range chosen by the user).

This version will spit out the equivalent width in both eV and mA, and 
give negative numbers for absorption, positive for emission.

Since ISIS doesn't store pieces of the fit (since it doesn't 
distinguish between additive or multiplicative models), it has to 
reevaluate the model at least once, preferably twice (to get yourself 
back to what the model was before, for plotting purposes, for example). 
  Thus, if you have especially time-consuming models, you might want to 
think before applying it.  (That's why I made the reevaluation at the 
end an option, not a necessity.)

-Mike

public define msg(str_array)
{
    () = printf("\n");
    foreach(str_array)
    {
       variable str = ();
       () = printf(" %s\n", str);
    }
    () = printf("\n");
    return;
}

public define eqw()
{
    variable indx, par, pmin=0, lo=NULL, hi=NULL, eva, eve;
    variable fv_hold = Fit_Verbose, noeval=0;

    switch(_NARGS)
    { case 2:
       (indx, par) = ();
    }
    { case 3:
       (indx, par, pmin) = ();
    }
    { case 4:
       (indx, par, pmin, lo) = ();
    }
    { case 5:
       (indx, par, pmin, lo, hi) = ();
    }
    { case 6:
       (indx, par, pmin, lo, hi, noeval) = ();
    }
    {
       variable str = [
       "(ew_ma, ew_ev) = eqw(indx, par, [pmin,] [lo, hi], [noeval]); "," 
",
       "  Calculate the equivalent width for a model component, in 
milli-angstrom and eV.",
       "  ",
       "  Inputs:",
       "   indx    : Data set index to which the model is applied",
       "   par     : Number for, *or* string with the name of, the 
*normalization*",
       "             parameter for the model component for which the EW 
will be calculated",
       " ",
       "  Optional Inputs:",
       "   pmin    : For parameter par, what value it should be set to 
when calculating",
       "             the continuum flux without the line (default = 0, 
but this allows ",
       "             other parameter toggles to be used).",
       "   lo, hi  : *Wavelength* ranges over which to restrict the",
       "             evaluation of the equivalent width. (Otherwise, 
uses",
       "             the full wavelength range of the arf; see 
get_model_flux();)"," ",
       "   noeval  : If = 1, don't do the final eval_counts (parameters 
will be reset",
       "             to initial values, but not evaluated - i.e., plots 
won't show",
       "             the fit with the line, and may otherwise be 
screwy)"," ",
       "  Outputs:",
       "   ew_ma, ew_ev : Equivalent widths in milli-angstrom and eV"];
       msg(str);
       return;
    }

    if( pmin == NULL )
    {
       pmin = 0;
    }

    Fit_Verbose = -1;
    () = eval_counts;

    variable a = get_model_flux(indx);

    variable phold = get_par_info(par);

    set_par(par,pmin,0,0,phold.max);
    () = eval_counts;
    variable b = get_model_flux(indx);

    set_par(par,phold.value,phold.freeze,phold.min,phold.max);
    if(noeval !=0)
    {
       () = eval_counts;
    }

    variable iw, iwa, iwe;
    if( lo != NULL and hi != NULL )
    {
       iw = where(a.bin_lo >= lo and a.bin_hi <= hi);
    }
    else
    {
       iw = [0:length(a.bin_lo)-1];
    }

    variable diff = a.value - b.value;
    variable bdena = b.value/(a.bin_hi-a.bin_lo);
    variable bdene = 
reverse(reverse(b.value)/(_A(a.bin_lo)-_A(a.bin_hi)));
    variable fdena = diff/(a.bin_hi-a.bin_lo);
    variable fdene = reverse(reverse(diff)/(_A(a.bin_lo) - 
_A(a.bin_hi)));
    variable flux = sum( diff[iw] );

    iwa = min(where( abs(fdena[iw]) == max(abs(fdena[iw])) ));
    iwe = min(where( abs(fdene[iw]) == max(abs(fdene[iw])) ));

    eva = flux/bdena[iw[iwa]]*1000.;
    eve = flux/bdene[iw[iwa]]*1000.;

    Fit_Verbose = fv_hold;

    return eva, eve;
}

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