function [a_sw,bb_sw]=water_iops_PF_TScorr(wave,T,S)

%function to obtain seawater absorption and backscattering spectra
%pure water absorption from Pope and Frye
%pure water backscattering from Morel
%corrected for in situ temperature and salinity conditions

if min(wave)<380ormax(wave)>730
    'input wavelengths  are beyond acceptable range'
else
%Pope and Frye pure water absorption 380-730 nm, then Smith and Baker 730-800 nm
    wl_water=[380	382.5	385	387.5	390	392.5	395	397.5	400	402.5	405	407.5	410	412.5	415	417.5	420	422.5	425	427.5	430	432.5	435	437.5	440	442.5	445	447.5	450	452.5	455	457.5	460	462.5	465	467.5	470	472.5	475	477.5	480	482.5	485	487.5	490	492.5	495	497.5	500	502.5	505	507.5	510	512.5	515	517.5	520	522.5	525	527.5	530	532.5	535	537.5	540	542.5	545	547.5	550	552.5	555	557.5	560	562.5	565	567.5	570	572.5	575	577.5	580	582.5	585	587.5	590	592.5	595	597.5	600	602.5	605	607.5	610	612.5	615	617.5	620	622.5	625	627.5	630	632.5	635	637.5	640	642.5	645	647.5	650	652.5	655	657.5	660	662.5	665	667.5	670	672.5	675	677.5	680	682.5	685	687.5	690	692.5	695	697.5	700	702.5	705	707.5	710	712.5	715	717.5	720	722.5	725	727.5	730  732.5  735.0  737.5  740.0  742.5  745.0  747.5  750.0  752.5  755.0  757.5  760.0  762.5  765.0  767.5  770.0  772.5  775.0  777.5 780.0  782.5  785.0  787.5  790.0  792.5  795.0  797.5  800.0];
    a_water=[0.01137 0.010044 0.00941 0.00917 0.00851 0.00829 0.00813 0.00775 0.00663 0.00579 0.0053 0.00503 0.00473 0.00452 0.00444 0.00442 0.00454 0.00474 0.00478 0.00482 0.00495 0.00504 0.0053	0.0058 0.00635	0.00696	0.00751	0.0083 0.00922	0.00969 0.00962 0.00957 0.00979 0.01005 0.01011 0.0102 0.0106 0.0109 0.0114 0.0121 0.0127 0.0131 0.0136 0.0144 0.015 0.0162 0.0173 0.0191 0.0204 0.0228 0.0256 0.028 0.0325 0.0372 0.0396 0.0399 0.0409 0.0416 0.0417 0.0428 0.0434 0.0447 0.0452 0.0466 0.0474 0.0489 0.0511 0.0537 0.0565 0.0593 0.0596 0.0606 0.0619 0.064 0.0642 0.0672 0.0695 0.0733 0.0772 0.0836 0.0896 0.0989 0.11 0.122 0.1351 0.1516 0.1672 0.1925 0.2224 0.247 0.2577 0.2629 0.2644 0.2665 0.2678 0.2707 0.2755 0.281 0.2834 0.2904 0.2916 0.2995 0.3012 0.3077 0.3108 0.322 0.325 0.335 0.34 0.358 0.371 0.393 0.41 0.424 0.429 0.436 0.439 0.448 0.448 0.461 0.465 0.478 0.486 0.502 0.516 0.538 0.559 0.592 0.624	0.663 0.704 0.756 0.827 0.914 1.007 1.119 1.231 1.356 1.489 1.678 1.7845 1.9333 2.0822 2.2311 2.3800 2.4025 2.4250 2.4475 2.4700 2.4900 2.5100 2.5300 2.5500 2.5400 2.5300 2.5200 2.5100 2.4725 2.4350 2.3975 2.3600 2.3100 2.2600 2.2100 2.1600 2.1375 2.1150 2.0925 2.0700];
    a_pw=interp1(wl_water,a_water,wave,'linear');

%Morel water backscattering
    %wl_b=[380	390	400	410	420	430	440	450	460	470	480	490	500	510	520	530	540	550	560	570	580	590	600	610	620	630	640	650	660	670	680	690	700 750];
    %b_water=[0.0073	0.0066	0.0058	0.0052	0.0047	0.0043	0.0039	0.0035	0.0032	0.0029	0.0027	0.0024	0.0022	0.002	0.0018	0.0017	0.0016	0.0015	0.0013	0.0013	0.0012	0.0011	0.00101	0.00094	0.00088	0.00082	0.00076	0.00071	0.00067	0.00063	0.00059	0.00055	0.00052 0.0005];
    %%choose pure water scattering function (divide by two for back-scattering):
    %bb_pw=0.5*interp1(wl_b,b_water,wl,'linear');

%log fit water backscattering
    bb_logfit=.0037000*(380^4.3)./(wave.^4.3);

%salinity correct:
if S>0,bb_sw=((1+0.01*S)*bb_logfit);,else bb_sw=bb_logfit;,end

%temp and salinity correction for water absorption (need to know at what T it was measured): 
if S==0,S=35;,end
if T==0,T=22;,end
T_pope=22.0;

%parameters for temp and salinity callibration (From Pegau et al Applied optics 1997):
M=[.18 .17 .52 1.4 4.6 2.1 4.3 9.6 1.6 34 18 42];
sig=[18 15 14 20 17.5 15 17 22 6 18 20 25];
lamda_c=[453 485 517 558 610 638 661 697 740 744 775 795];
M_T=[0.0045 0.002 0.0045 0.002 0.0045 -0.004 0.002 -0.001 0.0045 0.0062 -0.001 -0.001];

%computing the correction per degree C
for i=1:length(wave)
phi_T(i)=sum(M_T.*M./sig.*exp(-(wave(i)-lamda_c).^2/2./sig.^2));
end

%salinity correction based on Pegau and Zaneveld 1997:
wls=[412 440 488 510 532 555 650 676 715 750];
phi_S_PZ=[0.00012 -0.00002 -0.00002 -0.00002 -0.00003 -0.00003 0 -0.00002 -0.00027 0.00064];

%interpolate to compute salinity correction per psu
phi_S=interp1(wls,phi_S_PZ,wave,'spline');

%temperature and salinity corrections:
a_sw=(a_pw+phi_T'*(T-T_pope)+phi_S*S);
end