source: ETALON/reconstruction/f18/SRf.m

function [RSP,ISP,alp]=SRf(SP,SS)
%RECOVERY PART
w0=SP(1);
%Make LF extrapolation
DSP=diff(SS)./diff(SP);
rs=DSP(1);%must be negative
r=SS(1);
alp=rs/r/w0-log(r)/(w0*w0);
bet=2*log(r)/w0-rs/r;
x=1:(w0-1);
LFpart=exp(alp.*x.*x+bet.*x-alp-bet);
%Make interpolation
ISP = pchip(SP,SS,w0:max(SP));
%Make HF extrapolation

w=1:(2^15-length(ISP)-length(LFpart));
wmax=max(SP);
rosh=mean(DSP((end-4):end));
alp=(rosh/SS(end))*wmax;
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)];