1.zip_worldlyf_在matlab中计算水声简正波程序的代码_水声_简正_简正波 matlab

function [cp, cg, evf] = amodes (Z, c, frq, nm, zr) % % function [cp, cg, evf] = amodes (Z, c, frq, nm, zr) % % Calculation of Acoustic Modes % % Calculates acoustic modes modes, P, from % Pzz+ [w^2/c^2 - ev^2]*P=0, subject to P(-D)=0 and P(0)=0 % (I am at the moment a little unsure of the boundary conditions!) % (And how to implement them in this matrix formalism.) % % Uses Numerov's method following the suggestion of Jensen, Kuperman, % Porter, and Schmitt, Computational Ocean Acoustics, p. 282-286. % Googlebooks for a preview: % http://books.google.com/books?id=QHtx4zYPbzMC % % The approach here first calculates the components of the matrices required % for Numerov formulation of the differential equation. The eigenvalue % （特征值，本征值） % equation that results is: A*phi=(k^2)*B*phi, which is a standard eigenvalue % equation. The eigenvectors（特征向量） phi are found with the call to eigs(A,B,...). % % Roughly, B is tridiagonal（三对角） with 10/12 on the main diagonal（对角线） and 1/12 on % the +/-1 diagonals. % % INPUT: % % Z,c are input depth and sound speed. These must already be interpolated % to an evenly-spaced grid to achieve a desired accuracy. 5-m spacing seems % to work o.k.; this depends on frequency, of course. % % frq and nm are frequency and number of modes desired. % % zr is a vector of depths at which the mode shapes are desired. % This might correspond to the depths of hydrophones in a simulated % reception. amodes interpolates the modes it finds to the zr depths. % % Some timing on a 2.4 GHz Q6600 cpu: % At 100 Hz, using % 1000 sound speeds, 50 modes in 4.3 s; d(cp)=6e-4, d(cg)=1e-3 % 2500 sound speeds, 50 modes in 38.2 s; d(cp)=2e-5, d(cg)=1e-3 % 5000 sound speeds, 50 modes in 259.0 s; d(cp)=2e-5, d(cg)=1e-3 % d(cp) is difference in phase speed relative to chebmode, % Same for d(cg) and group speed. These taken to be some measure % of the accuracy of the calculations. % % OUTPUT: % % cp and cg are mode phase and group speed. % evf is a matrix of nm modes at the depths zr. % % Uses eigs to calculate only nm modes. The use of eig is possible, but % this calculates length(Z) modes, most of which were tossed out（丢弃，淘汰! % % B. Dushaw February-April 2009 % allocate for evf evf=zeros(length(zr),nm,length(frq)); % Grid spacing h=abs(Z(2)-Z(1)); % Frequency to radians/s w=2*pi*frq; % Calculate slowness. sp=1./c; % Main loop over frequency for kkk=1:length(frq), sp2=(w(kkk)^2)*sp.*sp; % frequency over sound speed squared % Diagonal elements dj=(-2/(h^2) + (10/12)*sp2)/(h*1024); % Upper off diagonal eja=(1/(h^2) + (1/12)*sp2(2:end))/(h*1024); % Lower off diagonal ejb=(1/(h^2) + (1/12)*sp2(1:(end-1)))/(h*1024); % Construct the matrix A - A is not symmetric. A1=diag(dj); A2=diag(eja,1); A3=diag(ejb,-1); A=A1+A2+A3; A(1,:)=[]; % Unsure about the boundary conditions A(:,1)=[]; % In their simpler example Jensen et al say to null out the first row and column. % Seems to work better with this as well - better accuracy, and modes are zero % at the surface and bottom like they are supposed to be. A(end,:)=[]; A(:,end)=[]; clear A1 A2 A3 % Construct the matrix B - B is symmetric. B1=(10/12)*eye(length(A)); B2=(1/12)*diag(ones(1,length(A)-1),1); B3=(1/12)*diag(ones(1,length(A)-1),-1); B=(B1+B2+B3)./(h*1024); clear B1 B2 B3 if 1==1, % Use eigs to just get the nm largest real eigenvalues % Cuts computation time by a factor of 3-4 over eig. % When only a small number of modes are requested, one needs to have % more Lanczos vectors for convergence. P=max(2*nm,75); P=min(P,length(A)); % P can be at most length(A) % tol is by default "eps" which is 2.2e-16 on my computer. Increasing this % tolerance makes things run much faster, without noticeably affecting the % accuracy of the modes or eigenvalues that I can tell. tol=1e-10; %tol=eps % Options for eigs: % 'tol'= numerical tolerance, % 'maxit'=maximum number of iterations, % 'p'=number of Lanczos vectors, % 'disp'=verbosity of calculation. % (there are other options) opts=struct('tol',tol,'maxit',500,'p',P,'disp',0); % All of the compute time is in this call. [V,kr,flag]=eigs(sparse(A),sparse(B),nm,'LR',opts); % Convergence=flag % If convergence is not zero, you have problems. kr=diag(kr)'; % eigenvalues - KH in the literature. % Arrange evalues and evectors in descending order % These may be in order already. [kr,I]=sort(kr,'descend'); V=V(:,I); else % Use eig to get all eigenvalues. % Solves for ALL modes; very computationally expensive. [V,kr]=eig(A,B); kr=diag(kr)'; % eigenvalues - KH in the literature. [kr,I]=sort(diag(kr),'descend'); V=V(:,I); kr=kr(1:nm); V=V(:,1:nm); end % kr are the wavenumbers kr=sqrt(kr); [n m]=size(V); V=[zeros(1,m);V;zeros(1,m)]; % mode is zero at the surface and bottom. % I don't understand how this works, but it % seems to. if 1==2, % calculate local vertical wavenumber for kk=1:nm, KV=sqrt(w(kkk)*w(kkk)*sp.*sp - kr(kk)*kr(kk)); LV=2*pi./KV; plot(LV/5,Z) title(['Mode Number: ' num2str(kk)]) grid on axis([0 35 -5000 0]) pause end end if 1==2, % Now the test...how well is the eigenvalue problem solved? V1=V; V1(1,:)=[]; for kk=1:4, Q=A*V1-kr(kk)*kr(kk)*B*V1; plot(Q(:,kk),Z(2:end)) hold on end pause % In my tests, Q came out to be at most 2e-16 end % Interpolate modes to 5-m spacing for better integration accuracy. s=sign(Z(end)); hs=5*s; % 5-m intervals seems to work fine; 10-m works ok as well, looks like. ZP=Z(1):hs:Z(end); hs=abs(hs); % Normalize the modes F=interp1(Z,V,ZP,'spline'); % always spline, never cubic! F=F.*F; N=hs*trapz(F); N=sqrt(N); N=ones(size(Z))*N; V=V./N; % Set the phase of modes at the top if 1==2, V=V.*(ones(length(V),1)*sign(V(3,:))); else for kook=1:nm, I=find(abs(V(:,kook))>1e-3); I=I(1); V(:,kook)=V(:,kook)*sign(V(I,kook)); end end % Phase speed cp(:,kkk)=w(kkk)./kr; % Group speed per Jensen et al, or Munk et al 1995 .... % V's are normalized just above, so we're good to go here. Q=sp*ones(1,nm); Q=V.*Q; Q=interp1(Z,Q,ZP,'spline'); % always spline, never cubic! Q=Q.*Q; icg=hs*(w(kkk)./kr).*trapz(Q); cg(:,kkk)=1./icg; % Interpolate the modes to the requested instrument depths. evf(:,:,kkk)=interp1(Z,V,zr,'spline'); end % end loop for frequencies.