find_sig_interface.m
by
Webmaster Legos
—
last modified
Jul 05, 2012 12:12 PM
find_sig_interface.m
—
Objective-C source code,
2 kB (2871 bytes)
File contents
function [psigint]=find_sig_interface(sigint,sigipref,... pref,temp,sali,Pres,maxdp,Pdep,gsecs) % KEY: find pressure of the given sigma interfaces % USAGE : % % % % % DESCRIPTION : % % % INPUT: % sigint(nlayer) sigma interface value % sigipref(nlayer) reference indices pressure for each interface % sigpref(npref) reference pressures % % % OUTPUT: % psigint(npair,nlayer) % % AUTHOR : A.Ganachaud (ganacho@gulf.mit.edu) , Jan 97 % % SIDE EFFECTS : % % SEE ALSO : % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CALLER: % CALLEE: getsigpres %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %FOR EACH REFERENCE PRESSURE disp('looking for sigma interface depths ...') npair=size(temp,2); nlayint=length(sigint); psigint=NaN*ones(npair,nlayint); for ipref=1:max(sigipref) disp(sprintf('reference pressure = %i db',pref(ipref))) %get the sigma indices for this reference pressure gipref=find(sigipref==ipref); %get sigmas relative to that pressure reference over the whole section sig_curpref = -1000+sw_pden(sali,temp,Pres,pref(ipref)); if nargin>=9 if any(strcmp(fieldnames(gsecs),'sp')); disp('SPLINE FILTERING (CSAPS)...') id1=min(find(Pres>gsecs.sp.pint(1))); int1=min(id1,maxdp); sig_curpref1=sig_curpref; for is=1:size(sig_curpref,2) sig_curpref1(1:int1(is),is)=csaps(Pres(1:int1(is)),... sig_curpref(1:int1(is),is),gsecs.sp.scoef(1),Pres(1:int1(is))); end %SECOND SPLINE id2=min(find(Pres>gsecs.sp.pint(2))); int2=min(id2,maxdp); gis=find(int2>int1+1); for i=1:length(gis) is=gis(i); sig_curpref1(int1(is)+1:int2(is),is)=csaps(Pres(int1(is)+1:int2(is)),... sig_curpref(int1(is)+1:int2(is),is),gsecs.sp.scoef(2),... Pres(int1(is)+1:int2(is))); end %THIRD SPLINE int3=maxdp; gis=find(int3>int2+1); for i=1:length(gis) is=gis(i); sig_curpref1(int2(is)+1:int3(is),is)=csaps(Pres(int2(is)+1:int3(is)),... sig_curpref(int2(is)+1:int3(is),is),gsecs.sp.scoef(3),... Pres(int2(is)+1:int3(is))); end %REMOVE SPIKES AT JUNCTIONS for is=1:size(sig_curpref1,2) imx=maxdp(is); sig_curpref1(1:imx,is)=csaps(Pres(1:imx),sig_curpref1(1:imx,is),... gsecs.sp.junccoef,Pres(1:imx)); end if 0 is=5;plot(sig_curpref(:,is),-Pres,'+',sig_curpref1(:,is),-Pres,'-') clf;vc=[26:.5:28.5,28.55:0.05:29.05,29.06:0.01:31]; Nd=900;Np=10; [c,h]=contour(1:Np,-Pres(1:Nd),sig_curpref1(1:Nd,:),vc,'k'); clabel(c,h);hold on;plot(1:Np,-Pres(maxdp)) end end end %interpolates the pressure to guess the sigma position psigint(:,gipref) = ... getsigpres(Pres,Pdep,maxdp,sig_curpref,sigint(gipref)); end %ipref=1:max(sigipref) disp('interface depths found')