source: ETALON/CLIO/Analysis/Data2May/KK/SRf.m @ 721

function [RSP,ISP,alp]=SRf(SP,SS)
%RECOVERY PART
w0=SP(1);
%Make LF extrapolation
DSP=diff(SS)./diff(SP);
rs=sum(DSP(1:3));%must be negative
r=SS(1);
alp=rs/r/w0-log(r)/(w0*w0);
bet=2*log(r)/w0-rs/r;
x=0:(w0-1);
% syms xx
% w0
% rs
% vpasolve(-2*w0*xx*exp(-xx*w0^2)==rs, xx)
%  LFpart=exp(alp.*x.*x+bet.*x);
LFpart=exp(-abs(log(r)/(w0^2)).*x.^2);
 LFpart( LFpart>1)= 1; 
%  LFpart=ones(size(x));
% LFpart=exp(-1.8262e-05*x.^2);
%  LFpart=1;
%Make interpolation
% SS=SS.*exp(-1.8262e-05*w0.^2)
ISP = pchip(SP,SS,w0:max(SP));
%Make HF extrapolation
w=1:(2^12-length(ISP)-length(LFpart));
wmax=max(SP);
rosh=DSP(end);

  alp=(rosh/ISP(end))*wmax;
        c=0;
        while (alp>0)
            c=c+1;
          rosh=  (SS(end)-SS(end-c))/(SP(end)-SP(end-c));
          alp=(rosh/ISP(end))*wmax;
        end
        if (alp>-2) alp=-2;end
HFpart=ISP(end)* (1+w/wmax).^(alp);
HFpart(HFpart<eps)=eps;

%Connect all together
ISP=[LFpart,ISP,HFpart];
%Mirror image of the spectrum
RSP=[ISP,fliplr(ISP)];