function [z,w,dz,ddz,Q] = pd(geom,ndisc,axlim)% PD.M Computes quadrature weights, discretized coordinates, Lagrange-%      interpolation polynomial-based differentiation arrays, and%      discretzed Jacobi polynomials for a 1-dimensional domain. Call%      as%%      [z,w,dz,ddz,Q] = pd(geom,ndisc,axlim);%%      INPUT Variables: %          geom  : geometry - 'slab', 'cyln', 'sphe', 'peri'%          ndisc : number of discretization points, including endpoints%          axlim : [lower upper] bounds on discretized interval%%      OUTPUT Variables:%          z     : (ndisc by 1) set of discretization points in%                  the specified interval;%          w     : (ndisc by 1) quadrature weight array; note that%                  the quadrature array contains the factor z.^a;%          dz    : (ndisc by ndisc) first order derivative array;%          ddz   : (ndisc by ndisc) Laplacian operator array;%          Q     : discretized Jacobi polynomials in column format.%%      EXAMPLE: Generate the location array, quadrature weights,%               differentiation arrays, and Jacobi polynomials%               for a cylindrical physical domain. Plot the first%               5 polynomials and demonstrate their orthogonality.%%               [r,w,dr,ddr,Q] = pd('disk',20);%               plot(r,Q(:,1:4))%               wip( Q(:,1:5),Q(:,1:5),w.*r.*(1-r) )% Written by Raymond A. Adomaitis, Hsiao-Yung Chang, and Yi-hung Lin % as part of MWRtools, version 2.% Copyright (c) by Raymond A. Adomaitis, 1998-2002if nargin < 3, axlim = [0 1]; endswitch lower(geom)   case 'slab', a = 0;   case 'cyln', a = 1;   case 'sphe', a = 2;   case 'peri', % not relevant   otherwise, disp('unknown geometry or not yet supported')   endif strcmp(lower(geom),'peri')% For periodic domains, axlim are set to [0, 2*pi - delT]   if ~mod(ndisc,2)      error('Odd number of points required w/ periodic coords')      end   z = 2*pi*[0:ndisc-1]'/ndisc;   w = ones(size(z))/ndisc;   Q   = dFs(z,ndisc);   dQ  = zeros(size(z));   ddQ = zeros(size(z));   for i = 1:(ndisc-1)/2      dQ(:,2*i)    = -i*sin(i*z);      dQ(:,2*i+1)  =  i*cos(i*z);      ddQ(:,2*i)   = -i^2*cos(i*z);      ddQ(:,2*i+1) = -i^2*sin(i*z);      end   dz  = dQ*inv(Q);   ddz = ddQ*inv(Q);   return   end% Non-periodic domainz        = qpts(a,ndisc); w        = qweights(z,a,ndisc); [dz ddz] = lDM(z,ndisc);Q        = dJs(z,a+1,ndisc-1,2);% Adjust for axis limitszBar = z;z0  = axlim(1);z1  = axlim(2);z   = z0 + (z1-z0)*zBar;dz  = dz/(z1-z0);ddz = ddz/(z1-z0)^2;w   = (z1-z0)*w;% Laplacian operator for a > 0 requires asymptotic value when z0 = 0if a ~= 0 & z0 > 0   ddz = a*dz./z(:,ones(ndisc,1)) + ddz;   endif a ~= 0 & z0 == 0   v = [2:ndisc];   ddz(v,:) =  a*dz(v,:)./z(v,ones(ndisc,1)) + ddz(v,:);   ddz(1,:) = (a+1)*ddz(1,:);   end% The quadrature weights, scaled for non-unit intervalsswitch acase 0   w = w;case 1   w0 = qweights(zBar,0,ndisc);    w = z0*w0 + w;   w = 2*pi*(z1-z0)*w;case 2   w0 = qweights(zBar,0,ndisc);    w1 = qweights(zBar,1,ndisc);    w = z0^2*w0 + 2*z0*(z1-z0)*w1 + (z1-z0)*w;   w = 4*pi*(z1-z0)*w;   end% ------------- quadrature points ------------------ %function [x] = qpts(a,N);   x  = [0; rJ(a+1,N-2); 1];   % ------------- quadrature weights ------------------ %function [w] = qweights(x,a,N);   dx = x(:,ones(1,N))' - x(:,ones(1,N)); % difference between all pointsp  = [ones(1,N); cumprod(dx)];dp = zeros(1,N);for i = 1 : N      dp = dx(i,:).*dp + p(i,:);   ende   = [1 meshgrid(a+1,ones(N-1,1))'];K   = [a+1 meshgrid((a+1)^2,ones(N-1,1))'; meshgrid(e,ones(N-1,1))];dp2 = (dp'*(1./dp)).^2;w = 1./sum(dp2.*K,2);% ------------- differentiation arrays ------------------ %function [A,B] = lDM(x,N);p  = [ones(N,1) zeros(N,3)];xt = x(:,ones(2));for j = 1 : N   p = (xt-x(j)).*p + [zeros(N,1) p(:,1) 2*p(:,2) 3*p(:,3)];   enddx = x(:,ones(1,N)) - x(:,ones(1,N))' + eye(N)*2;[p1 p2] = meshgrid(p(:,2));A = (triu(p2,1)+tril(p2,-1)+diag(p(:,3)))./p1./dx;B = 2*A.*(ndgrid(diag(A))-1./dx);B = triu(B,1) + tril(B,-1) + diag(p(:,4)./p(:,2)/3);% ------------- roots of a Jacobi polynomial ------------------ %function [r] = rJ(b,N);                      x = (cos(pi*[N*3-1:-1:0]'/(N*3-1))+1)/2;      r = rdf( x,dJs(x,b,N,1) );                     update = r;j = [1:length(r)]';while ~isempty(j)   update(j) = -dJs(r(j),b,N);   r(j)      = r(j) + update(j);   j         = find(abs(update)>sqrt(eps));   end