
function [ap_TSalScatCorr, deltaT] = 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

%     tmp = xlsread('Sullivan_etal_2006_instrumentspecific.xls');

% tmp = [0.0001;0.0001;0.0001;0.0001;0;0;0;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0002;0.0003;0.0003;0.0004;0.0005;0.0006;0.0006;0.0007;0.0008;0.0009;0.001;0.001;0.001;0.001;0.001;0.0009;0.0009;0.0008;0.0007;0.0006;0.0006;0.0005;0.0004;0.0003;0.0003;0.0002;0.0001;0.0001;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0002;0.0002;0.0002;0.0001;0.0001;0.0001;0;0;-0.0001;-0.0001;-0.0001;-0.0001;-0.0001;-0.0001;0;0;0.0001;0.0002;0.0003;0.0005;0.0007;0.0009;0.0013;0.0017;0.0021;0.0026;0.0032;0.0038;0.0045;0.0054;0.0063;0.0073;0.0083;0.0094;0.0104;0.0113;0.0121;0.0128;0.0133;0.0136;0.0136;0.0133;0.0129;0.0124;0.0116;0.0107;];
% 
%     psiT = interp1(tmp(:,1),tmp(:,2), wavel);
psiT = [0.000100000000000000,5.99999999999994e-05,0,5.00000000000000e-05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.49999999999989e-05,0.000100000000000000,0.000100000000000000,0.000100000000000000,5.00000000000114e-06,0,0,0,0,0,0,0,0,4.49999999999989e-05,0.000100000000000000,0.000300000000000000,0.000419999999999999,0.000600000000000000,0.000785000000000002,0.00100000000000000,0.00100000000000000,0.000960000000000002,0.000805000000000001,0.000600000000000000,0.000430000000000001,0.000280000000000001,0.000100000000000000,0,0,0,0,0.000100000000000000,0.000200000000000000,0.000139999999999998,9.49999999999989e-05,-1.99999999999989e-05,-0.000100000000000000,-0.000100000000000000,-4.49999999999989e-05,7.50000000000000e-05,0.000250000000000000,0.000610000000000002,0.00110000000000000,0.00184000000000001,0.00275000000000000,0.00380000000000000,0.00508499999999998,0.00674999999999999,0.00835500000000001,0.0103000000000000,0.0116200000000000,0.0128000000000000,0.0134350000000000,0.0136000000000000,0.0132200000000000,0.0125000000000000,0.0115550000000000];
 
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>=730, 1,'first');

    ap_TSalScatCorr = nan(size(ap_uncorr));
    deltaT = nan(size(ap_uncorr,1),1);
    fiterr = nan(size(ap_uncorr,1),1);
   
    for k = 1:nspect

        if all(isfinite(ap_uncorr(k,:)))
            % guess for paramters (beamc at lambda0, beamc slope)
            delT = 0;

            % minimization routine
            [deltaT(k), fiterr(k)] = fminsearch(@f_TS, 0, opts, ap_uncorr(k,:), cp(k,:), psiT, NIR, ref);

            bp = cp(k,:) - ap_uncorr(k,:);
            ap_TSalScatCorr(k,:) = ap_uncorr(k,:) - psiT.*deltaT(k) - ...
                    ((ap_uncorr(k,ref)-psiT(ref).*deltaT(k))./bp(ref)).*bp;
            
            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, cp, psiT, NIR, ref);

    bp = cp - ap;

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

return
    