%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
%% QPO (Lorentzian) %% 
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%

% Define a s-lang qpo that works better. Properly averages over bin,
% and thus yields a smoother fit.

define qpo_fit(lo,hi,par)
{
   variable l,r,q,f,al,ah;

   % Go from Angstrom to keV (which we pretend is Fourier Hz),
   % Hence it's F[(al-f)/f] - F[(ah-f)/f] below...

   al = _A(lo);
   ah = _A(hi);

   q = par[1];
   f = par[2];

   % In this formulation, par[0] is rms amplitude

   r = par[0]/(0.5 - atan(-2.*q)/PI);

   l = r^2/(al-ah)/PI * 
       ( atan(2.*q*(al-f)/f) - atan(2.*q*(ah-f)/f) );

   l = reverse(l);

   return l;
}

add_slang_function("qpo",["norm [rms]","Q [f/FWHM]","f [Hz]"]);

define qpo_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (1,0,0.01,1.e3);
   }
   {case 2:
      return (1,0,0.,1.e4);
   }
}

set_param_default_hook("qpo","qpo_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Zero Freq. Centered Lorentzian  %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

define zfc_fit(lo,hi,par)
{
   variable a,f,al,ah,l;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),
   % Hence it's atan(al/f) - atan(ah/f) below...

   al = _A(lo);
   ah = _A(hi);
   f = par[1];
   
   % In this formulation, par[0] is RMS

   a = par[0]^2 * PI/2./f;

   l = a*f/(al-ah) *  ( atan(al/f) - atan(ah/f) );

   l = reverse(l);

   return l;
}

add_slang_function("zfc",["norm [rms]","f [Hz]"]);

define zfc_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (1,0,1.e-6,1.e3);
   }
}

set_param_default_hook("zfc","zfc_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counts/Bin Powerlaw  %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

define plaw_fit(lo,hi,par)
{
   variable a,s,al,ah,l;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),

   al = _A(lo);
   ah = _A(hi);

   a = par[0];
   s = par[1];
   
   if(s != -1)
   {
      l = a/(s+1) * (al^(s+1) - ah^(s+1));
   }
   else
   {
      l = a * ( log(al) - log(ah) );
   }

   l = reverse(l/(al - ah));

   return l;
}

add_slang_function("plaw",["norm","slope"]);

define plaw_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (-2,0,-5,5);
   }
}

set_param_default_hook("plaw","plaw_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Broken Counts/Bin Powerlaw  %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

define bkn_plaw_fit(lo,hi,par)
{
   variable a,sa,b,sb,al,ah,aa,l;
   variable iwl,iwh,liwl,liwh;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),

   al = _A(lo);
   ah = _A(hi);
   aa = (al+ah)/2.;
   l = @al;

   a = par[0];
   sa = par[1];
   b = par[2];
   sb = par[3];

   iwl = where(aa<= b);
   iwh = where(aa > b);

   liwl = length(iwl);
   liwh = length(iwh);
   
   if(liwl != 0)
   {
      l[iwl] = a * aa[iwl]^sa;
   }

   if(liwh != 0)
   {
      l[iwh] = a * b^sa * ( aa[iwh]/b )^sb;
   }

   l = reverse(l);

   return l;
}

add_slang_function("bkn_plaw",["norm","slope_1","f_break [Hz]","slope_2"]);

define bkn_plaw_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (0,0,-5,5);
   }
   {case 2:
      return (1,0,0.,1.e4);
   }
   {case 3:
      return (-2,0,-5,5);
   }
}

set_param_default_hook("bkn_plaw","bkn_plaw_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Double Broken Counts/Bin Powerlaw  %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

define dbkn_plaw_fit(lo,hi,par)
{
   variable a,sa,ba,sb,bb,sc,al,ah,aa,l;
   variable iwl,iwm,iwh,liwl,liwm,liwh;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),

   al = _A(lo);
   ah = _A(hi);
   aa = (al+ah)/2.;
   l = @al;

   a = par[0];
   sa = par[1];
   ba = par[2];
   sb = par[3];
   bb = par[4];
   sc = par[5];

   iwl = where(aa <= ba);
   iwm = where(aa > ba and aa <= bb);
   iwh = where(aa > bb);

   liwl = length(iwl);
   liwm = length(iwm);
   liwh = length(iwh);
   
   if(liwl != 0)
   {
      l[iwl] = a * aa[iwl]^sa;
   }

   if(liwm != 0)
   {
      l[iwm] = a * ba^sa * ( aa[iwm]/ba )^sb;
   }

   if(liwh != 0)
   {
      l[iwh] = a * ba^sa * (bb/ba)^sb * (aa[iwh]/bb)^sc;
   }

   l = reverse(l);

   return l;
}

add_slang_function("dbkn_plaw",["norm","slope_1","f_brk_1 [Hz]",
                   "slope_2","f_brk_2 [Hz]","slope_3"]);

define dbkn_plaw_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (0,0,-5,5);
   }
   {case 2:
      return (1,0,0.,1.e4);
   }
   {case 3:
      return (-1,0,-5,5);
   }
   {case 4:
      return (10,0,0.,1.e4);
   }
   {case 5:
      return (-2,0,-5,5);
   }
}

set_param_default_hook("dbkn_plaw","dbkn_plaw_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counts/Bin Gaussian Normalized to RMS  %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#ifexists erf

define rms_gauss_fit(lo,hi,par)
{
   variable r,f,sig,al,ah,l;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),

   al = _A(lo);
   ah = _A(hi);

   r = par[0];
   f = par[1];
   sig = sqrt(2.)*par[2];

   l = r^2/2 * ( erf( (al-f)/sig ) - erf( (ah-f)/sig ) );
   l = reverse(l/(al-ah));

   return l;
}

add_slang_function("rms_gauss",["norm [rms]","f_0 [Hz]","sigma [Hz]"]);

define rms_gauss_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (1,0,0,1.e4);
   }
   {case 2:
      return (0.2,0,1.e-4,1.e3);
   }
}

set_param_default_hook("rms_gauss","rms_gauss_defaults");

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counts/bin Sine Wave %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

define sinwave_fit(lo,hi,par)
{
   variable a,f,phs,al,ah,l;
 
   % Go from Angstrom to keV (which we pretend is Fourier Hz),

   al = _A(lo);
   ah = _A(hi);

   a = par[0];
   f = par[1];
   phs = par[2];

   l = a * ( -cos(2*PI*(f*al-phs)) + cos(2*PI*(f*ah-phs)) )/2/PI/f/(al-ah);
   l = reverse(l);

   return l;
}

add_slang_function("sinwave",["norm","freq","phase"]);

define sinwave_defaults(i)
{
   switch(i)
   {case 0:
      return (0.1,0,0.,1.e8);
   }
   {case 1:
      return (1,0,0,1.e4);
   }
   {case 2:
      return (0,0,0.,1.);
   }
}

set_param_default_hook("sinwave","sinwave_defaults");