define apec() {  % T, abund, norm, binlo, binhi
   %
   %  Uses the built-in isis routine model_spectrum.
   %
   % flux = model_spectrum (bin_lo, bin_hi); 
   %
   variable T, abund, norm, binlo, binhi, inpbinlo, inpbinhi;
   variable i, index;
   
   if ((_NARGS != 4) and (_NARGS != 5)) {
      message("Call as:");
      message("apec(T, abund, norm, bin_lo, bin_hi) where");
      message("T = temperature (in keV)");
      message("abund = abundance (in solar units)");
      message("norm = emission measure / 4 pi D^2");
      message("bin_lo and bin_hi are in Angstroms");
      _pop_n(_NARGS);
      return NULL;
   }
   
   if (_NARGS == 4) { 
      %message("Warning: Code meant to run with integrate ON, not off.");
      (T, abund, norm, inpbinlo) = ();
      % 
      %  Check to see if backward
      %
      if (inpbinlo[1] < inpbinlo[0]) { 
	 index = [length(inpbinlo)-1:0:-1];
      } else {
	 index = [0:length(inpbinlo)-1:1];
      }
      binlo = inpbinlo[index];
      binhi = @binlo;
      for (i=0;i<length(binlo)-1;i++) binhi[i] = binlo[i+1];
      binhi[length(binlo)-1] = binlo[length(binlo)-1] + 0.001;
   }
   if (_NARGS == 5) {
      (T, abund, norm, inpbinlo, inpbinhi) = ();
      if (inpbinhi[0] < inpbinlo[0]) { 
	 index = [length(inpbinlo)-1:0:-1];
      } else {
	 index = [0:length(inpbinlo)-1:1];
      }
      binlo = inpbinlo[index];
      binhi = inpbinhi[index];
   }
   
   variable state = default_plasma_state();
   
   state.norm = norm;
   state.temperature = T/8.617e-8; 
   state.metal_abund = abund;
   
   define_model([state]);
   
   variable isis_result = model_spectrum(binlo, binhi);
   variable result = isis_result[index];

   if (Float_Type == _typeof(inpbinlo)) result = typecast(result, Float_Type);
   
   return result;
}
   
() = sherpa_register_model("apec",["kT","abund","norm"], 1,
                           [1.0 , 1.0, 1.0],  % default values
                           [0.01, 0.0, 0.0],  % Minimum values
                           [80. ,100.,1e10],  % Maximum values
                           [1   ,   0,   1]); % Float temp and norm