Can someone please convert this fortran program(.f90) to matlab (.m) for me. I would really appreciate
3 views (last 30 days)
Show older comments
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 Comments
Jan
on 19 Aug 2012
@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
on 20 Aug 2012
Formatting tip: tab is not recognized by the forum formatter as being whitespace.
Answers (1)
Ben Barrowes
on 23 Aug 2012
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 Comments
See Also
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!