
function [ap_TSalScatCorr,cp_corr, wl] = ACS_ResidTempScatCorr(ap_uncorr, cp, wavel)

wavel = wavel(:)';

[nspect, nwavel] = size(ap_uncorr);

if size(ap_uncorr)~=size(cp)
    error('ap,cp size mismatch')
end
if length(wavel)~=nwavel
    error('Invalid lambda')
end
for i=1:nspect
    I=find(isfinite(ap_uncorr(i,:))==1);
    if isempty(I)
        i
    else
        break
    end
end
wavel=wavel(I);
ap_uncorr=ap_uncorr(1:nspect,I);
cp=cp(1:nspect,I);

tmp = xlsread('Sullivan_etal_2006_instrumentspecific.xls');
psiT = interp1(tmp(:,1),tmp(:,2), wavel);

opts = optimset('fminsearch');
%opts = optimset(opts,'NonlEqnAlgorithm', 'gn');
opts = optimset(opts,'MaxIter',20000000);
opts = optimset(opts,'MaxFunEvals',20000);
opts = optimset(opts,'TolX',1e-8);
opts = optimset(opts,'TolFun',1e-8);

NIR = find(wavel>=710 & wavel<=750);  % spectral srange for optimization (710 to 750nm)
ref = find(wavel>=715, 1,'first');
VIS = find(wavel<=700 & wavel>=450);
ref2 = find(wavel>=650, 1,'first');
II = find(wavel<=700);

ap_TSalScatCorr = nan(size(ap_uncorr));
cp_corr = nan(size(cp));
deltaT = nan(size(ap_uncorr,1),1);
fiterr = nan(size(ap_uncorr,1),1);
fiterr1 = nan(size(ap_uncorr,1),1);

wl=wavel;

for k = 1:nspect
    
    if isfinite(ap_uncorr(k,:))==1
        % guess for paramters (beamc at lambda0, beamc slope)
        delT = 0;
        deltaT(k)=0;
        % minimization routine
        if isfinite(cp(k,ref2))
            [deltaT(k), fiterr(k)] = fminsearch(@f_TS, 0, opts, ap_uncorr(k,:), psiT, cp(k,:), NIR, ref);
            [x1, fiterr(k)] = fminsearch(@least_squares_cp, [cp(k,ref2) 1], opts, cp(k,VIS), wavel(VIS));
            [deltaTc(k), fiterr(k)] = fminsearch(@f_TSc, 0, opts, cp(k,:)-x1(1).*(650./wl).^x1(2), psiT, NIR);
            cp_corr(k,:)=cp(k,:)-psiT(ref).*deltaTc(k);
             bp = cp_corr(k,:) - ap_uncorr(k,:);
           [x1, fiterr1(k)] = fminsearch(@least_squares_cp, [cp(k,ref2) 1], opts, cp(k,II), wavel(II));
            x1
            fiterr1(k)
            if (x1(2)>0.2  & fiterr1(k) < 0.004)
                k
                ap_TSalScatCorr(k,:) = ap_uncorr(k,:) - psiT.*deltaT(k) - ...
                    ((ap_uncorr(k,ref)-psiT(ref).*deltaT(k))./bp(ref)).*bp;
                
            else
                [deltaTa(k), fiterr(k)] = fminsearch(@f_TSc, 0, opts, ap_uncorr(k,:), psiT, NIR);
                deltaTa(k)=0;
                ap_TSalScatCorr(k,:) = ap_uncorr(k,:) - psiT.*deltaTa(k)-(ap_uncorr(k,ref)-psiT(ref).*deltaTa(k));
                cp_corr(k,:)=NaN;
            end
            
        elseif isfinite(cp(k,ref2)) == 0
                [deltaTa(k), fiterr(k)] = fminsearch(@f_TSc, 0, opts, ap_uncorr(k,:), psiT, NIR);
                deltaTa(k)=0;
                ap_TSalScatCorr(k,:) = ap_uncorr(k,:) - psiT.*deltaTa(k)-(ap_uncorr(k,ref)-psiT(ref).*deltaTa(k));            
                cp_corr(k,:)=NaN;
        end
        figure(4)
        subplot(1,2,1)
        if isfinite(cp_corr(k,:))
            plot(wl,cp(k,:),'-',wl(II),x1(1).*(650./wl(II)).^x1(2),'r');
        end
        subplot(1,2,2)
        plot(wl,ap_TSalScatCorr(k,:),'-');
        %             pause
        disp(['  Residual T,S,Scat correction: k=' ...
            num2str(k) ' deltaT=' num2str(deltaT(k)) ' err=' num2str(fiterr(k))]);
        
    end
    
end




return






function costf = f_TS(delT, ap, psiT, cp, NIR, ref);

    bp = cp - ap;

    costf = sum(abs( ap(NIR) - psiT(NIR).*delT - ((ap(ref)-psiT(ref).*delT)./bp(ref)).*bp(NIR)))   ;

return

function costf = f_TSc(delT, cp, psiT, NIR);

    costf = sum(abs(cp(NIR) - psiT(NIR).*delT));

return
    