Code covered by the BSD License

# f2matlab

### Ben Barrowes (view profile)

16 Jun 2004 (Updated )

Converts basic fortran90 source code to MATLAB source code.

mcpbdn(varargin)
```function mcpbdn(varargin)

%From the book 'Computation of Special Functions.'
% by Shanjie Zhang, and Jianming Jin, New York : Wiley, c1996.
%       =============================================================
%       Purpose: This program computes parabolic cylinder functions
%                Dn(z) for an integer order and a complex argument
%                using subroutine CPBDN
%       Input :  x --- Real part of z
%                y --- Imaginary part of z
%                n --- Order of Dn(z)
%       Output:  CPB(|n|) --- Dn(z)
%                CPD(|n|) --- Dn'(z)
%       Example:
%                z = 5.0+ 5.0 i
%     n     Re[Dn(z)]      Im[Dn(z)]      Re[Dn'(z)]     Im[Dn'(z)]
%   -----------------------------------------------------------------
%     0   .99779828D+00  .66321897D-01 -.23286910D+01 -.26603004D+01
%     1   .46573819D+01  .53206009D+01  .26558457D+01 -.24878635D+02
%     2  -.43138931D+01  .49823592D+02  .14465848D+03 -.10313305D+03
%     3  -.28000219D+03  .21690729D+03  .12293320D+04  .30720802D+03
%     4  -.24716057D+04 -.46494526D+03  .38966424D+04  .82090067D+04
%    -1   .10813809D+00 -.90921592D-01 -.50014908D+00 -.23280660D-01
%    -2   .24998820D-02 -.19760577D-01 -.52486940D-01  .47769856D-01
%    -3  -.15821033D-02 -.23090595D-02 -.68249161D-03  .10032670D-01
%    -4  -.37829961D-03 -.10158757D-03  .89032322D-03  .11093416D-02
%       =============================================================
clear global; clear functions;

persistent cpb cpd ;
n=[];z=[];
if isempty(cpb), cpb=zeros(1,100+1); end;
if isempty(cpd), cpd=zeros(1,100+1); end;
writef(1,['%s \n'],'Please enter n, x and y ');
% WRITE(*,20)N,X,Y
n=5;
x=5;
y=5;
z=complex(x,y);
n0=abs(n);
[n,z,cpb,cpd]=cpbdn(n,z,cpb,cpd);
writef(1,['%0.15g \n']);
if(n >= 0)
writef(1,['%s %s \n'],'  n     Re[Dn(z)]       Im[Dn(z)]       ','Re[Dn''(Z)]      IM[DN''(Z)]');
else;
writef(1,['%s %s \n'],' -n     Re[Dn(z)]       Im[Dn(z)]       ','Re[Dn''(Z)]      IM[DN''(Z)]');
writef(1,['%s %s \n'], ' -n     Re[Dn(z)]       Im[Dn(z)]       ','Re[Dn''(Z)]      IM[DN''(Z)]');
end;
writef(1,['%s %s \n'],'-------------------------------------------','-------------------------');
for  i=0:n0;
writef(1,[repmat(' ',1,1),'%3i',repmat('%16.8f',1,4) ' \n'],i,cpb(i+1+1),cpd(i+1+1));
end;  i=n0+1;
%format(1x,i3,4d16.8);
%format(1X,'N =',i3,',   z =x+iy :',f6.2,'+',f6.2,' i');
end %program mcpbdn
function [n,z,cpb,cpd]=cpbdn(n,z,cpb,cpd);
%       ==================================================
%       Purpose: Compute the parabolic cylinder functions
%                 Dn(z) and Dn'(z) for a complex argument
%       Input:   z --- Complex argument of Dn(z)
%                n --- Order of Dn(z)  ( n=0,�,�,��)
%       Output:  CPB(|n|) --- Dn(z)
%                CPD(|n|) --- Dn'(z)
%       Routines called:
%            (1) CPDSA for computing Dn(z) for a small |z|
%            (2) CPDLA for computing Dn(z) for a large |z|
%       ==================================================
z1=[];cf1=[];n0=[];cfa=[];n1=[];cfb=[];
pi=3.141592653589793d0;
x=real(z);
a0=abs(z);
c0=complex(0.0d0,0.0d0);
ca0=exp(-0.25d0.*z.*z);
if(n >= 0)
cf0=ca0;
cf1=z.*ca0;
cpb(0+1+1)=cf0;
cpb(1+1+1)=cf1;
for  k=2:n;
cf=z.*cf1-(k-1.0d0).*cf0;
cpb(k+1+1)=cf;
cf0=cf1;
cf1=cf;
end;  k=n+1;
else;
n0=-n;
if(x <= 0.0||abs(z) == 0.0)
cf0=ca0;
cpb(0+1+1)=cf0;
z1=-z;
if(a0 <= 7.0)
[dumvar1,z1,cf1]=cpdsa(-1,z1,cf1);
else;
[dumvar1,z1,cf1]=cpdla(-1,z1,cf1);
end;
cf1=sqrt(2.0d0.*pi)./ca0-cf1;
cpb(1+1+1)=cf1;
for  k=2:n0;
cf=(-z.*cf1+cf0)./(k-1.0d0);
cpb(k+1+1)=cf;
cf0=cf1;
cf1=cf;
end;  k=n0+1;
else;
if(a0 <= 3.0)
[dumvar1,z,cfa]=cpdsa(-n0,z,cfa);
cpb(n0+1+1)=cfa;
n1=n0+1;
[dumvar1,z,cfb]=cpdsa(-n1,z,cfb);
cpb(n1+1+1)=cfb;
nm1=n0-1;
for  k=nm1:-1:0;
cf=z.*cfa+(k+1.0d0).*cfb;
cpb(k+1+1)=cf;
cfb=cfa;
cfa=cf;
end;  k=0-1;
else;
m=100+abs(n);
cfa=c0;
cfb=complex(1.0d-30,0.0d0);
for  k=m:-1:0;
cf=z.*cfb+(k+1.0d0).*cfa;
if(k <= n0)
cpb(k+1+1)=cf;
end;
cfa=cfb;
cfb=cf;
end;  k=0-1;
cs0=ca0./cf;
for  k=0:n0;
cpb(k+1+1)=cs0.*cpb(k+1+1);
end;  k=n0+1;
end;
end;
end;
cpd(0+1+1)=-0.5d0.*z.*cpb(0+1+1);
if(n >= 0)
for  k=1:n;
cpd(k+1+1)=-0.5d0.*z.*cpb(k+1+1)+k.*cpb(k-1+1+1);
end;  k=n+1;
else;
for  k=1:n0;
cpd(k+1+1)=0.5d0.*z.*cpb(k+1+1)-cpb(k-1+1+1);
end;  k=n0+1;
end;
return;
end %subroutine cpbdn
function [n,z,cdn]=cpdsa(n,z,cdn);
%       ===========================================================
%       Purpose: Compute complex parabolic cylinder function Dn(z)
%                for small argument
%       Input:   z   --- complex argument of D(z)
%                n   --- Order of D(z) (n = 0,-1,-2,��
%       Output:  CDN --- Dn(z)
%       Routine called: GAIH for computing �x), x=n/2 (n=1,2,...)
%       ===========================================================
va0=[];ga0=[];xn=[];g1=[];vt=[];g0=[];vm=[];gm=[];
eps=1.0d-15;
pi=3.141592653589793d0;
sq2=sqrt(2.0d0);
ca0=exp(-.25d0.*z.*z);
va0=0.5d0.*(1.0d0-n);
if(n == 0.0)
cdn=ca0;
else;
if(abs(z) == 0.0)
if(va0 <= 0.0&&va0 == fix(va0))
cdn=0.0d0;
else;
[va0,ga0]=gaih(va0,ga0);
pd=sqrt(pi)./(2.0d0.^(-.5d0.*n).*ga0);
cdn=complex(pd,0.0d0);
end;
else;
xn=-n;
[xn,g1]=gaih(xn,g1);
cb0=2.0d0.^(-0.5d0.*n-1.0d0).*ca0./g1;
vt=-.5d0.*n;
[vt,g0]=gaih(vt,g0);
cdn=complex(g0,0.0d0);
cr=complex(1.0d0,0.0d0);
for  m=1:250;
vm=.5d0.*(m-n);
[vm,gm]=gaih(vm,gm);
cr=-cr.*sq2.*z./m;
cdw=gm.*cr;
cdn=cdn+cdw;
if(abs(cdw) < abs(cdn).*eps)
break;
end;
end;
cdn=cb0.*cdn;
end;
end;
return;
end %subroutine cpdsa
function [n,z,cdn]=cpdla(n,z,cdn);
%       ===========================================================
%       Purpose: Compute complex parabolic cylinder function Dn(z)
%                for large argument
%       Input:   z   --- Complex argument of Dn(z)
%                n   --- Order of Dn(z) (n = 0,�,�,��
%       Output:  CDN --- Dn(z)
%       ===========================================================

cb0=z.^n.*exp(-.25d0.*z.*z);
cr=complex(1.0d0,0.0d0);
cdn=complex(1.0d0,0.0d0);
for  k=1:16;
cr=-0.5d0.*cr.*(2.0.*k-n-1.0).*(2.0.*k-n-2.0)./(k.*z.*z);
cdn=cdn+cr;
if(abs(cr) < abs(cdn).*1.0d-12)
break;
end;
end;
cdn=cb0.*cdn;
return;
end %subroutine cpdla
function [x,ga]=gaih(x,ga);
%       =====================================================
%       Purpose: Compute gamma function �x)
%       Input :  x  --- Argument of �x), x = n/2, n=1,2,��!       Output:  GA --- �x)
%       =====================================================

pi=3.141592653589793d0;
if(x == fix(x)&&x > 0.0)
ga=1.0d0;
m1=fix(x-1.0);
for  k=2:m1;
ga=ga.*k;
end;  k=m1+1;
elseif(x+.5d0 == fix(x+.5d0)&&x > 0.0) ;
m=fix(x);
ga=sqrt(pi);
for  k=1:m;
ga=0.5d0.*ga.*(2.0d0.*k-1.0d0);
end;  k=m+1;
end;
return;
end %subroutine gaih

function out=writef(fid,varargin)
% function out=writef(fid,varargin)
%  Catches fortran stdout (6) and reroutes in to Matlab's stdout (1)
%  Catches fortran stderr (0) and reroutes in to Matlab's stderr (2)
if isnumeric(fid)
if fid==6,      out=fprintf(1,varargin{:});
elseif fid==0,  out=fprintf(2,varargin{:});
elseif isempty(fid) %% treat empty array like a string array [sethg 2008-03-03]
out=sprintf(varargin{:});
if nargin>2 %set the calling var to out
if ~isempty(inputname(1)), assignin('caller',inputname(1),out); end
end
else,           out=fprintf(fid,varargin{:});
end
elseif ischar(fid)
out=sprintf(varargin{:});
if nargin>2 %set the calling var to out
if ~isempty(inputname(1)), assignin('caller',inputname(1),out); end
end
else,            out=fprintf(fid,varargin{:});
end
end```