PCdisp_源程序_管中导波频散曲线_PCdisp_导波频散曲线_频散曲线

% PCKFCURVES Traces the dispersion curves of a waveguide in the K-F space % % A lot of work has been put into optimizing this routine because tracing % the curves takes a lot of time, but remember that "you're bound to be % unhappy if you try to optimize everything" (Donald Knuth) % % PCDISP: Pochhammer-Chree dispersion in cylindrical waveguides % Written by Fernando Seco Granja (fernando.seco@csic.es) % Consejo Superior de Investigaciones Cient�ficas (CSIC) % http://www.car.upm-csic.es/lopsi/people/fernando.seco/pcdisp/ % % Last update: 3-10-2011 clear,clc % ********** PARTS OF THE SPECTRUM TO BE PLOTTED ******** flag_realroots=1; flag_imagroots=1; flag_complexroots=1; % ********** WAVEGUIDE PARAMETERS ********** pcwaveguide % Circumferential order to be plotted (n = -1 for torsional modes T(0,m), % n = 0 for longitudinal L(0,m), and n >= 1 for flexural, F(n,m)) n = 0; % *********** SOLVER PARAMETERS ********* % Solutions will be computed in the 3D space with 0<Re{k},Im{k}<=kmax, % and 0<f<=fmax fmax=1000e3; % Maximum frequency (Hz) kmax=2000; % Maximum wavenumber (1/m) % Tolerance error in the determination of the roots, in normalized units drtol=1e-6; % Maximum distance between consecutive points along a branch drmax=.05; drmin=5*drtol; dr_factor=1.33; % should be strictly bigger than 1 % Other solver parameters dthmax=pi/4; Nroots=20; % Number of points between consecutive roots tested for possible sign % changes in the bisection root finder. It should be increased if roots are % missing because they're very close npcheck=150; % ********* HERE WE GO! ********** % Normalized coordinates for 2D plot (x -> k, y -> f) xmax=kmax/sc_k; ymax=fmax/sc_freq; clc disp('PCKFCURVES running') disp(' ') figure(1),clf,hold on axis(1.01*[-kmax kmax 0 fmax/1e3]) xlabel('-Im\{k\} (1/m) Re\{k\}') ylabel('f (kHz)') handle_branch=[]; indbranch=0; xbranch=[]; ybranch=[]; npbranch=[]; flagbranch=[]; % ****** TRACING OF BRANCHES WITH REAL AND IMAGINARY WAVENUMBERS ****** tvect=[]; if(flag_realroots), tvect=[1]; end if(flag_imagroots), tvect=[tvect j]; end for t=tvect, if(t==1), color='b';xsign=1;name='real'; else, color='r';xsign=-1;name='imaginary'; end % First we collect the roots of the frequency equation along the edges % of the wavenumber and frequency boxes. These are the starting and % ending points for the branches, and are called "docking" points here disp(sprintf('Determining docking points for branches with %s wavenumber...',name)) dr=drmax/3; np=3; % number of pre-computed points at the beginning/end of each branch numleft=[];numright=[];numtop=[];numbottom=[];numcircle=[]; for ind=1:np, % Left edge of the k-f space (k=0) dummy=pcsolvebisection('k-f',t*dr*ind,[4/3*drmax ymax],drtol,-1,npcheck,n,waveguidedescription); numleft(ind)=length(dummy); yleft(ind,1:numleft(ind))=dummy; xleft(ind,1:numleft(ind))=dr*ind; % Right edge of the k-f space (k=kmax) dummy=pcsolvebisection('k-f',t*(xmax-dr*(ind-1)),[dr ymax],drtol,-1,npcheck,n,waveguidedescription); numright(ind)=length(dummy); yright(ind,1:numright(ind))=dummy; xright(ind,1:numright(ind))=xmax-dr*(ind-1); % Top edge of the k-f space (f=fmax) dummy=pcsolvebisection('f-k',ymax-dr*(ind-1),t*[0 xmax],drtol,-1,npcheck,n,waveguidedescription); numtop(ind)=length(dummy); xtop(ind,1:numtop(ind))=dummy; ytop(ind,1:numtop(ind))=ymax-dr*(ind-1); % Bottom edge of the k-f space (f=0) dummy=pcsolvebisection('f-k',dr*ind,t*[4/3*drmax xmax],drtol,-1,npcheck,n,waveguidedescription); numbottom(ind)=length(dummy); xbottom(ind,1:numbottom(ind))=dummy; ybottom(ind,1:numbottom(ind))=dr*ind; % Quarter circle around the origin (k=0,f=0) dummy=pcsolvebisection('0-0',t*dr*ind,[0 pi/2],drtol/(ind*dr),-1,npcheck,n,waveguidedescription); numcircle(ind)=length(dummy); xcircle(ind,1:numcircle(ind))=dr*ind*cos(dummy); ycircle(ind,1:numcircle(ind))=dr*ind*sin(dummy); end flag_plot_dockpoints=0; if(flag_plot_dockpoints), plot(xleft*sc_k*xsign,yleft*sc_freq/1e3,'Marker','o','Color',color) plot(xright*sc_k*xsign,yright*sc_freq/1e3,'Marker','x','Color',color) plot(xtop*sc_k*xsign,ytop*sc_freq/1e3,'Marker','+','Color',color) plot(xbottom*sc_k*xsign,ybottom*sc_freq/1e3,'Marker','*','Color',color) plot(xcircle*sc_k*xsign,ycircle*sc_freq/1e3,'Marker','s','Color',color) end % Checks that the starting points have been correctly computed if(any(numleft(2:np)~=numleft(1)) | any(numright(2:np)~=numright(1)) | ... any(numtop(2:np)~=numtop(1)) | any(numbottom(2:np)~=numbottom(1)) | ... any(numcircle(2:np)~=numcircle(1))), numleft,numright,numtop,numbottom,numcircle disp('Problems determining the starting points') stop else, numleft=numleft(1);numright=numright(1);numtop=numtop(1); numbottom=numbottom(1);numcircle=numcircle(1); end % Ocassionally the branch starting at the origin might have been % computed twice - we have to remove one instance if(numcircle==1), slope1=sum(ycircle(1:np,1))/sum(xcircle(1:np,1)); if(slope1<1/3), if(numbottom>0), xbottom=xbottom(1:np,2:numbottom); ybottom=ybottom(1:np,2:numbottom); numbottom=numbottom-1; end end if(slope1>3), if(numleft>0), xleft=xleft(1:np,2:numleft); yleft=yleft(1:np,2:numleft); numleft=numleft-1; end end end % Simple check if(t==1), if(numbottom>0 | numleft+numcircle~=numright+numtop), disp('There was a problem computing the docking points') numleft,numbottom,numright,numtop,numcircle stop end end % Arrays with the docking points xdock=[xleft(1:np,1:numleft) xright(1:np,1:numright) xtop(1:np,1:numtop) xbottom(1:np,1:numbottom) xcircle(1:np,1:numcircle)]'; ydock=[yleft(1:np,1:numleft) yright(1:np,1:numright) ytop(1:np,1:numtop) ybottom(1:np,1:numbottom) ycircle(1:np,1:numcircle)]'; Ndock=size(xdock,1); % Sort them [dummy,indord]=sort(xdock(:,1).^2+ydock(:,1).^2); xdock=xdock(indord,:);ydock=ydock(indord,:); % We plot them and get a handle for them in the figure figure(1), handle_dock=plot(xdock(:,1)*xsign*sc_k,ydock(:,1)*sc_freq/1e3,... 'Marker','o','Color',color,'LineStyle','none'); % Total number of branches Nbranches=Ndock/2; if(floor(Nbranches)~=Nbranches), disp('The number of docking points should be even!') stop else, disp(sprintf('Tracing %g branches with %s wavenumber',Nbranches,name)) end % Loop over all branches for indbranchtmp=1:Nbranches, numeval=0; xp=xdock(1,1:3); yp=ydock(1,1:3); % We eliminate the starting points from the arrays xdock=xdock(2:end,:);ydock=ydock(2:end,:); Ndock=Ndock-1; indbranch=indbranch+1; flagbranch(indbranch)=t; figure(1),handle_branch(indbranch)=plot(xp*xsign*sc_k,yp*sc_freq/1e3,'Marker','.','Color',color); flag_branch_complete=0; indp=3; drlast=[norm([xp(3)-xp(2) yp(3)-yp(2)])]; drlista=[drlast*dr_factor drlast drlast/dr_factor]; % When tracing the branch we keep track of the sign of the % determinant at the left side of the curve