tphfit_dke_yp

    Determination of the photon temperature from the energy spectrum,
   taking into account of an appropriate weight for each energy channel
   log(Ei*Ni/dE) = a+b*Ei -> Tph = -1/b
   See Numerical Recipes, W H. Press et al., Cambridge University Press, 2nd Edition, 1994

    Input:

        - Ei: Photon energy grid (keV) [n,m]
        - Ni: Number of counts (Counts) [n,m]
        - dNi: Uncertainty of the number of counts (Counts) [n,m]

    Output:

        - Tph: Estimated photon temperature (keV) [1,m]
        - eT: Estimated absolute error on the photon temperature (keV) [1,m]
        - A: Amplitude of the energy spectrum (Counts) [1,m]
        - eA: Estimated absolute error on the amplitude of the energy spectrum (Counts) [1,m]
       - Q: Goodness-of-fit of the data to the model (Q > 0.1 for the fit to be believable) [1,m]

by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>

0001 function [Tph,eTph,A,eA,Q] = tphfit_dke_yp(Ei,Ni,dNi)
0002 %
0003 %    Determination of the photon temperature from the energy spectrum,
0004 %   taking into account of an appropriate weight for each energy channel
0005 %   log(Ei*Ni/dE) = a+b*Ei -> Tph = -1/b
0006 %   See Numerical Recipes, W H. Press et al., Cambridge University Press, 2nd Edition, 1994
0007 %
0008 %    Input:
0009 %
0010 %        - Ei: Photon energy grid (keV) [n,m]
0011 %        - Ni: Number of counts (Counts) [n,m]
0012 %        - dNi: Uncertainty of the number of counts (Counts) [n,m]
0013 %
0014 %    Output:
0015 %
0016 %        - Tph: Estimated photon temperature (keV) [1,m]
0017 %        - eT: Estimated absolute error on the photon temperature (keV) [1,m]
0018 %        - A: Amplitude of the energy spectrum (Counts) [1,m]
0019 %        - eA: Estimated absolute error on the amplitude of the energy spectrum (Counts) [1,m]
0020 %       - Q: Goodness-of-fit of the data to the model (Q > 0.1 for the fit to be believable) [1,m]
0021 %
0022 %by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>
0023 %
0024 if nargin < 3,
0025     error('Wrong number of input arguments for tphfit_dke_yp');
0026 end
0027 %
0028 sNi = size(Ni);
0029 sEi = size(Ei);
0030 sdNi = size(dNi);
0031 if sEi(1)~=sNi(1) | sEi(1)~=sdNi(1) ,
0032     error('Size of data are not consistent !!!');
0033 end
0034 %
0035 Tph = zeros(1,sNi(2));
0036 eTph = zeros(1,sNi(2));
0037 A = zeros(1,sNi(2));
0038 eA = zeros(1,sNi(2));
0039 Q = zeros(1,sNi(2));
0040 chi2 = zeros(1,sNi(2));
0041 %
0042 dE = Ei(2,1)-Ei(1,1);
0043 yi = log(Ei.*Ni/dE);
0044 pi = (Ni./dNi).^2;% weigth
0045 xpiyi = sum(pi.*yi);
0046 xpiEiyi = sum(pi.*Ei.*yi);
0047 %
0048 xpi = sum(pi);
0049 xpiEi = sum(pi.*Ei);
0050 xpiEi2 = sum(pi.*Ei.*Ei);
0051 %
0052 %    Determination of the least square fit for each energy spectrum
0053 %
0054 for i =1:length(xpi),
0055     S = [xpi(i),xpiEi(i);xpiEi(i),xpiEi2(i)];
0056     Y = [xpiyi(i);xpiEiyi(i)];
0057     X = S\Y;
0058     E = sqrt([xpiEi2(i);xpi(i)]/det(S));
0059 %
0060     Tph(i) = -1.0/X(2);
0061     eTph(i) = E(2)*Tph(i)*Tph(i);    
0062     A(i) = exp(X(1));
0063     eA(i) = A(i)*E(1);
0064 %
0065 %    goodness-of-fit of the data to the model (Q > 0.1 for the fit to be believable)
0066 %
0067     chi2(i) = sum(((yi(:,i) - X(1) - X(2)*Ei(:,i))./pi(:,i)).^2);
0068     if chi2(i) >= 0,
0069         Q(i) = gammainc((sEi(1)-2)/2,chi2(i)/2);
0070     else
0071         Q(i) = NaN;
0072     end
0073 %
0074     clear X Y S
0075 end
0076 
0077 
0078 
0079 
0080