Can someone please convert this fortran program(.f90) to matlab (.m) for me. I would really appreciate
조회 수: 3 (최근 30일)
이전 댓글 표시
real a,b,p,z,r,h,c,cprim,Dprim,F ,alpha,Fprim,Gprim,Hi ,Hprim
double precision k ,kprim,L,Lprim,M,Mprim , Ar(200),At(200),Br(200),Bt(200),Sigmar(200),Sigmat(200)
dimension r(200), h(200),alpha(200),deltaT(200),E(200),bnew(200)
open(1,file='Contrainte radial.dat')
open(2,file='Contrainte tangentia.dat')
write(1,*)'variable= "r(n)" "Sigmar(n)"'
write(2,*)'variable= "r(n)" "Sigmat(n)"'
write(1,*)'zone T=zone2'
write(*,*)'donner les valeurs de a b z'
read(*,*) a, b, z
p=(b-a)/(z-1)
r(1)=a
Sigmar(1)=0
write(,)'r= ',r(1) write(,)'Sigmar= ',Sigmar(1)
write(1,*) r(1),Sigmar(1) write(2,*) r(1), Sigmat(1)
do n=1,z r(n+1)=r(n)+p
h(n)=-0.002*r(n)**5+0.041668*r(n)**4-0.28432*r(n)**3+0.68273*r(n)**2-0.66228*r(n)+4.5901
bnew(n)=0.00015*r(n)**3-0.0009*r(n)**2+0.0018*r(n)+0.33
E(n)=-2.192e4*r(n)**3+1.334e5*r(n)**2-2.908e5*r(n)+2.366e7
alpha(n)=2.718e-9*r(n)**3-1.631e-8*r(n)**2+3.476e-8*r(n)+9.129e-6
deltaT(n)=1.689*r(n)**3-10.110*r(n)**2+21.630*r(n)+588
enddo
do n=2,z C=r(n)*h(n)
Cprim=bnew(n)/E(n) + ((1+ bnew(n))*(r(n)-r(n-1)))/(2*E(n)*r(n))
D=0.5*(r(n)-r(n-1))*h(n)
Dprim=1/E(n) + ((1+ bnew(n))*(r(n)-r(n-1)))/(2*E(n)*r(n))
F=r(n-1)*h(n-1)
Fprim= bnew(n-1)/E(n-1) - ((1+ bnew(n-1))*(r(n)-r(n-1)))/(2*E(n-1)*r(n-1))
G=0.5*(r(n)-r(n-1))*h(n-1)
Gprim=1/E(n-1) - ((1+ bnew(n-1))*(r(n)-r(n-1)))/(2*E(n-1)*r(n-1))
Hi=0.5*1005.7*(r(n)-r(n-1))*(h(n)*r(n)*r(n)+h(n-1)*r(n-1)*r(n-1)) Hprim=(alpha(n)*deltaT(n))-(alpha(n-1)*deltaT(n-1)) K=(Fprim*D-F*Dprim)/(Cprim*D-C*Dprim) Kprim=(Fprim*C-F*Cprim)/(Cprim*D-C*Dprim) L= -(Gprim*D+G*Dprim)/(Cprim*D-C*Dprim) Lprim = -(Gprim*C+G*Cprim)/(Cprim*D-C*Dprim) M=(Hprim*D+Hi*Dprim)/(Cprim*D-C*Dprim) Mprim= (Cprim*Hi+C*Hprim)/(Cprim*D-C*Dprim) Ar(1)=0 At(1)=1 Br(1)=0 Bt(1)=0 Sigmar(1)=0
Ar(n)=K*Ar(n-1)+L*At(n-1)
At(n)=Kprim*Ar(n-1)+Lprim*At(n-1)
Br(n)=K*Br(n-1)+L*Bt(n-1) + M
Bt(n)=Kprim*Br(n-1)+Lprim*Bt(n-1) + Mprim
enddo sigmar(z)=8500
Sigmat(1)=(sigmar(z)-Br(z))/Ar(z) write(,)'Sigmat= ',Sigmat(1) do n=2,z
Sigmar(n)= Ar(n)*Sigmat(1) + Br(n)
Sigmat(n)= At(n)*Sigmat(1) + Bt(n)
write(,)'____________________________________' write(,)'n= ', n write(,)'r= ', r(n) write(,)'C= ', c write(,)'Cprim= ', cprim write(,)'Dprim= ', Dprim WRITE(,)'F= ', F write(,)'Fprim= ', Fprim write(,)'Gprim= ', Gprim write(,)'G = ', G write(,)'Hi= ', Hi write(,)'Hprim= ', Hprim write(,)'k= ', k write(,)'kprim= ', kprim write(,)'L= ', L write(,)'Lprim= ', Lprim write(,)'M= ', M write(,)'Mprim= ', Mprim write(,)'Ar= ',Ar(n) write(,)'At= ',At(n) write(,)'Br= ',Br(n) write(,)'Bt= ',Bt(n) write(,)'r= ',r(n) write(,)'Sigmar= ',Sigmar(n) write(,)'Sigmat= ',Sigmat(n) write(1,*) r(n),Sigmar(n) write(2,*) r(n),Sigmat(n)
enddo end
댓글 수: 3
Jan
2012년 8월 19일
@Basil: It is unlikely that a voluntary contributor in this forum will do the work for you. Did you have success with installing f2matlab?
Walter Roberson
2012년 8월 20일
Formatting tip: tab is not recognized by the forum formatter as being whitespace.
답변 (1개)
Ben Barrowes
2012년 8월 23일
Putting it through f2matlab, this gets you most of the way there:
function basil(varargin)
clear all; %clear functions;
global unit2fid; if ~isempty(unit2fid), unit2fid=[]; end
persistent a alpha ar at b bnew br bt c cprim deltat dprim e f fprim gprim h hi hprim k kprim l lprim m mprim p r sigmar sigmat z ;
if isempty(a), a=0; end;
if isempty(b), b=0; end;
if isempty(p), p=0; end;
if isempty(z), z=0; end;
if isempty(r), r=zeros(1,200); end;
if isempty(h), h=zeros(1,200); end;
if isempty(c), c=0; end;
if isempty(cprim), cprim=0; end;
if isempty(dprim), dprim=0; end;
if isempty(f), f=0; end;
if isempty(alpha), alpha=zeros(1,200); end;
if isempty(fprim), fprim=0; end;
if isempty(gprim), gprim=0; end;
if isempty(hi), hi=0; end;
if isempty(hprim), hprim=0; end;
if isempty(k), k=0; end;
if isempty(kprim), kprim=0; end;
if isempty(l), l=0; end;
if isempty(lprim), lprim=0; end;
if isempty(m), m=0; end;
if isempty(mprim), mprim=0; end;
if isempty(ar), ar=zeros(1,200); end;
if isempty(at), at=zeros(1,200); end;
if isempty(br), br=zeros(1,200); end;
if isempty(bt), bt=zeros(1,200); end;
if isempty(sigmar), sigmar=zeros(1,200); end;
if isempty(sigmat), sigmat=zeros(1,200); end;
if isempty(deltat), deltat=zeros(1,200); end;
if isempty(e), e=zeros(1,200); end;
if isempty(bnew), bnew=zeros(1,200); end;
thismlfid=fopen('Contrainte radial.dat','r+');
unit2fid=[unit2fid;1,thismlfid];
thismlfid=fopen('Contrainte tangentia.dat','r+');
unit2fid=[unit2fid;2,thismlfid];
writef(1,['%s '],'variable= "r(n)" "Sigmar(n)"');
writef(2,['%s '],'variable= "r(n)" "Sigmat(n)"');
writef(1,['%s '],'zone T=zone2');
writef(1,['%s '],'donner les valeurs de a b z');
a=input('','s');
b=input('','s');
z=input('','s');
p=(b-a)./(z-1);
r(1)=a;
sigmar(1)=0;
writef(1,['%s %0.15g '],'r= ',r(1));
writef(1,['%s %0.15g '],'Sigmar= ',sigmar(1));
writef(1,['%0.15g %0.15g '], r(1),sigmar(1));
writef(2,['%0.15g %0.15g '], r(1), sigmat(1));
for n=1:z;
r(n+1)=r(n)+p;
h(n)=-0.002.*r(n).^5+0.041668.*r(n).^4-0.28432.*r(n).^3+0.68273.*r(n).^2-0.66228.*r(n)+4.5901;
bnew(n)=0.00015.*r(n).^3-0.0009.*r(n).^2+0.0018.*r(n)+0.33;
e(n)=-2.192e4.*r(n).^3+1.334e5.*r(n).^2-2.908e5.*r(n)+2.366e7;
alpha(n)=2.718e-9.*r(n).^3-1.631e-8.*r(n).^2+3.476e-8.*r(n)+9.129e-6;
deltat(n)=1.689.*r(n).^3-10.110.*r(n).^2+21.630.*r(n)+588;
end; n=z+1;
for n=2:z;
c=r(n).*h(n);
cprim=bnew(n)./e(n) +((1+ bnew(n)).*(r(n)-r(n-1)))./(2.*e(n).*r(n));
d=0.5.*(r(n)-r(n-1)).*h(n);
dprim=1./e(n) +((1+ bnew(n)).*(r(n)-r(n-1)))./(2.*e(n).*r(n));
f=r(n-1).*h(n-1);
fprim= bnew(n-1)./e(n-1) -((1+ bnew(n-1)).*(r(n)-r(n-1)))./(2.*e(n-1).*r(n-1));
g=0.5.*(r(n)-r(n-1)).*h(n-1);
gprim=1./e(n-1) -((1+ bnew(n-1)).*(r(n)-r(n-1)))./(2.*e(n-1).*r(n-1));
hi=0.5.*1005.7.*(r(n)-r(n-1)).*(h(n).*r(n).*r(n)+h(n-1).*r(n-1).*r(n-1));
hprim=(alpha(n).*deltat(n))-(alpha(n-1).*deltat(n-1));
k=(fprim.*d-f.*dprim)./(cprim.*d-c.*dprim);
kprim=(fprim.*c-f.*cprim)./(cprim.*d-c.*dprim);
l=-(gprim.*d+g.*dprim)./(cprim.*d-c.*dprim);
lprim = -(gprim.*c+g.*cprim)./(cprim.*d-c.*dprim);
m=(hprim.*d+hi.*dprim)./(cprim.*d-c.*dprim);
mprim=(cprim.*hi+c.*hprim)./(cprim.*d-c.*dprim);
ar(1)=0;
at(1)=1;
br(1)=0;
bt(1)=0;
sigmar(1)=0;
ar(n)=k.*ar(n-1)+l.*at(n-1);
at(n)=kprim.*ar(n-1)+lprim.*at(n-1);
br(n)=k.*br(n-1)+l.*bt(n-1) + m;
bt(n)=kprim.*br(n-1)+lprim.*bt(n-1) + mprim;
end; n=z+1;
sigmar(z)=8500;
sigmat(1)=(sigmar(z)-br(z))./ar(z);
writef(1,['%s %0.15g '],'Sigmat= ',sigmat(1));
for n=2:z;
sigmar(n)= ar(n).*sigmat(1) + br(n);
sigmat(n)= at(n).*sigmat(1) + bt(n);
writef(1,['%s '],'____________________________________');
writef(1,['%s %0.15g '],'n= ', n);
writef(1,['%s %0.15g '],'r= ', r(n));
writef(1,['%s %0.15g '],'C= ', c);
writef(1,['%s %0.15g '],'Cprim= ', cprim);
writef(1,['%s %0.15g '],'Dprim= ', dprim);
writef(1,['%s %0.15g '],'F= ', f);
writef(1,['%s %0.15g '],'Fprim= ', fprim);
writef(1,['%s %0.15g '],'Gprim= ', gprim);
writef(1,['%s %0.15g '],'G = ', g);
writef(1,['%s %0.15g '],'Hi= ', hi);
writef(1,['%s %0.15g '],'Hprim= ', hprim);
writef(1,['%s %0.15g '],'k= ', k);
writef(1,['%s %0.15g '],'kprim= ', kprim);
writef(1,['%s %0.15g '],'L= ', l);
writef(1,['%s %0.15g '],'Lprim= ', lprim);
writef(1,['%s %0.15g '],'M= ', m);
writef(1,['%s %0.15g '],'Mprim= ', mprim);
writef(1,['%s %0.15g '],'Ar= ',ar(n));
writef(1,['%s %0.15g '],'At= ',at(n));
writef(1,['%s %0.15g '],'Br= ',br(n));
writef(1,['%s %0.15g '],'Bt= ',bt(n));
writef(1,['%s %0.15g '],'r= ',r(n));
writef(1,['%s %0.15g '],'Sigmar= ',sigmar(n));
writef(1,['%s %0.15g '],'Sigmat= ',sigmat(n));
writef(1,['%0.15g %0.15g '], r(n),sigmar(n));
writef(2,['%0.15g %0.15g '], r(n),sigmat(n));
end; n=z+1;
end %program basil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%extra functions as needed by the translation %%%%%%%%%%%
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|fid==1, out=fprintf(1,varargin{:});fprintf(1,char(10));
elseif fid==0, out=fprintf(2,varargin{:});fprintf(2,char(10));
elseif isempty(fid) %%treat empty array like a string array
out=sprintf(varargin{:});
if nargin>2 %set the calling var to out
if ~isempty(inputname(1)), assignin('caller',inputname(1),out); end
end
else
%translate unitIn to fidIn
global unit2fid
fidRow=find(unit2fid(:,1)==fid,1,'last');
if isempty(fidRow), error(['unknown fid in readFmt',]); end
fidIn=unit2fid(fidRow,2);
out=fprintf(fidIn,varargin{:});fprintf(fidIn,char(10));
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
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Startup and Shutdown에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!